text:
# 1 Diagrammatic structure of the various factors that enter the process-independent resummation formula ().

The properties of the transverse-momentum distributions of systems of high mass $`Q`$ produced at high-energy hadron colliders are important for QCD studies and for physics studies beyond the Standard Model (see, e.g., Refs. ). On the theoretical side the computation of these distributions is complicated by the presence of large logarithmic contributions of the type $`\mathrm{ln}Q^2/q_{}^2`$, which spoil the convergence of fixed-order perturbative calculations in the region of small transverse momenta $`q_{}`$. The logarithmically-enhanced terms have to be evaluated at higher perturbative orders, and possibly resummed to all orders in the QCD coupling $`\alpha _\mathrm{S}`$. The all-order resummation formalism was developed in the eighties . The structure of the resummed distribution is given in terms of a transverse-momentum form factor and of process-dependent contributions.

The transverse-momentum form factor is sometimes supposed to be universal, that is, independent of the process. Performing a process-independent calculation at relative order $`\alpha _\mathrm{S}^2`$, two of us have recently shown  in an explicit way that the form factor depends on the process. In this note, we give a general interpretation of the results in Ref. . We point out that this process-dependent features persist to higher orders, at least in the final form in which the resummed transverse-momentum distribution is usually organized. We also present a new form of the resummed distribution, where all the logarithmic contributions are embodied in universal (process-independent) factors.

We consider the inclusive hard-scattering process

$$h_1(p_1)+h_2(p_2)F(Q^2,q_{}^2)+X,$$
(1)
where the triggered final-state system $`F`$ is produced by the collision of the two hadrons $`h_1`$ and $`h_2`$ with momenta $`p_1`$ and $`p_2`$, respectively. We denote by $`\sqrt{s}`$ the centre-of-mass energy of the colliding hadrons $`(s=(p_1+p_2)^22p_1p_2)`$. The final state $`F`$ is a generic system of non-strongly interacting particles We do not consider the production of strongly interacting particles (hadrons, jets, heavy quarks, …), because in this case the resummation formalism of small-$`q_{}`$ logarithms has not yet been fully developed., such as one or more vector bosons $`(\gamma ^{},W,Z,\mathrm{})`$, Higgs particles ($`H`$) and so forth.

To simplify the discussion we limit ourselves to the case in which only the total invariant mass $`Q`$ and transverse momentum $`q_{}`$ (with respect to the direction of the colliding hadrons) of the system $`F`$ are measured. The extension to more general kinematic configurations, in which, for instance, the rapidities are measured, is straightforward (see the discussion below Eq. (25)). We also assume that at the parton level the system $`F`$ is produced with vanishing $`q_{}`$ (i.e. with no accompanying final-state radiation) in the leading-order (LO) approximation, so that the corresponding cross section is $`d\sigma ^{(LO)}/dq_{}^2\delta (q_{}^2)`$. Since $`F`$ is colourless, the LO partonic subprocess is either $`q\overline{q}`$ annihilation, as in the case of $`\gamma ^{},W`$ and $`Z`$ production, or $`gg`$ fusion, as in the case of the production of a Higgs boson $`H`$.

The transverse-momentum cross section for the process in Eq. (1) can be written as

$$\frac{d\sigma _F}{dQ^2dq_{}^2}=\left[\frac{d\sigma _F}{dQ^2dq_{}^2}\right]_{\mathrm{res}.}+\left[\frac{d\sigma _F}{dQ^2dq_{}^2}\right]_{\mathrm{fin}.}.$$
(2)
Both terms on the right-hand side are obtained as convolutions of partonic cross sections and the parton distributions $`f_{a/h}(x,Q^2)`$ ($`a=q_f,\overline{q}_f,g`$ is the parton label) of the colliding hadronsThroughout the paper we always use parton densities as defined in the $`\overline{\mathrm{MS}}`$ factorization scheme and $`\alpha _\mathrm{S}(q^2)`$ is the QCD running coupling in the $`\overline{\mathrm{MS}}`$ renormalization scheme.. The distinction between the two terms is purely theoretical. The partonic cross section that enters in the resummed part (the first term on the right-hand side) contains all the logarithmically-enhanced contributions $`\alpha _\mathrm{S}^n\mathrm{ln}^mQ^2/q_{}^2`$. Thus, this part has to be evaluated by resumming the logarithmic terms to all orders in perturbation theory. On the contrary, the partonic cross section in the second term on the right-hand side is finite order-by-order in perturbation theory when $`q_{}0`$. It can thus be computed by truncating the perturbative expansion at a given fixed order in $`\alpha _\mathrm{S}`$.

The finite component of the transverse-momentum cross section is obviously process-dependent, and we have nothing to add on it in this paper. In the following we discuss the structure of the resummed part.

The resummed component is§§§As discussed at the end of the paper, this expression can be generalized to include the dependence on the renormalization and factorization scales $`\mu _R`$ and $`\mu _F`$.

$`\left[{\displaystyle \frac{d\sigma _F}{dQ^2dq_{}^2}}\right]_{\mathrm{res}.}`$ $`=`$ $`{\displaystyle \underset{a,b}{}}{\displaystyle _0^1}𝑑x_1{\displaystyle _0^1}𝑑x_2{\displaystyle _0^{\mathrm{}}}𝑑b{\displaystyle \frac{b}{2}}J_0(bq_{})f_{a/h_1}(x_1,b_0^2/b^2)f_{b/h_2}(x_2,b_0^2/b^2)`$ (3)
$``$ $`W_{ab}^F(x_1x_2s;Q,b).`$
The Bessel function $`J_0(bq_{})`$ and the coefficient $`b_0=2e^{\gamma _E}`$ ($`\gamma _E=0.5772\mathrm{}`$ is the Euler number) have a kinematical origin. To correctly take into account the kinematics constraint of transverse-momentum conservation, the resummation procedure has to be carried out in the impact-parameter space. The transverse-momentum cross section (3) is then obtained by performing the inverse Fourier (Bessel) transformation with respect to the impact parameter $`b`$. The factor $`W_{ab}^F`$ is the perturbative and process-dependent partonic cross section that embodies the all-order resummation of the large logarithms $`\mathrm{ln}Q^2b^2`$ (the limit $`q_{}Q`$ corresponds to $`Qb1`$, because $`b`$ is the variable conjugate to $`q_{}`$).

The resummed partonic cross section is usually (see, e.g., the list of references in Sections 5.1 and 5.3 of Ref. ) written in the following form:

$`W_{ab}^F(s;Q,b)`$ $`=`$ $`{\displaystyle \underset{c}{}}{\displaystyle _0^1}𝑑z_1{\displaystyle _0^1}𝑑z_2C_{ca}^F(\alpha _\mathrm{S}(b_0^2/b^2),z_1)C_{\overline{c}b}^F(\alpha _\mathrm{S}(b_0^2/b^2),z_2)\delta (Q^2z_1z_2s)`$ (4)
$``$ $`\sigma _{c\overline{c}}^{(LO)F}(Q^2)S_c^F(Q,b).`$
Here, $`\sigma _{c\overline{c}}^{(LO)F}`$ is the cross section (integrated over $`q_{}`$) for the LO partonic subprocess $`c+\overline{c}F`$, where $`c,\overline{c}=q,\overline{q}`$ (the quark $`q_f`$ and the antiquark $`\overline{q}_f^{}`$ can possibly have different flavours $`f,f^{}`$) or $`c,\overline{c}=g,g`$. The expression $`\sigma _{c\overline{c}}^{(LO)F}`$ can include an overall factor $`\alpha _\mathrm{S}^p(Q^2)`$, as in the case of $`g+gH`$ through a triangular quark loop where $`p=2`$. The term $`S_c^F(Q,b)`$ is the quark $`(c=q)`$ or gluon $`(c=g)`$ Sudakov form factor. The resummation of the logarithmic contributions is achieved by exponentiation , that is by showing  that the form factor can be expressed as

$$S_c(Q,b)=\mathrm{exp}\left\{_{b_0^2/b^2}^{Q^2}\frac{dq^2}{q^2}\left[A_c(\alpha _\mathrm{S}(q^2))\mathrm{ln}\frac{Q^2}{q^2}+B_c(\alpha _\mathrm{S}(q^2))\right]\right\},$$
(5)
with $`c=q`$ or $`g`$. The functions $`A_c(\alpha _\mathrm{S}),B_c(\alpha _\mathrm{S})`$, as well as the coefficient functions $`C_{ab}(\alpha _\mathrm{S},z)`$ in Eq. (4), contain no $`\mathrm{ln}Q^2b^2`$ terms and are perturbatively computable according to the power expansionsNote that in Refs.  the perturbative coefficients are normalized to powers of $`\alpha _\mathrm{S}/2\pi `$ rather than $`\alpha _\mathrm{S}/\pi `$.

$`A_c(\alpha _\mathrm{S})`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\alpha _\mathrm{S}}{\pi }}\right)^nA_c^{(n)},`$ (6)
$`B_c(\alpha _\mathrm{S})`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\alpha _\mathrm{S}}{\pi }}\right)^nB_c^{(n)},`$ (7)
$`C_{ab}(\alpha _\mathrm{S},z)`$ $`=`$ $`\delta _{ab}\delta (1z)+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\alpha _\mathrm{S}}{\pi }}\right)^nC_{ab}^{(n)}(z).`$ (8)
The knowledge of the coefficients $`A^{(1)}`$ leads to the resummation of the leading logarithmic (LL) contributions. Analogously, the coefficients $`\{A^{(2)},B^{(1)},C^{(1)}\}`$ give the next-to-leading logarithmic (NLL) terms, the coefficients $`\{A^{(3)},B^{(2)},C^{(2)}\}`$ give the next-to-next-to-leading logarithmic (NNLL) terms, and so forth. The coefficients $`A^{(1)},A^{(2)},B^{(1)}`$ are known both for the quark  and for the gluon  form factors

$`A_q^{(1)}=C_F,A_g^{(1)}=C_A,`$
$`A_q^{(2)}={\displaystyle \frac{1}{2}}C_FK,A_g^{(2)}={\displaystyle \frac{1}{2}}C_AK,`$ (9)
$`B_q^{(1)}={\displaystyle \frac{3}{2}}C_F,B_g^{(1)}={\displaystyle \frac{1}{6}}(11C_A2N_f),`$
where

$$K=\left(\frac{67}{18}\frac{\pi ^2}{6}\right)C_A\frac{5}{9}N_f.$$
(10)
The best studied example among the processes in Eq. (1) is lepton-pair Drell–Yan (DY) production through the LO partonic subprocess $`q+\overline{q}\gamma ^{}(\mathrm{or}W,Z)l+l^{}`$. In this case the NNLL coefficient $`B^{(2)}`$ was computed by Davies and Stirling :

$$B_q^{(2)DY}=C_F^2\left(\frac{\pi ^2}{4}\frac{3}{16}3\zeta _3\right)+C_FC_A\left(\frac{11}{36}\pi ^2\frac{193}{48}+\frac{3}{2}\zeta _3\right)+C_FN_f\left(\frac{17}{24}\frac{\pi ^2}{18}\right),$$
(11)
where $`\zeta _n`$ is the Riemann $`\zeta `$-function $`(\zeta _3=1.202\mathrm{})`$. The coefficient $`B^{(2)}`$ for Higgs boson production has recently been computed by two of us . The corresponding LO partonic subprocess is gluon fusion, $`g+gH`$, through a massive-quark loop. In the limit of infinite quark mass the value of $`B^{(2)}`$ is

$$B_g^{(2)H}=C_A^2\left(\frac{23}{24}+\frac{11}{18}\pi ^2\frac{3}{2}\zeta _3\right)+\frac{1}{2}C_FN_fC_AN_f\left(\frac{1}{12}+\frac{\pi ^2}{9}\right)\frac{11}{8}C_FC_A.$$
(12)
The main issue that we want to discuss in this paper regards the process dependence of the various factors in the resummation formula (4). As denoted by the superscripts $`F`$ in the various terms on the right-hand side, the coefficient functions $`C_{ab}^F`$ depend on the process. This is confirmed by the calculations of the coefficients $`C_{ab}^{(1)F}`$, performed in the literature for several processes A general expression for the coefficients $`C_{ab}^{(1)F}`$ in terms of the one-loop matrix element of the corresponding process is given in Eq. (17) of Ref. .  . The form factor $`S_c(Q,b)`$ that enters Eq. (4) is (often) supposed to be universal (this is the reason why it is named quark or gluon form factor rather than DY, $`\gamma \gamma `$, $`WZ`$, $`H,\mathrm{}`$ form factor). However, this is not the case: the form factor $`S_c^F(Q,b)`$ in Eq. (4) is process-dependent. In the following, we first present a universal (process-independent) version of the resummation formula (4) and we sketch its physical origin. We then clarify the relation between Eq. (4) and our process-independent version.

The process-independent resummation formula is

$`W_{ab}^F(s;Q,b)`$ $`=`$ $`{\displaystyle \underset{c}{}}{\displaystyle _0^1}𝑑z_1{\displaystyle _0^1}𝑑z_2C_{ca}(\alpha _\mathrm{S}(b_0^2/b^2),z_1)C_{\overline{c}b}(\alpha _\mathrm{S}(b_0^2/b^2),z_2)\delta (Q^2z_1z_2s)`$ (13)
$``$ $`\sigma _{c\overline{c}}^F(Q^2,\alpha _\mathrm{S}(Q^2))S_c(Q,b).`$
It formally differs from Eq. (4) by the replacement $`\sigma _{c\overline{c}}^{(LO)F}(Q^2)\sigma _{c\overline{c}}^F(Q^2,\alpha _\mathrm{S}(Q^2))`$. While $`\sigma _{c\overline{c}}^{(LO)F}(Q^2)`$ is the cross section for the LO partonic subprocess, $`\sigma _{c\overline{c}}^F(Q^2,\alpha _\mathrm{S}(Q^2))`$ includes higher-order QCD corrections to it, according to

$$\sigma _{c\overline{c}}^F(Q^2,\alpha _\mathrm{S}(Q^2))=\sigma _{c\overline{c}}^{(LO)F}(Q^2)H_c^F(\alpha _\mathrm{S}(Q^2)),$$
(14)
where the function $`H_c^F(\alpha _\mathrm{S})`$ has a perturbative expansion similar to Eqs. (6)–(8):

$$H_c^F(\alpha _\mathrm{S})=1+\underset{n=1}{\overset{\mathrm{}}{}}\left(\frac{\alpha _\mathrm{S}}{\pi }\right)^nH_c^{(n)F}.$$
(15)
Note that the function $`H_c^F(\alpha _\mathrm{S})`$ depends on the process. Nonetheless, its introduction is sufficient to transform the process-dependent form factor $`S_c^F`$ and coefficient functions $`C_{ca}^F`$ of Eq. (4) into the process-independent form factor $`S_c`$ and coefficient functions $`C_{ca}`$ of Eq. (13).

The resummation formula in Eq. (13), which can be derived by the customary resummation methods  , has a simple physical origin. When the final-state system $`F`$ is kinematically constrained to have a small transverse momentum, the emission of accompanying radiation is strongly inhibited, so that only soft and collinear partons (i.e. partons with low transverse momenta $`q_t`$) can be radiated in the final state (Fig. 1). The process-dependent factor $`H_c^F(\alpha _\mathrm{S}(Q^2))`$ embodies hard contributions produced by virtual corrections at transverse-momentum scales $`q_tQ`$. The form factor $`S_c(Q,b)`$ contains real and virtual contributions due to soft (the function $`A_c(\alpha _\mathrm{S})`$ in Eq. (5)) and flavour-conserving collinear (the function $`B_c(\alpha _\mathrm{S})`$ in Eq. (5)) radiation at scales $`Q>q_t>1/b`$. At very low momentum scales, $`q_t<1/b`$, real and virtual soft-gluon corrections cancel because the cross section is infrared safe, and only real and virtual contributions due to collinear radiation remain (the coefficient functions $`C_{ab}(\alpha _\mathrm{S}(b_0^2/b^2),z)`$). Note that $`S_c(Q,b)`$ and $`C_{ab}(\alpha _\mathrm{S}(b_0^2/b^2),z)`$ are process-independent and only depend on the flavour and colour charges of the QCD partons.

The two versions (4) and (13) of the resummation formula can formally be related as follows. We use the renormalization-group identity

$$g_1(\alpha _\mathrm{S}(Q^2))=\mathrm{exp}\left\{_{b_0^2/b^2}^{Q^2}\frac{dq^2}{q^2}g_2(\alpha _\mathrm{S}(q^2))\right\}g_1(\alpha _\mathrm{S}(b_0^2/b^2)),$$
(16)
which is valid when

$$g_2(\alpha _\mathrm{S})=\beta (\alpha _\mathrm{S})\frac{d\mathrm{ln}g_1(\alpha _\mathrm{S})}{d\mathrm{ln}\alpha _\mathrm{S}},$$
(17)
where $`\beta (\alpha _\mathrm{S})`$ is the QCD $`\beta `$-function

$$\frac{d\mathrm{ln}\alpha _\mathrm{S}(q^2)}{d\mathrm{ln}q^2}=\beta (\alpha _\mathrm{S}(q^2)),$$
(18)
$$\beta (\alpha _\mathrm{S})=\beta _0\frac{\alpha _\mathrm{S}}{\pi }\beta _1\left(\frac{\alpha _\mathrm{S}}{\pi }\right)^2+\mathrm{},12\beta _0=11C_A2N_f.$$
(19)
Then, setting $`g_1(\alpha _\mathrm{S}(Q^2))=H_c^F(\alpha _\mathrm{S}(Q^2))`$ and inserting the right-hand side of Eq. (16) in Eq. (13), we immediately obtain Eq. (4). More precisely, the process-independent resummation formula in Eq. (13) implies the customary version in Eq. (4), provided the perturbative function $`B_c(\alpha _\mathrm{S})`$ in the form factor $`S_c^F`$ (see Eq. (5)) and the coefficient functions $`C_{ab}^F`$ are related to their process-independent analogues by the following all-order relations

$`C_{ab}^F(\alpha _\mathrm{S},z)`$ $`=`$ $`\left[H_a^F(\alpha _\mathrm{S})\right]^{1/2}C_{ab}(\alpha _\mathrm{S},z),`$ (20)
$`B_c^F(\alpha _\mathrm{S})`$ $`=`$ $`B_c(\alpha _\mathrm{S})\beta (\alpha _\mathrm{S}){\displaystyle \frac{d\mathrm{ln}H_c^F(\alpha _\mathrm{S})}{d\mathrm{ln}\alpha _\mathrm{S}}}.`$ (21)
While the perturbative function $`A_c(\alpha _\mathrm{S})`$ and the first-order coefficient $`B_c^{(1)}`$ of the function $`B_c(\alpha _\mathrm{S})`$ are process-independent, the result in Eqs. (21) and (20) shows that the coefficients $`B_c^{(2)},B_c^{(3)},\mathrm{}`$ and $`C_{ab}^{(1)}(z),C_{ab}^{(2)}(z),\mathrm{}`$ in Eq. (4) do depend on the process. The process dependence of the first few coefficients is explicitly given by

$`C_{ab}^{(1)F}(z)`$ $`=`$ $`C_{ab}^{(1)}(z)+\delta _{ab}\delta (1z){\displaystyle \frac{1}{2}}H_a^{(1)F},`$ (22)
$`C_{ab}^{(2)F}(z)`$ $`=`$ $`C_{ab}^{(2)}(z)+{\displaystyle \frac{1}{2}}H_a^{(1)F}C_{ab}^{(1)}(z)+\delta _{ab}\delta (1z){\displaystyle \frac{1}{2}}\left[H_a^{(2)F}{\displaystyle \frac{1}{4}}\left(H_a^{(1)F}\right)^2\right],`$ (23)
$`B_c^{(2)F}`$ $`=`$ $`B_c^{(2)}+\beta _0H_c^{(1)F},`$ (24)
$`B_c^{(3)F}`$ $`=`$ $`B_c^{(3)}+\beta _1H_c^{(1)F}+2\beta _0\left[H_c^{(2)F}{\displaystyle \frac{1}{2}}\left(H_c^{(1)F}\right)^2\right].`$ (25)
Note that the process dependence of the resummation formula (4) is not simply embodied in pure numerical coefficients. For example, considering the Higgs production mechanism $`g+gH`$ through a massive-quark loop, the corresponding coefficient $`H_g^{(1)H}`$ and, hence, the form-factor coefficient $`B_g^{(2)H}`$ depend on the mass of the quark in the triangular loop. Moreover, the results discussed so far can straightforwardly be extended to more general configurations, in which several kinematical variables (and not only the invariant mass $`Q`$) of the final-state system $`F`$ are measured. For example, when $`F`$ is a pair of vector bosons, we can consider the corresponding transverse-momentum cross section at fixed rapidities of the vector bosons. The extension simply amounts to including the dependence on the kinematics of the system $`F`$ in the cross section factor $`\sigma _{c\overline{c}}^F`$ (and, hence, in $`\sigma _{c\overline{c}}^{(LO)F}`$ and $`H_c^F`$) of the process-independent resummation formula (13). In this general case, Eqs. (21) and (20) imply that the form factor $`S_c^F`$ and the coefficient functions $`C_{ab}^F`$ in Eq. (4) depend on the kinematics of the system $`F`$ in a non-trivial manner.

The process-independent resummation formula (13) can alternatively be used to relate the factors on the right-hand side of Eq. (4) for different processes. Considering Eqs. (21) and (20) for two different processes $`F`$ and $`\stackrel{~}{F}`$, we obtain

$`C_{cb}^F(\alpha _\mathrm{S},z)`$ $`=`$ $`\left[H_c^{F\stackrel{~}{F}}(\alpha _\mathrm{S})\right]^{1/2}C_{cb}^{\stackrel{~}{F}}(\alpha _\mathrm{S},z),`$ (26)
$`B_c^F(\alpha _\mathrm{S})`$ $`=`$ $`B_c^{\stackrel{~}{F}}(\alpha _\mathrm{S})\beta (\alpha _\mathrm{S}){\displaystyle \frac{d\mathrm{ln}H_c^{F\stackrel{~}{F}}(\alpha _\mathrm{S})}{d\mathrm{ln}\alpha _\mathrm{S}}},`$ (27)
where $`H^{F\stackrel{~}{F}}=H^F/H^{\stackrel{~}{F}}`$. The perturbative expansion of these equations obviously leads to relations that are similar to those in Eqs. (22)–(24):

$`C_{cb}^{(1)F}(z)`$ $`=`$ $`C_{cb}^{(1)\stackrel{~}{F}}(z)+\delta _{cb}\delta (1z){\displaystyle \frac{1}{2}}H_c^{(1)F\stackrel{~}{F}},`$ (28)
$`C_{cb}^{(2)F}(z)`$ $`=`$ $`C_{cb}^{(2)\stackrel{~}{F}}(z)+{\displaystyle \frac{1}{2}}H_c^{(1)F\stackrel{~}{F}}C_{cb}^{(1)\stackrel{~}{F}}(z)+\delta _{cb}\delta (1z){\displaystyle \frac{1}{2}}\left[H_c^{(2)F\stackrel{~}{F}}{\displaystyle \frac{1}{4}}\left(H_c^{(1)F\stackrel{~}{F}}\right)^2\right],`$ (29)
$`B_c^{(2)F}`$ $`=`$ $`B_c^{(2)\stackrel{~}{F}}+\beta _0H_c^{(1)F\stackrel{~}{F}},`$ (30)
$`B_c^{(3)F}`$ $`=`$ $`B_c^{(3)\stackrel{~}{F}}+\beta _1H_c^{(1)F\stackrel{~}{F}}+2\beta _0\left[H_c^{(2)F\stackrel{~}{F}}{\displaystyle \frac{1}{2}}\left(H_c^{(1)F\stackrel{~}{F}}\right)^2\right].`$ (31)
Since the function $`H^{F\stackrel{~}{F}}`$ does not appear in Eq. (4), within the framework of the process-dependent resummation formulae, Eqs. (26)–(31) have to be regarded as non-linear relations between the form factors and coefficient functions of different processes. For instance, having computed the first-order coefficient $`C_{ab}^{(1)}`$ for two processes, one can check Eq. (28) and extract $`H_c^{(1)F\stackrel{~}{F}}`$. Then, the second-order coefficients $`B_c^{(2)}`$ of these processes must be related by Eq. (30). This constraint can either be checked or used to compute $`B_c^{(2)\stackrel{~}{F}}`$ from $`B_c^{(2)F}`$.

The relation in Eq. (28) (or Eq. (22)) can straightforwardly be checked by comparing the known first-order coefficients $`C_{ab}^{(1)F}(z)`$ for the production of DY lepton pairs , diphotons  and $`ZZ`$ pairs . In particular, we have

$$H_q^{(1)FDY}=H_q^{(1)F}H_q^{(1)DY}=\frac{1}{2}C_F\left(𝒱^F\pi ^2+8\right),F=\gamma \gamma ,ZZ,$$
(32)
where $`𝒱^{\gamma \gamma }`$ and $`𝒱^{ZZ}`$ are given in Eq. (11) of Ref.  and in Eq. (8)\**\**\**As pointed out in Ref. , the actual expression in Eq. (8) of Ref.  is not correct. of Ref. , respectively. The relation (30) (or (24)) for the NNLL coefficient $`B^{(2)F}`$ has explicitly been derived by the process-independent calculation at $`𝒪(\alpha _\mathrm{S}^2)`$ performed in Ref. .

The reasoning used to obtain Eq. (4) from Eq. (13) can also be used to show that Eq. (13) is invariant under the transformation

$`H_c^F(\alpha _\mathrm{S}(Q^2))`$ $``$ $`H_c^F(\alpha _\mathrm{S}(Q^2))\left[g(\alpha _\mathrm{S}(Q^2))\right]^1,`$
$`S_c(Q,b)`$ $``$ $`S_c(Q,b)\mathrm{exp}\left\{{\displaystyle _{b_0^2/b^2}^{Q^2}}{\displaystyle \frac{dq^2}{q^2}}\beta (\alpha _\mathrm{S}(q^2)){\displaystyle \frac{d\mathrm{ln}g(\alpha _\mathrm{S}(q^2))}{d\mathrm{ln}\alpha _\mathrm{S}(q^2)}}\right\},`$ (33)
$`C_{ab}(\alpha _\mathrm{S}(b_0^2/b^2),z)`$ $``$ $`C_{ab}(\alpha _\mathrm{S}(b_0^2/b^2),z)\left[g(\alpha _\mathrm{S}(b_0^2/b^2))\right]^{1/2},`$
where $`g(\alpha _\mathrm{S})=1+𝒪(\alpha _\mathrm{S})`$ is an arbitrary perturbative function. This renormalization-group symmetry of the resummation formula (13) implies that its factors, $`H_c^F,S_c`$ (more precisely, the function $`B_c`$) and $`C_{ab}`$, although they are process-independent, are not unambiguously computable (defined) order by order in perturbation theory. The same conclusion can be reached by expanding Eq. (13) in powers of $`\alpha _\mathrm{S}`$. At any given order in $`\alpha _\mathrm{S}`$, the coefficients of the various powers of $`\mathrm{ln}Q^2b^2`$ depend on the unknowns $`H_c^{(n)F},C_{ab}^{(n)},B_c^{(n)}`$, but the number of such coefficients is less that the number of unknowns. For instance, at the first relative order in $`\alpha _\mathrm{S}`$, the coefficient of $`\mathrm{ln}^2Q^2b^2`$ determines $`A^{(1)}`$, the coefficient of $`\mathrm{ln}Q^2b^2`$ determines $`B^{(1)}`$, while the coefficient of the constant term is a linear combination of $`H^{(1)F}`$ and $`C_{ab}^{(1)}`$.

The perturbative ambiguity of the decomposition of the right-hand side of Eq. (13) in the factors $`H_c^F,S_c`$ and $`C_{ab}`$ is a consequence of the fact that the transverse-momentum cross section is not a collinear-safe quantity. The effect of collinear radiation at low transverse-momentum scales, $`q_t<1/b`$, is divergent in perturbation theory. The (arbitrary) regularization procedure of these divergences introduces some ambiguity in the definition of the coefficient functions $`C_{ab}`$. Then, the ambiguity propagates to $`S_c`$ and $`H_c^F`$ through collinear evolution (see Eq. (16)). In this way of thinking, the ambiguity is similar to that encountered in the definition of the parton densities. As the parton densities have to be defined by fixing a factorization scheme (e.g. the $`\overline{\mathrm{MS}}`$ scheme or the DIS scheme), the factors on the right-hand side of Eq. (13) has to be defined by choosing a ‘resummation scheme’. Note that the choice of a ‘resummation scheme’ amounts to defining $`H_c^F`$ (or $`C_{ab}`$) for a single process***More precisely, $`H_c^F`$ has to be defined for two processes: one process that is controlled, at LO, by $`q\overline{q}`$ annihilation and another process that is controlled, at LO, by $`gg`$ fusion.. Having done that, the process-dependent factor $`H_c^F`$ and the universal factors $`S_c`$ and $`C_{ab}`$ in Eq. (13) are unambiguously determined for any other process in Eq. (1).

We can suggest three examples of ‘resummation schemes’, that is three prescriptions to use the process-independent factorization formula (13) in practice. A first ‘short-cut’ (because the DY transverse-momentum distribution is best studied) prescription, which we name ‘DY resummation scheme’, is obtained by setting $`H_q^{DY}(\alpha _\mathrm{S})1`$ in the case of the DY cross section (integrated over the rapidity of the vector boson). This unambiguously fixes the process-independent form factor $`S_c`$ and coefficients functions $`C_{ab}`$ as those determined from the DY process. In particular, the coefficient $`B_q^{(2)}`$ of the quark form factor is that in Eq. (11), and the quark coefficient functions $`C_{qb}^{(1)}(z)`$ $`(b=q,\overline{q},g)`$ are those computed in Refs.  for the DY process. Any other $`q\overline{q}`$-initiated cross section is then obtained by using Eq. (13) with a perturbatively computable hard function $`H_q^F(\alpha _\mathrm{S})`$. This prescription can be extended to $`gg`$-initiated cross sections by setting $`H_g^H(\alpha _\mathrm{S})1`$ in the case of Higgs production via $`gg`$ fusion. Thus, in the case of an infinitely-massive quark in the loop (the definition of the scheme will be different by keeping the mass of the quark finite), the coefficient $`B_g^{(2)}`$ of the gluon form factor is that given in Eq. (12) and the gluon coefficient functions $`C_{gb}^{(1)}(z)`$ $`(b=q,\overline{q},g)`$ are those computed in Ref. .

A second prescription, which is in some sense more physical, exploits the fact that the first moment (with respect to $`Q^2/s`$) of the flavour non-singlet (NS) hadronic (and partonic) cross sections defines collinear-safe quantities. This property follows from fermion-number conservationIn the case of the Altarelli–Parisi evolution of the parton densities, fermion-number conservation implies the vanishing of the first moment of the Altarelli–Parisi probabilities.. Thus, we can define a ‘NS resummation scheme’ by fixing the overall normalization of the NS coefficient function $`C_{ab}(\alpha _\mathrm{S},z)`$ in Eq. (13) in such a way that its first moment (with respect to $`z`$) vanishes. In this scheme (the subscript NS labels the scheme choice)

$$_0^1𝑑z\left[C_{qq,NS}(\alpha _\mathrm{S},z)C_{q\overline{q},NS}(\alpha _\mathrm{S},z)\right]1,$$
(34)
and thus, using the DY coefficients functions from Refs.  and Eq. (22), we obtain

$$H_{q,NS}^{(1)DY}=C_F\left(\frac{\pi ^2}{2}\frac{7}{2}\right).$$
(35)
Inserting Eq. (35) in Eq. (24), and using the DY coefficient in Eq. (11), we then obtain the corresponding coefficient of the quark form factor:

$`B_{q,NS}^{(2)}`$ $`=`$ $`B_q^{(2)DY}\beta _0H_{q,NS}^{(1)DY}`$ (36)
$`=`$ $`C_F^2\left({\displaystyle \frac{\pi ^2}{4}}{\displaystyle \frac{3}{16}}3\zeta _3\right)+C_FC_A\left({\displaystyle \frac{11}{72}}\pi ^2{\displaystyle \frac{13}{16}}+{\displaystyle \frac{3}{2}}\zeta _3\right)+C_FN_f\left({\displaystyle \frac{1}{8}}+{\displaystyle \frac{\pi ^2}{36}}\right).`$
The reason why this scheme can be considered ‘more physical’ is that the first moments of the NS cross sections are collinear-safe and, hence, free from ambiguities related to the regularization procedure of collinear singularities. As is known , the process-independent form factor coefficients in Eq. (9) measure the physical intensity of soft ($`A_c^{(1)}`$ and $`A_c^{(2)}`$) and collinear $`(B_c^{(1)})`$ radiation. Since the coefficient $`B_{q,NS}^{(2)}`$ in Eq. (36) appears in the form factor of an infrared and collinear safe quantity, it should measure the intensity of the collinear radiation from quarks at the second order in $`\alpha _\mathrm{S}`$ and it should enter in the form factor of other infrared and collinear safe quantities. For example, an analogous transverse-momentum form factor controls the back-to-back region of the energy–energy correlation (EEC) in $`e^+e^{}`$ annihilation . The coefficient in Eq. (36) should directly be related to the corresponding coefficient $`B^{(2)}`$ of the EEC, modulo (possible) effects due to the difference between the space-like kinematics of the transverse-momentum cross sections and the time-like kinematics of $`e^+e^{}`$ annihilation.

The generalization of the ‘NS resummation scheme’ to $`gg`$-fusion processes can be obtained, for instance, by fixing the overall normalization of the gluon coefficient function $`C_{gg}(\alpha _\mathrm{S},z)`$ in Eq. (13) in such a way that its first moment vanishes. Note, however, that in the gluon (or, more generally, flavour-singlet) channel there is no anologue of the fermion-number conservation rule. Flavour-singlet transverse-momentum cross sections are collinear unsafe quantities. Thus, a ‘physical’ (collinear-safe) interpretation of the higher-order coefficients of the function $`B_g(\alpha _\mathrm{S})`$ of the gluon form factor is less straightforward than in the quark channel.

The coefficient $`B_c^{(1)}`$ in Eq. (9) coincides with the coefficient of the end-point contribution (the term proportional to $`\delta (1z)`$) to the LO Altarelli–Parisi probability $`P_{cc}^{(1)}(z)`$. As can be noticed from Eqs. (22) and (23) in Ref. , the coefficient $`B_{q,NS}^{(2)}`$ in Eq. (36) (and the coefficients $`B_q^{(2)DY}`$ and $`B_g^{(2)H}`$ in Eqs. (11) and (12)) does not coincide with that of the end-point contribution to the NLO Altarelli–Parisi probability $`P_{cc}^{(2)\overline{\mathrm{MS}}}(z)`$ in the $`\overline{\mathrm{MS}}`$ factorization scheme. This should not be surprising, because beyond LO the Altarelli–Parisi probabilities are not collinear-safe quantities: they depend on the regularization procedure of the collinear singularities and on the factorization scheme. This interpretation is confirmed by the fact that the various coefficients $`B_c^{(2)}`$ differ from one another by terms proportional to the first coefficient $`\beta _0`$ of the $`\beta `$-function, as should be expected from factorization-scheme dependence. Using Eqs. (20)–(25), it is obviously possible to introduce an ‘$`\overline{\mathrm{MS}}`$ resummation scheme’ such that the function $`B_c(\alpha _\mathrm{S})`$ of the process-independent form factor in Eqs. (5) and (13) coincides by definition with the perturbative function $`B_{c,\overline{\mathrm{MS}}}(\alpha _\mathrm{S})`$ that controls the end-point contribution to the all-order Altarelli–Parisi probabilities in the $`\overline{\mathrm{MS}}`$ factorization scheme.

Note that our discussion on the relation between Eqs. (4) and (13) and on the resummation-scheme ambiguity does not imply that the two equations are equivalent. In fact, the process-independent resummation framework contains more information than its customary process-dependent version. To evaluate resummed cross section within the framework of Eq. (4), we should compute two functions ($`S_c^F`$ and $`C_{ab}^F`$) for each process $`F`$. Once a ‘resummation scheme’ has been defined, a similar evaluation by using Eq. (13) requires the computation of two universal functions ($`S_c`$ and $`C_{ab}`$) and of a single additional function $`(H_c^F)`$ for each process.

Of course, we can still continue to use Eq. (4). In this case, we can exploit the additional information contained in Eq. (13) by using Eqs. (26) and (27) (or their perturbative expansions in Eqs. (28)–(31)) to relate form factors and coefficient functions for different processes.

Independently of the resummation formula that is actually used, once the form factor coefficient $`B_c^{(2)}`$ has been computed for a single process by performing a calculation at relative order $`\alpha _\mathrm{S}^2`$, no further $`𝒪(\alpha _\mathrm{S}^2)`$-calculation is necessary for its universal implementation in all processes. The implementation can simply be performed by computing the coefficient functions at $`𝒪(\alpha _\mathrm{S})`$ and then using Eqs. (22) and (24) (or Eqs. (28) and (30)). Analogously, using Eqs. (23) and (25) (or Eqs. (29) and (31)), the universal implementation of the form factor coefficient $`B_c^{(3)}`$ requires its $`𝒪(\alpha _\mathrm{S}^3)`$ computation in a single process and the evaluation of the coefficient functions at $`𝒪(\alpha _\mathrm{S}^2)`$ for the various processes. This procedure extends to any higher orders.

At present, owing to the results in Eqs. (11) and (12), the NNLL coefficient $`B_c^{(2)}`$ can be included in all calculations of transverse-momentum cross sections in hadron collisions. Its inclusion cannot be regarded as the extension of these resummed calculations to full NNLL accuracy, because the corresponding coefficient $`A_c^{(3)}`$ is not yet known. Nonetheless, the knowledge of $`B_c^{(2)}`$ can certainly be used to improve the matching with fixed-order calculations at high $`q_{}`$. The matching can be performed by supplementing the resummed cross section with the finite contribution on the right-hand side of Eq. (2). The finite contribution is obtained, from the complete calculation at a fixed order in $`\alpha _\mathrm{S}`$, by subtraction of the corresponding perturbative terms already included in the resummed component. If $`B_c^{(2)}`$ is included in this component, the remaining finite component at relative order $`\alpha _\mathrm{S}^2`$ contains no $`\mathrm{ln}Q^2/q_{}^2`$ terms. Thus, the latter can be inserted in Eq. (2) uniformly with respect to $`q_{}`$, that is without including logarithmic terms that spoil the convergence of the matched result in the small-$`q_{}`$ region.

In this paper we have not discussed how non-perturbative effects (see Refs.  and references therein) affect the transverse-momentum cross sections. Although the process-independent resummation formula (13) can be recast in the form of Eq. (4) at the perturbative level, it is not evident whether the two formulations are equally suitable to deal with non-perturbative contributions. In particular, since the form factor $`S_c^F`$ and the coefficient functions $`C_{ab}^F`$ in Eq. (4) are process-dependentWe remind the reader that the process dependence can kinematically be quite non-trivial when the transverse-momentum cross section in not fully integrated over the kinematical variables of the final-state system $`F`$., a simple and uniform implementation of process-independent non-perturbative effects (such as those due to the initial intrinsic $`k_t`$ of the partons in the colliding hadrons) in Eq. (4) can be very involved. The process-independent version (13) of the resummed cross section can help to consistently introduce and implement non-perturbative effects in transverse-momentum distributions.

We conclude by discussing the dependence of the resummation formulae on the renormalization and factorization scales $`\mu _R`$ and $`\mu _F`$. This dependence is often parametrized by introducing some arbitrary coefficients, $`C_1,C_2,C_3,`$ as suggested in Refs. . This is a perfectly sensible and reasonable procedure to try to estimate the effect of higher-order corrections, but it does not exactly correspond to the procedure that is usually followed in fixed-order perturbative calculations. In the case of soft-gluon resummed calculations for event shapes and jet rates in $`e^+e^{}`$ annihilation and for threshold contributions in hadronic cross sections, a different procedure was introduced in Refs.  and . The latter procedure closely matches the use of $`\mu _R`$ and $`\mu _F`$ in fixed-order calculations, and it can directly be applied (see, e.g., Refs.  and several other references quoted in Section 5 of Ref. ) to compare the scale dependence of resummed and fixed-order predictions. Here, we would like to note that the procedure of Refs.  can be introduced also in the case of the transverse-momentum distributions§§§A similar observation can be found in Ref. ..

To this purpose, we first perform the replacement $`f_{a/h}(x,b_0^2/b^2)f_{a/h}(x,\mu _F^2)`$ in Eq. (3) and we rewrite the resummed component of the transverse-momentum cross section as

$`\left[{\displaystyle \frac{d\sigma _F}{dQ^2dq_{}^2}}\right]_{\mathrm{res}.}`$ $`=`$ $`{\displaystyle \underset{a,b}{}}{\displaystyle _0^1}𝑑x_1{\displaystyle _0^1}𝑑x_2f_{a/h_1}(x_1,\mu _F^2)f_{b/h_2}(x_2,\mu _F^2)`$ (37)
$``$ $`{\displaystyle _0^{\mathrm{}}}𝑑b{\displaystyle \frac{b}{2}}J_0(bq_{})W_{ab}^F(x_1x_2s;Q,b,\mu _F).`$
To obtain Eq. (37) we have simply used the scale dependence of the parton densities as given by

$$f_{a/h}(x,b_0^2/b^2)=\underset{b}{}_x^1\frac{dz}{z}U_{ab}(z;b_0^2/b^2,\mu _F^2)f_{b/h}(x/z,\mu _F^2),$$
(38)
where $`U_{ab}(z;b_0^2/b^2,\mu _F^2)`$ is the customary evolution operator matrix obtained by solving the Altarelli-Parisi evolution equations to the required logarithmic accuracy. The resummed partonic cross section in Eq. (37) is thus obtained from that in Eq. (13) (or (4)) by replacing the coefficient functions $`C_{ca}(\alpha _\mathrm{S}(b_0^2/b^2),x)`$ (or $`C_{ca}^F(\alpha _\mathrm{S}(b_0^2/b^2),x)`$) with the following convolution

$$C_{ca}(\alpha _\mathrm{S}(b_0^2/b^2),x;\mu _F^2)=\underset{b}{}_x^1\frac{dz}{z}C_{cb}(\alpha _\mathrm{S}(b_0^2/b^2),z)U_{ba}(z;b_0^2/b^2,\mu _F^2).$$
(39)
Note that, unlike the parton densities in Eq. (3), those in Eq. (37) do not depend on the impact parameter $`b`$, but only on the factorization scale $`\mu _F`$. Expression (37) can thus be useful in practical calculations, because it avoids the numerical integration of the parton densities with respect to the impact parameter. The only factor of the resummed cross section that has to be integrated over $`b`$ is now the perturbative component $`W_{ab}^F(x_1x_2s;Q,b,\mu _F)`$.

Then, to proceed as in Refs. , we simply observe that we can choose $`\mu _F`$ and $`\mu _R`$ to be of the order of the hard scale $`Q`$. Thus, $`W_{ab}^F(x_1x_2s;Q,b,\mu _F)`$ contains only an additional large parameter, $`L\mathrm{ln}Q^2b^2`$, and it can be written as

$`W^F(Q,b,\mu _F)`$ $``$ $`\mathrm{exp}\left\{Lf_1\left(\alpha _\mathrm{S}L{\displaystyle \frac{}{}}\right)+f_2(\alpha _\mathrm{S}L;{\displaystyle \frac{Q}{\mu _R}},{\displaystyle \frac{Q}{\mu _F}})+\alpha _\mathrm{S}f_3(\alpha _\mathrm{S}L;{\displaystyle \frac{Q}{\mu _R}},{\displaystyle \frac{Q}{\mu _F}})+\mathrm{}\right\}`$ (40)
$``$ $`\sigma ^F(Q^2,\alpha _\mathrm{S};Q^2/\mu _R^2),`$
where $`\alpha _\mathrm{S}=\alpha _\mathrm{S}(\mu _R^2)`$. The notation on the right-hand side is symbolical, because we have understood the dependence on the momentum fractions The convolutions with respect to the momentum fractions can as usual be diagonalized by taking $`N`$-moments and working in Mellin space. $`z_1,z_2`$ and on the parton indices $`a,b=q,\overline{q},g`$. In particular, the exponential should be understood as an exponential matrix that depends on the parton indices.

The form in Eq. (40) is straightforwardly obtained by using the exponentiated result in Eq. (5) for the form factor and the customary exponential form for the evolution operator $`U_{ab}(z;b_0^2/b^2,\mu _F^2)`$. In the exponent, the function $`Lf_1`$ resums all the LL contributions $`\alpha _\mathrm{S}^nL^{n+1}`$, the function $`f_2`$ resums the NLL terms $`\alpha _\mathrm{S}^nL^n`$, while $`\alpha _\mathrm{S}f_3`$ contains the NNLL terms $`\alpha _\mathrm{S}^nL^{n1}`$, and so forth. The functions $`f_i`$ are normalized as $`f_i(\alpha _\mathrm{S}L=0)=0`$ and can easily be obtained in analytic form.

Although it is symbolical, the form of Eq. (40) is sufficient to make our point. Indeed, it is exactly in the same form as the expressions that are obtained by performing soft-gluon resummation for $`e^+e^{}`$ event shapes and hadronic cross sections near threshold. Recasting the resummed transverse-momentum distributions in the form of Eq. (40), the theoretical accuracy of the resummed calculation can be investigated as in fixed-order calculations, by varying $`\mu _F`$ and $`\mu _R`$ around the value $`Q`$ of the typical hard scale. Future studies along these lines (for instance, comparisons of the scale dependence of resummed and fixed-order calculations) could be useful to increase our confidence in the quantitative reliability of resummed predictions for transverse-momentum cross sections.


format: text

text:
# Dirac-Brueckner Approach to Hyperon Interactions and Hypernuclei

## I Introduction

Hypernuclear studies are the natural extension of isospin dynamics in non-strange nuclei towards a more general theory of flavor dynamics in a baryonic environment. Modern hypernuclear theories are using non-relativistic and relativistic microscopic descriptions . Since from a QCD point of view hypernuclei as also isospin nuclei are deep in the non-perturbative low energy-momentum regime, such a description in terms of mesons and baryons should be adequate. Relativistic mean-field (RMF) theories of Walecka-type  have been applied successfully  with empirically adjusted meson-hyperon vertices. SU(3)f-symmetric field theories incorporating chirality  or accounting for the quark structure of hadrons  have been formulated and applied to hypernuclei. In a SU(3)f approach nucleon-hyperon and hyperon-hyperon interactions in free space  and in a nuclear environment  have been calculated.

In this contribution, hypernuclei are described in the Density Dependent Relativistic Hadron (DDRH) theory. Introduced originally as an effective field theory for isospin nuclei  it was extended to hypernuclei recently . In DDRH theory the in-medium modifications of meson-baryon vertices are incorporated by functionals of the fermion field operators. The functional dependence of the vertices on density is derived from infinite matter Dirac-Brueckner Hartree-Fock (DBHF) calculations. For isospin nuclei, a practically parameter-free model Lagrangian is obtained once a free space interaction is chosen. Lorentz-invariance, thermodynamical consistency and covariance of the field equations are retained.

The extension of DDRH theory to strange baryons is reviewed in sect.II. As the central theoretical result nucleon and hyperon vertices are found to be related by scaling laws. A reduced model, appropriate for relativistic Hartree calculations of single $`\mathrm{\Lambda }`$ nuclei, is introduced. A semi-microscopic derivation of the N$`\mathrm{\Lambda }`$ interaction is discussed where the $`\sigma `$ coupling is taken from a theoretical N$`\mathrm{\Lambda }`$ T-matrix  while the $`\omega `$ coupling is determined empirically. DDRH mean-field results for hypernuclei are presented in sect.III and compared to RMF calculations. On a global scale spectroscopic data are described satisfactorily well, including the reduced spin-orbit splitting in $`\mathrm{\Lambda }`$ nuclei. The paper closes with a summary, conclusions and an outlook to work in progress in sect.IV.

## II Density dependent hadron field theory with hyperons

### A The DDRH model Lagrangian

In hypernuclear models derived from a symmetry-broken SU(3)f Lagrangian  neither of the $`0^{}`$ $`\pi `$ and K meson fields contributes directly to a structure calculation, except in the u-channel through antisymmetrization. From the $`1^{}`$ vector meson octet condensed isoscalar $`\omega `$ and isovector $`\rho `$ meson fields will evolve. In a system with a large fraction of hyperons also condensed octet $`K^{}`$ and singlet $`\mathrm{\Phi }`$ mesons fields can appear. A shortcoming of a pure SU(3)f approach is the missing of $`0^+`$ scalar mesons and, hence, the absence of a binding mean-field. A satisfactory description of the $`0^+`$ meson channels, e.g. in terms of dynamical two-meson correlations , is an unsolved question.

In view of these problems we use an effective Lagrangian including the degrees of freedom which are relevant for the nuclear structure problem. This is achieved by extending the original DDRH proton-neutron Lagrangian to the $`1/2^+`$ S=-1 ($`\mathrm{\Lambda },\mathrm{\Sigma }^{\pm ,0}`$) and S=-2 ($`\mathrm{\Xi }^{,0}`$) hyperon multiplets . In the meson sector the isoscalar $`\sigma `$, $`\sigma _s`$ ($``$ scalar $`s\overline{s}`$ condensate), $`\omega `$ and $`\varphi `$ meson, the isovector $`\rho `$ meson and the photon $`\gamma `$ are included. This leads to the isospin-symmetric Lagrangian

$``$ $`=`$ $`_B+_M+_{int}`$ (1)
$`_B`$ $`=`$ $`\overline{\mathrm{\Psi }}_F\left[i\gamma _\mu ^\mu \widehat{M}\right]\mathrm{\Psi }_F`$ (2)
$`_M`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=\sigma ,\sigma _s}{}}\left(_\mu \mathrm{\Phi }_i^\mu \mathrm{\Phi }_im_{\mathrm{\Phi }_i}^2\mathrm{\Phi }_i^2\right){\displaystyle \frac{1}{2}}{\displaystyle \underset{\kappa =\omega ,\varphi ,\rho ,\gamma }{}}\left({\displaystyle \frac{1}{2}}F^{(\kappa )^2}m_\kappa ^2A^{(\kappa )^2}\right)`$ (3)
$`_{int}`$ $`=`$ $`\overline{\mathrm{\Psi }}_F\widehat{\mathrm{\Gamma }}_\sigma \mathrm{\Psi }_F\mathrm{\Phi }_\sigma \overline{\mathrm{\Psi }}_F\widehat{\mathrm{\Gamma }}_\omega \gamma _\mu \mathrm{\Psi }_FA_\omega ^\mu {\displaystyle \frac{1}{2}}\overline{\mathrm{\Psi }}_F\stackrel{}{\tau }\widehat{\mathrm{\Gamma }}_\rho \gamma _\mu \mathrm{\Psi }_F\stackrel{}{A}_\rho ^\mu `$ (5)
$`+\overline{\mathrm{\Psi }}_F\widehat{\mathrm{\Gamma }}_{\sigma _s}\mathrm{\Psi }_F\mathrm{\Phi }_{\sigma _s}\overline{\mathrm{\Psi }}_F\widehat{\mathrm{\Gamma }}_\varphi \gamma _\mu \mathrm{\Psi }_FA_\varphi ^\mu e\overline{\mathrm{\Psi }}_F\widehat{Q}\gamma _\mu \mathrm{\Psi }_FA_\gamma ^\mu ,.`$
Here, $`_B`$ and $`_M`$ are the free baryonic and mesonic Lagrangians, respectively. Baryons are described by the flavor spinor $`\mathrm{\Psi }_F`$

$$\mathrm{\Psi }_F=(\mathrm{\Psi }_N,\mathrm{\Psi }_\mathrm{\Lambda },\mathrm{\Psi }_\mathrm{\Sigma },\mathrm{\Psi }_\mathrm{\Xi })^T.$$
(6)
The diagonal matrix $`\widehat{M}`$ contains the free-space baryon masses and $`\widehat{Q}`$ is the electric charge operator. The meson-baryon interactions are contained in $`_{int}`$. The usual field strength tensor of either the vector mesons ($`\kappa =\omega ,\varphi ,\rho `$) or the photon ($`\kappa =\gamma `$) is denoted by $`F_{\mu \nu }^{(\kappa )}`$. Contractions of the field strength tensors are abbreviated as $`F^2=F_{\mu \nu }F^{\mu \nu }`$ etc..

An important difference of the DDRH Lagrangian, eq.(1), to standard RMF approaches  is the description of medium effects and non-linearities in terms of density dependent meson-baryon vertex functionals $`\widehat{\mathrm{\Gamma }}_{\alpha B}=\widehat{\mathrm{\Gamma }}_{\alpha B}(\overline{\mathrm{\Psi }}_F\mathrm{\Psi }_F)`$. The vertices are taken as Lorentz invariant functionals of the fermion field operators $`\mathrm{\Psi }_F`$ and since the baryon fields are treated as quantum fields, even in the mean-field limit a well defined class of quantum fluctuations with non-vanishing ground state expectation values is taken into account . Dynamically, the vertices contribute to the Dirac equations as rearrangement self-energies describing the static polarization of the medium . In bulk quantities, as for example total binding energies, the DDRH rearrangement self-energies are cancelled exactly but contribute to single particle quantities like separation energies, wave functions and density matrices .

A solvable model is obtained in the Hartree mean-field approximation. The operator-valued vertex functionals become c-number functions of the baryon densities . The meson fields are obeying classical field equations, while the baryons are treated as quantum fields by solving the Dirac equation with static but density dependent self-energies, including rearrangement contributions . For single $`\mathrm{\Lambda }`$ nuclei the hidden-strangeness $`\sigma _s`$ and $`\varphi `$ fields are neglected because they are of order $`𝒪`$(1/A)

### B Dirac-Brueckner Approach to In-medium Hyperon Interactions

It is obvious that the structure of the vertex functionals must be derived in a separate step. Provided that data are available the required information could be obtained from phenomenology, as in the approach of ref.  for isospin nuclei. A derivation on theoretical grounds has the advantage of providing \- at least in principle - a deeper insight into the structure of interactions and, especially, the origin of medium-dependencies and inter-relations between the flavour sectors. Hence, we follow the original DDRH approach  and derive the vertex functionals from Dirac-Brueckner calculations. Full scale SU(3)f DB calculations, however, are a considerable task whose solution is still pending. But the success of non-relativistic Brueckner calculations  clearly indicate the promising potential of such a microscopic approach.

The extension of Dirac-Brueckner calculations to the SU(3)f multiplet have been outlined in . Aiming at applications of the DB-results in relativistic mean-field calculations it is sufficient to have accurate knowledge on baryon self-energies rather than on the full momentum structure of in-medium interactions. A central result of  is the observation that Dirac-Brueckner interactions can be expressed in terms of medium-renormalized meson-exchange interactions with density dependent vertex factors,

$$\mathrm{\Gamma }_{\alpha B}(k_F^B)g_\alpha ^Bs_\alpha ^B(k_F^B)$$
(7)
given by the bare coupling constants $`g_\alpha ^B`$ and the density dependent renormalization factor $`s_\alpha ^B(k_F^B)`$ for baryons of type B interacting through the exchange of mesons $`\alpha `$.

With our choice of momentum independent, global vertices, the baryon self-energies are obtained as

$$\mathrm{\Sigma }_\alpha ^B(k|k_F^N,k_F^Y)=\mathrm{\Gamma }_{\alpha B}(k_F^B)\varphi _\alpha (k|k_F^N,k_F^Y,\mathrm{\Gamma })$$
(8)
where $`\mathrm{\Phi }_\alpha `$ denotes a condensed meson field. From this equation the link to the DDRH Lagrangian and RMF theory is evident: Evaluating eq.(1) in DHF approximation self-energies of the same structure are obtained.

A self-contained model is obtained by introducing an in-medium ”renormalization” scheme. For that purpose it is of advantage to consider symmetric hyper-matter, i.e. $`k_F^N=k_F^Y=k_F`$ where one finds the relation

$$\mathrm{\Gamma }_{\alpha Y}(k_F^Y)=\mathrm{\Gamma }_{\alpha N}(k_F^Y)\frac{\mathrm{\Sigma }_\alpha ^Y(k|k_F)}{\mathrm{\Sigma }_\alpha ^N(k|k_F)}|_{k=k_F,k_F^N=k_F^Y}.$$
(9)
which is exact in Hartree approximation. The significance of this result becomes apparent by expanding DB self-energies diagramatically with respect to the bare coupling constants $`g_\alpha ^B`$ . Because the leading order contributions are given by tadpole diagrams one derives easily that the nucleon and hyperon self-energies are related by scaling laws

$$R_\alpha ^Y=\frac{\mathrm{\Sigma }_\alpha ^Y}{\mathrm{\Sigma }_\alpha ^N}\frac{g_\alpha ^Y}{g_\alpha ^N}(1+𝒪(1\frac{M_N}{M_Y}))+\mathrm{}$$
(10)
where the realistic case $`g_\alpha ^Y<g_\alpha ^N`$ is considered. The scaling factors $`R_\alpha ^Y`$ are expected to be state-independent, universal constants whose values are close to the ratios of the bare coupling constants. For asymmetric hypermatter with a hyperon fraction $`\zeta _Y=\frac{\rho ^Y}{\rho ^N}1`$ a corresponding diagrammatic analysis shows that asymmetry terms are in fact suppressed because the asymmetry correction is of leading second order $`𝒪((\frac{g^Y}{g^N}\zeta _Y)^2)`$. Thus, even in a finite nucleus where $`\zeta _Y`$ may vary over the nuclear volume, we expect $`R_\alpha ^Y=const.`$ to a very good approximation.

These results agree remarkably well with the conclusions drawn from the analyis of single hypernuclei in purely phenomenological models. In the present context, eqs.(9) and (10) are of particular interest because they allow to extend the DDRH approach in a theoretically meaningful way to hypernuclei using the results available already from investigations of systems without strangeness. Below, the nucleon (Hartree) scalar and vector vertex functions $`\mathrm{\Gamma }_{\sigma ,\omega N}(k_F)`$ of  will be used as reference values. The hyperon scaling factors $`R_\alpha ^Y`$ are treated as phenomenological constants to be determined empirically. The $`\chi ^2`$ distribution from fitting DDRH RMF-calculations to existing hypernuclear data deduced from $`(\pi ,K)`$  experiments is shown in Fig. 1. The joined distribution of the scalar and vector scaling factors is characterized by a steep valley along the diagonal without a clearly pronounced minimum.

In order to stay as close as possible to the microscopic DDRH picture we use the theoretically derived value Rσ = 0.490 obtained from free space N$`\mathrm{\Lambda }`$ scattering with the Jülich potential . Since theoretical values for the $`\omega `$ vertex are not available $`R_\omega `$ is taken from Fig. 1 leading to Rω=0.553 with the above value of Rσ . The scaling factors are in surprisingly good agreement with RMF results  and are consistent with bounds on hyperon–nucleon couplings extracted from neutron star models .

## III Spectroscopy of single $`\mathrm{\Lambda }`$ hypernuclei

Relativistic DDRH Hartree theory and applications to isospin nuclei were discussed in great detail in ref.  and the references therein. Here, we present DDRH results only for single $`\mathrm{\Lambda }`$ hypernuclei. The vertices are taken from DB calculations with the Bonn A NN potential . The model parameters are compiled in Tab. I.

### A $`\mathrm{\Lambda }`$ Single particle states and spin-orbit splitting

Hyperon single particle spectra for $`S`$=-1 hypernuclei can be seen as a very clean fingerprint of mean-field dynamics, since they are only weakly affected by many-body effects. The bulk structure contains information on the mean-field and, indirectly, the nucleonic density distributions. Of particular theoretical interest are spin-orbit splittings. At high energy resolution the fine structure of the spectra provides information also on dynamical correlations beyond static mean-field dynamics.

DDRH $`\mathrm{\Lambda }`$ single particle levels for light to heavy nuclei are shown in Fig. 2. Two major differences between neutron and $`\mathrm{\Lambda }`$ spectra are observed :

1. $`\mathrm{\Lambda }`$ and neutron single particle spectra are overall related by a constant shift and an additional compression because the $`\mathrm{\Lambda }`$ central potential has a depth of only about -30 MeV, compared to -70 MeV for the neutrons.
2. The spin-orbit splitting of the $`\mathrm{\Lambda }`$ states is reduced further being less than what is expected from the overall reduction of the potential strength.

The DDRH calculations reproduce the experimentally observed very small spin-orbit splitting in $`\mathrm{\Lambda }`$ hypernuclei, e.g. , reasonably well even without an explicit dynamical suppression of spin-orbit interactions as e.g. a $`\mathrm{\Lambda }\omega `$ tensor coupling which, for example, is used in the QMC model . $`\mathrm{\Lambda }`$ and neutron spin-orbit potentials are compared in Fig. 3.

The reason for the small spin-orbit splitting is understood by considering the evolution with increasing mass number. From Fig. 4 it is seen that the splitting drops for higher masses. Such a behaviour is to be expected since the spin-orbit potential is a finite size effect and will vanish in the nuclear matter limit. A corresponding mass dependence is also found in pure isospin nuclei. The splitting also drops in the low mass region – now for the reason that the spin-orbit doublets approach the continuum threshold and get compressed before one of them or both become unbound. There is a remarkable similarity to the situation found in weakly bound neutron-rich exotic nuclei . In both cases weak binding is an important reason for the reduction of spin-orbit effects. The wave functions of weakly bound states \- either for the special dripline neutron states or $`\mathrm{\Lambda }`$ orbits in general - have a reduced overlap with the spin-orbit potential which remains well localized in the nuclear surface. As seen from Tab. II the binding energy effect leads to substantially larger spatial extensions of the $`\mathrm{\Lambda }`$ states.

### B Comparison to data and phenomenological RMF calculations

In Fig. 5 the DDRH single particle spectra are compared to spectroscopic data from $`(\pi ^+,K^+)`$ reactions. States in intermediate to high mass nuclei are described fairly well by the model, while for masses below about $`{}_{}{}^{28}Si_\mathrm{\Lambda }`$ deviations of up to 2.5 MeV arise, possibly indicating the limits of an approach using nuclear matter vertices in local density approximation. We seem to miss sytematically a decrease of the repulsive vector interaction for low mass single $`\mathrm{\Lambda }`$ hypernuclei. This tendency already becomes apparent going from $`{}_{}{}^{51}V_{\mathrm{\Lambda }}^{}`$ to $`{}_{}{}^{28}Si_\mathrm{\Lambda }`$. Actually, the dashed line for $`R_\omega =0.542`$ in Fig. 5 shows that a slight reduction of the vector repulsion improves considerably the description of the light mass data. This 2% variation of $`R_\omega `$ is largely within the uncertainties of the model parameters.

Especially for the heavier nuclei the microscopic DDRH results are at least of comparable quality as the phenomenological descriptions. This we consider as a remarkable success of the model because the coupling functionals were not especially adjusted to data except for the overall adjustment of the vector scaling factor $`R_\omega `$. In heavy nuclei the density dependence of the $`\mathrm{\Lambda }`$ vertices act mainly as an overall reduction factor because the $`\mathrm{\Lambda }`$ vector density is relatively small and varies only weakly across the nuclear volume. The essential density dependent effect is related to the dynamical redistribution of the surrounding nucleons by the rearrangement self-energies. Dynamically, it corresponds to a modification of the $`\mathrm{\Lambda }`$ core potential due to static polarization in the nucleonic sector. Obviously, this effect is not accounted for by conventional RMF models.

Finally, DDRH and RMF results for $`\mathrm{\Lambda }`$ particle-neutron hole configurations observed in ($`K^{},\pi ^{}`$) reactions  are given in Tab. III.

## IV Summary, conclusions and outlook

The DDRH theory introduced previously for isospin nuclei was extended to hypernuclei by including the full set of SU(3)f octet baryons. Interactions were described by a model Lagrangian including strangeness-neutral scalar and vector meson fields of $`q\overline{q}`$ (q=u,d) and $`s\overline{s}`$ quark character. The medium dependence of interactions was described by meson-baryon vertices chosen as functionals of the baryon field operators. The DDRH vertices are chosen to cancel Dirac-Brueckner ground state correlations. Hence, the approach corresponds to a resummation of ladder diagrams into the vertices under the constraint that infinite matter ground state self-energies and total binding energies are reproduced. As the central theoretical result it was found that the structure of Dirac-Brueckner interactions strongly indicates that the ratio of nucleon and hyperon in-medium vertices should be determined already by the ratio of the corresponding free space coupling constants being affected only weakly by the background medium.

Dynamical scaling of nucleon-hyperon vertices was tested in DDRH mean-field calculations for single $`\mathrm{\Lambda }`$ hypernuclei. Calculations over the full range of known single $`\mathrm{\Lambda }`$ nuclei led to a very satisfactory description of $`\mathrm{\Lambda }`$ separation energies. The deviations from the overall agreement for masses below A$``$16 are probably related to the enhancement of surface effects in light nuclei which are not described properly by static RMF calculations with DB vertices obtained in the local density approximation. In a recent non-relativistic calculation indeed sizable contributions of hyperon polarization self-energies especially in light nuclei  were found.

The results are encouraging and we conclude that DDRH theory is in fact an appropriate basis for a microscopic treatment of hypernuclei. The present formulation and applications are first steps on the way to a more general theory of in-medium SU(3)f flavor dynamics. Future progress on dynamical scaling and other theoretical aspects of the approach is depending on the availability of Dirac-Brueckner calculations for the full baryon octet including also the complete pseudoscalar $`0^{}`$ and vector $`1^+`$ meson multiplets.

As work in progress the production of hypernuclei in hadronic reactions is presently investigated. Initial and final state interactions of the incident and outgoing mesons in ($`\pi ^+,K^+`$) and ($`K^{},\pi ^{}`$) reactions are described in a relativistic eikonal approach. Good agreement with differential and total cross section data is obtained . A Lagrangian model, including s-channel production through nucleon resonances and t-channel production by e.g. $`K^{}`$ exchange is under investigation. The nuclear structure results discussed here are entering into the (coherent) production amplitude. Electro-production of hypernuclei will be used as an independent and important test for the production vertex and the nuclear structure input.

## Acknowledgements

This work was supported in part by DFG under contract Le439/4-3, Graduiertenkolleg Theoretische und Experimentelle Schwerionenphysik, GSI and BMBF. Discussions with C. Greiner are gratefully acknowledged.


text:
# Gas Flow and Star Formation in the ‘Antennae’ Galaxies

## References

Barnes J.E., 1988, ApJ, 331, 699

Jog C.J., Solomon P.M., 1992, ApJ, 387, 152

Rigopoulou D., et al., 1999, AJ, 118, 2625

Toomre, A., & Toomre, J. 1972, ApJ, 405, 142

Whitmore B.C., Zhang Q., Leitherer C., Fall S.M., 1999, 118, 1551

Mengel S., et al., 2000, ApJ, submitted

Mihos J.C., Hernquist L., 1996, ApJ, 464, 641

Scoville N.Z., et al., 2000, AJ, 119, 991

Soifer B.T., et al., 1987, ApJ, 320, 238


text:
# Solar Oscillations and Convection: II. Excitation of Radial Oscillations

## 1 Introduction

Two ideas for the source of $`p`$-mode excitation have been pursued – overstability (as in pulsating stars) and turbulent Reynolds stresses (as in jet noise) (Biermann, 1948; Schwarzschild, 1948; Lighthill, 1952; Stein, 1967, 1968; Crighton, 1975; Ando & Osaki, 1977; Goldreich & Keeley, 1977; Goldreich & Kumar, 1990; Balmforth, 1992; Musielak et al., 1994). We show here that it is the $`PdV`$ work of stochastic, non-adiabatic pressure fluctuations which is the primary mode excitation mechanism (Stein & Nordlund, 1991; Bogdan et al., 1993; Goldreich et al., 1994). Near the surface, the non-adiabatic gas pressure (i.e. entropy) fluctuations dominate. They are produced by radiative cooling at the solar surface which is not locally and instantaneously exactly balanced by convective heat deposition. Below the surface, non-adiabatic turbulent pressure (Reynolds stress) fluctuations dominate. They are produced by the turbulent convective motions.

In a previous paper (Nordlund & Stein, 2000) we derived an exact expression for the stochastic excitation rate of the radial solar $`p`$-modes by the $`PdV`$ work of non-adiabatic gas and turbulent pressure fluctuations on the mode compression. We now use realistic numerical simulations of near surface solar convection (depth about 2.5 Mm) to evaluate this expression (Stein & Nordlund, 1998). Because the largest entropy and pressure fluctuations occur near the solar surface and because modes with frequencies in the 3 – 4 mHz range, where the excitation rate is largest, are confined near the solar surface, these near surface simulations include the primary excitation and damping processes.

## 2 Mode Excitation: Formalism

The rate of energy input to the modes can be calculated starting with the kinetic energy equation for the modes (Nordlund & Stein, 2000). Neglecting the viscous stresses,

$`\rho {\displaystyle \frac{D}{Dt}}\left({\displaystyle \frac{1}{2}}u_z^2\right)`$ $`=`$ $`{\displaystyle \frac{}{z}}\left[u_z\left(P_g+P_t\right)\right]+\rho u_zg`$ (1)
$`+`$ $`\left(P_g+P_t\right){\displaystyle \frac{u_z}{z}}.`$
After integrating this equation over depth, the flux divergence term only gives contributions from the end points and is negligible. The buoyancy term is small because mass is conserved so there is no net mass flux. The last term is the $`PdV`$ work,

$$W=𝑑t𝑑z\delta P\frac{\dot{\xi }}{z}.$$
(2)
There are several contributions to this work. The displacement, $`\xi `$, has contributions from the modes as well as the random convective motions. The pressure, $`\delta P`$, has coherent contributions from the modes and incoherent contributions from the random convective motions. Both coherent and incoherent contributions can be further divided into adiabatic and non-adiabatic terms. The dominant driving comes from the interaction of the non-adiabatic, incoherent pressure fluctuations,

$$\delta P^{\mathrm{nad}}=\left(\delta \mathrm{ln}P\mathrm{\Gamma }_1\delta \mathrm{ln}\rho \right)P,$$
(3)
with the coherent mode displacement,

$$\frac{D\mathrm{ln}\rho }{Dt}=\frac{\dot{\xi }}{z}.$$
(4)
This is a stochastic process, so the pressure fluctuations occur with random phases with respect to the modes. Therefore one must average over all possible relative phases between them. The resulting rate of energy input to the modes is (Nordlund & Stein, 2000)

$$\frac{\mathrm{\Delta }E_\omega }{\mathrm{\Delta }t}=\frac{\omega ^2\left|_r𝑑r\delta P_\omega ^{}\frac{\xi _\omega }{z}\right|^2}{8\mathrm{\Delta }\nu E_\omega }.$$
(5)
Here, $`\delta P_\omega `$ is the time Fourier transform of the non-adiabatic total pressure. $`\mathrm{\Delta }\nu =`$ 1/(total time interval) is the frequency resolution with which $`\delta P_\omega `$ is computed. $`\xi _\omega `$ is the mode displacement eigenfunction, which is typically chosen to be real for an adiabatic mode. In that case, taking the complex conjugate of the pressure is not necessary, but we retain it for generality. The mode energy is

$$E_\omega =\frac{1}{2}\omega ^2_r𝑑r\rho \xi _\omega ^2\left(\frac{r}{R}\right)^2=M_\omega V_\omega ^2(R).$$
(6)
Here $`M_\omega `$ is the mode mass and $`V_\omega (R)`$ is the mode velocity amplitude at the surface. Eqn 5 is similar to the expression of Goldreich et al. (1994) (eqn. 26), but involves no approximations. Having the numerical simulation data, we can evaluate this expression exactly without having to make approximations in order to evaluate it analytically 111A detailed description of how eqn 5 is evaluated is given in the appendix. Its validity is tested by applying it to the three radial modes of the simulation domain (sec. 4).

## 3 Simulations

We simulate a small portion of the solar photosphere and the upper layers of the convection zone, a region extending $`6\times 6`$ Mm horizontally and from the temperature minimum at $`0.5`$ Mm down to $`2.59`$ Mm below the visible surface. We solve the equations of mass, momentum and energy conservation in the form:

$`{\displaystyle \frac{\mathrm{ln}\rho }{t}}`$ $`=`$ $`𝐮\mathrm{ln}\rho 𝐮,`$ (7)
$`{\displaystyle \frac{𝐮}{t}}`$ $`=`$ $`𝐮𝐮+𝐠{\displaystyle \frac{P}{\rho }}\mathrm{ln}P+{\displaystyle \frac{1}{\rho }}\underset{¯}{\underset{¯}{\sigma }},`$ (8)
$`{\displaystyle \frac{e}{t}}`$ $`=`$ $`𝐮e{\displaystyle \frac{P}{\rho }}𝐮`$ (9)
$`+Q_{rad}+Q_{visc},`$
where $`\underset{¯}{\underset{¯}{\sigma }}`$ is the viscous stress tensor, $`Q_{rad}`$ is the radiative heating and $`Q_{visc}`$ is the viscous dissipation.

We use a non-staggered grid of either $`125^2`$ cells horizontally $`\times 82`$ cells vertically or $`63^3`$ cells. The spatial derivatives are calculated using third order splines and the time advance is a third order leapfrog scheme (Hyman, 1979; Nordlund & Stein, 1990). The code is stabilized by a hyper-viscosity which removes short wavelength noise without damping the longer wavelengths.

A large fraction of the internal energy is in the form of ionization energy near the solar surface, so we use a realistic equation of state (including the effects of ionization and excitation of hydrogen and other abundant elements and the formation and ionization of $`H_2`$ molecules). The pressure is found by interpolation in a table of $`P(\mathrm{ln}\rho ,e)`$ and its derivatives which is calculated with the Uppsala stellar atmosphere package (Gustafsson et al., 1975).

Radiative energy exchange is critical in determining the structure of the upper convection zone. Near the surface of the sun, the energy flow changes from almost exclusively convective below the surface to radiative above the surface. The interaction between convection and radiation near the solar surface determines what we observe and produces the entropy fluctuations that lead to the buoyancy work which drives the convection and that give rise to the non-adiabatic pressure fluctuations which excite the p-mode oscillations. Hence, the interaction between convection and radiation is crucial for both the diagnostics and the dynamics of convection. Since the top of the convection zone occurs near the level where the continuum optical depth is one, neither the optically thin nor the diffusion approximations give reasonable results. We therefore include 3D, LTE, non-gray radiation transfer in our model.

We simulate only a small region near the solar surface and must therefore impose boundary conditions inside the convectively unstable region. What happens outside our computational domain in principle influences the convective motions inside. However, at the top boundary, the density is very low relative to the rest of the volume and hence whatever happens there has very little influence on the convective motions. At the bottom, the incoming fluid is to a very good approximation isentropic and featureless, and hence carries little information. The unknown influence of external regions should therefore be small. This assertion is indeed confirmed by experiments with boundaries located at different depths.

The horizontal directions are taken to be periodic. In the vertical direction, we have a transmitting boundary at the temperature minimum (Nordlund & Stein, 1990). This is achieved by a larger than normal zone at the top boundary. Across this zone we make the vertical derivative of the density hydrostatic, set the vertical derivative of the velocity zero and hold the internal energy at the top fiducial layer constant in time and space. At the bottom of the computational domain, outgoing fluid goes out with whatever properties it has. For incoming fluid, we adjust the pressure such that the net mass flux through the bottom boundary vanishes. (This ensures that there is no boundary work done on vertical oscillation modes.) The pressure on the bottom boundary thus varies in time but is uniform over the horizontal plane. We damp fluctuations of the horizontal and vertical velocity of the incoming fluid, using a long time constant. Finally, we adjust the density and energy of the incoming fluid, at constant pressure, to fix its entropy (in both space and time).

The ability to do a direct numerical simulation with the wide range of length scales matching the dimensionless parameters – Reynolds number, Rayleigh number and Prandtl number – of the solar convection zone is beyond the speed and memory capabilities of current computers. Thus, our simulations are of the type called “Large Eddy Simulations”. It is, however, possible to resolve the surface thermal boundary layer of the convection zone and this we have done. Indeed, this is required to achieve results that agree quantitatively with solar observations (Stein & Nordlund, 2000).

The picture of convection that has emerged from the simulations is the following: Convection is driven by radiative cooling in the thin thermal boundary layer at the solar surface. It consists of cool, low entropy, filamentary, turbulent, downdrafts that plunge through a warm, entropy neutral, smooth, diverging, laminar upflow. Upflows must diverge as they ascend into lower density layers in order to conserve mass. This divergence smooths out any variations in their properties that might arise. Only a small fraction of the fluid at depth reaches the surface to be cooled and form the cores of the downdrafts. Most fluid turns over within a scale height and is entrained by the downdrafts. These low entropy downflows are the site of most of the buoyancy work that drives the convection. (See Stein & Nordlund (1998) for more details.)

We have made numerous comparisons between the predictions of the simulations and solar observations:

* The depth of the convection zone depends on the opacities that determine the temperature versus pressure (and hence entropy) stratification of the surface layers, the spectral line blocking, the convective efficiency of the superadiabatic layers immediately below the surface that determines the transition to the asymptotic adiabat, and the equation of state that determines the further run of temperature versus pressure through the convection zone. Excellent agreement is obtained between the depth of the convection zone predicted by our numerical simulations and that inferred from helioseismology (Rosenthal et al., 1999).
* The cavity for high frequency modes is enlarged by turbulent pressure support and 3D non-linear opacity effects which increase the average temperature required to produce a given effective temperature in an inhomogeneous compared to a homogeneous atmosphere. The $`p`$-mode eigenfrequencies calculated from the mean simulation atmosphere are significantly closer to the observed mode frequencies than those for standard spherically symmetric, mixing length models (Rosenthal et al., 1998, 1999).
* The simulation granulation size spectrum and the distribution of emergent intensities, when smoothed by the point spread function appropriate for the Swedish Vacuum Solar Telescope on La Palma (which produces the best solar images available today), closely match the observations (Stein & Nordlund, 1998).
* The width of photospheric iron lines (whose thermal speed is small) is a signature of the convective velocities. The net Doppler blue shift and asymmetry of spectral lines is a signature of the correlation between velocity and temperature fluctuations. Both these signatures agree closely with observations (Asplund et al., 2000).

This gives us confidence that the simulations are properly modeling the crucial properties of near surface solar convection.

## 4 Mode Excitation: Simulation

Three radial modes exist in our simulation and we first apply eqn (5) to these modes and compare the rate of work on the modes it predicts with the actual excitation rate $`dE_\omega /dt=\mathrm{\Gamma }E_\omega `$ determined from the mode’s energies and widths in the simulation. The modes can be clearly seen in the spectrum of the horizontally averaged, depth integrated kinetic energy (fig 1). Their properties are given in table 1. The lowest mode (with no zero crossings inside the computational domain) has a frequency of 2.57 mHz and a FWHM of 19 $`\mu `$Hz. Hence a simulation significantly longer than 14.5 hrs is required to resolve this mode. We use a simulation of 43 hours, with a spatial resolution of $`63^3`$. Snapshots were saved at 30 s intervals.

The eigenfunctions of the velocity are calculated by taking the time Fourier transform of the velocity. To get the real eigenmode, the transform of the velocity at the frequency of the mode is divided by its most common phase among all depths. To reduce the noise, we average the result over a frequency band, approximately equal to the FWHM of the mode, centered on the mode. The eigenfunctions are essentially the same as the solar model modes of Christensen-Dalsgaard using his spherically symmetric model S (Christensen-Dalsgaard et al., 1996) at the mode frequencies (fig 2), when the modes are normalized by the square root of the mode energy, or which is equivalent by their amplitude at the surface. Hence, we choose to use the solar model eigenfunctions, because they are much denser (35 radial modes below the acoustic cutoff frequency instead of three) and because they are slightly smoother.

Another way of looking at the modes is via their kinetic energy. Figure 3 is an image of the logarithm of the kinetic energy as a function of depth and frequency. The three modes of the computational domain are clearly seen. Notice also that there is a continuous pattern of nodes and antinodes, with the nodes descending in depth as the frequency increases and new nodes starting at the surface. This pattern extends even into the propagating region above the acoustic cutoff frequency of about 5.3 mHz.

The actual mode energy decay rate of the simulation modes, which on average is equal to their excitation rate, is obtained by multiplying the energy of each mode by its decay rate $`\mathrm{\Gamma }=2\pi \mathrm{\Delta }\nu _{\mathrm{FWHM}}`$, which is its radian line width and is obtained from the fit to the modes (fig 1).

The total mode energy is the sum of the mode energies minus a fit to the background over the frequency range where the modes are significantly above the background multiplied by the area (36 Mm2) of the simulation domain (fig 1).

Eqn (5) is used to predict the mode excitation rate. The non-adiabatic total (gas + turbulent) pressure fluctuations are calculated directly from the simulation by first averaging the gas pressure, turbulent pressure and density over horizontal planes at each saved snapshot. These are then interpolated to the Lagrangian frame at each time. The non-adiabatic pressure fluctuations are calculated as in eqns A1 and A2. The oscillation modes of the domain are essentially adiabatic and do not affect the non-adiabatic pressure fluctuations as can be seen from its featureless, noisy spectrum (Fig 4). Thus the non-adiabatic pressure work is due primarily to the convection.

The mode compression is calculated from the Christensen-Dalsgaard modes (since these are essentially identical with the simulation modes at the resonant frequencies) normalized by the square root of their energy in the simulation domain. Because the simulation domain is shallow, while, especially at low frequency, the solar modes have substantial amplitude below the computational domain, the mode mass in the computational box is significantly smaller than the actual mode mass of the solar modes (Fig 5).

Finally, we compare the actual mode energy decay rate with the predicted rate of work by convection on the modes given by eqn (5) in figure 6. The squares are the actual mode energy decay rates for the three resonant modes of the box. The solid line is the running mean of the predicted work and the dashed line is a smooth two power law fit to the predicted work. The agreement is very good, but not perfect. There are significant deviations from the power law fit in the neighborhood of the modes and these are reflected in the actual mode decay rate. This is a phenomena that still remains to be explained. Notice also that the decrease in work toward low frequencies is much slower than for the Sun (Fig 7). The reason is the near constancy of the mode mass within the simulation domain at these low frequencies compared with the steeply rising mode mass with decreasing frequency on the Sun. This application of our excitation rate formula to the modes that are excited in our simulation verifies that the formula is correct and can be applied to the Sun.

## 5 Mode Excitation: Sun

### 5.1 Excitation Spectrum

The excitation rate for solar modes as a function of frequency is shown in Fig. 7. To obtain these results we used the shorter but higher resolution $`125^2\times 82`$ simulation because it has more high frequency power. The magnitude and frequency dependence we find for our 6 Mm square box is very close to the observed values for the entire sun. This means that the pressure fluctuations must be uncorrelated on larger scales, so there is no extra driving contribution.

What produces this frequency dependence? The low frequency behavior is controlled by the nature of the eigenmodes and the high frequency behavior is controlled by the nature of convection. The work integral (Eqn. 5) contains a factor pertaining to the modes: the radial gradient of the mode displacement normalized by the square root of the mode energy (which makes it independent of the mode amplitude). This is small at low frequencies and increases with frequency approximately as $`\nu ^{3.5}`$ (Fig 8). Part of this dependence is due to the radial gradient of the displacement. The mode dispersion relation is $`\omega =\sqrt{(gk)}`$, so $`k=\omega ^2/g`$, which accounts for two powers of the frequency. The remainder of the frequency dependence is due to the mode mass, which decreases with increasing frequency because higher frequency modes are more concentrated toward the surface (Fig. 5).

The mode excitation decreases at high frequency because the pressure fluctuation power decreases with increasing frequency roughly as $`\nu ^4`$ (Fig. 9). Convective motions, whose power decreases at small scales and high frequencies, produce the gas and turbulent pressure fluctuations. This then leads to a similar high frequency decline in the pressure spectrum.

Total pressure fluctuations are small at low frequency because the atmosphere is nearly in hydrostatic equilibrium. However, the turbulent pressure fluctuations are largest at low frequency because they are a convective effect and convection is a longer time scale phenomena. Hence, the gas pressure fluctuations must also be large and out of phase with the turbulent pressure fluctuations at low frequency in order to produce small total pressure fluctuations (Fig. 10). The turbulent pressure is small compared to the gas pressure ($``$ 15%), but the magnitude of fluctuations are comparable for the gas and turbulent pressures, so they can indeed cancel each other.

### 5.2 Excitation Location

At what depth does the driving occur? Consider the integrand of the work integral at different frequencies (Figs. 11 and 12). At low frequencies the integrand amplitude is similar over an extended depth range, but it is small. Where the integrand is large, in the region of peak driving around 3-4 mHz, the integrand becomes concentrated very close to the surface and most driving occurs between the surface and 500 km depth. At still higher frequencies the integrand becomes small again and even more concentrated near the surface.

Where in space does this driving occur? The warm granules have only small non-adiabatic pressure fluctuations. Large, negative fluctuations are concentrated in the downdrafts (Figs. 13, 14).

The maximum mode driving occurs in the frequency range of 3 – 4 mHz, by non-adiabatic pressure fluctuations in the same frequency range. By filtering the time sequence of these fluctuations we see (Fig 15) that in the peak driving range also the driving occurs predominantly in the intergranule lanes and near the edges of granules. The high frequency power near granule edges is due to the motion of the granule boundaries over the one hour time interval on which the filtering was performed. This is, in part, a result of changes in granule size as they evolve. No direct correlation of non-adiabatic pressure fluctuations in the range of 2-5 mHz with velocity is seen. Keep in mind, however, that a correlation plot does not reveal correlations of events that happen in the neighborhood of one another.

### 5.3 Excitation Source

What is the source of the non-adiabatic pressure fluctuations? Is it entropy fluctuations or Reynolds stresses? Both play a role, but the primary source of mode driving is turbulent pressure fluctuations (Reynolds stresses) (Fig 16). This is surprising, since the non-adiabatic gas pressure power is larger than the turbulent pressure power near the surface (Fig 9). The non-adiabatic gas and turbulent pressures contribute comparably to the work near the surface, but the contribution of the turbulent pressure extends deeper and provides the dominant contribution to the total work (Fig 17).

The gas pressure fluctuations are significantly larger than the turbulent pressure and have maxima at the frequencies of the three radial modes of the simulation. The non-adiabatic pressure fluctuation power, however, varies smoothly across these frequencies indicating that it is primarily due to stochastic convective processes (Fig 9).

There is a tight correlation between the non-adiabatic pressure fluctuations at the surface and the emergent intensity (Fig 18), which indicates that it is the fluctuating cooling at the surface that is the main source of stochastic mode excitation there. Indeed, the source of entropy fluctuations is the cooling of fluid that approaches optical depth unity (Stein & Nordlund, 1998). This correlated noise is believed to be responsible for the difference in asymmetry of the modal power spectra observed in velocity and intensity (Nigam et al., 1998; Kumar & Basu, 1999). Our discussion in terms of non-adiabatic pressure fluctuations is equivalent to the discussion in terms of entropy fluctuations by Goldreich et al. (1994).

We can use the energy equation to determine the processes producing the non-adiabatic gas pressure fluctuations (Stein & Nordlund, 1991).

$$\delta P_{\mathrm{nad}}=\delta P\frac{\mathrm{\Gamma }_1P}{\rho }\delta \rho =\left(\frac{P}{s}\right)_\rho \delta s=(\mathrm{\Gamma }_31)\delta Q,$$
(10)
and

$$\frac{DQ}{Dt}=\frac{}{z}\left(F_{\mathrm{rad}}+F_{\mathrm{conv}}\right)+\stackrel{˘}{𝐮}P+Q_{\mathrm{diss}},$$
(11)
where $`\stackrel{˘}{𝐮}=𝐮\rho 𝐮/\rho `$ is the convective component of the velocity. The divergence of the radiative and convective fluxes are the dominant terms, but they nearly cancel each other since energy transport shifts from convective to radiative near the surface. It is their slight instantaneous imbalance locally that leads to the non-adiabatic pressure fluctuations (Figs 19 and 20).

The separate contributions of the gas and turbulent pressure to the total non-adiabatic pressure spectrum are shown in figs. 21 and 22, with the same color scale in both. Near the surface, there is clearly more power in the gas pressure. However, in the peak driving range of 3-4 mHz, there is more power in the non-adiabatic turbulent pressure below the surface.

Another reason for the dominance of turbulent pressure in the mode driving is revealed in figs. 23 and 24. The left panels in each show the non-adiabatic gas and turbulent pressures (scaled by $`\rho ^{1/2}`$) as a function of time and depth. The turbulent pressure varies more slowly with depth and has a longer time scale than the gas pressure. The middle panels show one realization of the mode density or compression (scaled by $`\rho ^{1/2}`$). The right panels show the work integrand, $`\delta P_{\mathrm{nad}\omega }\xi _\omega /zE_\omega ^{1/2}`$. The gas pressure contribution to the integrand is more concentrated at the surface and alternates sign with time and depth more rapidly than the turbulent pressure contribution. The greater coherence of the turbulent pressure allows it to make a greater contribution to the net work.

The very largest pressure perturbations are associated with the sudden initiation of downdrafts and produce waves that propagate up into the chromosphere and steepen into shocks (Skartlien, 1998; Skartlien et al., 2000). These events may correspond to the large individual acoustic events observed by Rimmele et al. (1995) and Goode et al. (1998). They may be the tail of the distribution of the stochastic driving process. A comparison of observed acoustic events and those in the simulation has been made by Goode et al. (1998).

## 6 Summary

How does this all fit together? Figs 2527 show the time variation of the horizontally averaged non-adiabatic pressure fluctuations, the density (compression) in the modes and the product of the two which is the integrand of the work integral (eqn 5) as a function of depth and time for modes with frequencies 1,3, and 5 mHz.

The dominant contribution to the work comes from the turbulent pressure (Reynolds Stress). Near the surface non-adiabatic gas pressure fluctuations, produced by an instantaneous imbalance between the divergences of the radiative and convective fluxes, also contribute. Their divergence is individually large, but since energy transport switches between convection and radiation in the surface layers, they are nearly equal and opposite. At each instant they do not exactly cancel, which leads to heating and cooling and hence entropy fluctuations. Excitation is small at low frequencies due to mode properties – the mode compression decreases and the mode mass increases at low frequency. In addition, at low frequency the mode amplitude is small near the surface where the non-adiabatic pressure fluctuations are large. (At very low frequencies, driving occurs deeper than the simulation domain.) Excitation is small at high frequencies due to the pressure fluctuation spectrum – pressure fluctuations become small at high frequencies because they are due to convective motions which have a longer time scale than the dominant $`p`$-mode periods. Large non-adiabatic pressure fluctuations occur primarily in the intergranular lanes (due to turbulence and surface radiative cooling) and near the edges of granules (due to granule expansion and collapse).

## Acknowledgments

We thank Dali Georgobiani for making the 43 hour simulation. This work was supported in part by NASA grants NAG 5-4031 and NAG 5-8053, NSF grants AST 9521785 and AST 9819799, and the Danish Research Foundation, through its establishment of the Theoretical Astrophysics Center. The calculations were performed at the National Center for Supercomputer Applications, which is supported by the National Science Foundation, at Michigan State University and at UNI$``$C, Denmark. This valuable support is greatly appreciated.

## Appendix A Evaluation of the Mode Excitation Rate (eqn 5)

The mode excitation rate (eqn 5) can be evaluated in two different ways: in Fourier space and in real space. We describe both methods. The variables that appear in eqn 5 for the mode excitation are the non-adiabatic total pressure fluctuation, the derivative of the mode displacement and the mode energy. These quantities are calculated as follows for both methods of evaluating the excitation.

The non-adiabatic pressure fluctuations are obtained from the simulation results. First the gas pressure, turbulent pressure ($`P_t=\rho u_z^2`$) and density are averaged over horizontal planes and saved at 10 s or 30 s intervals. These are then interpolated to the Lagrangian frame, that moves with the vertical (radial oscillation) motions. The Lagrangian frame is determined by calculating the mass column density at each time and interpolating the variables to the time average mass column density. Next, the non-adiabatic total (gas plus turbulent) pressure is calculated from

$$P^{\mathrm{nad}}=\left(\mathrm{ln}(P_{\mathrm{gas}}+P_{\mathrm{turb}})\mathrm{\Gamma }_1\mathrm{ln}\rho \right)(P_{\mathrm{gas}}+P_{\mathrm{turb}}).$$
(A1)
Finally, the fluctuation of the non-adiabatic total pressure about its time average,

$$\delta P^{\mathrm{nad}}(z,t)=P^{\mathrm{nad}}(z,t)P^{\mathrm{nad}}(z,t)_t,$$
(A2)
is determined.

The mode displacement $`\xi _\nu (z)`$ for the radial mode of frequency $`\nu `$ is obtained from the eigenmode calculations of Christensen-Dalsgaard, using his spherically symmetric model S (Christensen-Dalsgaard et al., 1996). Alternatively, one could extract modes directly from the numerical simulations, using Fourier decomposition. A drawback with this method is, apart from slightly noisier mode structure, that there are only three resonant modes in the shallow simulation box—the 35 radial modes from standard solar models provide much better frequency coverage. As we have shown (fig 2), the modes from the simulation are essentially the same as the solar modes at the resonant frequencies of the computational domain, where the solar modes have a node at the depth of the bottom of the computational domain.

The mode energy is calculated from the displacement according to eqn. 6. The mode displacements are interpolated to the simulation grid at the same total pressure as in the Christensen-Dalsgaard model and the derivative of the displacement is calculated. Since the derivative of the mode displacement always appears normalized by the square root of the mode energy, the mode amplitude normalization cancels.

In Fourier space, the work integral which appears in the energy input rate to the modes (eqn 5) is evaluated by taking the time Fourier transform of the total non-adiabatic pressure fluctuations at each depth, $`\delta P^{\mathrm{nad}}(z,t)`$. For each frequency and depth, this is multiplied by the spatial derivative of the mode displacement at that depth interpolated to the frequencies of the Fourier transform (and normalized by the square root of the mode energy), $`\omega (\xi _\omega /r)E_\omega ^{1/2}`$. The spatial dependence of the modes varies slowly and continuously with frequency, so such interpolation is possible. This product is integrated over the depth of the simulation domain for each frequency. The energy input rate to the modes is the square of the absolute value of this work integral, divided by the frequency interval for the Fourier transform, $`\mathrm{\Delta }\nu `$, (which equals multiplying by the time interval of the simulation), multiplied by the area of the simulation (36 $`\mathrm{Mm}^2`$) and divided by eight (Eqn 5).

In real space, the integrand of the work integral for each frequency, $`\nu `$, is calculated as the non-adiabatic total pressure fluctuation at each depth and time, $`\delta P^{\mathrm{nad}}(z,t)`$, multiplied by the normalized derivative of the mode displacement for that frequency at each depth, $`\omega (\xi _\omega (z)/z)/E_\omega ^{1/2}`$, multiplied by a phase function $`\mathrm{cos}\left(\varphi +\omega t\right)`$ for each time, $`t`$, a snapshot was saved in the simulation (10 s or 30 s intervals). For each depth, this is summed over all saved snapshots in the longest interval that is an integral number of mode periods for the given mode frequency $`\nu `$, and divided by the number of snapshots summed over to get the average value. This integrand is then integrated over depth. This is done for two values of the phase, $`\varphi =0\mathrm{and}\pi /2`$ (since the phases between the modes and the pressure fluctuations are random, we average these two orthogonal cases). The value of the integral for each phase is squared and these two values are summed. The energy input rate to the mode at this frequency is this average multiplied by the time interval integrated over (for that frequency), multiplied by the area of the simulation (36 $`\mathrm{Mm}^2`$), and divided by eight.