text:
# Proper Motion Measurements with the VLA: I. Wide-field Imaging and Pulse Gating Techniques

## 1 Introduction

Interferometric proper motions have always been hindered by the low flux densities of many pulsars. This restriction results in a bias in which only the closest pulsars are observed. Although these pulsars move quickly on angular scales, their absolute motion must be small for them to still be near the plane of the Galaxy. The fastest pulsars, however, have traveled the greatest distances and may be more than a kpc above the Galactic plane. Previous interferometric proper motion projects have been unable to include these distant pulsars due to their low flux densities. The study by Fomalont et al., begun in 1984, was limited to pulsars with $`S>2`$ mJy so that an adequate signal to noise could be achieved (Fomalont et al., 1992). Few faint pulsars have proper motion measurements as a result of these selection effects.

This study of proper motions of distant, high $`z`$ pulsars using the Very Large Array (VLA) of the National Radio Astronomy Observatory111The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. began in 1992. We have chosen a group of 28 pulsars, 15 of which have $`|z|>`$ 400 pc, $`d4000`$ pc, and $`|b|<30^{}`$. This selection criteria minimizes the uncertainty in the determination of the proper motion in the $`z`$ direction which enables the study of proper motions of pulsars relative to the Galactic plane. Upon completion, this project will double the number of proper motions measured for distant, high $`z`$ pulsars. The flux density restriction is overcome by gating the VLA. In order to have the pulsar in the same field as stationary sources, we must create wide-field images with position accuracies of a few mas. This requirement forces us to account for small effects such as annual aberration and Lorentz contraction over the fields.

This paper focuses on the method of data reduction and error analysis used to calculate the proper motions. Three pulsars (B0919+06, B1237+25 and B1937–26) are used as “case studies” demonstrating the accuracy of the method over a range of declinations and flux densities. Proper motions for all 28 pulsars will be presented by Brisken et al. in Paper II.

## 2 Data and Observational Techniques

Five epochs of data centered on 1992.96, 1994.23, 1995.52, 1998.32 and 1999.47 have been gathered at the VLA, amounting to more than 120 hours of observations. Every pulsar has been observed twice; 25 pulsars have three or more epochs of data.

### 2.1 Observational Setup

Observations were made in the A array at the VLA with the correlator in the 2 AD mode. A maximum baseline of 36 km provides the highest possible resolution (1.1′′). With this resolution, large sources are resolved out leaving mostly distant, extra-galactic sources in the images. These sources should be stationary over time and are used as reference sources against which to measure the motion of the pulsar.

The data were gathered at 20 cm (1452.4 MHz) as a compromise between the reduction in pulsar flux density at higher frequencies and decreased resolution at lower frequencies. A 25 MHz bandwidth was divided into fifteen channels ($`\mathrm{\Delta }\nu =1.56`$ MHz) to provide the large usable field of view ($`36^{}`$) necessary for wide-field imaging (see section 2.2). Observations were made in two circular polarizations to provide the gating capability for pulsars without losing the ability to image weak continuum sources (see section 2.4). Thirteen channels were used to produce final images (see section 3).

Each pulsar was typically observed four times for 15 minutes during an observing run. Each segment was bracketed by five minutes of observation of a bright nearby VLA calibrator that was subsequently used for phase calibration and secondary amplitude calibration.

### 2.2 Delay Beam $`\&`$ Bandwidth Smearing

The size of the primary beam of the VLA is given by $`\theta _p\lambda /d`$ where $`d=25`$ m is the diameter of the dish. Thus $`\theta _p30^{}`$ at 20 cm defines the field of view contained inside the first null of the primary beam. In practice, observations are made over a band of non-zero width ($`\mathrm{\Delta }\nu 0`$). The range of frequencies in the band results in bandwidth smearing in which sources far from field center are smeared in the radial direction.

Bandwidth smearing is characterized by the delay beam, which is defined as $`\theta _d=2c/\mathrm{\Delta }\nu D`$ where $`D`$ is the maximum baseline (36 km for the A array). To maximize the usable field of view, observations are set up so that $`\theta _d\theta _p`$. Since a single continuum channel ($`\mathrm{\Delta }\nu =25`$ MHz) has a usable field of view of $``$2, observations must be made in the spectral line mode with narrow channels. By dividing the band into fifteen channels ($`\mathrm{\Delta }\nu =1.56`$ MHz) a usable field of view of $``$36 is attained. Bandwidth smearing is only an important effect for sources more than $``$18 from field center. The added uncertainty in position of reference sources more than 20 from field center is taken into account in the proper motion fitting program (see section 6)

### 2.3 Gating

The VLA correlator was gated using the Princeton Mark III Pulsar Timing Machine (Stinebring et al., 1992). This computer allows real-time adjustment of the pulsar gate enabling the correlator to record data only when the pulsar is “on.” All the data during the off-pulse are discarded and only on-pulse portions are retained.

Gating can increase the signal-to-noise ratio (SNR) of a pulsar by up to a factor of five. The gate is matched to the duty cycle of the pulsar to maximize the improvement in the SNR. This improvement is approximately proportional to the inverse square-root of the duty cycle. Because the VLA has only a single gating circuit for all channels, the width of the gate must account for dispersion smearing across the full 25 MHz bandwidth and typically ranges between 5 and 15% of the period for pulsars in this study. For B0919+06, gating increased the flux density from 6 to 120 mJy while the noise went from 0.1 to 1.5 mJy. B1237+25 showed an increase in flux density of 5 to 55 mJy and an increase in noise of 0.2 to 0.7 mJy. Finally, the flux density of B1937–26 increased from 1 to 15 mJy while the noise increased from 0.2 to 0.6 mJy when the VLA correlator was gated. Therefore, gating increased the SNR of B0919+06 by a factor of $``$1.3 while B1235+25 and B1937–26 both showed an increase in SNR of a factor of $``$4. In general, gating is not applied to pulsars that are brighter than typical reference sources ($`10`$ mJy). However, technique testing (such as LL, RR alignment) was done using gated, bright pulsars while hardware availability prevented gating on some weak pulsars.

For our observations, the VLA is gated in the right circular polarization (RR) in order to increase the SNR of the pulsar. The left circular polarization (LL) is not gated, to permit optimal detection of reference sources in the field. Since the position of the pulsar is compared to those of the reference sources, positions in the right and left polarizations must agree to a few mas. Positions of strong pulsars and reference sources show a typical agreement between polarizations of 3 to 5 mas.

## 3 Calibration

The processing of the data begins with the removal of discrepant data points using the procedure TVFLG in $`AIPS`$ (Astronomical Image Processing System). During the first observation of a pulsar, the first 3 to 5 minutes are typically used to determine the optimum position of the gate. Since the ungated data have significantly lower SNR than the gated data, they must be removed from the data set.

Calibrators for each pulsar are positioned at the phase center and are unaffected by bandwidth smearing. Therefore, basic calibrations are done in the continuum data set, which is a sum of the inner 75% of the observing band, and the resulting tables are copied to the spectral line data. 3C286 (J1331+305), 3C48 (J0137+331) and 3C147 (J0542+498) were used as the flux density and bandpass calibrators for the observations. Channels on the edge of the band pass filter are degraded by a roll-off in the band which results in reduced sensitivity in these channels. The first and last channels were removed from the data leaving thirteen channels to be used for final images. Calibration and editing information was applied when splitting the uv data into separate data sets for each pulsar.

For the clean algorithm to converge properly, the brightness and structure of all sources must remain constant, in both time and frequency, over the entire integration. The regular amplitude calibration and removal of data taken without the gate ensure that the time constancy is achieved, and the bandpass calibration flattens the spectrum. Although pulsars have considerably steeper spectra than do most of the reference sources and calibrators, their mean flux density changes by less than 5% over our band and produces no noticeable cleaning artifacts. Scintillation in frequency and time violate the source constancy that is assumed in the clean algorithm, resulting in cleaning artifacts and an increase in the background noise. Since the noise near each source is used to determine the uncertainty of its position fit, this source of uncertainty is automatically included in the analysis.

### 3.1 UVFIX

Before imaging, UVFIX is run on the single source data sets. This $`AIPS`$ program recovers the correct $`u,v`$ and $`w`$ coordinates of a source that are only approximated by the VLA correlator. The omission of these terms results in a movement in the tangential direction of up to 60 mas for a source 10 from the field center (Fomalont et al., 1992). UVFIX also corrects for Lorentz contraction of the field due to the motion of the Earth along the line of sight. The magnitude of this effect depends on the time of year in which the observations were made and can result in a maximum radial movement of 60 mas for a source 10 from field center (Fomalont et al., 1992).

L. Kogan at NRAO-Socorro revised UVFIX to account for these two effects by recalculating the correct $`u,v`$ and $`w`$ coordinates as well as moving to a stationary frame relative to the Sun. The results of the updated UVFIX program were tested by comparing reference source position agreement between epochs before and after UVFIX was run. As seen in Figure Proper Motion Measurements with the VLA: I. Wide-field Imaging and Pulse Gating Techniques, UVFIX removed the rotation in the field near B0919+06 improving reference source alignment by a factor of three. This new version of UVFIX was incorporated into AIPS in early 1998.

## 4  Imaging Techniques

Images are made using the Clark “clean” algorithm in $`AIPS`$. A pixel size of $`0.15^{\prime \prime }`$ provides more than 5 pixels/beam. Clean boxes are placed around regions containing sources to ensure that sidelobes are removed from these areas. Separate clean boxes are used for the pulsar and each reference source. A robust weighting of zero, which indicates a compromise between natural and uniform weighting, is used for all imaging.

### 4.1 Wide-field Imaging

A wide-field image ($`60^{}\times 60^{}`$) is made in the LL polarization with ungated data for the detection of reference source candidates. Positions of bright sources are recorded and new images are made with separate fields for each source. Ideally, only point sources are retained as reference sources. In cases where there are only a few point sources we choose slightly extended sources as additional reference sources. The larger error in the position estimates of these extended sources results in additional uncertainty in the final proper motion. All sources, including extended and weak sources are included in the imaging process to remove sidelobes from other fields. Due to the wide field used to image the reference sources and the non-coplanarity of the VLA, the uv coverage at different points within the field of view can differ by an appreciable amount (Perley, 1999). To account for this, the pulsar and all reference sources are imaged using uv values appropriate for their part of the sky.

Multiple reference sources have been found for all 28 pulsars. There are typically about eight good reference sources for each pulsar. The quality of each reference source is determined by measuring its proper motion relative to the other reference sources. Extended reference sources are generally not used and point-like reference sources with large motions are also omitted. Sources used as reference sources for the final proper motion calculations have their proper motions listed in Table 1. For the three pulsars considered in this paper, the maximum distance of any detected source from field center is 27.6. However, no source further than 22.3 was used for the proper motion calculations.

### 4.2 Self-calibration

If there is enough flux density in compact sources inside the inner third of the primary beam then self-calibration can improve the signal-to-noise ratio of the detections. In self-calibration, one assumes that the image is degraded by antenna based gain and phase errors which vary in time and prevent perfect calibration. By re-calibrating the data using an initial set of images as a model, corrections to these phase and amplitude errors can be calculated as a function of time.

An accurate model of the flux density distribution across the field is necessary for self-calibration to be successful. Self-calibration is performed in the left hand polarization since the left hand data typically contain more flux density than the gated right hand polarization. In addition, pulsar scintillation is much more prominent in the gated polarization and could introduce additional amplitude errors in the calibration. Solutions from the left hand data are applied to both polarizations. Self-calibration can result in a shift in the field of as much as 20 mas when the initial model is not complete. Since the pulsar and reference sources are equally affected by this shift, the proper motion calculations are not compromised. Self calibration and imaging loops were iterated between 2 and 4 times. Solution intervals of between two and five minutes were used.

## 5 Position Determination

Positions of the pulsar and reference sources are found from final images using the AIPS task JMFIT. This task fits a gaussian profile to a point source with a width based on the size of the synthesized beam. For extended sources, the program solves for the width of the source as well.

JMFIT also reports an uncertainty in the position estimate. Since this uncertainty directly affects our confidence in the proper motions, tests were run on the output of JMFIT to ensure that it gives a reasonable error estimate. A continuum point source in 41 line-free channels of an HI image of NGC4688 was used for this test. These data were chosen because a bright point source and a large number of channels were available. Using identical parameters, JMFIT was run on each channel individually and the position and error output were recorded. The uncertainties reported by JMFIT agreed with the standard deviation of the 41 positions reported for the point source.

## 6 Proper Motions

The pulsar positions and proper motions are determined through a global least squares fit. Positions of reference sources and shifts in the coordinate system between epochs are also fit. The beam shape is used with the uncertainty reported by JMFIT to produce a correctly oriented elliptical gaussian uncertainty. An additional uncertainty of $`FR\mathrm{\Delta }\nu /\nu `$ is added in quadrature in the radial direction to account for bandwidth smearing where $`R`$ is the source’s distance from the phase center, $`\mathrm{\Delta }\nu /\nu `$ is the single channel fractional bandwidth, which in our case is close to 0.001, and $`F`$ is an empirically determined constant. By measuring the scatter in position measurements of point-like reference sources far from the field center in the same way that JMFIT was tested, we have determined that $`F=0.08`$ for our data.

A small systematic coordinate offset remains for each epoch, even after UVFIX is applied. This effect ($``$30 mas for a source 10 from field center) is approximately half as large as the correction made by UVFIX and cannot be adequately modeled by a simple rotation and dilation. Although the source of this offset is not fully understood, it is effectively removed by fitting for a six parameter general linear transformation, $`\stackrel{}{X}^{}=𝐀\stackrel{}{X}+\stackrel{}{B}`$, between the coordinate systems of each epoch and the first epoch. $`𝐀`$ is a two by two matrix whose elements typically deviate from the identity matrix by a few times $`10^5`$ and $`\stackrel{}{B}`$ is a coordinate frame shift between the two epochs. The amount of correction can be characterized by a dimensionless number equal to the RMS value of the matrix elements. If we define $`N_{epochs}`$ as the number of epochs in which the pulsar was observed and $`\delta 𝐀_{ij}=𝐀_{ij}\delta _{ij}`$, where $`\delta _{ij}`$ is the identity matrix, then

$$\mathrm{RMS}=\sqrt{\frac{1}{4(N_{epochs}1)}\underset{e=2}{\overset{N_{epochs}}{}}(\delta 𝐀_{e;11}^2+\delta 𝐀_{e;12}^2+\delta 𝐀_{e;21}^2+\delta 𝐀_{e;22}^2)}.$$
(1)
The RMS correction for 28 pulsars is plotted against source declination in Figure Proper Motion Measurements with the VLA: I. Wide-field Imaging and Pulse Gating Techniques. Since most of the observations were made near transit the elevation is approximately $`90\left|declat\right|`$, where $`lat`$ is the lattitude of the VLA, about $`+34^{}`$. The required correction seems to follow the secant of the elevation suggesting that the atmosphere may cause this effect.

A Monte-Carlo bootstrap method was used to better determine the final error ellipses. In this test, data are randomly resubstituted and the solution is fit thousands of times. The solutions for $`\mu _\alpha `$ and $`\mu _\delta `$ are plotted in a scatter plot. The uncertainties are estimated from the width of the distribution with an orientation based on the shape of the distribution in right ascension and declination. The uncertainties obtained in this manner agree well with the least squares errors. For a complete description of this method, see Numerical Recipes in C (Press et al., 1992).

## 7 Case Studies

B0919+06, B1237+25 and B1937–26 demonstrate the proper motion calculation as well as the accuracy of the error analysis. B0919+06 and B1237+25 were chosen because they both have previously published proper motions and can be used to test our results against those obtained from other methods. B0919+06 is especially interesting because it has a VLBA determination of the proper motion accurate to $`<1`$ mas yr-1 (Chatterjee et al., in preparation). B1937–26 is included because it is a weak pulsar ($``$1 mJy) and is also located at a low declination. At –26, the synthesized beam of the VLA is 2.5′′ in declination compared to only 1.1′′ in right ascension. This elongated beam makes accurate measurement of $`\mu _\delta `$ difficult. In addition, calibrations for low elevation sources depend strongly on pointing direction. Depending on the distance from pulsar to phase calibrator, this can make self-calibration difficult. The calculated proper motions for all three pulsars can be seen in Table 2.

### 7.1 B0919+06

B0919+06 provides the most stringent test of the accuracy of our method. We initially compared our measurement of the proper motion ($`\mu _\alpha =`$18.8 $`\pm `$ 0.9 mas yr-1, $`\mu _\delta =`$86.4 $`\pm `$ 0.7 mas yr-1) to that published by Fomalont et al. (1999) ($`\mu _\alpha =`$17.7 $`\pm `$ 0.3 mas yr-1, $`\mu _\delta =`$79.2 $`\pm `$ 0.5 mas yr-1). There was an obvious disagreement between the two results, especially in $`\mu _\delta `$. A recent analysis of more extensive VLBA data (Fomalont et al. data plus additional VLBA observations from October, 1998) by Chatterjee, Fomalont et al. (in preparation) has led to a revision of the initial VLBA result. The revised VLBA values are: $`\mu _\alpha =`$18.4 $`\pm `$ 0.2 mas yr-1, $`\mu _\delta =`$86.7 $`\pm `$ 0.3 mas yr-1 with a parallax of $`\pi =`$1.15 $`\pm `$ 0.25 mas (Chatterjee et al., in preparation). We agree with this result to less than 1$`\sigma `$ in $`\mu _\alpha `$ and $`\mu _\delta `$. It is important to note that the errors from our new technique are comparable to those of the updated VLBA observations. The new measurement also agrees to within 1$`\sigma `$ of the proper motion measurement by Harrison, Lyne & Anderson (1993) of $`\mu _\alpha =`$13 $`\pm `$ 29 mas yr-1, $`\mu _\delta =`$64 $`\pm `$ 37 mas yr-1. Figure Proper Motion Measurements with the VLA: I. Wide-field Imaging and Pulse Gating Techniquesa shows the new proper motion measurement along with the Chatterjee et al. (in preparation) and Harrison, Lyne & Anderson (1993) results.

### 7.2 B1237+25

The proper motion measured for B1237+25 ($`\mu _\alpha =`$–104.5 $`\pm `$ 1.1 mas yr-1, $`\mu _\delta =`$49.4 $`\pm `$ 1.4 mas yr-1) also agrees with both previous measurements (see Figure Proper Motion Measurements with the VLA: I. Wide-field Imaging and Pulse Gating Techniquesb). Our measurement deviates by $`<`$1$`\sigma `$ from the Fomalont et al. (1992) result of $`\mu _\alpha =`$–113 $`\pm `$ 13 mas yr-1, $`\mu _\delta =`$43 $`\pm `$ 14 mas yr-1. It also agrees reasonably well with the Lyne, Anderson & Salter (1982) measurement of $`\mu _\alpha =`$–106 $`\pm `$ 4 mas yr-1, $`\mu _\delta =`$42 $`\pm `$ 3 mas yr-1.

### 7.3 B1937–26

A proper motion of $`\mu _\alpha =`$12.1 $`\pm `$ 2.4 mas yr-1, $`\mu _\delta =`$–9.9 $`\pm `$ 3.8 mas yr-1 was measured for B1937–26 (see Figure Proper Motion Measurements with the VLA: I. Wide-field Imaging and Pulse Gating Techniquesc). Due to the low flux density of this pulsar ($``$1 mJy), B1937–26 has no previous measurement of its proper motion. Despite the elongated beam at low declination, the value for $`\mu _\delta `$ has a small error. The motion of the pulsar is significant when compared to the small motion of its reference sources (see Table 1).

Without gating, we would not be able to determine the proper motion of this pulsar to such high accuracy. The proper motion of B1937-26 was recalculated using only the ungated, LL polarization for the positions of the pulsar and the reference sources. The resulting proper motion was $`\mu _\alpha =`$20 $`\pm `$ 8 mas yr-1, $`\mu _\delta =`$–15 $`\pm `$ 15 mas yr-1. The increase in error due to the low signal to noise of the pulsar makes the measurement much less significant than the gated case.

The successful proper motion measurements for these three pulsars show the accuracy of this technique. B1937-26 is also limited by low flux density. A precise proper motion measurement would have been impossible without gating. These results confirm our ability to produce proper motions of weak pulsars that are accurate to less than five mas yr-1.

## 8 Conclusion

We have developed a new technique to calculate proper motions of weak pulsars at the VLA. By gating the VLA in one polarization, distant, high $`z`$ pulsars with flux densities as low as $``$1 mJy can be detected with a high signal-to-noise ratio. Previous work by Lyne, Anderson & Salter (1982) and Harrison, Lyne & Anderson (1993) also utilized the advantages of in-beam reference sources and gating. These observations were limited by a small number of baselines (one baseline in Lyne, Anderson & Salter (1982) and two baselines in Harrison, Lyne & Anderson (1993)) and the positions determined using fringe rate mapping rather than imaging could not correct for extended structure in reference sources. In addition, their data was degraded by the ionosphere which is worse at 408 MHz than in our observations at 1452.4 MHz. The 351 baselines and larger total collecting area of the VLA enable detection of much weaker sources and improved imaging. The addition of wide field imaging and a detailed understanding of systematic effects at the VLA also makes this technique more useful. Wide-field imaging enables sources located up to $`20^{}`$ from the field center to be used as reference sources. By imaging the reference sources, we can exclude sources which show extended structure. These sources are very important in the calculation of plate solutions as well as the pulsars’ motions. Proper motions of the pulsars are fit simultaneously to all reference sources in the field.

Both B0919+06 and B1237+25 agree well with previously published proper motions confirming the accuracy of this technique. The measurement of the proper motion of B1937–26 shows the real success of the method. Ungated, this pulsar has a flux density of only $``$1 mJy at 20 cm. This weak flux density excluded it from all previous proper motion projects. By gating the VLA, proper motions of faint pulsars can now be obtained with accuracies of a few mas yr-1 in less than ten years. These accuracies are comparable to those obtained with VLBA observations. Although the larger synthesized beam of the VLA reduces positional accuracies, the VLA has the advantage of more straightforward data reduction and more thoroughly understood systematic effects. The presence of reference sources in the field of view of the VLA images also makes this a favorable technique since VLBA observations generally require out-of-beam reference sources. Therefore, proper motions of faint pulsars can now be obtained at the VLA in just a few years. Proper motions for all 28 pulsars in this study will be reported in Paper II by Brisken et al.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research was funded in part by the National Science Foundation including the REU program at NRAO-Socorro and a graduate student fellowship. Additional funding was provided by Caltech and Princeton Universities. S. E. Thorsett is a Sloan Research Fellow. We would like to thank Leonia Kogan for help with the new version of UVFIX, Liese van Zee for help in acquiring the HI data and Ed Fomalont for advice.

CAPTIONS


format: text

text:
# Gravity-driven instability in a spherical Hele-Shaw cell

## I Introduction

The Saffman-Taylor instability  has been the object of extensive study during the last four decades . It arises at the interface separating two viscous fluids constrained to flow in the narrow gap between closed spaced parallel plates, a device known as Hele-Shaw cell. The cell thickness is smaller than any other length scale in the problem, so the flow is effectively two-dimensional. The instability arises either from a pressure gradient advancing the less viscous fluid against the more viscous one, or by gravity acting on the density difference between the fluids. The action of such driving-force mechanisms leads to the celebrated viscous fingering patterns .

Most Saffman-Taylor investigations analyze flow between flat Hele-Shaw cells. In a previous work  we started studying the Saffman-Taylor problem on curved surfaces by considering flow in a spherical Hele-Shaw cell (figure 1). The interfacial instability was produced by a nonzero flow injection rate $`Q`$, and gravitational effects were completely negleted. We examined the effect of cell curvature on the shape of the patterns and showed that positive spatial curvature inhibits finger tip-splitting.

In the present paper we focus on the influence of gravity in a spherical Hele-Shaw cell. We consider both $`Q=0`$ and $`Q>0`$ at fixed cell curvature. The unperturbed domain shape is a polar cap of some size, presenting an initially nearly circular boundary, which is gravitationally unstable and evolves at constant area ($`Q=0`$), or slowly increasing area ($`Q>0`$), without change of topology.

The study of viscous flow in a nonplanar Hele-Shaw cell is of interest for both scientific and practical reasons. On the scientific level, the influence of spatial curvature on hydrodynamic flow is a matter of fundamental interest. It also provides a simple mathematical model to describe more general situations involving the filling of a thin cavity between two walls of a given shape with fluid. On the practical level, it may have applications in a number of industrial, manufacturing processes, ranging through pressure moulding of molten metals and polymer materials , and formation of coating defects in drying paint thin films .

The gravity-driven instability on a spherical Hele-Shaw cell also allows one to gain insight into the properties of the dynamically similar, but more complex, Rayleigh-Taylor instability  for flow on substrates of arbitrary shapes. One familiar example of this type of instability is the formation of fingering patterns when chocolate syrup drains, due to the action of gravity, from the top of a scoop of ice cream. Despite the apparent simplicity of this example, it is in fact a rather complicated three-dimensional problem , much less amenable to analytic treatment than its two-dimensional Hele-Shaw counterpart.

The geometrically constrained spherical Hele-Shaw cell forces the flow to become essentially two-dimensional, and the interface one-dimensional. High viscosity eliminates inertial terms from the equations of motion, making the problem simpler yet. In contrast, the conventional Rayleigh-Taylor problem is inertially driven and three-dimensional effects become important .

The outline of the work is the following: section II derives the nonlinear equation of motion including both injection and gravitational driving. Section III discusses the resulting motion. Section III A considers linear stability analysis for purely gravitational driving with $`Q=0`$. The growth is purely exponential, and the linear growth rate depends on the geodesic distance from the sphere’s north pole. Larger distance causes both faster growth and more unstable modes, though the characteristic wavelength is nearly constant. The nonvanishing injection case $`Q>0`$ is studied in section III B. Injection of a high viscosity fluid tends to stabilize the interface. For $`Q>0`$ linear growth is non-exponential due to evolution of the linear growth rate with distance. Section III C studies the coupling of a small number of modes. It is shown that for gravity-driven flow fingers have tendency to finger tip-sharpening.

## II Stability analysis and mode coupling

Consider two immiscible, incompressible, viscous fluids, flowing in a spherical Hele-Shaw cell of thickness $`b`$ (see figure 1). The effetively two-dimensional flow takes place on the surface of a sphere endowed with the metric

$$ds^2=d\rho ^2+a^2\mathrm{sin}^2\left(\frac{\rho }{a}\right)d\phi ^2,$$
(1)
where $`a`$ is the radius of curvature of the sphere, $`0\phi <2\pi `$ denotes the azimuthal angle measured on the sphere and $`0\rho \pi a`$ is the geodesic distance from the sphere’s north pole. The polar angle $`\theta =\rho /a`$. Denote the viscosities and densities of the upper and lower fluids, respectively as $`\eta _1`$, $`\varrho _1`$ and $`\eta _2`$, $`\varrho _2`$. Consider the case $`\varrho _1>\varrho _2`$ and $`\eta _1>\eta _2`$, and examine flow in the northern hemisphere $`\theta <\pi /2`$. Between the two fluids there exists a surface tension $`\sigma `$. The flows are assumed to be irrotational, except at the interface. Fluid 1 is injected into fluid 2 through an inlet located at the sphere’s north pole, at a given flow rate $`Q`$, which is the area covered per unit time. Fluid 2 is simultaneously withdrawn, at the same rate, through an outlet placed at the south pole. The acceleration of gravity is constant, represented by g, and points from north to south pole.

During the flow, the fluid-fluid interface has a perturbed shape described as $`\rho =R(t)+\zeta (\phi ,t)`$. The interface perturbation amplitude is represented by $`\zeta (\phi ,t)`$, and $`R`$ denotes the time-dependent unperturbed radius

$$R(t)=a\mathrm{arccos}\left(C_0\frac{Qt}{2\pi a^2}\right),$$
(2)
where $`C_0=\mathrm{cos}(R_0/a)`$, and $`R_0`$ is the unperturbed radius at $`t=0`$. The unperturbed shape is a polar cap of geodesic radius $`\rho =R`$, surface area $`𝒜=4\pi a^2\mathrm{sin}^2(R/2a)`$ and circumference $`=2\pi a\mathrm{sin}(R/a)`$.

We express the net perturbation $`\zeta (\phi ,t)`$ as a Fourier series

$$\zeta (\phi ,t)=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\zeta _n(t)\mathrm{exp}(in\phi ),$$
(3)
where $`\zeta _n(t)`$ denotes the complex Fourier mode amplitudes, and $`n`$=0, $`\pm 1`$, $`\pm 2`$, $`\mathrm{}`$ is the discrete azimuthal wave number. The area of the perturbed shape is kept independent of the perturbation $`\zeta `$ by expressing the zeroth Fourier mode in (3) as $`\zeta _0(t)=(1/2a)\mathrm{cot}(R/a)\underset{n0}{}|\zeta _n(t)|^2`$.

We consider a generalized version of Darcy’s law , adjusted to describe flow between concentric spheres

$$𝐯_j=\frac{b^2}{12\eta _j}\left[p_j\varrho _jg\mathrm{sin}\left(\frac{\rho }{a}\right)\widehat{\rho }\right],$$
(4)
where $`𝐯_j=𝐯_j(\rho ,\phi )`$ and $`p_j=p_j(\rho ,\phi )`$ are, respectively, the velocity and pressure in fluids $`j=1`$ and $`2`$. The gradient in equation (4) is associated with the metric (1) and is obtained from the corresponding three-dimensional expression for the gradient in spherical coordinates $`(r,\theta ,\phi )`$, by keeping $`r=a`$ and noting that $`\theta =\rho /a`$. In contrast to gravity-driven flows in flat Hele-Shaw cells, the gravity term in (4) is not constant, but depends on the radial distance $`\rho `$. This is a manifestation of the cell spatial curvature together with its embedding in three dimensional space.

We can exploit the irrotational flow condition to define the velocity potential $`𝐯_j=\varphi _j`$. Using the velocity potential, we evaluate equation (4) for each of the fluids on the interface, subtract the resulting equations from each other, and divide by the sum of the two fluids’ viscosities to get

$$A\left(\frac{\varphi _1|_{}+\varphi _2|_{}}{2}\right)\left(\frac{\varphi _1|_{}\varphi _2|_{}}{2}\right)=\alpha (\kappa )|_{}\gamma a\mathrm{cos}\left(\frac{}{a}\right),$$
(5)
where A=$`(\eta _2\eta _1)/(\eta _2+\eta _1)`$ is the viscosity contrast, $`\alpha =b^2\sigma /[12(\eta _1+\eta _2)]`$ contains the surface tension, and $`\gamma =b^2g(\varrho _1\varrho _2)/[12(\eta _1+\eta _2)]`$ is a measure of gravitational force. To obtain (5) we used the pressure boundary condition $`p_2p_1=\sigma \kappa `$ at the interface $`\rho =`$, where $`\kappa `$ is the interfacial curvature .

Following steps similar to those performed in Ref. , we define Fourier expansions for the velocity potentials, which obey Laplace’s equation. We express $`\varphi _j`$ in terms of the perturbation amplitudes $`\zeta _n`$ by considering the kinematic boundary condition for flow on a sphere. As in the flat cell case, this condition refers to the continuity of the normal velocity across the fluid-fluid interface . Substituting these relations into equation (5), and Fourier transforming, yields the mode coupling equation of the Saffman-Taylor problem in a spherical Hele-Shaw cell, taking into account both injection and gravity

$$\dot{\zeta }_n=\lambda (n)\zeta _n+\underset{n^{}0}{}\left[F(n,n^{})\zeta _n^{}\zeta _{nn^{}}+G(n,n^{})\dot{\zeta }_n^{}\zeta _{nn^{}}\right],$$
(6)
where

$$\lambda (n)=\left[\frac{Q}{2\pi a^2S^2}\left(A|n|C\right)\frac{\alpha }{a^3S^3}|n|(n^21)+|n|\frac{\gamma }{a}\right]$$
(7)
is the linear growth rate, and

$$F(n,n^{})=\frac{|n|}{aS}\left\{\frac{QAC}{2\pi a^2S^2}\left[\frac{1}{2}sgn(nn^{})\right]\frac{\alpha C}{a^3S^3}\left[1\frac{n^{}}{2}(3n^{}+n)\right]+\frac{\gamma C}{2a}\right\}+\frac{Q}{4\pi a^3S},$$
(8)
$$G(n,n^{})=\frac{1}{aS}\left\{A\right|n|[1sgn(nn^{})]C\}$$
(9)
are the second-order mode coupling terms, with $`S=\mathrm{sin}(R/a)`$ and $`C=\mathrm{cos}(R/a)`$. The overdot denotes total time derivative and the sign function $`sgn(nn^{})=1`$ if $`(nn^{})>0`$ and $`sgn(nn^{})=1`$ if $`(nn^{})<0`$.

## III Discussion

### A Linear stability analysis with $`Q=0`$

First concentrate on the purely gravity-driven case, $`Q=0`$. In this situation, we consider flow in a closed cell obeying mass conservation. To drive the interface gravitationally, one could first allow the system to form a stable, unperturbed spherical cap at the south pole, and then invert the spherical cell to put the denser fluid on top in the unstable position at the north pole.

We begin by investigating the dispersion relation (7). Notice that the gravity term in Eq. (7) contains no explicit dependence on radial distance. This finding is somewhat surprising, since the gravity term in Darcy’s law (4) clearly presents a dependence on $`\rho `$ . Physical intuition suggests a factor of $`S`$ should multiply the gravity term in (7), making it vanish at the poles and become maximal at the equator. This apparent missing factor reappears if we rewrite the linear growth rate in terms of the wave number $`k=|n|/aS`$. The variables $`n`$ and $`k`$ are both useful: for example, $`n`$ occurs in integer values and counts the number of fingers. On the other hand, $`k`$ determines the characteristic wavelength of a perturbation.

To better understand the physical information behind the description of the linear stage in terms of $`n`$ and $`k`$, we plot the linear growth rate at a sequence of radial distances $`R=a\theta `$ in figures 2a and 2b. We use typical experimental parameters given in a recent experimental work in rotating, flat Hele-Shaw cells : fluid 1 is a silicone oil ($`\eta _10.5`$ $`\mathrm{g}`$/$`\mathrm{cm}`$ $`\mathrm{s}`$, $`\varrho _11.0`$ $`\mathrm{g}`$/$`\mathrm{cm}^3`$) and fluid 2 is air ($`\eta _20`$, $`\varrho _20`$). The thickness of the cell $`b=0.1`$ $`\mathrm{cm}`$ and the surface tension $`\sigma =20.7`$ $`\mathrm{dyne}/\mathrm{cm}`$. We set the radius of the sphere $`a=5`$ $`\mathrm{cm}`$ and acceleration of gravity $`g=980`$ $`\mathrm{cm}`$/$`\mathrm{s}^2`$.

Comparing figures 2a and 2b, it is clear that the fastest growing mode (number of fingers $`n`$ developing at a maximum growth rate) moves to large $`n`$ for large radial distance (and small $`n`$ for small radial distance) in order to keep the corresponding fastest growing $`k`$, and thus the apparent wavelength, nearly constant. The same effect occurs in radial flow in flat space . In addition to the shift to larger $`n`$, the growth rate of the fastest growing mode increases with radial distance because the net gravitational force depends on $`R`$. This differs from the case of radial flow in flat space , where the interface velocity falls off as $`1/R`$ causing the growth rate of the fastest growing mode to decrease.

Consider the purely linear contribution, which appears as the first term on the right hand side of equation (6). The condition $`Q=0`$ simplifies the theoretical description: $`R(t)=R_0=`$ constant, and consequently the linear growth rate $`\lambda (n)`$ is time independent. This implies that the actual relaxation or growth of mode $`n`$ is purely exponential

$$\zeta _n^{(Q=0)}(t)=\zeta _n(0)\mathrm{exp}[\lambda (n)t],$$
(10)
where $`\zeta _n(0)`$ is the initial perturbation amplitude. To see the overall effect of Eq. (10), we plot evolved interfaces using the same experimental parameters as those used in figures 2a and 2b. It is convenient to rewrite the net perturbation (3) in terms of cosine and sine modes $`\zeta (\theta ,t)=\zeta _0+\underset{n>0}{}\left[a_n(t)\mathrm{cos}(n\theta )+b_n(t)\mathrm{sin}(n\theta )\right]`$, where $`a_n=\zeta _n+\zeta _n`$ and $`b_n=i\left(\zeta _n\zeta _n\right)`$ are real-valued. We take into account modes $`n`$ ranging from $`n=1`$ up to $`20`$. Figure 3 depicts the evolution of the interfaces, for a random choice of phases, at time $`t=5`$ $`\mathrm{s}`$. We evolve from two distinct initial radii (a) $`R_0=a\pi /8`$, and (b) $`R_0=a\pi /4`$. In both cases $`|\zeta _n(0)|=0.05`$ $`\mathrm{cm}`$ and we use the same randomly chosen phases. It is evident that for larger radial distances (or equivalently, larger polar angles) we have both faster growth and more unstable modes, though the characteristic wavelength is not very different, in agreement with the predictions made from figures 2a and 2b.

### B Linear stability analysis with $`Q>0`$

Now consider the nonvanishing injection case. Taking $`Q>0`$ introduces two important new effects. First, inspecting equation (7), we see that $`Q`$ multiplies a term linear in $`n`$ that must be added to the growth rate $`\lambda `$ obtained with $`Q=0`$. If the inner fluid is high viscosity, so that $`A<0`$, then this new term in $`\lambda `$ is negative. All modes grow more slowly with $`Q>0`$, and some may even become stable. Thus, $`Q>0`$ diminishes the strength of the instability. Essentially, the flow fills in the gaps between the fingers.

The second effect is more subtle. Since the unperturbed interface radius $`R(t)`$ is now time dependent, the linear growth rate $`\lambda `$ evolves with time as the radial distance steadily grows. At any instant the interface evolves exponentially with an instantaneous growth rate $`\lambda `$ depending on the current $`R(t)`$. For sufficiently small $`Q`$ the instantaneous $`\lambda `$ at $`R(t)`$ nearly equals the $`Q=0`$ value for the same $`R(t)`$. We call this a quasistatic approximation, and within this approximation it is clear that an increasing number of modes become unstable as time progresses due to the steady increase of $`R(t)`$. The evolution of $`\lambda `$ in the quasistatic approximation is given by the series of curves shown in figure 2.

If we wish to study the linear growth without making the quasistatic approximation, we may integrate the equation of motion (6) exactly, keeping only terms of first order in $`\zeta `$ on the right-hand side. Because $`\lambda `$ is time dependent, through the variation of $`R(t)`$, the exact solution of the linear equation of motion becomes non-exponential and can be written as

$$\zeta _n^{(Q>0)}(t)=\zeta _n(0)\frac{S_0}{S}\left[\frac{\mathrm{tan}(R/2a)}{\mathrm{tan}(R_0/2a)}\right]^{A|n|}\mathrm{exp}\left\{\frac{2\pi |n|}{Q}\left[\frac{\alpha (n^21)}{a}\left(\frac{C}{S}\frac{C_0}{S_0}\right)+\gamma a\left(C_0C\right)\right]\right\}$$
(11)
where $`S_0=\mathrm{sin}(R_0/a)`$, and $`R`$, $`C`$ and $`S`$ all depend on time.

### C Mode coupling analysis

We use the mode coupling equation (6) to investigate the effect of the nonlinear terms in the interface evolution. We are interested in studying how gravity influences the shape of the fingering structures, focusing on finger-tip behavior. Finger tip-splitting and tip-sharpening phenomena are related to the influence of a fundamental mode $`n`$ on the growth of its harmonic $`2n`$. Without loss of generality we may choose the phase of the fundamental mode so that $`a_n>0`$ and $`b_n=0`$. We replace the time derivative terms $`\dot{a}_n`$ and $`\dot{b}_n`$ by $`\lambda (n)a_n`$ and $`\lambda (n)b_n`$, respectively, for consistent second order expressions. Under these circumstances the equation of motion for the harmonic cosine mode becomes

$$\dot{a}_{2n}=\lambda (2n)a_{2n}+\frac{1}{2}T(2n,n)a_n^2$$
(12)
where the finger-tip function is

$$T(2n,n)=\frac{2n}{aS}\left\{\frac{Q}{2\pi a^2S^2}\left[\frac{(C^2+1)}{4n}AC\right]+\frac{3\alpha C(2n^21)}{2a^3S^3}\right\}.$$
(13)
The corresponding equation for the sine modes $`b_{2n}`$ is not as interesting as (12), since growth of $`b_{2n}`$ is uninfluenced by $`a_n`$.

The sign of $`T(2n,n)`$ dictates whether finger tip-splitting or finger tip-sharpening is favored by the dynamics . If $`T(2n,n)<0`$, at second order the result is a driving term of order $`a_n^2`$ forcing growth of $`a_{2n}<0`$. With this particular phase of the harmonic forced by the dynamics, the $`n`$ outwards-pointing fingers of the fundamental mode $`n`$ tend to split. In contrast, if $`T(2n,n)>0`$ growth of $`a_{2n}>0`$ would be favored, leading to outwards-pointing finger tip-sharpening.

A noteworthy point about equation (13) is that it shows no dependence whatsoever on gravity. In the evaluation of $`T(2n,n)`$ from equation (6) we found that the term involving gravity in $`F(2n,n)`$ exactly cancels against the term involving gravity in $`\lambda (n)G(2n,n)`$. The second order term driving tip-splitting in equation (12) is therefore independent of the force of gravity, though gravity does generally influence mode coupling at second and higher orders.

Inspecting Eq. (13) we find that, since $`A=1`$ and $`C>0`$, the finger-tip function $`T(2n,n)>0`$ for both $`Q=0`$ and $`Q>0`$. Eq. (13) predicts that gravity-driven flow on a sphere leads to patterns showing enhanced finger tip-narrowing. Informal studies of chocolate syrup fingers on ice cream scoops support this claim.

ACKNOWLEDGMENTS

J.A.M., F.P., and F.M. thank CNPq and FINEP for financial support. Work of M.W. was supported in part by the National Science Foundation grant No. DMR-9732567.

FIGURE CAPTIONS

FIG. 1: Schematic configuration of flow in a spherical Hele-Shaw cell.

FIG. 2: Linear growth rate as a function of (a) mode number $`n`$ and (b) wave number $`k`$ at a sequence of radial distances (1) $`R=a\pi /8`$, (2) $`R=a\pi /4`$, (3) $`R=a3\pi /8`$ and (4) $`R=a\pi /2`$.

FIG. 3: Interface evolution according to Eq. (10) at $`t=5`$ $`\mathrm{s}`$, for (a) $`R_0=a\pi /8`$ and (b) $`R_0=a\pi /4`$. Other parameters are given in the text.


text:
# The Structure of the Virtual Photon1footnote 11footnote 1Talk given at the International Conference on the Structure and Interactions of the Photon (PHOTON 2000), Ambleside, UK, August 26^{𝑡ℎ}-31^{𝑠𝑡}, 2000.

## INTRODUCTION

It has been long established that real photons have a partonic structure from measurements of the photon structure function $`F_2^\gamma `$ in $`e\gamma `$ interactions eg  and observation of resolved photon processes in $`\gamma p`$ interactions gp  and of single- and double-resolved processes in $`\gamma \gamma `$ interactions gg . Thus, it is natural to expect that virtual photons also have a partonic structure SaS ; DG ; theory ; predictions . QCD predicts that the parton densities of virtual photons become logarithmically suppressed as the virtuality of the probed photon $`P^2`$ increases for fixed $`\mu ^2`$, the scale of the interaction; $`\mu ^2`$ is usually taken as $`Q^2`$ ($`Q^2`$ is the virtuality of the probing photon) in deep inelastic scattering (DIS) $`e\gamma `$ and the jet transverse energy ($`E_T^{jet}`$) in jet production.

The virtual photon structure function $`F_2^\gamma ^{}`$ at leading order (LO) in QCD is given by

$$F_2^\gamma ^{}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)=\underset{q}{}x_\gamma ^{}e_q^2[f_{q/\gamma ^{}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)+f_{\overline{q}/\gamma ^{}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)],$$

where $`x_\gamma ^{}`$ is the fraction of the virtual photon momentum taken by the interacting parton, $`\mu _{F_\gamma ^{}}^2`$ is the virtual photon fragmentation scale and $`e_q`$ is the quark charge. The virtual photon quark densities $`f_{q/\gamma ^{}}`$ contain two terms, as in the case of the real photon,

$$f_{i/\gamma ^{}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)=f_{i/\gamma ^{}}^{\mathrm{had}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)+f_{i/\gamma ^{}}^{\mathrm{anom}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2),$$

one term associated to the non-perturbative hadronic component ($`f^{\mathrm{had}}`$) and a term $`f^{\mathrm{anom}}`$, unique to the photon, which expresses the direct coupling of the photon to a $`q\overline{q}`$ pair, calculable in perturbative QCD (pQCD). The hadronic component is the one expected to decrease as $`P^2`$ increases. The first measurements of the structure function of the virtual photon were done by PLUTO pluto . Only recently new measurements from LEP have become available l3 .

At HERA, the virtual photon structure is studied in jet production mediated by virtual photons. In LO QCD, two processes are expected to contribute to the jet cross section: resolved processes (figure 1a) in which the virtual photon interacts with the proton via its hadronic component giving two jets in the final state, and the direct processes (figure 1b) in which the virtual photon interacts as a point-like particle with a parton from the proton.

The dijet production cross sections in LO QCD for direct and resolved processes are given by

$$\sigma _D^{LO,ep}=𝑑\mathrm{\Omega }f_{\gamma ^{}/e}(y,P^2)f_{j/p}(x_p,\mu _{F_p}^2)𝑑\sigma (\gamma ^{}j\mathrm{jet}\mathrm{jet})$$

and

$$\sigma _R^{LO,ep}=𝑑\mathrm{\Omega }f_{\gamma ^{}/e}(y,P^2)f_{i/\gamma ^{}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)f_{j/p}(x_p,\mu _{F_p}^2)𝑑\sigma (ij\mathrm{jet}\mathrm{jet}),$$

where the integrals are performed over the phase space represented by “$`d\mathrm{\Omega }`$”; $`f_{\gamma ^{}/e}(y,P^2)`$ is the flux of virtual photons in the positron and $`y`$ is the fraction of the positron energy taken by the virtual photon; $`f_{j/p}(x_p,\mu _{F_p}^2)`$ are the parton densities in the proton and $`x_p`$ is the fraction of the proton momentum taken by parton $`j`$; $`\mu _{F_p}^2`$ is the proton fragmentation scale; and $`d\sigma (\gamma ^{}(i)j\mathrm{jet}\mathrm{jet})`$ is the subprocess cross section, calculable in pQCD. In the case of resolved processes, there is an additional ingredient: the parton densities in the virtual photon $`f_{i/\gamma ^{}}(x_\gamma ^{},P^2,\mu _{F_\gamma ^{}}^2)`$, for which up to the present there is only very little experimental information. This allows the study of the virtual photon structure by measuring jet cross sections.

The contribution to the dijet cross section from resolved processes is expected to decrease relative to the contribution from the direct processes as $`P^2\mu ^2`$; this means that the partonic content of the virtual photon becomes suppressed as $`P^2`$ increases. Experimentally, resolved and direct processes are separated by using the variable $`x_\gamma ^{OBS}`$,

$$x_\gamma ^{OBS}=\frac{1}{2yE_e}(E_T^{jet1}e^{\eta ^{jet1}}+E_T^{jet2}e^{\eta ^{jet2}}),$$
(1)
where $`\eta ^{jet1,2}`$ is the jet pseudorapidity and $`E_e`$ is the positron beam energy. $`x_\gamma ^{OBS}`$ measures the fraction of the virtual photon energy invested in the production of the dijet system and it is well defined at all orders in pQCD. For resolved (direct) processes, $`x_\gamma ^{OBS}`$ takes low (high) values.

The measurements presented here were performed with the ZEUS detector at HERA. During 1995 to 1997 HERA operated with positrons of energy $`E_e=27.5`$ GeV and protons of energy $`E_p=820`$ GeV.

## DIJET CROSS SECTIONS IN THE LABORATORY FRAME

The $`x_\gamma ^{OBS}`$ distribution displays a high sensitivity to the virtual photon structure. Figure 2 shows the $`x_\gamma ^{OBS}`$ distribution for a sample of dijet events in three regions of $`P^2`$. In all the $`P^2`$ regions studied, a significant resolved contribution is needed in the Monte Carlo to describe the data for $`x_\gamma ^{OBS}<0.75`$. Direct processes alone (the shaded area in figure 2) cannot describe the data lab .

Dijet cross sections as a function of $`x_\gamma ^{OBS}`$ have been measured for jets found using the longitudinally invariant $`k_T`$ cluster algorithm in the inclusive mode kt . The measurements have been performed in the kinematic region given by $`0.2<y<0.55`$ and in three regions of $`P^2`$: $`P^20`$ (the photoproduction regime in which the scattered positron is lost in the beam pipe, figure 3a), $`0.1<P^2<0.55`$ GeV2 (the intermediate $`P^2`$ region, in which the scattered positron is detected in the beam pipe calorimeter, figure 3b) and $`1.5<P^2<4.5`$ GeV2 (the low $`P^2`$ DIS regime, in which the scattered positron is detected in the uranium-scintillator calorimeter, figure 3c). The events are required to have at least two jets with $`E_T^{jet}>5.5`$ GeV and $`1.125<\eta ^{jet}<2.2`$. The data show that the shape of the measured cross section changes markedly with increasing $`P^2`$: the cross section for low $`x_\gamma ^{OBS}`$ values decreases faster than at high $`x_\gamma ^{OBS}`$ values.

The measurements have been compared to several Monte Carlo models using various parametrisations of the photon parton densities (PDFs). The GRV grv  and WHIT2 whit2  photon PDFs have been extracted for real photons and have no $`P^2`$ dependence. On the other hand, the SaS 1D SaS  photon PDFs consist of two components with different $`P^2`$ dependence: the hadronic part decreases as $`(m_\rho ^2/(m_\rho ^2+P^2))^2`$ with increasing $`P^2`$, while the $`P^2`$ dependence of the anomalous component goes like $`\mathrm{log}(\mu ^2/P^2)`$. The predictions of HERWIG herwig  using the GRV PDFs do not describe the data, as expected since they have no $`P^2`$ dependence. The prediction based on WHIT2 provides a reasonable description of the data in the photoproduction and in the intermediate $`P^2`$ regimes. The predictions using SaS 1D agree with the data at high $`x_\gamma ^{OBS}`$ in all regimes. The DIS Monte Carlo LEPTO lepto  is compared to the data in the DIS regime. It underestimates the data at low $`x_\gamma ^{OBS}`$, i.e. parton shower effects are not the sole contribution to this region of $`x_\gamma ^{OBS}`$ and $`P^2`$.

The low $`E_T^{jet}`$ region is affected by contributions from the presence of a possible underlying event, which may mask the $`P^2`$ evolution of the virtual photon. Therefore, measurements of the dijet cross sections have been performed at higher transverse energies. For these measurements, two jets were required with $`E_T^{jet,1(2)}>7.5(6.5)`$ GeV and $`1.125<\eta ^{jet1,2}<1.875`$. The cross sections are shown in figure 4. At low $`x_\gamma ^{OBS}`$ the cross sections are also observed to decrease faster with $`P^2`$ than at high $`x_\gamma ^{OBS}`$. The predictions from HERWIG based on the GRV PDFs are in good agreement with the data in the two lowest $`P^2`$ ranges, but fail at higher $`P^2`$. On the other hand, the predictions based on SaS 1D agree with the data at high $`P^2`$ but underestimate the data at low and intermediate $`P^2`$ at low $`x_\gamma ^{OBS}`$. The predictions from LEPTO at high $`P^2`$ underestimate the data.

The $`P^2`$ evolution of the virtual photon structure was studied further by measuring the ratio of the dijet cross section for $`x_\gamma ^{OBS}<0.75`$ (enriched in resolved processes) to the dijet cross section for $`x_\gamma ^{OBS}>0.75`$ (enriched in direct processes). The ratio (see figure 5) falls steeply with $`P^2`$ which may be interpreted as the suppression of the resolved photon component as the virtuality of the photon increases. The prediction of HERWIG based on the GRV PDFs is constant, as expected from a photon structure without $`P^2`$ dependence. The prediction of SaS 1D decreases with $`P^2`$ but lies below the data at low $`P^2`$ and the prediction from LEPTO shows that the contribution to the ratio from parton shower effects alone is not enough to explain the $`P^2`$ dependence of the data. QCD predictions have been calculated using the program JetViP jetvip  and different parametrisations of the photon PDFs. The predictions show sensitivity to the choice of PDF but lie well below the data. Hadronisation corrections ($`2030\%`$) are insufficient to explain the discrepancy.

## DIJET CROSS SECTIONS IN THE $`\gamma ^{}p`$ FRAME

For $`P^20`$, the definition of $`x_\gamma ^{OBS}`$ (equation 1) is still valid. However, the photon remnant may have high transverse energy in the laboratory frame and be considered as a jet. The solution to avoid this problem is to transform to the $`\gamma ^{}p`$ frame. In such a frame the photon remnant has very low transverse energy and will not be mistaken for a jet emanating from the hard interaction.

Figure 6a shows the measurements of the dijet cross section as a function of $`x_\gamma ^{OBS}`$ in the $`\gamma ^{}p`$ frame gpf . The measurements have been performed in the kinematic region given by $`0.2<y<0.55`$ and $`0.1<P^2<10^4`$ GeV2. The events are required to have at least two jets with $`E_{T,\gamma ^{}p}^{jet1(2)}>7.5(6.5)`$ GeV and $`3<\eta _{\gamma ^{}p}^{jet1,2}<0`$. Also for these measurements the shape of the measured cross section changes with increasing $`P^2`$: the cross section for low $`x_\gamma ^{OBS}`$ values falls more rapidly with increasing $`P^2`$ than at high $`x_\gamma ^{OBS}`$ values.

The predictions of HERWIG based on the SaS 1D PDFs (see figure 6a) give a good description of the data in the high $`P^2`$ region but fail in the intermediate $`P^2`$ region. A resolved photon component is needed to describe the data up to $`P^249`$ GeV2. Above this value, the HERWIG prediction is dominated by the direct component and describes the data well. Therefore, for $`P^2(E_T^{jet})^2`$, direct processes alone are able to describe the data.

The ratio of the cross section for $`x_\gamma ^{OBS}<0.75`$ and $`x_\gamma ^{OBS}>0.75`$ also falls steeply in the $`\gamma ^{}p`$ frame with increasing $`P^2`$ (figure 6b). The prediction of SaS 1D shows a decrease with $`P^2`$ but lies below the data in the whole $`P^2`$ range. The ratio of the data to the SaS 1D prediction (inset in figure 6b) has a constant value of $`1.3`$. This indicates that the resolved component suppression included in these PDFs is in agreement with the data, but they underestimate the fraction of resolved component by $`30\%`$.

## SUMMARY AND CONCLUSIONS

Dijet cross sections as a function of $`x_\gamma ^{OBS}`$ have been measured in the laboratory and $`\gamma ^{}p`$ frames in the kinematic region given by $`0.2<y<0.55`$ and $`0P^2<10^4`$ GeV2 for jets found using the longitudinally invariant $`k_T`$ cluster algorithm in the inclusive mode. The $`x_\gamma ^{OBS}`$ dependence of the measured cross sections changes with increasing $`P^2`$: the cross section for low $`x_\gamma ^{OBS}`$ values decreases faster than at high $`x_\gamma ^{OBS}`$ values. The predictions of HERWIG based on the SaS 1D photon PDFs describe the data for $`1.5<P^2<10^4`$ GeV2 but fail in the region $`0.1<P^2<0.55`$ GeV2. A resolved photon component is needed to describe the data up to $`P^249`$ GeV2, i.e. where the virtual photon is probed at a scale comparable to the hard interaction scale $`(E_T^{jet})^2(7\mathrm{GeV})^2`$. The ratio of the dijet cross section $`\sigma (x_\gamma ^{OBS}<0.75)/\sigma (x_\gamma ^{OBS}>0.75)`$ decreases with $`P^2`$; the predicted $`P^2`$ dependence of the ratio agrees with the data, but it underestimates the fraction of the resolved component by $`30\%`$. This result can be interpreted in terms of a resolved photon component that is suppressed as the photon virtuality increases.

Acknowledgements. I would like to thank my colleagues from ZEUS for their help in preparing this report and the organisers of the conference for providing a warm atmosphere and hospitality.


text:
# Invariant Polynomial Functions on 𝑘 qudits

## 1 Introduction

In quantum mechanics, a combination of states in Hilbert spaces $`H_1`$, .., $`H_k`$ leads to a state in the tensor product Hilbert space $`H_1\mathrm{}H_k`$. Such a state will be called here a tensor state. In this paper we take $`H_1=\mathrm{}=H_k=^n`$ where $`n>1`$. Then a tensor state is a joint state of $`k`$ qudits. It would be very interesting to classify tensor states in $`(^n)^k`$ up to the action of the product $`U(n)^k`$ of unitary groups of local symmetries. A natural approach to this is to study the algebra of invariant polynomials. This approach was developed by Rains \[R\], by Grassl, Rötteler and Beth \[G-R-B1\] \[G-R-B2\], by Linden and Popescu \[L-P\] and by Coffman, Kundu and Wootters \[C-K-W\]. These authors study the ring of invariant polynomials in the components of a tensor state in $`(^n)^k`$ and in their complex-conjugates. For $`k`$ qubits, explicit descriptions of invariants are given in \[G-R-B1\], \[G-R-B2\], \[L-P\] and in \[C-K-W\].

In this paper the symmetry group we consider is the product $`G=SU(n)^k`$ of special unitary groups; one thinks of $`G`$ as the special group of local symmetries. We study the $`G`$-invariant polynomial functions $`Q`$ on the tensor states in $`(^n)^k`$ (we discuss in §2 how this is relevant to the description of the $`G`$-orbits). We consider polynomials in the entries of a tensor state, in other words, holomorphic polynomials.

Let $`_{n,k,d}`$ be the space of homogeneous degree $`d`$ polynomial functions on tensor states in $`(^n)^k`$. Let $`_{n,k,d}^G`$ be the space of $`G`$-invariants in $`_{n,k,d}`$. See §2 for more discussion. We reduce the problem of computing $`_{n,k,d}^G`$ to a problem in the invariant theory of the symmetric group $`𝔖_d`$ (Proposition 2.1). In particular, $`_{n,k,d}^G`$ is non-zero only if $`d`$ is a multiple of $`n`$. So the “first”case is $`d=n`$; we examine this in §3. We find that if $`k`$ is odd then $`_{n,k,n}^G=0`$ while if $`k`$ is even then $`_{n,k,n}^G`$ is $`1`$-dimensional. In the latter case we write down (§3) explicitly the corresponding invariant polynomial $`P_{n,k}`$ in $`_{n,k,n}`$; we find $`P_{n,k}`$ is a natural generalization of the determinant of a square matrix.

For fixed $`n,d`$ the direct sum $`_k_{n,k,d}`$ is an associative algebra. We study the asymptotic behavior of $`dim_{n,k,d}^G`$ as $`k\mathrm{}`$ in §4. In §5, we specialize to the case of $`k`$-qubits, i.e. $`n=2`$. We compute the dimension of the space $`_{2,k,4}^G`$ of degree $`4`$ invariants as well as the dimension of the space of invariants in $`_{2,k,4}^G`$ under the natural action of $`𝔖_k`$. We show that $`_k_{2,k,4}^{𝔖_k\text{ }\times G}`$ is a polynomial algebra on $`2`$ generators. For $`k5`$ we describe the representation of $`𝔖_k`$ on $`_{2,k,4}^G`$. For $`k=4`$ we find some interesting relations with the results on classification of tensor states in $`(^2)^4`$ given in \[B\].

We thank Markus Grassl for his useful comments on the first version of this paper.

## 2 Polynomial invariants of tensor states

We will consider $`(^n)^k`$ as a space of contravariant tensor states $`u`$. Then (once we fix a basis of $`^n`$) $`u`$ is given by $`n^k`$ components $`u^{p_1p_2\mathrm{}p_k}`$. We consider the algebra $`_{n,k}`$ of polynomial functions on $`(^n)^k`$. So $`_{n,k}`$ is the polynomial algebra $`[x_{p_1p_2\mathrm{}p_k}]`$ in the $`n^k`$ coordinate functions $`x_{p_1p_2\mathrm{}p_k}`$. We have a natural algebra grading $`_{n,k}=_{d=0}^{\mathrm{}}_{n,k,d}`$ where $`_{n,k,d}`$ is the space of homogeneous degree $`d`$ polynomial functions.

A function in $`_{n,k,d}`$ amounts to a symmetric degree $`d`$ covariant tensor $`Q`$ in $`(^n)^k`$. So $`Q`$ has $`n^{dk}`$ components $`Q_{i_{11}\mathrm{}i_{dk}}`$ where we think of the indices $`i_{ab}`$ as being arranged in a rectangular array of $`d`$ rows and $`k`$ columns and $`Q_{i_{11}\mathrm{}i_{dk}}`$ is invariant under permutations of the rows of the array. Then $`Q`$ defines the function

$$uQ_{i_{11}\mathrm{}i_{dk}}u^{i_{11}i_{12}\mathrm{}i_{1k}}u^{i_{21}i_{22}\mathrm{}i_{2k}}\mathrm{}u^{i_{d1}i_{d2}\mathrm{}i_{dk}}$$
(2.1)
where we used the usual Einstein summation convention. In this way, $`_{n,k}`$ identifies with $`S^d((^n)^k)`$.

Now the group $`G=SU(n)^k`$ acts on our tensor states $`u`$ and tensors $`Q`$ as follows. Let the matrix $`g_{ij}`$ live in the $`m`$-th copy of $`SU(n)`$ and let $`g^{ij}`$ be the inverse matrix. Then $`g_{ij}`$ transforms $`u^{p_1p_2\mathrm{}p_k}`$ into $`g_{p_mq_m}u^{q_1q_2\mathrm{}\mathrm{}q_k}`$ and $`Q_{i_{11}\mathrm{}i_{dk}}`$ into $`Q_{j_{11}\mathrm{}j_{dk}}g^{j_{1m}i_{1m}}g^{j_{2m}i_{2m}}\mathrm{}g^{j_{dm}i_{dm}}`$. The identification of $`_{n,k,d}`$ with $`S^d((^n)^k)`$ is $`G`$-equivariant.

We are interested in the algebra $`_{n,k}^G=_{d=0}^{\mathrm{}}_{n,k,d}^G`$ of $`G`$-invariants. We view this as a first step towards studying the orbits of $`G`$ on $`(^n)^k`$. One can first study the orbits of the complex group $`G_{}=SL(n,)^k`$ and then decompose the $`G_{}`$-orbits under the $`G`$-action. Note that a polynomial is $`G`$-invariant if and only if its is $`G_{}`$-invariant. The closed $`G_{}`$ orbits play a special role – they are the most degenerate orbits. Given any orbit $`Y`$, its closure contains a unique closed orbit $`Z`$; then points in $`Y`$ degenerate to points in $`Z`$. The $`G_{}`$-invariant functions separate the closed orbits; they take the same values on $`Y`$ and on $`Z`$. The set of closed orbits of $`G_{}`$ in $`(^n)^k`$ has the structure of an affine complex algebraic variety with $`_{n,k}^G`$ as its algebra of regular functions. Thus a complete description of $`_{n,k}^G`$ would lead to a precise knowledge of the closed $`G_{}`$-orbits.

Our approach is thus somewhat different from that of \[R\] \[G-R-B1\] \[G-R-B2\] \[L-P\] \[C-K-W\] who study the invariant functions on $`(^n)^k`$ which are polynomials in the $`x_{p_1\mathrm{}p_k}`$ and in their complex conjugates; these can also be described as the invariant polynomial functions on $`(^n)^k\overline{(^n)^k}`$.

At this point it is useful to examine the case $`k=2`$. We can identify $`(^n)^2`$ with the space $`M_n()`$ of square matrices and then $`G=SU(n)^2`$ acts on $`M_n()`$ by $`(g,h)u=guh^1`$. So $`_{n,k,d}^G`$ is the space of homogeneous degree $`d`$ polynomial functions $`Q`$ of an $`n`$ by $`n`$ matrix $`u`$ which are bi-$`SL(n,)`$-invariant, i.e. $`Q(guh^1)=Q(u)`$ for $`g,hSL(n,)`$. Then $`Q`$ is, up to scaling, the $`r`$th power of the determinant $`D`$ for some $`r`$. Hence $`d=rn`$. It follows that $`_{n,2}^G`$ is the polynomial algebra $`[D]`$. Thus the space of closed orbits for $`SL(n,)^2`$ identifies with $``$, where $`\lambda `$ corresponds to the unique closed orbit $`Z_\lambda `$ inside the set $`X_\lambda `$ of matrices of determinant $`\lambda `$. For $`\lambda 0`$, $`Z_\lambda =X_\lambda `$ while for $`\lambda =0`$, $`Z_0`$ reduces to the zero matrix.

We view $`S^d((^n)^k)`$ as the space of invariants for the symmetric group $`𝔖_d`$ acting on $`((^n)^k)^d`$. So

$$_{n,k,d}^G=(((^n)^k)^d)^{G\times 𝔖_d}=(((^n)^d)^k)^{G\times 𝔖_d}$$
(2.2)
Recall the Schur decomposition $`(^n)^d=_\alpha S^\alpha (^n)E_\alpha `$ where $`\alpha `$ ranges over partitions of $`d`$ with at most $`n`$ rows, $`S^\alpha (^n)`$ is the irreducible covariant representation of $`SU(n)`$ given by the Schur functor $`S^\alpha `$, and $`E_\alpha `$ is the corresponding irreducible representation of $`𝔖_d`$. We use the convention that $`E_\alpha `$ is the trivial representation if $`\alpha =[d]`$, while $`E_\alpha `$ is the sign representation if $`\alpha =[1^d]`$. Thus we have

$$((^n)^d)^k=\underset{|\alpha _1|=\mathrm{}=|\alpha _k|=d}{}S^{\alpha _1}(^n)\mathrm{}S^{\alpha _k}(^n)E_{\alpha _1}\mathrm{}E_{\alpha _k}$$
(2.3)
Now taking the invariants under $`G\times 𝔖_d`$ we get

$$_{n,k,d}^G=\underset{|\alpha _1|=\mathrm{}=|\alpha _k|=d}{}S^{\alpha _1}(^n)^{SU(n)}\mathrm{}S^{\alpha _k}(^n)^{SU(n)}\left(E_{\alpha _1}\mathrm{}E_{\alpha _k}\right)^{𝔖_d}$$
(2.4)
The representation $`S^{\alpha _j}(^n)`$, since it is irreducible, has no $`SU(n)`$-invariants except if $`S^{\alpha _j}(^n)=`$ is trivial. This happens if and only if $`\alpha _j`$ is a rectangular partition with all columns of length $`n`$. This proves:

###### Proposition 2.1.

If $`n`$ does not divide $`d`$, then $`_{n,k,d}^G=0`$. If $`d=nr`$, then $`_{n,k,d}^G`$ is isomorphic to $`(E_\pi ^k)^{𝔖_d}`$ where $`\pi =[r^n]`$.

The permutation action of $`𝔖_k`$ on $`(^n)^k`$ induces an action of $`𝔖_k`$ on $`_{n,k,d}^G`$.

###### Corollary 2.2.

The isomorphism of Proposition 2.1 intertwines the $`𝔖_k`$-action on $`_{n,k,d}^G`$ with the action of $`𝔖_k`$ on $`(E_\pi ^k)^{𝔖_d}`$ given by permuting the $`k`$ factors $`E_\pi `$.

## 3 The generalized determinant function

Given $`n`$ and $`k`$, we want to find the smallest positive value of $`d`$ such that $`_{n,k,d}^G0`$. By Proposition 2.1, the first candidate is $`d=n`$.

###### Corollary 3.1.

$`_{n,k,n}^G0`$ iff $`k`$ is even. In that case, $`_{n,k,n}^G`$ is one-dimensional and consists of the multiples of the function $`P_{n,k}`$ given by

$$P_{n,k}(u)=_{\sigma _2,\mathrm{},\sigma _k𝔖_n}ϵ(\sigma _2)\mathrm{}ϵ(\sigma _k)_{h=1}^nu^{hh_{\sigma _2}\mathrm{}h_{\sigma _k}}$$
(3.1)
where $`h_{\sigma _j}=\sigma _j(h)`$.

###### Proof.

By Proposition 2.1, we need to compute $`(E_\pi ^k)^{𝔖_d}`$. For $`d=n`$, $`\pi =[1^n]`$ and so $`E_\pi `$ is the sign representation of $`𝔖_n`$. Then $`(E_\pi ^k)`$ is one-dimensional and carries the trivial representation if $`k`$ is even, or the sign representation if $`k`$ is odd.

Now for $`k`$ even, we can easily compute a non-zero function $`P=P_{n,k}`$ in $`_{n,k,n}`$. For $`S^\pi (^n)`$ is the top exterior power $`^n^n`$. Thus $`P`$ is a non-zero element of the one-dimensional subspace $`(^n^n)^k`$ of $`((^n)^n)^k`$. The tensor components of $`P`$ are then given by $`P_{i_{11}\mathrm{}i_{nk}}=\frac{1}{n!}ϵ(\sigma _1)\mathrm{}ϵ(\sigma _k)`$ if for each $`j`$, the column $`i_{1j},\mathrm{},i_{nj}`$ is a permutation $`\sigma _j`$ of $`1,\mathrm{},n`$ and $`0`$ otherwise. Then we get

$$P_{n,k}(u)=\frac{1}{n!}_{\sigma _1,\mathrm{},\sigma _k𝔖_n}ϵ(\sigma _1)\mathrm{}ϵ(\sigma _k)_{h=1}^nu^{h_{\sigma _1}\mathrm{}h_{\sigma _k}}$$
(3.2)
where $`h_{\sigma _i}=\sigma _i(h)`$. The expression is very redundant, as each term appears $`n!`$ times. We remedy this by restricting the first permutation $`\sigma _1`$ to be $`1`$. This gives (3.1). ∎

$`P_{n,k}`$ is a *generalized determinant*; $`P_{n,k}`$ is invariant under the $`𝔖_k`$-action. For $`k=2`$, (3.1) reduces to the usual formula for the matrix determinant.

Recall that the rank $`s`$ of a tensor state $`u`$ in $`(^n)^k`$ is the smallest integer $`s`$ such that $`u`$ can be written as $`u=v_1+v_2+\mathrm{}+v_s`$, where the $`v_i`$ are decomposable tensor states $`v_i=w_{i1}w_{i2}\mathrm{}w_{ik}`$. There is a relation between the rank and the vanishing of $`P_{n,k}`$ as follows:

###### Corollary 3.2.

If the tensor state $`u`$ in $`(^n)^k`$ has rank less than $`n`$, then $`P_{n,k}(u)=0`$.

It is easy to find a tensor state $`u`$ of rank $`n`$ such that $`P_{n,k}(u)`$ is non-zero. For instance, $`P_{n,k}(u)=1`$ if $`u`$ has all components zero except $`u^{1\mathrm{}1}=\mathrm{}=u^{n\mathrm{}n}=1`$. For $`k=2`$, $`P_{n,k}(u)=0`$ implies $`u`$ has rank less than $`n`$. For bigger (even) $`k`$, this is false, if $`n`$ is large enough. This happens essentially because the rank of $`u`$ can be very large (at least $`\frac{n^k}{knk+1}`$). Thus $`P_{n,k}`$ gives only partial information about the rank.

## 4 Asymptotics as $`k\mathrm{}`$

Suppose we fix $`n`$ and $`d`$ where $`d=rn`$. Then there is a $`G`$-invariant associative graded algebra structure $`PQ`$ on the direct sum $`_k_{n,k,d}^G`$. Indeed, the product of tensors induces a $`(G\times 𝔖_d)`$-invariant map $`V^kV^lV^{(k+l)}`$ where $`V=(^n)^d`$. The induced multiplication on the spaces of $`(G\times 𝔖_d)`$-invariants gives the product on $`_k_{n,k,d}^G`$, where we use the identification in (2.2). This multiplication corresponds, under the isomorphism of Proposition 2.1, to the product map $`E_\pi ^kE_\pi ^lE_\pi ^{(k+l)}`$. This structure is very useful. For instance, if $`d=n`$, then $`P_{n,k}P_{n,l}=\frac{1}{n!}P_{n,k+l}`$. Thus the determinant $`P_{n,2}`$ determines $`P_{n,2m}`$ in that the $`m`$-fold product $`P_{n,2}\mathrm{}P_{n,2}`$ is equal to $`(n!)^{m+1}P_{n,2m}`$.

We will study the size of the algebra $`_k_{n,k,d}^G`$ by finding an asymptotic formula for the dimension of $`_{n,k,d}^G`$. We do this for $`r2`$. Indeed for $`r=1`$ we already know $`dim_{n,k,n}^G`$ is $`1`$ if $`k`$ is even or $`0`$ if $`k`$ is odd; we call this the *static* case. The asymptotics involve the number

$$p=dimE_\pi =d!_{m=0}^{n1}\frac{m!}{(m+r)!}$$
(4.1)
where $`\pi =[r^n]`$ as in Proposition 2.1. Our formula for $`p`$ is immediate from the hook formula for the dimension of an irreducible symmetric group representation.

###### Proposition 4.1.

Assume $`d=rn`$ with $`r2`$. Then $`dim_{n,k,d}^Gc{\displaystyle \frac{p^k}{d!}}`$ as $`k\mathrm{}`$, where $`c=1`$ with one exception: $`c=4`$ if $`n=2,d=4`$.

###### Proof.

Let $`s=dim_{n,k,d}^G=dim(E_\pi ^k)^{𝔖_d}`$. Then $`s=\frac{1}{d!}_{\sigma 𝔖_d}\chi (\sigma )^k`$ where $`\chi :𝔖_d`$ is the character of $`E_\pi `$. If $`\sigma `$ acts trivially on $`E_\pi `$, then $`\chi (\sigma )=p`$. If $`\sigma `$ acts non-trivially, we claim $`|\chi (\sigma )|<p`$. To show this, it suffices to show that $`\sigma `$ has at least two distinct eigenvalues on $`E_\pi `$; this is because $`\chi (\sigma )`$ is the sum of the $`p`$ eigenvalues of $`\sigma `$. Now the set $`𝔗_d`$ of $`\sigma 𝔖_d`$ which act on $`E_\pi `$ by a scalar is a normal subgroup of $`𝔖_d`$. So if $`d5`$, then $`𝔗_d`$ is $`\{1\}`$, the alternating group $`𝔄_d`$ or $`𝔖_d`$. We can easily rule out the latter two possibilities, so $`𝔗_d=\{1\}`$, which proves our claim. If $`d4`$, then (since $`r>1`$ and $`n>1`$), we have $`d=4`$, $`n=2`$ and $`\pi =[2,2]`$. Our claim is clear here since $`𝔖_4`$ acts on $`E_\pi `$ through the reflection representation of $`𝔖_3`$ on $`^2`$.

Therefore we have $`s=c\frac{p^k}{d!}+o(p^k)`$ as $`k\mathrm{}`$ where $`c`$ is cardinality of the kernel of $`𝔖_dAutE_\pi `$. Our work in the previous paragraph computes $`c`$. ∎

Proposition 4.1 implies that the algebra $`_k_{n,k,d}^G`$ is far from commutative, as it has roughly $`1/N`$ times the size of the tensor algebra $`_k(^p)^k`$. We note however that the $`𝔖_k`$-invariants in $`_k_{n,k,d}^G`$ form a commutative subalgebra, isomorphic to $`S(E_\pi )^{𝔖_d}`$.

## 5 Quartic invariants of $`k`$ qubits

The case $`n=2`$ is of particular interest, as here the qudits are qubits, and this is the case being most discussed in quantum computation. Here we can give some precise non-asymptotic results for the first non-static case, namely $`_{2,k,4}^G`$. We put $`E=E_\pi =E_{[2,2]}`$. The proof of Proposition 4.1 easily gives

###### Corollary 5.1.

We have $`dim_{2,k,4}^G=\frac{1}{3}(2^{k1}+(1)^k)`$.

The first few values of $`dim_{2,k,4}^G`$, starting at $`k=1`$, are $`0,1,1,3,5,11,21,43`$. For $`k=2`$ and $`k=3`$ the unique (up to scalar) invariants are, respectively, the squared determinant $`P_{2,2}^2`$ and the Cayley hyperdeterminant $`H_{2,3}`$ (see \[G-K-Z\]). We note that the hyperdeterminant is very closely related to the relative tangle of $`3`$ entangled qubits discussd in \[C-K-W\].

It would be useful to study $`_{2,k,4}^G`$ as a representation of $`𝔖_k`$, where $`𝔖_k`$ acts by permuting the $`k`$ qubits. The $`𝔖_k`$-invariants in $`_{2,k,4}^G`$ are the $`(𝔖_k\text{ }\times G)`$-invariants in $`_{2,k,4}`$. These $`(𝔖_k\text{ }\times G)`$-invariant polynomials are very significant as they separate the closed orbits of the extended symmetry group $`𝔖_k\text{ }\times SL(2,)^k`$ acting on $`(^2)^k`$. We can compute the dimension of the $`𝔖_k`$-invariants as follows:

###### Proposition 5.2.

The dimension of the space of $`𝔖_k\text{ }\times G`$-invariants in $`_{2,k,4}`$ is $`M_k=\left[\frac{k}{6}\right]+r_k`$ where $`r_k=0`$ if $`k1mod6`$, or $`r_k=1`$ otherwise. Furthermore the algebra $`_k_{2,k,4}^{𝔖_k\text{ }\times G}`$ is the polynomial algebra $`[P_{2,2}^2,H_{2,3}]`$.

###### Proof.

We have isomorphisms $`_{2,k,4}^{𝔖_k\text{ }\times G}(E^k)^{𝔖_k\times 𝔖_3}S^k(E)^{𝔖_3}`$ since the representation of $`𝔖_4`$ on $`E`$ factors through $`𝔖_3`$. Thus the algebra $`_k_{2,k,4}^{𝔖_k\text{ }\times G}`$ identifies with $`S(E)^{𝔖_3}`$. Now $`S(E)^{𝔖_3}`$ is the algebra of $`𝔖_3`$-invariant polynomial functions on traceless $`3\times 3`$ diagonal matrices, and so is a polynomial algebra on the functions $`ATr(A^2)`$ and $`ATr(A^3)`$. These invariants correspond (up to scaling) to $`P_{2,2}^2`$ and $`H_{2,3}`$. The formula for the dimension follows easily. ∎

For instance, we have: $`M_1=0`$, $`M_k=1`$ for $`2k5`$, and $`M_6=2`$. We remark that by replacing $`S(E)^{𝔖_3}`$ by $`(E)^{𝔖_3}`$, it is easy to prove that the sign representation of $`𝔖_k`$ does not occur in $`(E^k)^{𝔖_4}`$ for any $`k2`$.

We can determine the $`𝔖_k`$-representation on $`_{2,k,4}^G`$ for small $`k`$ by explicit trace computations. For $`k=2`$ and $`k=3`$ we have the trivial $`1`$-dimensional representation. For $`k=4`$, we find $`_{2,4,4}^G`$ is the direct sum $`E_{[4]}E_{[2,2]}`$. The trivial representation $`E_{[4]}`$ of $`𝔖_4`$ is spanned by $`P_{2,4}^2`$, while the $`2`$-dimensional representation $`E=E_{[2,2]}`$ is spanned by the determinants $`\mathrm{\Delta }(ijkl)`$ introduced in \[B\]. Here $`(ijkl)`$ is a permutation of $`(1234)`$. Given a tensor state $`u(^2)^4`$, we can view it as an element $`v`$ of $`^4^4`$, where the first (resp. second) $`^4`$ is the tensor product of the $`i`$-th and $`j`$-th copies of $`^2`$ (resp. of the $`k`$-th and $`l`$-th copies). Then $`\mathrm{\Delta }(ijkl)(u)`$ is the determinant of $`v`$. As shown in \[B\], the $`\mathrm{\Delta }(ijkl)`$ span the representation $`E`$ of $`𝔖_4`$. The significance of the $`\mathrm{\Delta }(ijkl)`$ is that their vanishing describes the closure of the set of tensor states in $`(^2)^4`$ of rank $`3`$. For $`k=5`$ the representation $`_{2,5,4}^G`$ of $`𝔖_5`$ is $`E_{[5]}E_{[2,1,1,1]}`$.