text:
# Instabilities and Spatio-temporal Chaos of Long-wave Hexagon Patterns in Rotating Marangoni Convection

## 1 Introduction

In the classic, buoyancy-driven Rayleigh-Bénard convection rotation of the system has been found to lead to very interesting new dynamic patterns. The origin of the dynamics is the Küppers-Lortz instability , which renders all straight roll patterns unstable with respect to a set of rolls rotated with respect to the original ones. It occurs above a critical rotation rate. In simple low-dimensional models the instability is found to a heteroclinic cycle connecting roll patterns rotated with respect to each other by $`120^o`$ . In experiments (e.g. ) as well as in numerical computations of the Navier-Stokes equations (e.g. ) and of model equations (e.g. ) more complex behavior is found. Patches of convection rolls of various orientations arise and the roll orientation changes through the propagation of the fronts separating them or by the nucleation of rolls with a new orientation.

The other classic convection system is surface-tension-driven Marangoni convection of a liquid layer in contact with another liquid or a gas layer (for a recent review see ). Two different instability mechanisms have been identified that drive this type of convection: a short-wave instability, which occurs independent of any deflections of the interface  and a long-wave instability, for which deformations of the interface are essential . Not too far from threshold, the short-wave instability typically leads to convection patterns with hexagonal planform (e.g. ), while the long-wave instability induces isolated dry spots or high spots . The stability of the hexagonal short-wave pattern with respect to sideband instabilities (the ’stability balloon’) has been determined in weakly nonlinear  and in fully nonlinear numerical approaches . Further above onset the short-wave hexagon pattern is often observed to undergo a transition to squares . If the mean surface tension of the fluid is not too large the short-wave convection mode has to be coupled to a weakly damped long-wave mode that corresponds to the local layer height of the fluid . A different type of long-wave instability has been identified theoretically for the case when the fluid layer is bounded above and below by poor thermal conductors . The mechanism for this instability is the same as for the short-wave instability. Due to the poor conductivity the critical wavenumber is, however, shifted to small values. Again one finds that near onset the planform is usually hexagonal. It turns out that experimentally this limit of long wavelengths is difficult to achieve. Theoretically, however, it is a very useful limit since it allows to derive simplified evolution equations systematically from the basic fluid equations without the restriction to small amplitudes, although restrictions on the control parameter remain . Moreover, compared to the weakly nonlinear Ginzburg-Landau equations the long-wave equations have the additional great advantage of preserving the isotropy of the system. This will play an important role in the present paper.

With the interesting behavior ensuing from the Küppers-Lortz instability, it is natural to ask whether rotation has also an interesting impact on hexagonal convection patterns. Within the general context of small-amplitude expansions and Ginzburg-Landau equations this question has been addressed previously . These analyses were aimed at a general investigation of hexagonal patterns in rotating systems and apply also to weakly non-Boussinesq Rayleigh-Bénard convection. It was found that in a rotating system a new state of oscillating hexagons can appear. This transition replaces then the transition from hexagons to rolls that is obtained in non-Boussinesq Rayleigh-Bénard convection. In Marangoni convection this transition does not occur and no oscillating hexagons of this type are expected. Instead, in the absence of rotation a transition to squares is found . The effect of rotation on this transition has not been studied yet. For the usual, steady hexagons the weakly nonlinear analysis indicates the possibility of persistent chaotic dynamics coexisting with the steady hexagons for the same parameter values . However, these dynamics cannot be captured by the Ginzburg-Landau equations employed. To make reliable statements about this chaotic state, which is characterized by a spatially isotropic Fourier spectrum, equations have to be studied that preserve the isotropy of the system. A first step in this direction was the study of a Swift-Hohenberg model, which was generalized to capture the breaking of the chiral symmetry that characterizes the rotation. A detailed analysis of the resulting equation, which was not derived from any physical system, revealed that steady hexagons of all wavenumbers can become unstable already for amplitudes below the transition to the oscillating hexagons, leaving a range in the control parameter over which no simple periodic solution is stable  and in which chaotic dynamics are found in numerical simulations.

Motivated by the previous results of weakly nonlinear approaches to hexagons in rotating systems, our goal for the present paper is to obtain a systematic description of Marangoni convection in a rotating system that is not limited to small amplitudes and that preserves the isotropy of the system. We therefore turn to the long-wave limit that is obtained for poorly conducting top and bottom boundaries and extend the derivation of  to the case of a rotating system. The poor heat conductivity is expressed in terms of a small Biot number $`Bi`$. While the true limit of long wavelengths ($`Bi1`$) may not be attainable in experiments  we expect to gain valuable insight into the system for moderate wavelengths by expanding in $`Bi`$. The form of the resulting long-wave equation is the same as that found for buoyancy-driven convection with poorly conducting boundaries  and for buoyancy-driven convection in binary mixtures with positive separation ratio . In the latter case the critical wavenumber goes to 0 at a critical value of the separation ratio. In both cases quadratic terms arise only if the vertical symmetry is broken by imposing different boundary conditions at the top and the bottom plate (free-slip vs. no-slip or fixed temperature vs. fixed heat flux). The focus of  is the effect of rotation on the system. In contrast to our work, however, there the competition between strictly periodic roll and square planforms is considered. The work is performed analytically in the weakly nonlinear regime, supplemented by a few simulations of a small system.

In sec.2 we introduce the basic fluid equations and associated boundary conditions. The general linear stability analysis that identifies the long-wave regime is performed in sec.3. The long-wave equation is derived in sec.4. The stability and dynamics of the hexagon patterns is then studied numerically in sec.5. The main results of the simulations are periodic and chaotic dynamics arising from secondary instabilities of the hexagons as well as from a resonant mode interaction.

## 2 Formulation of the problem

### 2.1 Basic equations

We consider a thin viscous and incompressible fluid layer heated from below with a free surface and rotating around a vertical axis. The mean surface tension is assumed to be large enough to suppress surface deformations even in the presence of rotation. The Navier-Stokes and the heat equation then read

$`\rho \left(_t+𝐯\right)𝐯`$ $`=`$ $`p+{\displaystyle \frac{\rho }{2}}|\mathrm{\Omega }𝐞_z\times 𝐫|^2+2\rho \mathrm{\Omega }𝐞_z\times 𝐯+\eta \mathrm{\Delta }𝐯\rho g𝐞_z,`$ (1)
$`𝐯`$ $`=`$ $`0,`$ (2)
$`\left(_t+𝐯\right)T`$ $`=`$ $`\kappa ^2T.`$ (3)
In this paper we focus on fluid layers that are sufficiently thin that buoyancy can be ignored. The density is then taken to be constant, $`\rho =\rho _0`$. The dynamic viscosity $`\eta `$ is related to the kinematic viscosity $`\nu `$ through $`via`$ $`\nu =\eta /\rho _0`$ and $`\kappa `$ is the thermal diffusivity. The equations in dimensionless form can be obtained if the space, time, angular velocity, temperature, pressure and velocity fields are respectively divided by the constants $`d`$, $`d^2/\kappa `$, $`\kappa /d^2`$, $`\mathrm{\Delta }T`$, $`\rho _0\kappa \nu /d^2`$ and $`\kappa /d`$. Here $`d`$ is the thickness of the fluid layer and $`\mathrm{\Delta }T`$ is the temperature difference across the liquid layer in the conductive state. They are written as follows,

$`{\displaystyle \frac{1}{𝒫}}\left(_t+𝐮\right)𝐮`$ $`=`$ $`p+{\displaystyle \frac{𝒫}{2}}|𝛀\times 𝐫|^2+2𝛀\times 𝐮+^2𝐮𝒢𝐞_z,`$ (4)
$`𝐮`$ $`=`$ $`0,`$ (5)
$`\left(_t+𝐮\right)\theta `$ $`=`$ $`^2\theta ,`$ (6)
where $`𝒫=\nu /\kappa `$ is the Prandtl number and $`𝒢=d^3\rho g/(\eta \nu )`$ is the Galileo number. In this paper we consider the limit of infinite Prandtl number.

We eliminate the pressure from the Navier-Stokes equations (4) by using the incompressibility condition (2). Taking the curl and the double curl of (1) we obtain

$`^4w2\mathrm{\Omega }_z\zeta `$ $`=`$ $`0,`$ (7)
$`^2\zeta +2\mathrm{\Omega }_zw`$ $`=`$ $`0,`$ (8)
where $`w`$ is the $`z`$-component of the velocity $`𝐮=(u,v,w)`$ and $`\zeta `$ is the $`z`$-component of the vorticity $`\times 𝐮`$. We also introduce the horizontal Laplacian $`\mathrm{\Delta }_H=(_{xx}^2+_{yy}^2)`$. Note that the rotation couples the vertical velocity and the vertical vorticity. Thus, while in the absence of rotation the vertical vorticity vanishes in many solutions of interest, with rotation the vertical vorticity is in general non-zero.

The equations for the $`x`$ and $`y`$-components of the velocity field are obtained considering that the velocity field $`𝐮=(u,v,w)`$ is expressed in terms of the poloidal and the toroidal stream functions $`\varphi `$ and $`\xi `$ as follows: $`𝐮=(_x_z\varphi +_y\xi ,_y_z\varphi _x\xi ,\mathrm{\Delta }_H\varphi )`$. Then the equations for $`u`$ and $`v`$ are:

$`\mathrm{\Delta }_Hu`$ $`=`$ $`_x_zw+_y\zeta ,`$ (9)
$`\mathrm{\Delta }_Hv`$ $`=`$ $`_y_zw_x\zeta .`$ (10)
These equations together with (7)-(8) are equivalent to (4)-(5) .

### 2.2 Boundary conditions

In order to obtain the velocity boundary condition at the top free surface we take into account thermocapillary effects and assume a linear dependence of the surface tension on the temperature, $`\sigma =\sigma _0\gamma (TT_0)`$. Furthermore, for sufficiently large mean surface tension and not too high rotation rates we can take the free surface to be flat. The boundary conditions at the upper free surface are then given by

$`\rho \nu \left({\displaystyle \frac{w}{x}}+{\displaystyle \frac{u}{z}}\right)`$ $`=`$ $`\gamma \left({\displaystyle \frac{T}{x}}\right),`$ (11)
$`\rho \nu \left({\displaystyle \frac{w}{y}}+{\displaystyle \frac{v}{z}}\right)`$ $`=`$ $`\gamma \left({\displaystyle \frac{T}{y}}\right).`$ (12)
After taking derivatives of (11) and (12) with respect to $`x`$ and $`y`$, respectively, and addition of both expressions, one is led to

$`\rho \nu \left({\displaystyle \frac{^2w}{z^2}}\right)`$ $`=`$ $`\gamma \mathrm{\Delta }_HT.`$ (13)
Introducing the Marangoni number $`M=\gamma d\mathrm{\Delta }T/(\rho \nu \kappa )`$ the boundary condition can be written as

$`{\displaystyle \frac{^2w}{z^2}}`$ $`=`$ $`M\mathrm{\Delta }_H\theta \text{ at }z=1.`$ (14)
From the Marangoni equations (11-12) and from the underformability of the surface the other boundary conditions at the top are

$`w=_z\zeta =0,`$ $`\mathrm{at}z=1.`$ (15)
At the bottom we take no-slip boundary conditions for the velocity ,

$`w=\zeta =_zw=0,`$ $`\mathrm{at}z=0.`$ (16)
The critical wavenumber at which convection sets in depends sensitively on the temperature boundary conditions. Long wavelengths are obtained if the boundaries are poor conductors. The Newton law at the bottom plate can be expressed as

$`\mathrm{k}_zT_{fb}`$ $`=`$ $`h_1(T_{fb}T_b)\text{ at }z=0,`$ (17)
where $`T_{fb}`$ is the temperature of the fluid at the bottom plate and $`T_b`$ is the temperature of the plate at some distance below this interface. The thermal conductivity of the liquid is given by k and $`h_1`$ is the heat exchange coefficient at the bottom. We non-dimensionalize the temperature field by setting $`\theta =(TT_{fbc})/\mathrm{\Delta }T`$, where $`T_{fbc}`$ is the temperature of the fluid in the conductive state at the bottom plate and $`\mathrm{\Delta }T`$ is again the temperature difference across the layer in the conductive state, i.e. $`\mathrm{\Delta }T=T_{fbc}T_{ftc}`$ with $`T_{ftc}`$ the upper temperature of the fluid in the conductive state. In terms of $`\theta `$ the boundary condition (17) is then given by

$`_z\theta `$ $`=`$ $`Bi_1\theta 1\text{ at }z=0`$ (18)
with the Biot number $`Bi_1=h_1d/\mathrm{k}`$. Here we have made use of the conservation of heat flux at the interface, $`h_1(T_{fbc}T_b)=\mathrm{k}/d(T_{fbc}T_{ftc})`$. In a similar way the temperature boundary condition at the top is found to be

$`_z\theta `$ $`=`$ $`Bi_2\theta 1Bi_2\text{ at }z=1.`$ (19)
These expressions agree with those used in  when the limit of a single layer is taken (in their notation this limit corresponds to $`\mathrm{\Theta }_1=\mathrm{\Theta }_2`$, $`\kappa =1`$ and $`a=0`$).

## 3 Linear stability analysis

In order to motivate the long-wavelength analysis we first discuss the linear stability of the conductive solution. The normal modes of the linearized system (6-8) together the boundary conditions (14-16), (18-19), take the form

$$(𝐮,\theta )=(\widehat{𝐮}(z),\widehat{\theta }(z))e^{i𝐪𝐫+\mathrm{\Lambda }t},$$

where $`\mathrm{\Lambda }`$ is the growth rate, $`𝐪`$ is the wavevector and r is the position vector $`(x,y)`$. Since we are interested only in thermocapillary effects the Marangoni number is the sole control parameter. At the critical value $`M_c`$ of the Marangoni number, determined by a vanishing real part of $`\mathrm{\Lambda }`$, the conductive state first becomes unstable. If the growth rate is $`\mathrm{\Lambda }=0`$ at $`M_c`$ the bifurcation is steady, while if $`\mathrm{\Lambda }=i\omega `$ for some nonzero real frequency $`\omega `$, the bifurcation is oscillatory. In this work we examine only the steady bifurcation. It is known that for sufficiently large rotation rates an oscillatory bifurcation may occur . In the absence of rotation, the linearized system (6-8) yields the critical Marangoni number $`M_c`$ as a function of the wave number $`|𝐪|`$ and the Biot numbers $`Bi_{1,2}`$

$`M=8q((q^2\mathrm{cosh}q+q\mathrm{cosh}q^2\mathrm{sinh}q)(Bi_1+Bi_2)`$ (20)
$`+\left(q\mathrm{sinh}q+\mathrm{cosh}q^3\mathrm{cosh}q\right)Bi_1Bi_2`$
$`+q^2\mathrm{cosh}q^3q^2\mathrm{cosh}qq^3\mathrm{sinh}q)/A,`$
where,

$`\mathrm{A}=\left(\mathrm{cosh}q^2\mathrm{sinh}q\mathrm{sinh}qq^3\mathrm{cosh}q\right)Bi_1+q^3\mathrm{cosh}q`$ (21)
$`2q^2\mathrm{sinh}q+q\mathrm{cosh}q^3q\mathrm{cosh}qq^4\mathrm{sinh}q).`$ (22)
As shown in Fig.1a the critical wavenumber $`q_c`$ vanishes for $`Bi_1=Bi_2=0`$, while for $`Bi_1>0`$ and $`Bi_2>0`$ one has $`q_c>0`$. A series expansion of (20) around $`Bi_1=Bi_2=0`$ shows that $`q_c=((Bi_1+Bi_2)405/27)^{1/4}`$. Thus, for poorly conducting boundary conditions the wavenumber of the growing pattern is still small enough to allow the description by a long-wave model. In Fig.1b the marginal curve is shown for $`Bi_1=Bi_2=0`$ and two values of the rotation rate, $`\mathrm{\Omega }=1`$ and $`\mathrm{\Omega }=5`$. We see that small values of the rotation rate do not modify the critical wavenumber. However for $`\mathrm{\Omega }>\mathrm{\Omega }_c3.7513`$ a minimum develops in the neutral curve even for $`Bi_1=Bi_2=0`$ and the critical wavenumber changes sharply from zero to a finite value. Thus, our study will be valid only for $`\mathrm{\Omega }<\mathrm{\Omega }_c`$. In this regime the critical Marangoni number at $`q=0`$ as a function of $`\mathrm{\Omega }`$ is given by

$`M_c(\mathrm{\Omega })|_{q=0}={\displaystyle \frac{4\mathrm{\Omega }^{3/2}\left(2\mathrm{cos}(\sqrt{\mathrm{\Omega }})\mathrm{sin}(\sqrt{\mathrm{\Omega }})\mathrm{sinh}(2\sqrt{\mathrm{\Omega }})\right)}{4\mathrm{sinh}(\sqrt{\mathrm{\Omega }})\mathrm{sin}(\sqrt{\mathrm{\Omega }})+\mathrm{cos}(2\sqrt{\mathrm{\Omega }})\mathrm{cosh}(2\sqrt{\mathrm{\Omega }})}}.`$ (23)
## 4 Derivation of the Long-Wave Equation

Here we derive the long-wavelength model for rotating Marangoni convection and consider the system (3)-(10) together with the boundary conditions (14)-(15). Expanding $`\theta `$ and $`𝐮`$ around the conductive solution, we rewrite the equations in terms of the deviations from the basic solution. To avoid excessive notation we keep the same notation $`\theta `$ and $`𝐮`$ for deviations of the temperature and velocity fields and obtain

$`^4w2\mathrm{\Omega }_z\zeta `$ $`=`$ $`0,`$ (24)
$`^2\zeta +2\mathrm{\Omega }_zw`$ $`=`$ $`0,`$ (25)
$`\mathrm{\Delta }_Hu+_x_zw+_y\zeta `$ $`=`$ $`0,`$ (26)
$`\mathrm{\Delta }_Hv+_y_zw_x\zeta `$ $`=`$ $`0,`$ (27)
$`\left(_t+𝐮\right)\theta `$ $`=`$ $`w+\mathrm{\Delta }\theta .`$ (28)
The boundary conditions are:

$`\theta _zBi_1\theta =w=\zeta =w_z=0,`$ $`\mathrm{at}z=0,`$ (29)
$`w_{zz}M\mathrm{\Delta }_H\theta =\theta _z+Bi_2\theta =w=\zeta _z=0,`$ $`\mathrm{at}z=1.`$ (30)
We define a small parameter $`ϵ`$ which measures distance from the critical Marangoni number by setting:

$`M=M_c(1+ϵ^2\mu ).`$ (31)
Here $`\mu `$ is a parameter which allows to introduce changes in the Marangoni number independently of $`ϵ`$. We assume that the heat loss is small,

$`Bi=ϵ^4b.`$ (32)
As discussed above, this implies that the dimensionless wavelength of the initial instability is long, of order $`ϵ^1`$. This allow us to introduce the slow scales,

$`X=ϵx,Y=ϵy.`$ (33)
The vertical coordinate $`z`$ is not rescaled. The evolution occurs on the slow time $`T`$,

$`t=ϵ^4T,`$ (34)
as suggested by the linear stability analysis . These scalings of space and time suggest the scaling of the fields involved. Specifically, the Marangoni condition (30) suggests that the vertical velocity be taken to be small while the deviation of the temperature from the linear profile can be large,

$`(w,\theta )=(ϵ^2W,\mathrm{\Theta }).`$ (35)
Equations (26) and (27) point out the scaling for the $`x`$ and $`y`$ component of the velocity,

$`(u,v)=(ϵU,ϵV).`$ (36)
From the definition of the $`z`$-component of the vorticity and considering the previous scaling we have,

$`\zeta =ϵ^2Z.`$ (37)
With these definitions the equations (24)-(28) become:

$`\mathrm{D}^4W+2ϵ^2\mathrm{D}^2\mathrm{\Delta }_HW+ϵ^4_H^42\mathrm{\Omega }\mathrm{D}Z`$ $`=`$ $`0,`$ (38)
$`\mathrm{D}^2Z+2\mathrm{\Omega }\mathrm{D}W+ϵ^2\mathrm{\Delta }_HZ`$ $`=`$ $`0,`$ (39)
$`\mathrm{\Delta }_HU+W_{Xz}+Z_Y`$ $`=`$ $`0,`$ (40)
$`\mathrm{\Delta }_HV+W_{Yz}Z_X`$ $`=`$ $`0,`$ (41)
$`ϵ^4\mathrm{\Theta }_T\mathrm{D}^4\mathrm{\Theta }ϵ^2\mathrm{\Delta }_H\mathrm{\Theta }ϵ^2W`$ $`=`$ $`ϵ^2U\mathrm{\Theta }_Xϵ^2W\mathrm{\Theta }_Yϵ^2W\mathrm{D}\mathrm{\Theta },`$ (42)
with D$`/z`$. The boundary conditions are,

$`W=\mathrm{D}W=Z=\mathrm{D}\mathrm{\Theta }b_1ϵ^4\mathrm{\Theta }=0,`$ $`\mathrm{on}z=0,`$ (43)
$`W=\mathrm{D}^2WM_c(1+ϵ^2\mu )\mathrm{\Delta }_H\mathrm{\Theta }=\mathrm{D}Z=\mathrm{D}\mathrm{\Theta }b_2ϵ^4\mathrm{\Theta }=0,`$ $`\mathrm{on}z=1.`$ (44)
Proceeding as in Refs. , we expand the temperature as follows,

$$\mathrm{\Theta }=\mathrm{\Theta }_0(X,Y,z,T)+ϵ^2\mathrm{\Theta }_1(X,Y,z,T)+\mathrm{},$$

and similarly for the other fields. Because in (38)-(44) only even powers of $`ϵ`$ appear, the expansion of the fields can also be restricted to even powers. If we insert these expansions in (38)-(44) we obtain at $`𝒪(ϵ^0)`$

$`\mathrm{\Theta }_0`$ $`=`$ $`F(X,Y,T),`$ (45)
$`W_0`$ $`=`$ $`M\mathrm{\Delta }_HFf_1(\mathrm{\Omega },z),`$ (46)
$`Z_0`$ $`=`$ $`M\mathrm{\Delta }_HFf_2(\mathrm{\Omega },z),`$ (47)
$`U_0`$ $`=`$ $`M(F_X\mathrm{D}f_2+F_Yf_1),`$ (48)
$`V_0`$ $`=`$ $`M(F_Y\mathrm{D}f_2F_Xf_1).`$ (49)
In contrast to the non-rotating case the vertical vorticity $`Z`$ does not vanish and at each order we must solve two coupled inhomogeneous differential equations instead of only one as in . In these expressions $`F(X,Y,T)`$ is the planform function whose governing equation we are looking for. The functions $`f_1`$ and $`f_2`$ are complicated functions of $`\mathrm{\Omega }`$ and $`z`$. In the limit $`\mathrm{\Omega }0`$ the function $`f_2`$ vanishes and $`f_1`$ is the same as in .

At $`𝒪(ϵ^2)`$ we obtain,

$$\mathrm{D}^2\mathrm{\Theta }_1=\mathrm{\Delta }_HFMf_1\mathrm{\Delta }_HFM\mathrm{D}f_1|_HF|^2.$$

The critical Marangoni number at $`k=0`$ is obtained from the solvability condition in agreement with the expression (23) obtained in the linear stability analysis. For the fields we obtain

$`\mathrm{\Theta }_1`$ $`=`$ $`f_3(\mathrm{\Omega },z)\mathrm{\Delta }_HF+f_4(\mathrm{\Omega },z)|_HF|^2+F_2(X,Y,T),`$ (50)
$`W_1`$ $`=`$ $`f_5(\mathrm{\Omega },z)\mathrm{\Delta }_HF+f_6(\mathrm{\Omega },z)\mathrm{\Delta }_H^2F+f_7(\mathrm{\Omega },z)\mathrm{\Delta }_H|_HF|^2+`$ (51)
$`f_8(\mathrm{\Omega },z)\mathrm{\Delta }_HF_2,`$
$`Z_1`$ $`=`$ $`f_9(\mathrm{\Omega },z)\mathrm{\Delta }_HF+f_{10}(\mathrm{\Omega },z)\mathrm{\Delta }_H^2F+f_{11}(\mathrm{\Omega },z)\mathrm{\Delta }_H|_HF|^2+`$ (52)
$`f_{12}(\mathrm{\Omega },z)\mathrm{\Delta }_HF_2,`$
$`U_1`$ $`=`$ $`f_9F_Yf_{10}\mathrm{\Delta }_HF_Yf_{11}(|_HF|^2)_Yf_{12}F_{2,Y}f_{5,z}F_X`$ (53)
$`f_{6,z}\mathrm{\Delta }_HF_Xf_{7,z}(|_HF|^2)_Xf_{8,z}F_{2,X},`$
$`V_1`$ $`=`$ $`f_9F_X+f_{10}\mathrm{\Delta }_HF_X+f_{11}(|_HF|^2)_X+f_{12}F_{2,X}f_{5,z}F_Y`$ (54)
$`f_{6,z}\mathrm{\Delta }_HF_Yf_{7,z}(|_HF|^2)_Yf_{8,z}F_{2,Y}.`$
Here $`F_2`$ is the $`O(ϵ^2)`$-correction to the planform function $`F`$. The functions $`f_3`$, $`f_4`$, $`f_5`$, $`f_6`$, $`f_7`$, $`f_8`$, $`f_9`$, $`f_{10}`$, $`f_{11}`$ and $`f_{12}`$ are again complicated functions of $`z`$ and $`\mathrm{\Omega }`$. The functions $`f_3`$ and $`f_4`$ are related to integrals of $`f_1`$ with respect to $`z`$. It turns out that after the integration these functions contain terms that are singular at $`\mathrm{\Omega }=0`$ and do not depend on $`z`$ in this limit. In order to cancel these singularities and to recover the expressions of  in the limit $`\mathrm{\Omega }0`$, the planform function $`F_2`$, which is undetermined at this order, must contain the same singular constants with the opposite sign. Since the integral expressions are too complicate to present here, this is best illustrated by the simple example

$`{\displaystyle \mathrm{sin}(\mathrm{\Omega }z)𝑑z}={\displaystyle \frac{1}{\mathrm{\Omega }}}\mathrm{cos}(\mathrm{\Omega }z)+C.`$ (55)
This integral is not finite in the limit $`\mathrm{\Omega }0`$. If a finite value is required for it, the arbitrary constant $`C`$ must be chosen as $`1/\mathrm{\Omega }`$. Analogously, in $`f_3`$ and $`f_4`$ constants appear that need to be chosen appropriately. Since they do not depend on $`z`$ they can be absorbed in $`F_2`$ and therefore do not affect the long-wave equation for $`F`$.

Finally, at $`𝒪(ϵ^4)`$ the solvability condition for D$`{}_{}{}^{2}\mathrm{\Theta }_{2}^{}`$ yields the required equation for the planform function $`F(X,Y,T)`$,

$`F_T`$ $`=`$ $`(b_1+b_2)F\mu \mathrm{\Delta }_HF\lambda \mathrm{\Delta }_H^2Fc_1_H|_HF|^2_HF`$ (56)
$`c_2_H\mathrm{\Delta }_HF_HFc_3\mathrm{\Delta }_H|_HF|^2c_4𝐞_z_HF\times _H\left\{\mathrm{\Delta }F\right\}`$
$`+c_5𝐞_z_HF\times _H\left\{(_HF)^2\right\}.`$
Here $`\mu `$ is the control parameter (cf. (31)) The constants $`\lambda `$, $`c_1`$, $`c_2`$, $`c_3`$, $`c_4`$ and $`c_5`$ are complicated functions of the rotation rate $`\mathrm{\Omega }`$. We plot them in Fig. 2. A series expansion of these functions around $`\mathrm{\Omega }=0`$, shows that in the limit $`\mathrm{\Omega }0`$, the same equation as in  is obtained:

$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{15}}{\displaystyle \frac{7}{1200}}\mathrm{\Omega }^2+𝒪\left(\mathrm{\Omega }^4\right),`$
$`c_1`$ $`=`$ $`{\displaystyle \frac{48}{35}}{\displaystyle \frac{104}{121275}}\mathrm{\Omega }^2+𝒪\left(\mathrm{\Omega }^4\right),`$
$`c_2`$ $`=`$ $`{\displaystyle \frac{1}{10}}+{\displaystyle \frac{1}{1575}}\mathrm{\Omega }^2+𝒪\left(\mathrm{\Omega }^4\right),`$
$`c_3`$ $`=`$ $`{\displaystyle \frac{3}{5}}+{\displaystyle \frac{1}{1575}}\mathrm{\Omega }^2+𝒪\left(\mathrm{\Omega }^4\right),`$
$`c_4`$ $`=`$ $`{\displaystyle \frac{233}{525}}\mathrm{\Omega }+{\displaystyle \frac{404}{72765}}\mathrm{\Omega }^3+𝒪\left(\mathrm{\Omega }^5\right),`$
$`c_5`$ $`=`$ $`{\displaystyle \frac{8}{35}}\mathrm{\Omega }{\displaystyle \frac{122}{1576575}}\mathrm{\Omega }^3+𝒪\left(\mathrm{\Omega }^5\right).`$
## 5 Numerical results

We now turn to the stability and dynamics of hexagon patterns described by (56). One of our goals is to identify parameter regimes in which regular hexagon patterns become unstable at all wavenumbers. In the weakly nonlinear case hexagons can be described by three coupled Ginzburg-Landau equations. Within this framework their stability and dynamics has been studied in detail without and with rotation . It has been found that for small amplitudes there is always a range of wavenumbers over which the hexagons are stable. In a simple Swift-Hohenberg model, situations were identified in which somewhat further above onset the stability balloon closes, leading to persistent chaotic dynamics . In this paper we therefore do not reduce (56) to small-amplitude Ginzburg-Landau equations. Instead determine the stability of the hexagon solutions of (56) and the resulting dynamics numerically without any restriction to small amplitudes.

The stability of the hexagon patterns is determined with a Fourier-Galerkin code on a hexagonal lattice. We then determine their stability with respect to side-band perturbations using a Floquet analysis with a Floquet parameter $`𝐩=p(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$. Due to the isotropy of the system and the discrete rotational symmetry of the hexagons we can restrict our analysis to perturbation wavevectors subtending an angle $`\theta `$ of up to $`\pm 30^o`$ relative to the hexagon modes. It turns out that accurate solutions of (56) require a large number of Fourier modes. One measure of the accuracy of a solution is obtained with the stability analysis by determining the magnitude of the eigenvalue corresponding to the translation mode given by the Floquet parameter $`𝐩=0`$. For the exact solution one has $`\sigma (𝐩)=0`$, while numerically truncated solutions will have $`\sigma (𝐩)0`$. The dependence of $`\sigma (𝐩)`$ as a function of the number $`n`$ of harmonics retained is illustrated in table 1, where $`n=1`$ corresponds to the basic hexagon solution.

Fig.3 shows the resulting stability limits in the absence of rotation ($`\mathrm{\Omega }=0`$) for $`b_1+b_2=0.025`$. They are given by steady long-wave instabilities with $`p1`$ (cf.). Numerical simulations for selected parameter values (marked by plusses and crosses in fig.3) confirm that the instabilities lead to the formation of defects which then change the wavevector of the pattern.

Rotation strongly modifies the stability limits as shown in fig.4 for $`\mathrm{\Omega }=3`$ and $`b_1+b_2=0.025`$. The most striking aspect is that the stability balloon closes for $`\mu >0.07`$. This closing is due to a steady short-wave instability (triangles) which supersedes the usual steady long-wave instability (open circles) for larger $`\mu `$. In the short-wave regime the largest growth rates are obtained for $`p0.6q`$ and $`\theta 25^o30^o`$. The appearance of a short-wave instability is somewhat similar to the result found in the weakly nonlinear analysis. There it has been found that in the presence of rotation a steady or oscillatory instability arises at finite modulation wavenumbers and preempts the long-wave oscillatory instability that emerges from the interaction of the two long-wave phase modes .

A further interesting feature of the stability limits for $`\mathrm{\Omega }=3`$ is the oscillatory instability that limits the stability balloon on the low-$`q`$ side (solid symbols). For low values of $`\mu `$ it is a long-wave side-band instability, i.e. its growth rate and frequency vanish for $`𝐩=0`$ (solid circles). For larger $`\mu `$, however, the largest growth rate occurs at $`𝐩=0`$ and the hexagon pattern is unstable to a homogeneous oscillatory mode (solid squares). In the weakly nonlinear analysis such a homogeneous oscillatory instability has only been identified in the context of the transition from hexagons to rolls . No such connection is apparent in the present system.

Anticipating interesting dynamics for values of $`\mu `$ above the closing of the stability balloon, we have studied its dependence on the rotation rate. As a function of $`\mathrm{\Omega }`$, fig.5 shows the maximal control parameter $`\mu `$ for which there are still stable steady hexagons (solid squares). To indicate the $`\mu `$-range of stable hexagons the critical control parameter $`\mu _c`$ is also given (dashed line). With decreasing $`\mathrm{\Omega }`$ that $`\mu `$-range grows rapidly. To resolve the stability limits for yet smaller values of $`\mathrm{\Omega }`$ a very large number of modes is needed; $`n=19`$ is not sufficient (cf. table 1).

The nonlinear evolution ensuing from the linear instabilities is investigated by simulating (56) numerically with a pseudo-spectral Fourier code using an integrating-factor, fourth-order Runge-Kutta scheme. It turns out that for larger values of $`\mu `$ the nonlinear gradient terms in (56) can under certain conditions lead to an apparently spurious blow-up of the numerical solution in which very high wavenumbers are excited <sup>3</sup><sup>3</sup>3Even though the numerical time-stepping routine may lead to a blow-up, no change in the behavior of the stationary states is found with the Galerkin code for those parameters. Apparently, in these cases the Fourier representation of the derivatives in the time-stepping code does not provide sufficient numerical dissipation. By introducing an artificial dissipation term these divergencies can be avoided. To obtain sufficiently smooth starting solutions we have at times used a filter that affects only wavenumbers above a certain high wavenumber and corresponds then to a linear 8th-order derivative term. However, all the simulations presented in this paper have been performed with that filter turned off.. However, if a sufficiently smooth initial condition was used no such numerical divergencies occurred. As already indicated earlier (cf. table 1) the proper resolution of the patterns requires a large number of Fourier modes. In all simulations we have used a rectangular computational domain of size $`L_X=2L`$ and $`L_Y=2L/\sqrt{3}`$ to accommodate undeformed hexagon patterns. Specifically, for $`L=N\mathrm{\hspace{0.17em}2}\pi /q`$ this domain allows $`N`$ periods of the pattern in the $`x`$\- and $`y`$-direction.

To test the stability results we have performed simulations near the stability limits depicted in figs.3,4. The plusses denote parameters for which a small modulation of the pattern did not destabilize it, while at the crosses the hexagons became unstable. These simulations were performed in systems of size $`L=42\pi /q`$. For this aspect ratio some of the instabilities are still slightly suppressed shifting the stability limits somewhat. The limits obtained in the simulations are, however, consistent with the stability results if the modulation wavenumbers $`𝐩`$ are restricted to those compatible with the computational domain used in the simulations.

We first address the nonlinear evolution of the homogeneous oscillatory instability (solid squares). Simulations for $`\mu =0.060`$ and $`q=0.93`$ show that this stability limit is due to a supercritical Hopf bifurcation, i.e. the spatially homogeneous oscillations saturate and lead to a state of oscillating hexagons. They are similar in appearance to the oscillating hexagons that arise in the Hopf bifurcation that replaces the instability of hexagons to the mixed mode, which is associated with the transition to rolls. As mentioned earlier, this transition to the mixed mode and to rolls has not been seen in the present system. A time sequence of snapshots of the oscillating hexagons is shown in fig.6<sup>4</sup><sup>4</sup>4A movie of the temporal evolution of the oscillating hexagons can be found at http://www.esam.northwestern.edu/riecke/research/Marangoni/marangoni.html..

Now we turn to the nonlinear evolution arising from the steady short-wave instability (open left-triangles in fig.4). In fig.7 and fig.8 a sequence of the pattern and the temporal evolution of two modes is shown, respectively in a system of size $`L_X=20.8`$ for $`\mu =0.07`$ (using $`128\times 128`$ modes). For this value of the control parameter no stable hexagon patterns exist in the infinite system. In this case the temporal evolution does not become periodic for the time period investigated.

In fig.9 typical snapshots for the evolution of the pattern in a system of size $`L_X=48`$ are shown for $`\mu =0.08`$ (using $`256\times 256`$ modes). The associated time series for two Fourier modes are shown in fig.10. Again, for these parameter values no steady stable hexagons exist in large systems resulting in persistent chaotic dynamics.<sup>5</sup><sup>5</sup>5A movie of the chaotic dynamics for $`\mu =0.075`$ and $`L_X=48`$ can be found at http://www.esam.northwestern.edu/riecke/research/Marangoni/marangoni.html. While the pattern is spatially quite disordered a tendency to form locally square arrangements can be observed. This presumably reflects the fact that in the absence of the quadratic terms in (56) even in the rotating system square patterns are preferred over roll patterns .

Finally, in Figs.11, 12 we show the temporal evolution in a quite small system of size $`L_X=17`$ for $`\mu =0.058`$ (using $`64\times 64`$ Fourier modes). For such a small system the stability limits as shown in figs.4 do not really apply. The initial condition was a regular hexagon pattern with wavenumber $`q=1.478`$. After a transient during which defects appeared and changed the number of cells the time dependence of the pattern becomes quite regular. The time series (cf. fig.12) is, however, not strictly periodic but shows small ($`1\%`$) irregular variations of the oscillation amplitude. A closer inspection of the neutral curve, which in terms of the system size and of the number $`N`$ of wavelengths contained in the system is given by

$`\mu (L_x)=b\left({\displaystyle \frac{L_X}{4\pi N}}\right)^2+\lambda \left({\displaystyle \frac{4\pi N}{L_X}}\right)^2,`$ (57)
shows that the parameters for this simulation are close to the codimension-two point at which hexagons with $`q=0.74`$ and hexagons with $`q=1.48`$ first arise. This is illustrated in fig.13 which gives the neutral curve (57) for various values of $`N`$. In this small system we associate the origin of the dynamics with the low-dimensional interaction of the two modes associated with the codimension-two point.

## 6 Conclusion

In this paper we have considered Marangoni convection in the presence of rotation. In order to derive systematically a set of reduced equations that - unlike the Ginzburg-Landau equations - preserve the isotropy of the system we have considered poorly conducting boundary conditions. Under these conditions convection sets in first at long wavelengths, which allows the derivation of a long-wave equation. Extending previous analyses , we have derived the long-wave equation in the presence of rotation and obtained the coefficients as a function of the rotation rate $`\mathrm{\Omega }`$. Using a Galerkin approach we have then determined the linear stability of steady hexagonal convection with respect to general side-band perturbations without and with rotation. We have employed direct numerical simulations of the long-wave equation to investigate the fate of these instabilities.

In the absence of rotation the band of stable steady hexagons in large systems is determined solely by a steady long-wave instability for the range of the control parameter considered. With rotation, however, steady long- and short-wave instabilities as well as oscillatory instabilities arise. The latter can be a long-wave side-band instability or an instability that is spatially homogeneous at its onset. Interestingly, the latter is supercritical and leads to stable, spatially homogeneous oscillations of the hexagons. This state resembles the oscillating-hexagon branch that arises in non-Boussinesq buoyancy-driven convection in connection with the transition from hexagon to roll convection . In the present system no transition to rolls has been identified and in experiments on Marangoni convection a transition from hexagons to squares is observed instead . We find temporally almost periodic and spatially more complex states that appear to arise from a resonant $`q:2q`$ mode interaction. This and other resonances could be studied in future work. Most interesting is the fact that for sufficiently large rotation rates the stability limits close for larger values of the control parameter. In the numerical simulations this is shown to lead to spatially and temporally chaotic patterns.

So far, there have been very few experimental studies of spatio-temporal chaos arising from a hexagonal planform . Therefore, it would be particularly interesting to see whether the complex dynamics of the disordered hexagonal patterns predicted by the long-wave model (56) can be obtained in experiments. It is yet to be seen to what extent in the current system the deformation of the free surface by the centrifugal force modifies the results presented here. Similarly, the derivation of the long-wave equation depends on the poor thermal conductivity of the horizontal boundaries. For experimentally realizable boundary conditions the wavelength of the convecting state may not be as large . As is often the case, however, the instabilities predicted by the long-wave equation (56) may persist qualitatively beyond the asymptotic long-wave limit. A study of this question will presumably require the derivation of the Ginzburg-Landau equations with nonlinear gradient terms or, more likely due to the finite amplitudes at the closing of the stability balloon, a numerical analysis of the full fluid equations (cf.).

AMM thanks the Gobierno de Navarra for a Postdoctoral Fellowship and the Universidad de Navarra for partial support. Part of this work was done during a visit by AMM at Northwestern University, which she would like to thank for its hospitality and support. We would also like to thank the Centro de Astrobiología and INTA for providing access to their computer systems where most of the numerical simulations were performed. HR gratefully acknowledges support by NASA (NAG3-2113), the Engineering Research Program of the Office of Basic Energy Sciences at the Department of Energy (DE-FG02-92ER14303), and by a grant from NSF (DMS-9804673).


format: text

text:
# Testing of Hamamatsu R5900-00-M64 Multi-Pixel PMTs for MINOS

## 1 Introduction

The MINOS (Main Injector Neutrino Oscillation Search) Experiment has been designed to identify neutrino oscillations and measure their parameters in the controlled environment of the NuMI neutrino beam at Fermi National Accelerator Laboratory (FNAL). The beam of predominantly muon neutrinos will pass through the Near Detector (ND) a calorimeter composed of planes of plastic scintillator and steel, at FNAL where the original flavour composition of the beam will be measured. The Far Detector (FD), a 5.4KT calorimeter of similar composition, situated in the Soudan mine in Northern Minnesota, will measure the flavour composition of the same beam $`735`$ km further away.

Both detectors will use extruded scintillator strips, read out via wavelength shifting (WLS) fibres glued into a central groove, to detect energy depositions of particles traversing the detector. The light from the WLS fibres is transmitted by clear acrylic fibres to 16-pixel Hamamatsu R5900-00-M16 photomultiplier tubes (PMTs) in the FD and to 64-pixel Hamamatsu R5900-00-M64 (M64) PMTs in the ND. Two different types of PMTs are required in order to take economical advantage of optical summing, possible in the low rate environment of the FD but not in the high-rate environment of the ND.

The results from testing fifteen M64s are reported in this paper.

The article is organised as follows: Section 2 describes the requirements on the tubes from the MINOS experiment while Section 3 gives technical details about the M64. Results from dark noise measurements are given in Section 4, while Section5 describes the details of the measurements of different PMT properties.

## 2 Required performance of PMTs for MINOS

The main requirements on the PMT characteristics are the following:

* good single photoelectron (p.e.) sensitivity;
* linearity up to 200 p.e.;
* single-channel gain above $`3\times 10^5`$ and up to $`10^6`$;
* small inter-pixel cross-talk;
* uniformity of response that does not significantly change in the presence of a magnetic field of up to $`5`$ gauss;
* tube dark count rate of less than $`3`$ kHz above a $`1/3`$ p.e. threshold.

In order to achieve the proposed physics goals, MINOS has to detect minimum ionising particles (MIPs) and has to measure electromagnetic showers. For a MIP, the energy released in the scintillator corresponds to about $`28`$ p.e. per strip traversed. On the other hand, electromagnetic showers can produce light signals as large as 200 p.e. It is therefore of the utmost importance that the photodetectors used have a wide dynamic range: they must be able to measure single p.e. signals precisely while still behaving as linear devices at light levels of up to 200 p.e.

The digitisation of the PMT signals is carried out roughly 1 m away from the tubes. Therefore, to avoid any effect due to electromagnetic pick-up noise, the PMT gain should be greater than $`3\times 10^5`$.

Due to the high beam intensity at the near detector, up to $`200`$ neutrino events may be present within the same $`11\mu `$s spill shortly after the primary protons hit the neutrino production target. To minimise complications in the event reconstruction algorithms, cross-talk to all neighbouring pixels of the same PMT should be below 10%.

The PMTs will be located close to the magnetised steel planes, where a fringe field of about $`5`$ gauss is expected. The PMT will be installed inside a steel housing, which will act as a shield and should reduce the field by at least $`80\%`$. However, it is important to verify that the tube response is not significantly affected by the remaining field and to have a safety margin.

All hits are continuously digitised by the front-end electronics and sent via readout processors to a “trigger farm”, which decides, depending on the event pattern, whether or not to save the data for further analysis. As the processing power and the bandwidth of the transfer are limited, it is important that the performance of MINOS is not affected by the presence of noise. The dark count rate of a single PMT channel should therefore be limited to about $`50`$ Hz. This is mainly due to bandwidth limitations of the DAQ system. The digitisation of all anode channels is started by the presence of a sizeable signal on the common last dynode. It is therefore important to study the noise properties of the anodes as well as of the last dynode.

## 3 Technical characteristics of M64 PMTs

Before describing in detail the procedure to evaluate the performance of the PMTs, it is useful to summarise the technical characteristics of the average M64 PMT as provided by Hamamatsu Photonics.

The Hamamatsu R5900-00-M64 is a $`12`$-dynode multi-anode PMT, with a transverse section of $`1\times 1`$ inches. The $`64`$ pixels (2 mm $`\times `$ 2 mm) are arranged on an $`8\times 8`$ square grid. These pixels are separated from each other by a distance of less than 0.3 mm, but there is a region of reduced collection efficiency extending to 0.3 mm.

Focusing electrodes run over the pixels which guide the photoelectrons into the dynode structure. Each pixel has two identical dynode chains associated with it which terminate in a single anode pad.

The cathode and the housing of the tube are at negative potential and a tapered bleeder circuit is used to apply voltage to the dynodes. The divider ratio adopted for the present measurements is the one recommended by Hamamatsu for optimal performance ($`3:2:2:1:1:1:1:1:1:1:1:2:5`$). Decoupling capacitors at the last stages help to keep the potential differences constant between the last dynodes and the anode in the presence of large instantaneous currents. The main parameters for the average PMT, as provided by Hamamatsu, are summarised in Table 1.

The only data available concerning the quantum efficiency (QE) of the PMTs are those given in the Hamamatsu data sheets. Fig. 1 shows QE as a function of the incident light wavelength $`(\lambda )`$ for $`15`$ tubes. The mean value of QE at $`\lambda =520`$ nm, the wavelength of interest for MINOS, is 13.5% with a spread of 0.9%.

For ten of the tubes, Hamamatsu measured the anode dark count rates at $`900`$ V and after $`30`$ minutes in the dark. The distribution of measured rates is shown in Fig. 2. Integrating over all the channels, the dark count rate is approximmately $`330`$ Hz. This is well below the maximum rate acceptable for MINOS.

## 4 Dark Noise

The light signal released by a minimum ionising particle in MINOS can be as low as two photoelectrons. It is therefore extremely important to study the PMT dark noise at very low threshold levels. Since the last dynode signal will be used in MINOS to trigger the anode readout, the dark count rate has to be measured for dynode thresholds corresponding to anode signals of approxiately $`1/3`$ photoelectron.

### 4.1 Set-up description

A block scheme of the experimental set-up used for these measurements  is shown in Fig. 3. The PMT under test is connected to its base and placed into a dark box. A removable plug is used to allow light into the otherwise light-tight box. The high voltage is supplied to the PMT by a Hewlett-Packard DC power supply (HP-E3631A) through a DC/DC converter circuit and monitored by a digital multi-meter (HP-34401A). Both are connected to a PC ($`500`$ MHz, Pentium III processor) through a GPIB interface.

The signals from the last dynode and from the anode pixel under study are sent via a $`50\mathrm{\Omega }`$ coaxial cable and an amplification stage (Phillips 779, nominal gain $`10.0\pm 0.2`$) to a four-channel Tektronix digital oscilloscope (TDS 540A). All the remaining anode channels are terminated through $`50\mathrm{\Omega }`$ to prevent charge build-up. The oscilloscope bandwidth is limited to 100 MHz to reduce noise.

The PC communicates with the oscilloscope via a LabVIEW based program. Once the triggering conditions have been set and the channel whose data will be acquired has been selected, the program reads in the waveform of the PMT pulse event by event and measures:

* the baseline level, to account for any DC offset by the amplifier;
* the signal amplitude.

The oscilloscope acquires 910 pulses at a rate in excess of 10 kHz before transferring the data to the PC. The oscilloscope is inhibited during the transfer which takes about one second.

To obtain an absolute measurement of the gain, although not crucial for the determination of the dark count rate, a voltage-to-charge conversion factor of $`(0.81\pm 0.16)\times 10^{11}\text{C}/\text{V}`$ has been measured by assuming shape universality for the PMT pulses.

### 4.2 Anode dark noise

Each of the tubes under test was kept in the dark and under high voltage ($`900`$ V) for at least six hours before performing the actual measurement. Low-light spectra and dark spectra were taken at high voltages between $`900`$ and $`1000`$ V for all pixels in one quadrant of the tube. The oscilloscope trigger threshold was set just above the anode pedestal peak which is at zero, as evidenced in Fig. 4 which shows the falling edge of the pedestal at very low pulse heights. The dark count rate was determined using an offline threshold of 1/3 of a photoelectron.

For all of the measured tubes, the low-light spectra of the pixels belonging to the inner $`6\times 6`$ matrix show a distinct single-photoelectron peak, well separated from the electronic noise. These pixels typically have relatively low dark count rates (a few Hz). Typical low-light and dark spectra of such pixels are shown in Fig. 4. Pixels along the edges of the PMT generally have higher noise rates, although in most cases still within acceptable limits. Of greater concern is the shape of their low-light spectra, since quite frequently the single-photoelectron peak is not very well distinguishable from the dynode noise contribution. Moreover a second peak, distinct from the noise but at pulse heights below the single-photoelectron peak, is often present (see Fig. 4). This is due to photons passing through the photocathode and hitting the first dynode where they are absorbed and a single photoelectron is emitted. This results in a peak in the distribution corresponding to an amplification of one dynode stage less than the main one photoeletron peak. The different behaviour between the inner and the outer pixels is expected to be due to different focusing conditions as well as to the difficulty in keeping the fabrication process uniform across the entire photocathode area.

The dependence of the dark noise on the PMT high voltage has also been studied. This is shown in Fig. 5. The dark count rate increases at higher values of high voltage. The single anode dark count rate at $`900`$ V is measured to be $`5.2`$ Hz, with a spread of $`3.0`$ Hz.

The temperature dependence of the dark count rate has been measured for one PMT. For both inner and outer pixels the gain does not vary substantially with temperature, while, as expected for bi-alkali photocathodes, the dark count rate increases exponentially with temperature (($`16\pm 3)\%`$ per C; see Fig. 6).

### 4.3 Dynode dark noise

The last dynode is common to all anodes, and so its single-photoelectron spectrum can be viewed as the sum of $`N_{pixel}`$ charge distributions, each of them corresponding to the single-photoelectron signal on one of the $`N_{pixel}`$ pixels. As the gain is different for each channel, in general no clear peak is visible when the spectrum is acquired by triggering on the dynode signal itself.

It is necessary to determine the noise rate at the dynode caused by signals above $`1/3`$ p.e. from the lowest gain pixel. To do this, the PMT was illuminated at a very low level. Initially, the spectrum from the lowest gain pixel was digitized and confirmed to be mostly one p.e. signals. Then, the pulses from this pixel were used to trigger the oscilloscope while the dynode signal was digitized. This enabled us to measure the contribution to the dynode spectrum from this channel and to determine the dynode threshold level which corresponded to $`1/3`$ p.e. on the lowest gain pixel as shown in Figure 7.

After light-sealing the dark box again, the PMT was kept under high voltage and in the dark for more than 24 hours. A dark spectrum of the last dynode was taken, triggering on the dynode signal itself. In order to reach the necessary sensitivity to the dark count rate, the trigger threshold had to be tuned somewhat below the $`1/3`$ p.e. level ascertained in the previous measurement, but high enough not to be dominated by electronic noise. Knowing the acquisition time, the integral dynode dark count rate as a function of the trigger threshold can be computed from the pulse height distribution. Knowing the relative positions of the anode and of the last dynode 1 p.e. peaks, the gain of the last amplification stage can be computed from simple charge conservation considerations. At $`950`$ V, the gain of the last dynode is $`3.5`$.

The integral dynode dark count rate as a function of the trigger threshold, which was obtained for both settings of high voltage, is shown in Fig. 7 for one of the PMTs. The count rate for a threshold of $`1/3`$ p.e. is equal to about $`275\text{Hz}`$ at $`950\text{V}`$ and to about $`330\text{Hz}`$ at $`1000\text{V}`$; as expected, the higher the gain, the higher the dark noise.

Following the procedure described above, the last dynode dark count rate at $`1/3`$ p.e. has been measured for all the available PMTs at $`950`$ V and after $`24`$ hrs in the dark. The resulting distribution for the 15 tubes is shown in Fig. 8. The mean value is about $`220`$ Hz, which confirms that the tubes under test are extremely quiet.

## 5 Results from PMT scanning

Most of the properties of the PMTs have been studied by measuring the phototube response to a point-like light source, which has been moved across the photocathode surface. This procedure is referred to as a “PMT scan”.

### 5.1 Set-up

Two different set-ups have been independently developed, and most phototubes have been tested in both setups. While in one set-up the PMT was moved to scan the light source across the photocathode, the other assembly used a moving plastic fibre. The results achieved with both setups are consistent. We therefore only describe the latter apparatus in detail. A schematic view of the setup with the movable fibre is shown in Fig. 9 .

The PMT under test is mounted into the scanning apparatus with the help of a tightly fitted plastic collar, which in turn is mounted onto an aluminium PMT holder. In order to prevent bad connections between the PMT and the base sockets, the base is screwed to the PMT holder from the back.

Light output from a programmably pulsed blue LED at a repetition rate of 16 kHz is used to illuminate the PMT. The blue light is carried by a clear acrylic fibre to an attenuation box. In this box the clear fibre is mounted at a fixed distance (90 degree angle) to a 1 m long WLS fibre in a similar LED illumination scheme as used in the MINOS scintillator modules. The WLS fibre carries the light to the PMT surface. For the measurements described here, the level was set to be between a few and 200 photoelectrons. To prevent possible transmission changes in the WLS fibre while moving it across the PMT cathode, a routing device was used to constrain the fibre to slide in an arc-shaped groove ($`25`$ cm diameter, 0.5 cm width). The fibre was connected to a non-magnetic armature that is in turn screwed to a motorised platform. This arrangement moves the fibre across the photocathode. The fibre itself is fitted into a black plastic cylinder which is gently pressed against the glass window by a spring.

The high voltage is supplied to the PMT as in the dark noise set up. The motorised platform is powered by a SmartDrive SA46 base unit, containing two D28 drives. The drives are controlled by a Trio Motion Technology Euro205 Motion Coordinator linked to the PC via a serial interface.

The 64 PMT anodes are connected to LeCroy 2249A ADC modules ($`10`$ bits, $`0.25`$ pC/count sensitivity); these are read out by a Jorway 73A CAMAC crate-controller coupled to the PC via a SCSI connection. Trigger signals from the LED pulser, coinciding with the light pulses, are used to gate the ADCs. The whole set-up is controlled by a LabVIEW 5.1 program and events are taken at a rate of $`1.52`$ kHz, i.e., not every light pulse is digitised.

The standard program of measurements for a PMT consists of the following steps:

1. an alignment run, to get the PMT centre coordinates with respect to the scanning table home position;
2. $`8\times 8`$ point scans (WLS fibre at the centre of the pixels) at different high voltages ($`U=`$ 850, 900, 950 V) and light levels (from a few to 100 photoelectrons);
3. determination of the lowest gain pixel and measurement of the dynode dark count rate at $`1/3`$ p.e. (see Section 4).

At the beginning of each scan a pedestal run is acquired, to determine the zero offset for each ADC channel. For each fibre position $`10,000`$ events were recorded, and the ADC information for all of the $`64`$ channels was stored as one-dimensional histograms.

### 5.2 Tools and Techniques

#### 5.2.1 Alignment

Every time a tube was installed in the PMT holder, an alignment run was performed before running the actual scan in order to determine the position of the PMT centre with respect to the home position of the scanning table. The home position is known to within $`10\mu `$m.

The alignment run consists of two scans, one along a row of the PMT $`(x`$-axis$`)`$ and the other along a column $`(y`$-axis$`)`$, with a step size of $`250\mu `$m $`(100`$ points per scan$`)`$ and with statistics of $`10,000`$ events per point. For each point, the pulse height in the pixel with the largest signal is plotted against the position of the WLS fibre: the profiles obtained in this way are similar to those in Fig. 10. The structure of the profile closely reflects the geometry of the PMT pixel: the peak amplitudes are proportional to the corresponding pixel gain. Each peak is fitted with a Gaussian function. The central position is calculated as the average of the eight fitted centres.

Several tests have been conducted in order to assess the reliability of the apparatus, especially concerning the reproducibility of its geometry and of the alignment procedure. Our studies show that the position of the PMT centre on the transverse plane is known to within $`12.5\mu `$m along the $`x`$ axis and to within $`23.0\mu `$ along the $`y`$ axis. Moreover, no appreciable effect from any rotation around the PMT symmetry axes has been observed.

#### 5.2.2 Photon Statistics

In order to make measurements of gain and efficiency, it is necessary to have an estimate of the amount of light impinging on the photocathode. It is possible to make such an estimate by assuming that the number of photoelectrons $`𝒩`$ emitted at the photocathode is described by Poisson statistics and that the amplification processes initiated by different photoelectrons along the dynode multiplicative system are independent. In this case the mean value ($`\mu `$) and the RMS ($`\sigma `$) of the pedestal-subtracted charge distribution of the PMT are given by:

$$\mu =𝒢𝒩\sigma =𝒢\sqrt{𝒩}$$
(1)
where $`𝒢`$ is the PMT gain. The number of photoelectrons and the PMT gain are then

$$𝒩=\left(\frac{\mu }{\sigma }\right)^2,𝒢=\frac{\mu }{𝒩}=\frac{\sigma ^2}{\mu }.$$
(2)
This method, referred to as the ‘photon-statistics method’, is the one used in the following paragraphs to calculate gain and number of photoelectrons for any analysed PMT. The method is known to be inaccurate in several respects:

* It does not take into account the fluctuations in the number of electrons produced at the first dynodes, i.e. the width of the one-photoelectron peak). Due to the sensitivity of our electronics we could not measure these fluctuations and did not correct for them.
* If the tube response is non-linear, $`𝒢`$ becomes a function of $`𝒩`$ and Eqs. (2) are no longer strictly valid. However, the measurements of efficiency and gain were all done at a level of approximately 10 p.e. where the tubes are in their linear regime.

It is still meaningful to use Eqs. (2) in order to separate gain and number of photoelectrons to compare different PMTs. We have verified that the LED pulser is a stable light source (Fig. 11); therefore the number of photoelectrons gives us a measure of the relative efficiency of the pixels within a tube. The efficiency here is a product of the collection efficiency (CE) and the quantum efficiency (QE). An absolute measure of the efficiency is not possible as the light source was not calibrated. The absolute value of the gain is estimated to be accurate to within approximately $`20\%`$.

### 5.3 Results

#### 5.3.1 Gain dependence on high-voltage

Scans at approximately 10 p.e. have been used in order to minimise possible effects arising from non-linearities, but still having sufficient sensitivity. The gains for all pixels of the tested PMTs have been measured at different high voltages (HV).

The plots in Fig. 12 show the average gain of all pixels for all of the 15 tubes as a function of HV ($`850,900,950`$ V). The gains of the pixels with the highest and lowest gains in each tube are also shown. The average gain increases from about $`1.3\times 10^6`$ at $`850`$ V to about $`3.9\times 10^6`$ at $`950`$ V. Although MINOS running will require gains of $`1.0\times 10^6`$, corresponding to high voltages of lower than 850V, the sensitivity of the electronics used in this test was not sufficient to evaluate the PMTs at that voltage.

The nominal gain spread ($`_𝒢=𝒢_{min}/𝒢_{max}`$, the gain ratio of the minimum and maximum gain pixels) for M64 PMTs is 1:3 ($`0.33`$). The average ratio measured is about $`0.4`$ at $`850`$ V and reaches $`0.46`$ at $`950`$ V with small differences between different tubes. The fact that $`_𝒢`$ decreases with decreasing HV is probably an artifact of saturation effects, which become more important at higher gain. Similar features can be recognised in Fig. 13, which shows the ratio $`_i^{(k)}=𝒢_i^{(k)}/𝒢_{max}^{(k)}`$ between each channel gain and the maximum gain of the same PMT, where $`i`$ runs over the PMT index and $`k`$ runs over the pixel number. The 15 entries at $`_i^{(k)}=1`$ in each histogram correspond to the 15 maximum gain pixels.

#### 5.3.2 Cross-talk

As all channels were read out constantly during the PMT scans, it was straightforward to measure the PMT cross-talk. The cross-talk $`𝒞`$ for a fibre position $`j`$ has been defined to be

$$𝒞=\frac{1}{Q_j}\underset{ij}{}Q_i,$$
(3)
where $`Q_j`$ is the charge of channel $`j`$ and the sum runs over a certain class of neighbouring channels. All charges are averaged over 10,000 events. The cross-talk into direct or diagonal channels as well as the total cross-talk has been measured. Most of the cross-talk is into direct neighbours and only much reduced cross-talk into diagonal pixels was observed. Results are summarised in Fig. 14, where the cross-talk for all the tested phototubes is displayed. The total cross-talk is about $`10\%`$ on average, with long tails extending to $`25\%`$. The cross-talk into non-neighbouring pixels sums up to about $`3\%`$.

As an example we have plotted the charge distributions for pixel 51 in Fig. 15. This measurement was performed with a modified base (voltage divider ratio $`3:2:2:1:1:1:1:1:1:1:1:1:1`$) and using a high bandwidth $`\times 10`$ amplifier to achieve the required sensitivity. Two spectra are compared: The solid spectrum shows the distribution when the pixel itself is pulsed with an average of 0.9 photoelectrons. The dashed distribution had been recorded when the neighbouring pixel 59 was illuminated with approximately 30 photoelectrons. From plots like this one can determine the kind of cross-talk present:

1. A distinct single photoelectron peak is visible when illuminating a neighbouring pixel. This suggests that some of the photons coming through the WLS fibre are not producing a photoelectron at the fibre location, but on adjacent pixels. This can be explained by the divergence of the light leaving the fibre and internal reflections in the PMT glass window.
2. The pedestal peak visible in Fig. 15 is different from the one observed when pulsing the pixel directly. It is shifted to higher values. This shift is proportional to the light level. It is an indication that a constant fraction of charge is seen at the neighbouring anode due, presumably, to charge leaking across to the neighbouring dynode structure during multiplication.

#### 5.3.3 Linearity studies

Photomultiplier tubes are in general reasonably linear devices: the charge at the anodes is proportional to the amount of light incident on the cathode. Nevertheless, when the incident light is too intense or when the bias voltage is too high, significant deviations from the ideal linear response are observed. This is mainly related to the dynode chain linearity characteristics, which are independent of the incident light wavelength and which, at fixed high-voltage values, are functions only of the peak current passing through the last dynode stages. Space charge effects are very important for PMTs operated in pulsed mode when intense light pulses can reach the photocathode and induce large current flows through the last dynode stages in addition to leaving charge between these dynodes. The magnitude of the consequent saturation effect is related to the internal tube geometry, the design of the resistive bleeder circuit, as well as to the electric field intensity and its spatial distribution.

The use of a tapered bleeder circuit (see Sec. 3) reduces such effects, as the large potential differences in the last few stages of the amplification chain enhance the voltage gradient where the electron space density becomes higher. Nonetheless, the M64 phototubes still exhibit a considerable non-linear response, especially when compared to M16 PMTs which is a consequence of the M64’s comparitively smaller pixel size. The different behaviour of the two types of phototubes is very clear in Fig. 16, where the relative outputs of the two photodetectors are shown as a function of the instantaneous PMT current when one pixel is illuminated (Hamamatsu data): while M16s are essentially linear for currents up to about $`1`$ mA, M64s depart from linearity already at about $`0.1`$ mA.

The photon-statistics method, described in Sec. 5.2.2, relies on the assumption that the PMT is a linear device. Assuming that the light source is stable, any deviation from linearity should have a direct effect on the quantities calculated in Eq. (2). This is what is observed in Fig. 17, where $`\mathrm{N}_{\mathrm{pe}}`$, the number of photoelectrons, is plotted as a function of HV. The spread in $`\mathrm{N}_{\mathrm{pe}}`$ at fixed HV reflects a general change in light level, when testing different sets of tested tubes. (We had to demount and remount the fibres in the attenuation box leading to a change in overall light level.) It appears that the calculated number of photoelectrons is a function of the applied voltage; but of course the light level cannot really depend on the voltage. To verify that this dependence is an artifact of the photon-statistics computation and not an effect correlated with the PMT quantum efficiency and collection efficiency, in Fig. 18 $`\mathrm{N}_{\mathrm{pe}}`$ is plotted as a function of QE at $`520`$ nm and of the calculated gain $`𝒢`$ separately: $`\mathrm{N}_{\mathrm{pe}}`$ seems to be independent of QE, while a clear dependence on $`𝒢`$ is observed. Simulations have confirmed that tube non-linearities can influence the calculations of $`\mathrm{N}_{\mathrm{pe}}`$ in the observed way when the light level is 30-50 photoelectrons as shown in Fig. 19

At light levels as low as 3-5 photoelectrons the PMT gain fluctuations (width of the single-photoelectron peak) become dominant and influence the calculation. This width varies with gain. The relative width of the single-photoelectron peak becomes narrower as the gain increases, leading to a higher calculated number of photoelectrons.

#### 5.3.4 Sensitivity to Magnetic Fields

To test the potential consequences of a stray magnetic field on the PMT performances, a Helmholtz coil was used to generate a homogeneous field of up to $`40`$ gauss. A full scan (a total of $`40\times 40`$ points, $`5\times 5`$ points per pixel) of several tubes was performed for each value of the magnetic field. No effect was observed for a field orientation perpendicular to the PMT longitudinal axis ($`z`$), either on the tube gain or on collection efficiency. However, if the field is parallel to the $`z`$-axis, significant effects on the collection efficiency of the outer pixels are observed at intensities of $`7`$ gauss and above. Cross-talk was observed not to be affected by the presence of a magnetic field, regardless of its orientation.

## 6 Conclusions

We have measured the properties of the Hamamatsu R5900-00-M64 multi-anode photomultiplier tube and find that these tubes are well suited for the MINOS experiment. We have found that the main properties of the tube are:

* The quantum efficiency at 520 nm is $`13\pm 0.5\%`$
* The average tube gain is 1.3 $`\times 10^6`$ at 850 V while the ratio of minimum to maximum pixel gain is better than 1:3.
* There are only small variations in the tube efficiency. The typical variation has been measured to be $`20\%`$.
* There are already sizable non-linearities at a tube gain of $`10^6`$.
* When placing a 1.2 mm acrylic fibre on the photocathode, the cross-talk to neighbouring pixels is on average $`10\%`$ which includes contributions from both optical and charge cross-talk
* Typical dark count rates at $`20^{}\text{C}`$ are less than 10 Hz per pixel for signals above 1/3 of a photoelectron.
* The tubes are to a large extent insensitive to magnetic fields of up to 5 gauss.

We would like to thank Hamamatsu for their cooperation and for providing extensive test data for their PMTs. We also thank all the technicians and engineers who helped in preparing these measurements. This work was supported in part by PPARC, the UK Particle Physics and Astronomy Research Council, and by DOE, the United States Department of Energy.


text:
# Systematics of Heavy Quark Production at RHIC

## 1 Introduction

It is important to have an accurate measure of the charm and bottom cross sections for several reasons. Heavy quark decays are expected to dominate the lepton pair continuum from the $`J/\psi (c\overline{c})`$ and $`\mathrm{{\rm Y}}(b\overline{b})`$ up to the mass of the $`Z^0`$ . Thus the Drell-Yan yield and any thermal dilepton production will essentially be hidden by the heavy quark decay contributions . The shape of the charm and bottom contributions to this continuum could be significantly altered by heavy quark energy loss . If the loss is large, it may be possible to extract a thermal dilepton yield if it cannot be determined by other means . Heavy quark production in a quark-gluon plasma has also been predicted . This additional yield can only be determined if the $`AA`$ rate can be accurately measured. Finally, the total charm rate would be a useful reference for $`J/\psi `$ production since enhancement of the $`J/\psi `$ to total charm ratio has been predicted in a number of models .

## 2 Baseline Rates in $`pp`$

We first discuss some new calculations of the $`Q\overline{Q}`$ total cross sections in $`pp`$ collisions with the most recent nucleon parton distribution functions. At leading order (LO) heavy quarks are produced by $`gg`$ fusion and $`q\overline{q}`$ annihilation while at next-to-leading order (NLO) $`qg`$ and $`\overline{q}g`$ scattering is also included. To any order, the partonic cross section may be expressed in terms of dimensionless scaling functions $`f_{ij}^{(k,l)}`$ that depend only on the variable $`\eta `$ ,

$`\widehat{\sigma }_{ij}(\widehat{s},m_Q^2,\mu ^2)={\displaystyle \frac{\alpha _s^2(\mu )}{m^2}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\left(4\pi \alpha _s(\mu )\right)^k{\displaystyle \underset{l=0}{\overset{k}{}}}f_{ij}^{(k,l)}(\eta )\mathrm{ln}^l\left({\displaystyle \frac{\mu ^2}{m_Q^2}}\right),`$ (1)
where $`\widehat{s}`$ is the partonic center of mass energy squared, $`m_Q`$ is the heavy quark mass, $`\mu `$ is the scale and $`\eta =\widehat{s}/4m_Q^21`$. The cross section is calculated as an expansion in powers of $`\alpha _s`$ with $`k=0`$ corresponding to the Born cross section at order $`𝒪(\alpha _s^2)`$. The first correction, $`k=1`$, corresponds to the NLO cross section at $`𝒪(\alpha _s^3)`$. It is only at this order and above that the dependence on renormalization scale, $`\mu _R`$, enters the calculation since when $`k=1`$ and $`l=1`$, the logarithm $`\mathrm{ln}(\mu ^2/m_Q^2)`$ appears. The dependence on the factorization scale, $`\mu _F`$, the argument of $`\alpha _s`$, appears already at LO. We assume that $`\mu _R=\mu _F=\mu `$. The next-to-next-to-leading order (NNLO) corrections to next-to-next-to-leading logarithm have been calculated near threshold  but the complete calculation only exists to NLO.

The total hadronic cross section is obtained by convoluting the total partonic cross section with the parton distribution functions (PDFs) of the initial hadrons,

$`\sigma _{pp}(s,m_Q^2)={\displaystyle \underset{i,j=q,\overline{q},g}{}}{\displaystyle _{\frac{4m_Q^2}{s}}^1}{\displaystyle \frac{d\tau }{\tau }}\delta (x_1x_2\tau )F_i^p(x_1,\mu ^2)F_j^p(x_2,\mu ^2)\widehat{\sigma }_{ij}(\tau ,m_Q^2,\mu ^2),`$ (2)
where the sum $`i`$ is over all massless partons and $`x_1`$ and $`x_2`$ are fractional momenta. The PDFs, denoted by $`F_i^p`$, are evaluated at scale $`\mu `$. All our calculations are fully NLO, applying NLO parton distribution functions and the two-loop $`\alpha _s`$ to both the $`𝒪(\alpha _s^2)`$ and $`𝒪(\alpha _s^3)`$ contributions, as is typically done .

To obtain the $`pp`$ cross sections at RHIC and LHC, we first compare the NLO cross sections to the available $`c\overline{c}`$ and $`b\overline{b}`$ production data by varying the mass, $`m_Q`$, and scale, $`\mu `$, to obtain the ‘best’ agreement with the data for several combinations of $`m_Q`$, $`\mu `$, and PDF. We use the recent MRST HO central gluon , CTEQ 5M , and GRV 98 HO  distributions. The results for the $`c\overline{c}`$ cross section in $`pp`$ interactions is shown in Fig. 1.

On the left-hand side, $`\mu =m_c`$ for $`1.2m_c1.8`$ GeV, while on the right-hand side, $`\mu =2m_c`$ for the same masses, all calculated with MRST HO. The scale is not decreased below $`m_c`$ because the minimum scale in the PDF is larger than $`m_c/2`$. The cross sections with $`\mu =m_c`$ are all larger than those with $`\mu =2m_c`$ for the same $`m_c`$ because $`\alpha _s(m_c)>\alpha _s(2m_c)`$ by virtue of the running of $`\alpha _s`$. Evolution of the PDFs with $`\mu `$ tends to go in the opposite direction. At higher scales the two effects tend to compensate and reduce the scale dependence but the charm quark mass is not large enough for this to occur.

The best agreement with $`\mu =m_c`$ is for $`m_c=1.4`$ GeV and $`m_c=1.2`$ GeV is the best choice for $`\mu =2m_c`$ for the MRST HO and CTEQ 5M distributions. The best agreement with GRV 98 HO is $`\mu =m_c=1.3`$ GeV while the results with $`\mu =2m_c`$ lie below the data for all $`m_c`$. All five results agree very well with each other for $`ppc\overline{c}`$, as shown on the left side of Fig. 2. There is more of a spread in the $`\pi ^{}pc\overline{c}`$ results, shown on the right side of Fig. 2.

This is because the $`\pi ^{}`$ PDFs are not very well known. The last evaluations, SMRS , Owens-$`\pi `$ , and GRV-$`\pi `$  were 10-15 years ago and do not reflect any of the latest information on the low $`x`$ behavior of the proton PDFs, e.g. the distributions are all flat as $`x0`$ with no low $`x`$ rise. These pion evaluations also depend on the behavior of the proton PDFs used in the original fit, including the value of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. Thus the pion and proton PDFs are generally incompatible. Note that the $`\pi ^{}p`$ cross sections are a bit lower than the data compared to the $`pp`$ cross sections, suggesting that lighter quark masses would tend to be favored for this data. The reason is because the low $`x`$ rise in the proton PDFs depletes the gluon density for $`x>0.02`$ relative to a constant at $`x0`$ for $`\mu =\mu _0`$ the initial scale of the PDF. The $`\pi ^{}p`$ data are in a relatively large $`x`$ region, $`0.1x=2\mu /\sqrt{s}0.3`$, where this difference is important.

We have tried to play the same game with the $`b\overline{b}`$ total cross sections but these have mostly been measured in $`\pi ^{}p`$ interactions. The typical $`x`$ values of $`b\overline{b}`$ production are even larger than those for $`c\overline{c}`$ but it is not clear that $`\pi ^{}pb\overline{b}`$ also favors lower masses. At the fixed target energies of $`b\overline{b}`$ production, $`q\overline{q}`$ annihilation dominates while $`gg`$ fusion is still most important for $`c\overline{c}`$ production . The valence-valence $`\overline{u}_\pi u_p`$ contribution is most important since valence distributions dominate at large $`x`$. For all three PDFs used, we find $`m_b=\mu =4.75`$ GeV, $`m_b=\mu /2=4.5`$ GeV, and $`m_b=2\mu =5`$ GeV most compatible with the sparse data. Attempts to measure the $`b\overline{b}`$ total cross section in fixed-target $`pp`$ interactions have been less successful. Hopefully the HERA-B experiment at DESY  will soon provide a new measurement.

Our calculations can then be extrapolated to RHIC and LHC energies. The result for $`c\overline{c}`$ is shown in Fig. 3.

Even though the cross sections agree within 30% at 40 GeV, by the Pb+Pb energy of the LHC they differ by a factor of 2.3. The spread in the $`b\overline{b}`$ cross sections is considerably smaller, $`2030`$% at the ion collider energies. Our results for $`pp`$ interactions at 40 GeV, 200 GeV, and 5.5 TeV are given in Table 1.

The $`AA`$ rates per event at $`b=0`$ with the same energies can be obtained by multiplying these cross sections by $`T_{AA}(b=0)`$, 29.3/mb for Au+Au and 30.4/mb for Pb+Pb. We find 8-13 $`c\overline{c}`$ pairs and $`0.05`$ $`b\overline{b}`$ pairs at RHIC with 97-225 $`c\overline{c}`$ pairs and $`5`$ $`b\overline{b}`$ pairs at LHC without nuclear shadowing. The shadowing effect is rather small for $`c\overline{c}`$ at RHIC and actually enhances the $`b\overline{b}`$ rate. The only important modification due to shadowing in the total cross section is on the $`c\overline{c}`$ rate at the LHC which is reduced to 67-150 pairs. Energy loss does not affect the total rate . As noted by Thews, this $`c\overline{c}`$ rate is large enough at RHIC for independent $`c`$ and $`\overline{c}`$ quarks to dynamically recombine to form $`J/\psi `$’s . The baseline rates of $`Q\overline{Q}`$ production are thus important for studying these effects.

## 3 Nuclear Gluon Distribution in $`pA`$

We now turn to a calculation of the nuclear gluon distribution in $`pA`$ interactions . We show that the dilepton continuum can be used to study nuclear shadowing and reproduces the input shadowing function well, in this case, the EKS98 parameterization . To simplify notation, we refer to generic heavy quarks, $`Q`$, and heavy-flavored mesons, $`H`$. The lepton pair production cross section is

$`{\displaystyle \frac{d\sigma ^{pAl\overline{l}+X}}{dM_{l\overline{l}}dy_{l\overline{l}}}}`$ $`=`$ $`{\displaystyle d^3\stackrel{}{p}_ld^3\stackrel{}{p}_{\overline{l}}d^3\stackrel{}{p}_Hd^3\stackrel{}{p}_{\overline{H}}\delta (M_{l\overline{l}}M(p_l,p_{\overline{l}}))\delta (y_{l\overline{l}}y(p_l,p_{\overline{l}}))}`$ (3)
$`\times {\displaystyle \frac{d\mathrm{\Gamma }^{Hl+X}(\stackrel{}{p}_H)}{d^3\stackrel{}{p}_l}}{\displaystyle \frac{d\mathrm{\Gamma }^{\overline{H}\overline{l}+X}(\stackrel{}{p}_{\overline{H}})}{d^3\stackrel{}{p}_{\overline{l}}}}{\displaystyle \frac{d\sigma ^{\mathrm{p}AH\overline{H}+X}}{d^3\stackrel{}{p}_Hd^3\stackrel{}{p}_{\overline{H}}}}`$
$`\times \theta (y_{\mathrm{min}}<y_l,y_{\overline{l}}<y_{\mathrm{max}})\theta (\varphi _{\mathrm{min}}<\varphi _l,\varphi _{\overline{l}}<\varphi _{\mathrm{max}})`$
where $`M(p_l,p_{\overline{l}})`$ and $`y(p_l,p_{\overline{l}})`$ are the invariant mass and rapidity of the $`l\overline{l}`$ pair. The decay rate, $`d\mathrm{\Gamma }^{Hl+X}(\stackrel{}{p}_H)/d^3\stackrel{}{p}_l`$, is the probability that meson $`H`$ with momentum $`\stackrel{}{p}_H`$ decays to a lepton $`l`$ with momentum $`\stackrel{}{p}_l`$. The $`\theta `$ functions define single lepton rapidity and azimuthal angle cuts used to simulate detector acceptances.

Using a fragmentation function $`D_Q^H`$ to describe quark fragmentation to mesons, the $`H\overline{H}`$ production cross section can be written as

$`{\displaystyle \frac{d\sigma ^{pAH\overline{H}+X}}{d^3\stackrel{}{p}_Hd^3\stackrel{}{p}_{\overline{H}}}}`$ $`=`$ $`{\displaystyle \frac{d^3\stackrel{}{p}_Q}{E_Q}\frac{d^3\stackrel{}{p}_{\overline{Q}}}{E_{\overline{Q}}}E_QE_{\overline{Q}}\frac{d\sigma ^{pAQ\overline{Q}+X}}{d^3\stackrel{}{p}_Qd^3\stackrel{}{p}_{\overline{Q}}}𝑑z_1𝑑z_2D_Q^H(z_1)D_{\overline{Q}}^{\overline{H}}(z_2)}`$ (4)
$`\times `$ $`\delta ^{(3)}(\stackrel{}{p}_Hz_1\stackrel{}{p}_Q)\delta ^{(3)}(\stackrel{}{p}_{\overline{H}}z_2\stackrel{}{p}_{\overline{Q}}).`$
Our calculations were done with two different fragmentation functions. We found that our results were independent of $`D_Q^H`$. The hadronic heavy quark production cross section per nucleon in $`pA`$ collisions can be factorized into the general form

$`{\displaystyle \frac{1}{A}}E_QE_{\overline{Q}}{\displaystyle \frac{d\sigma ^{pAQ\overline{Q}+X}}{d^3\stackrel{}{p}_Qd^3\stackrel{}{p}_{\overline{Q}}}}={\displaystyle \underset{i,j}{}}{\displaystyle 𝑑x_1𝑑x_2f_i^p(x_1,\mu ^2)f_j^A(x_2,\mu ^2)E_QE_{\overline{Q}}\frac{d\widehat{\sigma }^{ijQ\overline{Q}}}{d^3\stackrel{}{p}_Qd^3\stackrel{}{p}_{\overline{Q}}}}`$ (5)
where $`f_i^p=F_i^p/x`$ and $`f_i^A=F_i^A/x`$ with $`F_i^A=F_i^pR_i^A`$. The shadowing ratio $`R_i^A`$ is that of EKS98 . The partonic cross section is the differential of Eq. (1) at $`k=0`$. Note that the total lepton pair production cross section is equal to the total $`Q\overline{Q}`$ cross section multiplied by the square of the lepton branching ratio.

We compare the ratios of lepton pair cross sections with the input $`R_g^A`$ in Fig. 4. All the results are integrated over the rapidity intervals appropriate to the PHENIX and ALICE dilepton coverages. The ratio follows $`R_g^A`$ at all energies. The higher the energy, the better the agreement: at the LHC the two agree very well.

The ratio always lies below $`R_g^A`$ for two reasons. First, $`q\overline{q}`$ annihilation is included and quark shadowing is different than gluon shadowing. The $`q\overline{q}`$ contribution decreases with energy, leading to better agreement at the LHC. Second, the phase space integration smears the shadowing effect relative to $`R_g^A(x_2,\mu )`$. Note that the ratio deviates slightly more from $`R_g^A(x_2,\mu )`$ for $`e^+e^{}`$ than for $`\mu ^+\mu ^{}`$ because the curvature of $`R_g^A`$ with $`x`$ is stronger at larger values of $`x`$ and, due to the differences in rapidity coverage, the average values of $`x_2`$ are larger for $`e^+e^{}`$.

The average $`x_2`$ decreases with energy. We have $`0.14x_20.32`$ at the SPS where $`R_g^A`$ is decreasing. At RHIC, $`0.003x_20.012`$, where $`R_g^A`$ is increasing quite rapidly. Finally, at the LHC, $`3\times 10^5x_22\times 10^4`$ where $`R_g^A`$ is almost independent of $`x`$. The values of $`x_2`$ are typically larger for electron pairs at collider energies because the electron coverage is more central than the muon coverage.

The average $`\mu ^2`$ increases with energy and quark mass. For $`c\overline{c}`$ we have $`7.58\mu ^248.5`$ GeV2 at the SPS, $`9.46\mu ^2141`$ GeV2 at RHIC, and $`11.4\mu ^2577`$ GeV2 at the LHC. For $`b\overline{b}`$ production, $`32.0\mu ^254.3`$ GeV2 at RHIC and $`37.9\mu ^2156`$ GeV2 at LHC.

## 4 Heavy Quarks in $`AA`$

### 4.1 Effects of Energy Loss

Energy loss would best be determined by reconstruction of $`D`$ and $`B`$ meson decays and comparing with distributions expected from $`pp`$ and $`pA`$ extrapolations that do not consider energy loss. Whether energy loss is measurable in reconstructed $`D`$ and $`B`$ decays or not, the change in the dilepton continuum should surely be present if the loss is nonzero and will bias the interpretation of the dilepton continuum. So far, the amount of the energy lost by heavy quarks is unknown. While a number of calculations have been made of the collisional loss in a quark-gluon plasma , only recently has radiative loss been applied to heavy quarks . The radiative loss can be rather large, $`dE/dx5`$ GeV/fm for a 10 GeV heavy quark, and increasing with energy, but the collisional loss is smaller, $`dE/dx12`$ GeV/fm, and nearly independent of energy . We note that energy loss will suppress high $`p_T`$ and large invariant mass quark pairs as long as $`|dE/dx|p_T/R_A`$ .

It is important to note that energy loss does not reduce the number of $`Q\overline{Q}`$ pairs produced but only changes their momentum. However, an effective reduction in the observed heavy quark yield can be expected in a finite acceptance detector because fewer leptons from the subsequent decays of the heavy quarks will pass kinematic cuts such as a minimum lepton $`p_T`$.

If the loss or the $`p_T`$ cut is large, the Drell-Yan and thermal dileptons could emerge from under the reduced $`D\overline{D}`$ and $`B\overline{B}`$ decay contributions at large masses. Even without considering energy loss, Gallmeister et al. suggested that thermal dileptons could be detected by increasing the minimum lepton $`p_T`$ because, in the $`D`$ and $`B`$ rest frames, the maximum energy of the individual leptons is limited to 0.9 and 2.2 GeV respectively. The lepton $`p_T`$ from thermal production has no such limitation .

### 4.2 Quarkonium normalization

Heavy quark production in $`AA`$ collisions is also interesting because of the prominent effect it could have on quarkonium. Initial nucleon-nucleon collisions may not be the only source of quarkonium production. Regeneration of quarkonium in the plasma phase  could counter the effects of suppression, ultimately leading to enhanced quarkonium production. In the plasma phase, there are two basic approaches: statistical and dynamical coalescence. Both these approaches depend on being able to measure the quarkonium rate relative to total $`Q\overline{Q}`$ production. Thus the $`Q\overline{Q}`$ rate is preferable as a normalization of quarkonium production, particularly since both share the same production mechanisms and approximate $`x`$, $`\mu `$ values. However the final-state effects such as energy loss will make the total rate difficult to quantify without substantial detailed studies. These secondary production models should be testable already at RHIC where enhancements of factors of 2-3 are expected from coalescence .

Other processes besides heavy quark production have been suggested as references for quarkonium production. Using the $`Z^0`$ as a reference  as a reference would eliminate the uncertainty due to final-state effects on the $`Z^0l^+l^{}`$ decays but the different production mechanisms and masses leaves it less desirable. It has also been suggested that the $`\psi ^{}/J/\psi `$ and $`\mathrm{{\rm Y}}^{}/\mathrm{{\rm Y}}`$ ratios be studied as a function of $`p_T`$  since deviations from the $`pp`$ ratios should reflect quark-gluon plasma characteristics. The only drawback to such a mechanism is that strong suppression may result in poor statistics.

## Acknowledgements

I would like to thank K.J. Eskola and V.J. Kolhinen for collaboration on part of this work. This work was supported by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Nuclear Physics Division of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.