text:
# Has the nonlinear Meissner effect been observed?

## Abstract

We examine recent high-precision experimental data on the magnetic field, $`𝐇`$, dependence of the penetration depth $`\lambda (H)`$ in $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ (YBCO) for several field directions in the $`ab`$ plane. In a new theoretical analysis that incorporates the effects of orthorhombic symmetry, we show that the data at sufficiently high magnetic fields and low temperatures are in quantitative agreement with the theoretical predictions of the nonlinear Meissner effect.

It is widely accepted that the symmetry of the order parameter (OP) in high temperature superconductors (HTSC’s) is at least predominantly $`d`$-wave, vanishing at nodal lines approximately ninety degrees apart in a quasi two-dimensional Fermi surface (FS). Many details of the OP in these materials remain quite unclear, however. Are the nodes exactly at right angles? Are they true nodes or only very deep minima, “quasinodes”? Addressing these and similar questions is important for obtaining clues about the nature of the superconductivity. They are particularly difficult to answer for the bulk OP, which may well differ  from the more easily observed surface state. To probe the bulk OP it is best to use electromagnetic techniques, since electromagnetic fields penetrate the sample over a depth of the order of the London penetration length $`\lambda `$ which is, in these materials, several orders of magnitude larger than the coherence length $`\xi `$ which characterizes the range of typical surface probes. Indeed, one of the early key results in support of bulk $`d`$-wave superconductivity was the measurement of the linear temperature dependence of the penetration depth in a high-purity YBCO single crystal. However, such data indicated only the existence of nodal lines without the angular resolution needed to identify their position. Consequently, intensive efforts to precisely determine the structure of the bulk OP have continued.

It was first pointed out years ago that nodes in the OP yield distinctive and measurable nonlinear effects in the field and angular dependence of the penetration depth when the superconductor is in the Meissner state. In subsequent theoretical work more emphasis was placed on the existence, due to this nonlinear Meissner effect (NLME), of a component of the diamagnetic moment normal to the applied field, and on the torque associated with this transverse component. These phenomena were deemed to be easier to measure than the changes in $`\lambda `$ itself. It was shown in this context that the NLME can be used to perform node spectroscopy, that is, not just to infer the existence of nodes, but to locate their positions on the FS, and to determine whether they are true nodes or not. Thus, the NLME is potentially a very important tool for the study of the pairing state in HTSC ’s, as well as other materials in the ever increasing list of those for which the proposed OP leads to an energy gap with nodes. Yet, the experimental situation is rather confusing. The best experimental effort to measure the transverse diamagnetic moment in YBCO was inconclusive. Subsequently, results for the magnetic field dependent penetration depth to a precision of $`0.1`$Å were reported. The NLME should be observable in such a high precision experiment, more precise than existing transverse moment measurements. Unfortunately, no theory of the NLME contribution to the penetration depth for orthorhomobic structures such as YBCO was available when Ref.  was written. Only very recently have the necessary calculations been performed. This has resulted in contradictory claims as to whether observed results are in agreement or not with NLME theory. Thus, a certain amount of skepticism has developed as to the observability of the NLME.

In this paper we show that measurements of the field dependent penetration depth $`\lambda (\psi ,H)`$ as a function of the angle $`\psi `$ that an applied field $`𝐇`$ in the $`ab`$ plane forms with the $`a`$-axis, must be analyzed very carefully. The anisotropy of the linear penetration depth tensor has a drastic effect on the NLME for $`\lambda (\psi ,H)`$. This may have been overlooked because the anisotropy effects in the transverse moment are known  to be relatively minor. One must also eliminate several other factors that may mask the signal at low fields and which are very difficult to account for theoretically. Thus we reanalyze here the best data available for the penetration depth in YBCO. We find that, although some questions remain, the low temperature data are quantitatively in agreement with theoretical expectations for the NLME in this material.

We focus here on YBCO, the most experimentally studied  HTSC in this context. Hence, the relevant material parameters are well known, thus reducing the uncertainty in the fitting procedures. We perform our analysis primarily on the most complete available high resolution data of Ref.  which includes results for four different directions of the applied field in the $`ab`$ plane.

The angular and field dependent increase in the penetration depth due to the NLME for materials with orthorhombic anisotropy of the YBCO type was first calculated in Ref. . The details will not be repeated here. The sample is assumed to have its larger faces parallel to the $`ab`$ plane (this is the case for crystals grown by the usual methods) and thickness large compared with the penetration depth. One has for the quantity $`\mathrm{\Delta }\lambda (\psi ,H)\lambda (\psi ,H)\lambda (\psi ,0)`$:

$$\mathrm{\Delta }\lambda (\psi ,H)=\frac{1}{6}\frac{H}{H_0}\lambda 𝒴(\psi ).$$
(1)
Here $`\lambda `$ is the geometric mean of the two in-plane principal values, $`\lambda _a`$ and $`\lambda _b`$, of the zero field penetration depth tensor, $`H_0`$ is a characteristic field of order $`\mathrm{\Phi }_0/\pi ^2\lambda \xi `$ ($`\mathrm{\Phi }_0`$ is the flux quantum), and $`𝒴`$ carries the angular dependence. The orthorhombicity, very important in this case, is incorporated into $`𝒴`$ through two parameters: one is the ratio $`\mathrm{\Lambda }\lambda _a/\lambda _b`$, and the other is the angle $`\alpha `$ that the Fermi velocity at the node located in the first quadrant forms with the $`a`$ direction. Because of the orthorhombic distortion of the FS, this angle does not have to exactly equal $`\pi /4`$ even for a pure $`d_{x^2y^2}`$ state, while the quantity $`\mathrm{\Lambda }`$, for YBCO, considerably exceeds unity. Here we take the zero field quantities $`\lambda _a`$ and $`\lambda _b`$ fixed at their experimental values ($`1050`$Å and $`1575`$Å), giving $`\mathrm{\Lambda }=1.5`$. In this case, the full expressions for $`𝒴(\psi )`$ simplify somewhat and can be written as:

$`𝒴(\psi )`$ $`=`$ $`{\displaystyle \frac{18\mathrm{\Lambda }}{2+\mathrm{\Lambda }}}\mathrm{cos}^2\alpha \mathrm{sin}\alpha \mathrm{cos}\psi \mathrm{sin}^2\psi +{\displaystyle \frac{2}{\mathrm{\Lambda }^2(1+2\mathrm{\Lambda })}}`$ (3)
$`\times `$ $`\mathrm{sin}^3\alpha \mathrm{cos}^3\psi \left[1+2\mathrm{\Lambda }+(4\mathrm{\Lambda }1)\left({\displaystyle \frac{\mathrm{tan}\psi }{\mathrm{tan}\psi _1}}\right)^{\frac{3\mathrm{\Lambda }}{\mathrm{\Lambda }1}}\right]`$ (4)
$`+`$ $`{\displaystyle \frac{2\mathrm{\Lambda }^2(2\mathrm{\Lambda }^210\mathrm{\Lambda }1)}{(2+\mathrm{\Lambda })(1+2\mathrm{\Lambda })}}\mathrm{cos}^3\alpha \mathrm{sin}^3\psi \left({\displaystyle \frac{\mathrm{tan}\psi }{\mathrm{tan}\psi _1}}\right)^{\frac{3}{\mathrm{\Lambda }1}},`$ (5)
$`\psi `$ $``$ $`[0,\psi _1],`$ (6)
$`𝒴(\psi )`$ $`=`$ $`{\displaystyle \frac{18}{1+2\mathrm{\Lambda }}}\mathrm{sin}^2\alpha \mathrm{cos}\alpha \mathrm{cos}^2\psi \mathrm{sin}\psi +2\mathrm{\Lambda }^2\mathrm{cos}^3\alpha \mathrm{sin}^3\psi ,`$ (7)
$`\psi `$ $``$ $`[\psi _1,{\displaystyle \frac{\pi }{2}}].`$ (8)
where the angle $`\psi _1`$ is given by $`\psi _1\mathrm{arctan}(\mathrm{tan}\alpha /\mathrm{\Lambda })`$.

Because of the orthorhombicity, the angular dependence of $`\mathrm{\Delta }\lambda (\psi ,H)`$ is quite different for $`\mathrm{\Lambda }=1.5`$ than that found for the tetragonal case ($`\mathrm{\Lambda }=1`$, $`\alpha =\pi /4`$). This is unlike the situation for the transverse magnetic moment, which from symmetry considerations vanishes when $`𝐇`$ is along the $`a`$ or the $`b`$ axis. Thus, moderate orthorhombic anisotropy induces only a relatively minor distortion in the curve between $`\psi =0`$ and $`\psi =\pi /2`$ since these points are so to speak, anchored. This is not the case for $`\mathrm{\Delta }\lambda (\psi ,H)`$: when the field is applied along $`\psi =0`$ the currents flow over a region of thickness determined by $`\lambda _b`$ while if $`\psi =\pi /2`$ the relevant skin depth is $`\lambda _a`$. The effect is nonzero, and different, in either case. This difference is compounded by the nonlinearity and the apparent overall symmetry of $`𝒴`$ is $`\pi `$ rather than $`\pi /2`$ even for moderate orthorhombicity. Failure to take this into account leads to erroneous conclusions concerning the angular dependence of $`\mathrm{\Delta }\lambda (\psi ,H)`$.

To analyze data in terms of Eqs. (1) and (1), additional considerations are needed. These expressions, indicating that $`\mathrm{\Delta }\lambda `$ is proportional to $`H`$, are valid at low temperature. Here, “low” temperature is a field dependent concept. The characteristic temperature separating the high and low $`T`$ regimes is  $`T^{}(H)\mathrm{\Delta }_0(H/H_0)`$, where $`\mathrm{\Delta }_0`$ is the gap amplitude. At any finite $`T`$, the validity of the above equations will break down at sufficiently small $`H`$. Further, the effect of impurities is not included. For the clean samples used in experiments this should affect only the small field results. The same is true of possible nonlocal effects. If they are present at all in this geometry, they would affect results at fields below 20 gauss. Ideally, one would like to take into account all of these effects by modifying the above formulas. However, it is not feasible at present to include all of these factors simultaneously in a reliable manner. It is therefore best to perform the analysis in a consistent manner in terms of data in the higher range of fields available, where these additional effects are all weak.

In Fig. 1 we show best straight line fits to the 1.2 K data of Ref.  for $`𝐇`$ along the $`a`$ and $`b`$ directions. All data in the range $`H>60`$ gauss are included in the fit. The cutoff of 60 gauss was chosen as the point below which deviations from a straight line begin and it will be shown below to lead to a consistent interpretation. The straight line does not intercept the origin of the original plot, which has to be shifted downward. This is as expected, since the experimental $`\mathrm{\Delta }\lambda `$ includes the previously mentioned temperature, impurity, and possible nonlocal effects which increase this quantity with respect to the theoretical, clean, zero temperature, local value. The shift is small, of order 1 Å, indicating that the sample is clean and any such spurious effects are small. The two slopes of the lines obtained from these fits are the quantities $`𝒴(0)\lambda /H_0`$ and $`𝒴(\pi /2)\lambda /H_0`$ respectively, where from Eq.(1), $`𝒴(0)=(2/\mathrm{\Lambda }^2)\mathrm{sin}^3\alpha `$, and $`𝒴(\pi /2)=2\mathrm{\Lambda }^2\mathrm{cos}^3\alpha `$. Thus, the ratio of the two slopes depends only on $`\alpha `$ since we have used the independent experimental value $`\mathrm{\Lambda }=1.5`$. We then determine the $`\alpha `$ that fits this ratio and subsequently find the characteristic field $`H_0`$ from either one of the slopes. The results are very sensible: we obtain $`\alpha =\pi /4+\pi /17`$ and $`H_0=5660`$ gauss. The value for the angle between the Fermi velocity at the node and the $`a`$-axis exceeds $`\pi /4`$ by a small amount, as one would expect for a pure $`d_{x^2y^2}`$ pairing state and a tight binding FS with a slight orthorhombic distortion. The value of the characteristic field is consistent with expectations and also with our cutoff choice for the field: in the range of fitting we have $`H/H_0>0.01`$. This means that in this field range, the characteristic temperature $`T^{}(H)`$ introduced above is of order of 4 to 12 K. Hence the 1.2 K data included in the fitting are in the low temperature regime and the procedure is consistent.

Up to now, we have, however, fit two quantities with two parameters, although the reasonable values obtained for these parameters are encouraging. To go beyond, we now use the obtained values of $`\alpha `$ and $`H_0`$ to plot the predicted slope of the high field data at $`|\psi |=\pi /4`$ without any additional parameters. This is done in Fig. 2. Experimental results for fields applied in the $`\psi =\pm \pi /4`$ directions are included. These results ought to be identical (even with the orthorhombic distortion) and their small discrepancy reflects systematic errors in the experiment. Nevertheless, the fit is excellent in the high field range. We also plot, (inset) with these parameter values, the predicted angular dependence of $`\mathrm{\Delta }\lambda `$ for YBCO. One can see that, with the orthorhombicity, this angular dependence differs considerably from that obtained for a tetragonal system, also plotted for comparison. The actual curve is not symmetric about $`\pi /4`$ and its maximum is much less pronounced than that for the tetragonal case, which is characterized by a factor of $`\sqrt{2}`$ between maxima and minima. Because of these differences, the attempt made by the authors of Ref.  to reconcile the angular dependence of their data with the theory of the NLME in a tetragonal system, had to fail.

In Fig. 3 we compare the theoretical results with other more recent data on YBCO at $`T`$=1.4 K. The parameters used are exactly the same as previously obtained. No new fits were performed. Results for the two directions available (field applied along the two principal axes) are shown. This data is in a more restricted, lower field range, and it has considerably more scatter than that of Bidinosti et al. All that can be said with certainty is that it is consistent with the NLME theory with the same parameter values.

In summary, the main result of the analysis presented here is that the best low temperature, high field, data on the nonlinear penetration depth in YBCO is in quantitative agreement, in its magnitude and angular and and field dependence, with the NLME theoretical expectations. Other data are also consistent with theory. Failure to observe the NLME in the transverse moment seems to be attributable to the actual precision in that experiment being just slightly less than what was in fact required.

Two remarks must be added: first, the crucial influence of the orthorhombic anisotropy in the angular dependence of $`\mathrm{\Delta }\lambda `$, which becomes very different (see Fig. 2) from that found for tetragonal symmetry, must be emphasized. Second, one sees the need to finesse the temperature, impurity and possibly other problems associated with smaller fields by obtaining and analyzing data at the highest possible fields below that of first flux penetration. Fortunately, this field is in the range 200-400 gauss for typical YBCO crystals.

The question of the temperature dependence of the results is less clear and needs further discussion: results obtained at 7 K for the same sample mainly discussed here are not substantially different from those at 1.7 K. With the characteristic temperature $`T^{}(H)`$ in the range estimated above, it can be that the high field results are not yet affected by the temperature at 7 K while those at low fields are dominated by largely temperature-independent impurity effects. Indeed, it appears that a straight line fit to the 7K data at the highest fields (see Fig. 4 of Ref. ) has a larger (in absolute value) vertical axis intercept than that for the 1.7 K data, which would be consistent with this scenario. Nevertheless, the weak temperature dependence of the data will remain a puzzle so long as a rigorous calculation including impurities, temperature, and possibly, nonlocal and other effects is not feasible. It is possible that these effects combine to yield a temperature dependence weaker than what the naive theory would predict.

Finally, our analysis indicates that there is no significant $`is`$ admixture to the $`d_{x^2y^2}`$ gap, since such an admixture would lead to quasinodes and to a considerable reduction in $`\mathrm{\Delta }\lambda `$. Furthermore, the nearness of $`\alpha `$ to $`\pi /4`$ is consistent with the absence of a real $`s`$ component as well.

It would be desirable to perform measurements of $`\mathrm{\Delta }\lambda `$ in YBCO at additional values of the angle $`\psi `$ to verify in more detail if the curve in the inset of Fig. 2 is indeed closely followed. The comments made here on the proper way to analyze experimental data should also be taken into account in any attempts to use the NLME to elucidate the pairing states of other suspected unconventional superconducting materials for which it is estimated that the sensitivity of present penetration depth measurements is sufficient to probe the NLME.

We thank C.P. Bidinosti and A. Bhattacharya for very useful information about their respective measurements. This work was supported in part by Petroleum Research Fund, administered by the ACS (at Minnesota) and by DARPA and ONR (N00140010028) at Maryland.


format: text

text:
# 1 Introduction

## 1 Introduction

The study of the properties of $`r`$-modes in rotating compact stars and their relevance to relativistic astrophysics has received considerable attention since the discovery (Andersson 1998, Friedman & Morsink 1998) that these modes are unstable to the emission of gravitational radiation. The $`r`$mode instability provides an explanation for the spin-down of rapidly rotating young neutron stars to Crab-like spin-periods and for the spin-distribution of millisecond pulsars and accreting neutron stars. Additionally, it is considered to be a strong source of continuous gravitational radiation (see Friedman & Lockitch 1999 for a review). Moreover, if $`r`$-modes induce differential rotation, their interaction with the magnetic field in neutron stars can enhance the toroidal magnetic field of the star (Rezzolla, Lamb & Shapiro 2000).

Before the instability can have an effect on the spin evolution of a young neutron star, the $`r`$-mode grows to an amplitude where it is saturated by nonlinear effects. Motivated by the absence of studies of such nonlinear saturation we performed a numerical study of $`r`$-mode hydrodynamical evolutions in rapidly rotating relativistic stars. We tried to elucidate what is the maximum amplitude such modes can reach, before nonlinear saturation (via hydrodynamical coupling) sets in (Stergioulas & Font 2000). The present contribution highlights the main findings of our study. We note that the saturation is most likely to set in on a hydrodynamical timescale, although it cannot be excluded that weak hydrodynamical couplings saturate the $`r`$-mode amplitude on longer timescales (but shorter than the growth timescale due to gravitational radiation reaction). However, at present, those long timescales cannot be achieved accurately in 3-D simulations, even with the largest available supercomputers.

## 2 Numerical framework

For our study we use a numerical code based on the 3-D `CACTUS` code (Font et al 2000, Alcubierre et al 2000), into which we implemented the 3rd order PPM method (Colella & Woodward 1984) for the hydrodynamics, and initial data for equilibrium and perturbed rapidly rotating relativistic stars (Stergioulas & Friedman 1995). In Font, Stergioulas & Kokkotas (2000) it was shown that the PPM method is suitable for long-term evolutions of rotating relativistic stars. We focus on a representative, rapidly (and uniformly) rotating model with gravitational mass $`M=1.63M_{}`$, equatorial circumferential radius $`R=17.25`$km and spin period $`P=1.26`$ms. We use the $`N=1.0`$ relativistic polytropic equation of state. Our simulations employ $`116^3`$ Cartesian grid-points, yielding a resolution of $`0.31`$km per zone.

The excitation of $`r`$-modes is achieved by perturbing the initial stationary model, adding a specific perturbation $`\delta v^i`$ to the contravariant components of the equilibrium 3-velocity, $`v^i`$. As there is no exact linear eigenfunction available in the literature that would correspond to an $`l=m=2`$ $`r`$-mode eigenfunction for rapidly rotating relativistic stars, we use an approximate eigenfunction, valid in the slow-rotation $`O(\mathrm{\Omega })`$ limit to the first post-Newtonian order (Lockitch 1999).

We also note that since $`r`$-modes are basically fluid modes we only evolve the hydrodynamical variables, keeping all spacetime variables fixed at their initial, unperturbed, values. The computational requirements for an accurate, coupled spacetime and hydrodynamical evolution, by far exceed current supercomputing resources.

## 3 Discussion

Figure 1 (left panel) displays the evolution of the axial velocity in the equatorial plane ($`v^z`$ along the $`y`$-axis) at a radius of $`r0.75R`$. The evolution is a superposition of several modes, the $`l=m=2`$ $`r`$-mode being the dominant component. The chosen amplitude of the eigenfunction is $`\alpha =1.0`$ (in the Newtonian limit, our definition of the amplitude agrees with that of Owen et al 1998). The perturbed star is evolved for more than 25ms (26 $`l=m=2`$ $`r`$-mode periods), during which the amplitude of the oscillation decreases due to numerical viscosity only. With much larger amplitudes the evolution is still similar to that in Figure 1, with no sign of nonlinear saturation of the $`r`$-mode amplitude on a dynamical timescale. Therefore, unless nonlinear saturation sets in at timescales much longer than the dynamical one, unstable $`r`$-modes could be driven to large amplitudes (of order one) by gravitational radiation reaction. An absolute upper limit on the $`r`$-mode amplitude is set by causality, requiring $`\sqrt{v_iv^i}<c`$.

A Fourier transform of the time-evolution, as a function of the frequency in the inertial frame, is shown in the right panel of Figure 1. It reveals that the initial data we are using, excite mainly the $`l=m=2`$ $`r`$-mode ($`r_2`$), with a frequency of 1.03 kHz and, with smaller amplitudes, several inertial modes ($`i_3`$, $`i_4`$, $`i_5`$) and higher harmonics. The frequencies of the various modes are the same at any given point inside the star, which indicates that the spectrum is a sum of discrete modes, in agreement with the conclusions of Lockitch, Andersson & Friedman (2000).

We have also investigated the possible appearance in our evolutions of the kinematical drift reported by Rezolla et al (2000). In a rotating star with a poloidal magnetic field, this kinematical drift may wind up the magnetic field lines and limit the effect of the $`r`$-mode instability. Our results do indicate that the perturbed star is rotating slower near the surface (compared to the unperturbed star) and that this drift scales roughly as $`\alpha ^2`$ for amplitudes of order $`\alpha 1.0`$, although its magnitude is significantly smaller than that estimated by Rezzolla et al (2000).

Our finding that gravitational radiation could drive unstable $`r`$-modes to a large amplitude implies that $`r`$-modes can easily melt the crust in newly-born neutron stars (Lindblom, Owen & Ushomirsky 1999), leaving the initial conclusions about the $`r`$-mode instability being a strong source of gravitational waves (see Friedman & Lockitch 1999) essentially unchanged.


text:
# Thermal convection for large Prandtl numbers

## Abstract

The Rayleigh-Benard theory by Grossmann and Lohse \[J. Fluid Mech. 407, 27 (2000)\] is extended towards very large Prandtl numbers $`Pr`$. The Nusselt number $`Nu`$ is found here to be independent of $`Pr`$. However, for fixed Rayleigh numbers $`Ra>10^{10}`$ a maximum around $`Pr2`$ in the $`Nu(Pr)`$-dependence is predicted which is absent for lower $`Ra`$. We moreover offer the full functional dependences of $`Nu(Ra,Pr)`$ and $`Re(Ra,Pr)`$ within this extended theory, rather than only giving the limiting power laws as done in ref. . This enables us to more realistically describe the transitions between the various scaling regimes, including their widths.

In thermal convection, the control parameters are the Rayleigh number $`Ra`$ and the Prandtl number $`Pr`$. The system responds with the Nusselt number $`Nu`$ (the dimensionless heat flux) and the Reynolds number $`Re`$ (the dimensionless large scale velocity). The key question is to understand the dependences $`Nu(Ra,Pr)`$ and $`Re(Ra,Pr)`$. In experiments, traditionally the Prandtl number was more or less kept fixed . However, the recent experiments in the vicinity of the critical point of helium gas  and of SF6  or with various alcohols  allow to vary both $`Ra`$ and $`Pr`$ and thus to explore a larger domain of the $`RaPr`$ parameter space of Rayleigh-Benard (RB) convection, in particular that for $`Pr1`$. While the experiments of Steinberg’s group  suggest a decreasing Nusselt number with increasing $`Pr`$, namely $`Nu=0.22Ra^{0.3\pm 0.03}Pr^{0.2\pm 0.04}`$ in $`10^9Ra10^{14}`$ and $`1Pr93`$, the experiments of the Ahlers group suggest a saturation of $`Nu`$ with increasing $`Pr`$ for fixed $`Ra`$, at least up to $`Ra=10^{10}`$ . The same saturation (at fixed $`Ra=610^5`$) is found in the numerical simulations by Verzicco and Camussi  and Herring and Kerr .

The large $`Pr`$ regime of the latest experiments has not been covered by the recent theory on thermal convection by Grossmann and Lohse (GL, ), which otherwise does pretty well in accounting for various measurements. In particular, it explains the low $`Pr`$ measurements of Cioni et al.  ($`Pr=0.025`$), the low $`Pr`$ numerics which reveal $`NuPr^{0.14}`$ for fixed $`Ra`$ , and the above mentioned experiments by Niemela et al.  and Xu et al. .

In the present paper we extend the GL theory in a natural way to the regime of very large $`Pr`$, on which no statement has been made in the original paper . We find $`Nu`$ to be independent of $`Pr`$ in that regime. We in addition present the complete functional dependences $`Nu(Ra,Pr)`$ and $`Re(Ra,Pr)`$ within the GL theory, rather than only giving the limiting power laws and superpositions of those as was done in . This enables us to more realistically describe the transitions between the various scaling regimes found already in .

Approach: To make this paper selfcontained we very briefly recapitulate the key idea of the GL theory, which is to decompose in the volume averages of the energy dissipation rate $`ϵ_u`$ and the thermal dissipation rate $`ϵ_\theta `$ into their boundary layer (BL) and bulk contributions,

$`ϵ_u`$ $`=`$ $`ϵ_{u,BL}+ϵ_{u,bulk},`$ (1)
$`ϵ_\theta `$ $`=`$ $`ϵ_{\theta ,BL}+ϵ_{\theta ,bulk}.`$ (2)
For the left hand sides the exact relations $`ϵ_u=\frac{\nu ^3}{L^4}(Nu1)RaPr^2`$ and $`ϵ_\theta =\kappa \frac{\mathrm{\Delta }^2}{L^2}Nu`$ are used, where $`\nu `$ is the kinematic viscosity, $`\kappa `$ the thermal diffusivity, $`L`$ the height of the cell, and $`\mathrm{\Delta }`$ the temperature difference between the bottom and the top plates. Next, the local dissipation rates in the BL and in the bulk (right hand sides of eqs. (1) and (2)) are modelled as the corresponding energy input rates, i.e., in terms of $`U`$, $`\mathrm{\Delta }`$, and the widths $`\lambda _u`$ and $`\lambda _\theta `$ of the kinetic and thermal boundary layers, respectively. For the thickness of the thermal BL we assume $`\lambda _\theta =L/(2Nu)`$ and for that of the kinetic one $`\lambda _u=L/(4\sqrt{Re})`$ as it holds in Blasius type layers ; as for the prefactor $`1/4`$ cf. , Sec. 4.3. For very large $`Ra`$ the laminar BL will become turbulent and $`\lambda _u`$ will show a stronger $`Re`$ dependence. Note that whereas the thermal BLs only build up at the top and bottom wall, the kinetic BL occurs at all walls of the cell and therefore the contribution of $`ϵ_{u,BL}`$ to $`ϵ_u`$ is larger than a simple minded argument would suggest. The two eqs. (1) and (2) then allow to calculate the two dependent variables $`Nu`$ and $`Re`$ as functions of the two independent ones $`Ra`$ and $`Pr`$.

Input rate modeling: The modeling of the dissipation rates on the rhs of eqs. (1) and (2) is guided by the Boussinesq equations. Depending on whether the BL or the bulk contributions are dominant, one gets different expressions on the rhs of eqs. (1), (2) and thus different relations for $`Nu`$, $`Re`$ vs. $`Ra`$, $`Pr`$, defining different main regimes, see ref.  and fig. 1.

The rhs thermal dissipation rates depend on whether the kinetic BL of thickness $`\lambda _u`$ is within the thermal BL of thickness $`\lambda _\theta `$ ($`\lambda _u<\lambda _\theta `$, small $`Pr`$) or vice versa ($`\lambda _u>\lambda _\theta `$, large $`Pr`$). The line $`\lambda _u=\lambda _\theta `$, corresponding to $`Nu=2\sqrt{Re}`$, splits the phase diagram into a lower (small $`Pr`$) and an upper (large $`Pr`$) part, which we shall label by “$`\mathrm{}`$” and “$`u`$”.

We first consider $`\lambda _u<\lambda _\theta `$ (i.e., small $`Pr`$, regime “$`\mathrm{}`$”). Then (see )

$$ϵ_{u,bulk}\frac{U^3}{L}\frac{\nu ^3}{L^4}Re^3,$$
(3)
$$ϵ_{u,BL}\nu \frac{U^2}{\lambda _u^2}\frac{\lambda _u}{L}\frac{\nu ^3}{L^4}Re^{5/2},$$
(4)
$$ϵ_{\theta ,bulk}\frac{U\mathrm{\Delta }^2}{L}\kappa \frac{\mathrm{\Delta }^2}{L^2}PrRe,$$
(5)
$$ϵ_{\theta ,BL}\kappa \frac{\mathrm{\Delta }^2}{L^2}(RePr)^{1/2}.$$
(6)
The last expression is concluded  from the heat transfer equation $`u_x_x\theta +u_z_z\theta =\kappa _z^2\theta `$, which implies $`U/L\kappa /\lambda _\theta ^2`$ giving $`Re^{1/2}Pr^{1/2}Nu`$.

If now larger $`Pr`$ are considered, the kinetic boundary will eventually exceed the thermal one, $`\lambda _u>\lambda _\theta `$, upper range “$`u`$”. The relevant velocity at the edge between the thermal BL and the thermal bulk now is less than $`U`$, about $`U\lambda _\theta /\lambda _u`$. To describe the transition from $`\lambda _u`$ being smaller to being larger than $`\lambda _\theta `$ we introduce the function $`f(x)=(1+x^n)^{1/n}`$ of the variable $`x_\theta =\lambda _u/\lambda _\theta =c_\theta Nu/2\sqrt{Re}`$, $`f`$ being 1 in the lower range “$`\mathrm{}`$” (small $`Pr`$) and $`1/x_\theta `$ in “$`u`$” (large $`Pr`$), respectively. $`c_\theta `$ is about $`1`$. The relevant velocity then is $`Uf(x_\theta )`$. We take $`n=4`$ to characterize the sharpness of the transition.

This generalizes (5), (6) to

$$ϵ_{\theta ,bulk}\kappa \frac{\mathrm{\Delta }^2}{L^2}PrRef(c_\theta Nu/2\sqrt{Re}),$$
(7)
$$ϵ_{\theta ,BL}\kappa \frac{\mathrm{\Delta }^2}{L^2}\sqrt{PrRef(c_\theta Nu/2\sqrt{Re})},$$
(8)
while $`ϵ_{u,bulk}`$ and $`ϵ_{u,BL}`$ are still given by (3),(4). Introducing (7),(8),(3),(4) into (1),(2) leads to the $`Ra`$,$`Pr`$ dependences of $`Nu`$,$`Re`$ in the upper regime “$`u`$”. The pure power laws $`Nu(Ra,Pr)`$ and $`Re(Ra,Pr)`$ in both the lower$`\mathrm{}`$” ($`\lambda _u<\lambda _\theta `$, small $`Pr`$) and the upper$`u`$” ($`\lambda _u>\lambda _\theta `$, large $`Pr`$) regimes are summarized in table I.

Very large $`Pr`$ regime: We now extend the theory to very large Prandtl numbers. As long as $`Ra>Ra_c=1708`$ there still is wind, even for very large $`Pr`$, as $`Ra_c`$ is independent of $`Pr`$. However, to keep $`Ra`$ fixed, one has to increase the temperature difference $`\mathrm{\Delta }`$ with increasing $`\nu `$ or $`Pr`$. $`Pr1`$ will lead to a smaller and smaller large scale wind $`Re`$. The flow will eventually become laminar throughout the cell. $`\lambda _u`$ can no longer continue to increase according to $`\lambda _uRe^{1/2}`$ with decreasing $`Re`$, but will saturate to a constant value of order $`L`$. This is the central new point of the very large $`Pr`$ regime. In  we had assumed that this happens at about $`Re=50`$. This corresponds to a saturation at $`\lambda _u=L/4\sqrt{Re}=L/28`$. To model the smooth transition to the very large $`Pr`$ regime beyond the line $`Re=50`$, i.e., if $`\lambda _u=L/4\sqrt{Re}`$ approaches $`\lambda _u=L/28`$ we use the crossover function $`g(x)=x(1+x^n)^{1/n}`$ of the crossover variable $`x_L=28\lambda _u/L=7/\sqrt{Re}`$, and again $`n=4`$. The function $`g`$ increases linearly, $`g(x_L)=x_L`$, below the transition ($`x_L`$ small) and is 1 in the very large $`Pr`$ regime with $`Re50`$. In the above modelings for the local dissipations rates we have to replace each $`\lambda _u`$ by $`g(x_L)L/28`$.

The resulting formulae, given momentarily, will lead, depending on $`Ra`$, to three new regimes, valid for very large $`Pr`$, denoted as $`I_{\mathrm{}}^<`$, $`I_{\mathrm{}}^>`$, and $`III_{\mathrm{}}`$, see Fig. 1 and Table I.

While eq. (3) for $`ϵ_{u,bulk}`$ is assumed to still hold in the very large $`Pr`$ range (where $`ϵ_{u,bulk}`$ of course hardly contributes to $`ϵ_u`$ due to the large extension of the kinetic BLs), eqs. (4), (5), and (6) have to be generalized. First generalize (4) for $`ϵ_{u,BL}`$,

$$ϵ_{u,BL}\nu \frac{U^2}{g(x_L)L^2}\frac{\nu ^3}{L^4}\frac{Re^2}{g(7/\sqrt{Re})}.$$
(9)
Next $`ϵ_{\theta ,bulk}`$: Above the wind velocity $`U`$ in (5), which sets the time scale of the stirring, has already been generalized to $`Uf(\lambda _u/\lambda _\theta )`$. This equals $`U`$ itself in “$`\mathrm{}`$” and $`U\lambda _\theta /\lambda _u`$ in the “$`u`$” regimes. Now in addition the explicit $`\lambda _u`$ is to be replaced by $`g(x_L)L/28`$, i.e.,

$$ϵ_{\theta ,bulk}\kappa \frac{\mathrm{\Delta }^2}{L^2}PrRef\left(\frac{Nu}{14}g\left(\frac{7}{\sqrt{Re}}\right)\right).$$
(10)
This simplifies for large enough $`Ra`$ (therefore large $`f`$-argument) and very large $`Pr`$ (thus large $`g`$-argument) to

$$ϵ_{\theta ,bulk}\kappa \frac{\mathrm{\Delta }^2}{L^2}\frac{PrRe}{Nu}.$$
(11)
Inserting (9) and (11) into the rhs of (1) and (2) leads to the new power laws describing the heat flux and the wind velocity in the regime $`III_{\mathrm{}}`$ beyond $`III_u`$, cf. Figure 1,

$$NuRa^{1/3}Pr^0,ReRa^{2/3}Pr^1.$$
(12)
Finally $`ϵ_{\theta ,BL}`$: In the thermal boundary layer range beyond $`I_{\mathrm{}}`$, relevant for medium $`Ra`$, eq.(6) stays valid, because its derivation did not involve $`\lambda _u`$ and also $`f=1`$. The range $`I_{\mathrm{}}^<`$ has to be described by eqs.(1),(2),(9) (with $`g=1`$), and (6), resulting in the same power laws as in regime $`III_{\mathrm{}}`$, i.e., eqs. (12). For $`Pr`$ values above regime $`I_u`$ eq.(6) no longer holds. It originated from the heat transport equation. There we have to use now $`Uf(x_\theta )`$ instead of merely $`U`$. The balance from the heat transfer equation then reads $`Uf(x_\theta )/L\kappa /\lambda _\theta ^2`$. In the $`f`$-argument $`x_\theta =\lambda _u/\lambda _\theta `$ the kinetic BL width $`\lambda _u`$ and therefore the crossover function $`g`$ appears, leading to

$$PrRef\left(\frac{Nu}{14}g(\frac{7}{\sqrt{Re}})\right)Nu^2.$$
(13)
In the very large $`Pr`$ regime (where $`g(x_L)=1`$) above $`I_u`$ (where $`f(x_\theta )=x_\theta ^1`$) one obtains from (13) the relation $`(RePr)^{1/3}Nu`$, valid in $`I_{\mathrm{}}^>`$. Together with (1), (2), and (9) one derives the new scaling laws in the interior of $`I_{\mathrm{}}^>`$

$$NuRa^{1/5}Pr^0,ReRa^{3/5}Pr^1.$$
(14)
The scaling behavior $`NuRa^{1/5}`$ has earlier been suggested by Roberts . Note that in all three very large $`Pr`$ regimes $`Nu`$ does not depend on $`Pr`$. Furthermore the Constantin-Doering  upper bound $`NuconstRa^{1/3}(1+\mathrm{log}Ra)^{2/3}`$, holding in the limit $`Pr\mathrm{}`$, is strictly fulfilled.

$`Nu`$ and $`Re`$ in the whole parameter plane: Plugging now the general expressions for the local dissipation rates (3), (9), (10), and (13) into the balance eqs. (1) and (2) finally results in

$$NuRaPr^2=c_1\frac{Re^2}{g(7/\sqrt{Re})}+c_2Re^3,$$
(15)
$$Nu=c_3Re^{1/2}Pr^{1/2}\left[f\left(c_\theta \frac{Nu}{14}g(\frac{7}{\sqrt{Re}})\right)\right]^{1/2}+c_4PrRef\left(c_\theta \frac{Nu}{14}g(\frac{7}{\sqrt{Re}})\right).$$
(16)
Here we have added dimensionless prefactors where appropriate to complete the modeling of the dissipation rates. In ref.  $`c_1`$ through $`c_5`$ were already adopted to Chavanne et al.’s experimental data . The result was $`c_1=1028`$, $`c_2=9.38`$, $`c_3=1.42`$, $`c_4=0.0123`$, $`1/c_\theta =c_5=2.0`$, see eqs. (4.4) through (4.9) of . Note that the $`c_i`$ may depend on the aspect ratio and are not universal.

The set of eqs.(15) and (16) is the second main result of this paper. It allows to calculate $`Nu(Ra,Pr)`$ and $`Re(Ra,Pr)`$ in the whole $`RaPr`$ parameter space, including all crossovers from any regime to any neighboring one. Technically, eq. (15) is solved for $`Nu(Ra,Re,Pr)`$ which then is inserted into eq. (16), leading to one implicit equation for $`Re(Ra,Pr)`$.

All limiting, pure scaling regimes which can be derived from eqs. (15) and (16) are listed in table I which extends table II of ref. , now including the three new regimes $`I_{\mathrm{}}^<`$, $`I_{\mathrm{}}^>`$, and $`III_{\mathrm{}}`$. The corresponding phase diagram is shown in figure 1, completing that of  towards very large $`Pr`$.

Though in the phase diagram we have drawn lines to indicate transitions between the regimes, defined by either $`ϵ_{u,BL}=ϵ_{u,bulk}`$ or $`ϵ_{\theta ,BL}=ϵ_{\theta ,bulk}`$ or $`\lambda _u=\lambda _\theta `$ or $`\lambda _u=L/28`$, we note that the crossovers are nothing at all but sharp. All transitions are smeared out over broad ranges, the more, the more similar the scaling exponents of the neighboring regimes are.

Discussion of $`Nu(Pr)`$ and $`Nu(Ra)`$: The functions $`Nu(Pr)`$ (for fixed values of $`Ra`$) and $`Nu(Ra)/Ra^{1/4}`$ (for fixed values of $`Pr`$) resulting from eqs. (15) - (16) are shown in figure 2. Indeed, for $`Ra10^610^9`$ the Nusselt number saturates with increasing $`Pr`$ and the regime $`I_u`$ with $`NuPr^{1/12}`$ is suppressed, in full agreement with the experimental and numerical results . However, for larger $`Ra`$ beyond $`10^{11}`$ the theory predicts a maximum for the curve $`Nu(Pr)`$. The decreasing branches of the curve for $`Ra10^{12}`$ and $`Pr50`$ may be consistent with above mentioned experimental results by Ashkenazi and Steinberg . This observation may resolve the apparent descrepancy between the Ahlers et al. and the Steinberg et al. data: Simply different regimes in the $`RaPr`$ phase space are probed.

The importance of transition ranges is highlighted in fig. 2b. E.g., in regime $`II_l`$ with $`NuRa^{1/5}`$ (for fixed $`Pr`$) the corresponding scaling exponent $`1/51/4=1/20`$ only becomes observable for very small $`Pr10^310^4`$. At $`Pr=10^2`$ the roughly four decades of regime $`II_l`$ (see figure 1) between regime $`I_l`$ and $`IV_l`$, which both have a stronger $`Ra`$ dependence of $`Nu`$, are not sufficient to reveal the scaling exponent. And at $`Pr=10^1`$ the roughly 2.5 decades of regime $`II_l`$ are nowhere sufficient to lead to a local scaling exponent $`dNu/dRa`$ smaller than $`1/4`$!

Similarly, for $`Pr=1`$ only for very large $`Ra\text{ }\stackrel{>}{}\text{ }10^{15}`$ pure scaling $`NuRa^{1/3}`$ is revealed. For smaller $`Ra`$ regimes $`I_u`$ and $`I_l`$ with $`NuRa^{1/4}`$ strongly contribute, resulting in an effective local scaling exponent (increasing with $`Ra`$) in the range between $`0.28`$ (“$`2/7`$”) and $`0.31`$, just as observed in experiment .

Direct comparison to experiment: Finally, let us point out that the present approach simplifies the comparison with experimental data: In experiment, it is hard to vary either $`Ra`$ or $`Pr`$ over many decades and at the same time to keep the other variable fixed. So most measurements are done along curved lines in the phase space figure 1, mixing the $`Ra`$ and $`Pr`$ dependences. Now eqs. (15) and (16) allow to calculate $`Nu`$ and $`Re`$ along such a curve $`Pr(Ra)`$ given by the experimental restrictions, which connect $`Ra`$ and $`Pr`$. E.g., figure 3 shows $`Nu(Ra,Pr(Ra))`$ with $`Pr(Ra)`$ as resulting from the Niemela et al.  and the Chavanne et al.  experiments.

The present theory suggests that beyond $`Ra10^{12}`$ the Nusselt number in the Chavanne et al. experiment  should be smaller than in the Niemela et al. experiment , in contrast to what is found. This result is also in contrast to our former speculation of ref. , that the strong $`NuRa^{3/7}`$ dependence of regime $`III_u`$, to which the Chavanne et al. data are closer, would account for their experimental findings. What we had overseen is that the larger $`Pr`$ numbers of the Chavanne et al. experiment (for fixed $`Ra`$) and the $`NuPr^{1/7}`$ dependence of regime $`III_u`$ overcompensates the strong $`Ra`$ dependence of $`Nu`$ for the Chavanne et al. data close to regime $`III_u`$.

This present theory cannot resolve the paradox between the Chavanne et al. and the Niemela et al. data. Possibly, different temperature boundary conditions have been applied. Possibly, Chavanne et al. have already observed the transition from a laminar to a turbulent kinetic BL (dotted line in the phase diagram fig. 1).

Acknowledgement: This research work was prompted by the experiments of G. Ahlers . We thank him for sharing his results with us prior to publication and for various discussions.


text:
# XMM-Newton monitoring of X-ray variability in the quasar PKS~0558–504 Based on observations with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).

## 1 Introduction

Active Galactic Nuclei (AGN) are variable in every observable wave band. The X-ray flux exhibits variability on shorter time scales than any other energy band, indicating that the emission occurs in the innermost regions of the central engine. Therefore, a study of the X-ray variability is a powerful tool providing upper limits on the sizes of the emitting regions and the masses of central black holes and it allows one to probe the extreme physical processes operating in the inner parts of the accretion flow close to the event horizon. Although X-ray variability has been observed in AGN for more than two decades, its origin is still poorly understood. Narrow-Line Seyfert 1 galaxies (NLS1) often display rapid, large amplitude X-ray variability as well as extreme long term changes (Forster & Halpern forst (1996), Boller et al. (1997), Brandt et al. (1999)), and therefore they represent the ideal class for an X-ray temporal analysis.

PKS~0558–504 ($`z=0.137,m_\mathrm{B}=14.97`$) is one of the few radio-loud NLS1. In X-rays it is characterized by a steep spectrum \[the photon index $`\mathrm{\Gamma }`$ ranges between 2.2 and 3.1, obtained from different instruments and energy bands (a detailed spectral analysis of XMM-Newton observations has been performed in a separate paper, O’Brien et al. (2001))\], high luminosity \[$`(25)\times 10^{45}\mathrm{erg}\mathrm{s}^1`$, with $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1,q_0=0.5`$\], and strong variability. A Ginga observation (Remillard et al. (1991)) showed an increase of the X-ray flux by $`67\%`$ in 3 minutes, implying the presence of a relativistic flare. Recent ROSAT High Resolution Imager (HRI) observations (Gliozzi et al. (2000)) confirmed the presence of strong and persistent X-ray variability, which suggests the presence of a rotating black hole, and ruled out external contributions to the high luminosity and variability from other nearby sources.

In this paper we report the results of four XMM-Newton observations of PKS~0558–504 taken with the European Photon Imaging Camera (EPIC) PN, which is the best suited instrument on board the ESA satellite for timing analysis purposes, due to its high time resolution. In Sect.2 we describe the observations and data reduction. The long and short term X-ray variability are discussed in Sect.3 and Sect.4, respectively. Sect.5 contains the main conclusions and a summary.

## 2 Observations

PKS~0558–504 was observed with the EPIC PN camera during orbits 30, 32, 42 and 84 (corresponding to February 7 and 10, March 2 and May 24, respectively). The camera was operated in different observing modes and with various filters. In the Full Frame mode all the pixels of all CCDs are read out with time resolution of 70 ms and the full field of view ($``$ 30′) is covered. In the Partial Window modes (Large Window and Small Window) a two dimensional imaging readout over part of the CCDs array is performed with increased time resolution (48 ms and 6 ms, respectively). The primary reason for using filters is that the EPIC CCDs are not only sensitive to X-ray photons, but also to IR, visible and UV light. Therefore for objects with high optical to X-ray flux ratio, the X-ray signal might be contaminated by those photons. This is not the case of PKS~0558–504, for which different filters were used only for calibration purposes. Table 1 provides a summary of the individual observations. The data reduction and analysis were performed using the last release of the XMM-Newton SAS (Science Analysis Software) package. From the available observations, the last reprocessed data sets without any indication of timing problems and only single events (i.e. events where the charge released by a photon is contained in a single CCD pixel) have been analyzed. The main reason for using only singles is related to the energy dependence in the relative fraction of “multiple” events. To quantify their contribution, a test with all valid events has been performed. The main effect is an increase of the count rate by $`10\%`$ and variations of the hardness ratio by $`20\%`$, which would introduce additional uncertainties. Light curves were obtained by extracting photons from circular region of radius of $``$ 33″ around the source center and subtracting the background taken from a source-free region (which amounted to $`1\%`$). After confirming the lack of significant variability on time scales less than minutes, the photons were binned into 200 s intervals providing a good signal to noise ratio. PKS~0558–504 was chosen as a CAL/PV target with the purpose of investigating the photon pile-up in the EPIC cameras. With this aim we compared the spectra from the central pixels with those from the outer parts of the PSF and found no significant indication for photon pile-up.

## 3 Long term variability

In Fig. 1 we show the total light curve for PKS~0558–504 from February 7 to May 24 2000. Each data point is an average count rate over the individual observation in the range 0.2-10 keV. The indicated error of $`1\sigma `$, represents the dispersion around the mean and gives an idea of the variability within each observation. The error bar is particularly large for the first data point because we merged two exposures (with very different mean count rates) taken during orbit 30 nearly one hour apart. The total light curve is characterized by a moderate variability, if compared with the ROSAT light curve where the count rate varied by a factor 5 in less than 3 days (Gliozzi et al. (2000)). In the EPIC PN light curve the most prominent variation in the count rate is by a factor of 2 over a period of 20 days.

We sought evidence for spectral variability by computing the hardness ratios as a function of the total 0.2-10 keV flux. As hardness ratio we define the difference between the count rates in the $`1\mathrm{keV}\mathrm{E}10\mathrm{keV}`$ band and those in the $`0.2\mathrm{keV}\mathrm{E}<1\mathrm{keV}`$ band, divided by the 0.2-10 keV count rate. These ranges where chosen to provide a good signal to noise in both the soft and the hard light curves. The hardness ratios versus the total count rate with data from each observation binned in 2000 s intervals is shown in Fig. 2. The presence of a positive correlation was quantitatively tested by performing a linear least square fit, which confirmed the result at 5$`\sigma `$ confidence level. On long time scales the spectral variability is correlated with flux variations such that the spectrum becomes harder when the count rate increases.

Another interesting correlation was found between the variability in the hard band and the total count rate . The variability was studied using the excess variance (Nandra et al. (1997)), which is obtained by computing the variance of the overall light curve, subtracting the variance due to measurement error and dividing by the mean squared. Care must be taken in interpreting the excess variance of different observations as it is related to the length and the sampling of the time series (e.g. Leighly 1999a ). However this drawback is circumvented in our analysis by considering the ratio between the excess variance in the hard band and that in the soft range during each individual observation. Fig. 3 shows a clear correlation between the hard to soft excess variance ratio and the mean count rate in the 0.2-10 band: the contribution of the hard band to the variability increases and becomes dominant as the mean count rate increases. From a direct inspection of the individual values of the excess variance in the soft and hard band, this behaviour turns out to be due more to the actual increase of the hard variability than to a depletion of the soft variability level.

Based on ROSAT and ASCA data Fiore et al. (1998) and Leighly (1999b ) reported the presence of correlations between variability and X-ray luminosity and the steepness of the X-ray spectrum. A first analysis seems to indicate that the XMM data of PKS~0558-504 are consistent with the previous results, however a detailed study of such correlations is beyond the scope of this paper. A spectral study of the XMM observations of PKS~0558-504 can be found in O’Brien et al. (2001).

## 4 Short term variability

The PKS~0558–504 EPIC PN light curves for each individual observations, binned into 200 s intervals, are shown in Fig. 4. Strong and persistent X-ray flux variability is observed during each observation. The most prominent increase in the count rate is observed during orbit 30 (February 7), when the 0.2-10 keV flux increased by a factor 2 in less than 4 hours. On the other hand the most extreme temporal variation in the count rate occurred on February 10, when the count rate decayed by $`18\%`$ in 1500 seconds. This value is in agreement with that measured by the EPIC MOS during orbit 84 (O’Brien et al. (2001)), which led to the calculation of a lower bound for the radiative efficiency which exceeds the theoretical maximum for a non rotating black hole, under the assumption that the photon diffusion through a spherical mass of accreting matter is dominated by Thomson scattering.

A visual examination of the individual light curves indicates that all of them show a similar long term variability pattern with a common rise time of about 2 hours and, in some cases, small amplitude flares superimposed. The best method to quantify time variability without the problems encountered in the traditional Fourier analysis technique in case of unevenly sampled data is a structure function analysis (e.g. Simonetti et al. (1985), Hughes et al. (1992)). The first-order structure function is the mean deviation for data points separated by a time lag $`\tau `$, $`SF(\tau )=[F(t)F(t+\tau )]^2`$. One of the most useful features of the structure function is its ability to discern the range of time scales that contribute to the variations in the data set: the characteristic time scales of the variability are identified by the maxima and slope changes in the $`\tau SF`$ plane. For a stationary random process the structure function reaches a plateau state for lags longer than the longest correlation time scale. If a light curve contains cycles of period $`P`$, the $`SF`$ will rise to a maximum at $`\tau =P/2`$ (Smith et al. (1993)).

Fig. 5 shows that all the observations have either a relative or absolute maximum around 2 hours. This means that a common typical time scale (which probably reflects the similar rise time in the light curves) characterizes all observations. In addition, during orbits 30 and 84 a further characteristic time scale around 3 hours seems to be present. A first confirmation for the presence of quasi-periodicity in the temporal behaviour of PKS~0558–504 comes from a periodogram analysis, which yielded a strong signal at 2.3 hr. However, further longer observations are necessary for a firm confirmation.

In order to seek evidence for spectral variability and for its origin on short time scales we have split each of the broad band light curves into a hard and a soft light curve and plotted the hardness ratios versus the time. While during orbits 32 and 84 (when the flux was at intermediate values) no significant spectral variability was found, orbits 30 and 42 present an interesting and somewhat puzzling spectral and temporal behaviour shown in Fig. 6 and Fig. 7, respectively. The mean count rates and the excess variances in the soft, hard and broad bands are summarized in Table 2.

During the first exposure of orbit 30 the variability, quantified by the excess variance $`\sigma _{\mathrm{rms}}^2`$ is strongly dominated by the soft photons ($`4.6\times 10^2\mathrm{vs}\mathrm{.\; 4.0}\times 10^3`$, for the soft and the hard band respectively), with the soft flux steadily increasing and the hard fluctuating around a constant value. As a consequence, while the total flux increases the hardness ratio decreases in time, with a slope of $`(2.9\pm 0.7)\times 10^2`$/hr obtained from a linear least square fit. A totally different trend is observed during the second exposure, starting just one hour later. In this case the variability is dominated by the hard flux ($`\sigma _{\mathrm{rms}}^2:5.2\times 10^2\mathrm{and}2.7\times 10^2`$, for the hard and soft band, respectively), which increases faster than the soft emission and seems to peak before it. The hardness ratio is steadily increasing with the total flux and the time (slope: $`(1.2\pm 0.4)\times 10^2`$/hr). A similar spectral behaviour is observed during orbit 42: the variability is strongly dominated by the hard flux ($`\sigma _{\mathrm{rms}}^2:9.7\times 10^2\mathrm{vs}\mathrm{.\; 2.7}\times 10^2`$) and the spectrum becomes harder when the total count rate increases (slope: $`(1.3\pm 0.3)\times 10^2`$/hr).

On the basis of the five exposures available (two for orbit 30 and one each for orbits 32, 42 and 84), the spectral variability on short time scales of PKS~0558–504 can be summarized in the following way: when the average broad band count rate is at a low level (the mean value is around $`8\mathrm{cts}\mathrm{s}^1`$), the variability is dominated by the soft photons and the hardness ratio decreases as the total flux raises. Conversely, when the average count rate is above a certain level ($`12.6\mathrm{cts}\mathrm{s}^1`$ and $`16.3\mathrm{cts}\mathrm{s}^1`$ for the second exposure of orbit 30 and orbit 42, respectively), the variability is dominated by the hard band and the spectrum becomes harder when the count rate increases. At intermediate values of count rate (the mean is around $`10\mathrm{cts}\mathrm{s}^1`$ both for orbit 32 and 84), the variability in the two bands is very similar and no significant spectral variability is detected.

## 5 Summary and conclusions

We have presented the results of the EPIC PN observations of the NLS1 PKS~0558–504 taken during the Cal/PV phase. We found:

* During the XMM pointings PKS~0558–504 showed a moderate but persistent variability both on long (months) and short (minutes and hours) time scales.
* Strong correlations between hardness ratio and total flux have been found. When the mean flux is above a certain threshold, the spectrum becomes harder as the count rate increases. When the flux is below a certain level, the spectral behaviour is in the opposite sense. At intermediate count rates no significant spectral variability is observed.
* A further strong correlation between the broad band flux and the variability (expressed in terms of excess variance) in the hard energy band has been discovered: the contribution to the total variability of the hard band increases and becomes dominant as the mean count rate raises.
* In all the observations the variability is characterized by smooth similar modulations, with typical time scale $`2.2`$ hr, determined with the structure function analysis.
* The most extreme temporal variation of the count rate implies a radiative efficiency slightly larger than the theoretical maximum for accretion onto a Schwarzschild black hole. It is important to point out that this is a strictly lower limit, since during our analysis only single events (which represent $`70\%`$ of the total X-ray photons collected) have been used.

A detailed model of PKS~0558–504 which accounts for all timing and spectral features discovered is beyond the scope of this paper. Therefore, we limit ourselves to some general considerations. The most intriguing result is the tentative detection of a common characteristic time scale of $``$ 2 hours. However a confirmation of this characteristic time scale has to await for longer observations, to sample a higher number of cycles. A more detailed analysis to assess its significance and robustness (e.g. as in Boller et al. (2001)) is deferred to a follow-up publication. This time scale might be associated to periodic phenomena occurring in the accretion disk, whose contribution seems to be dominant also on the basis of the spectral analysis of EPIC MOS data (O’Brien et al. (2001)). Two possible explanations which need to be quantitatively tested are: i) the presence of a hot spot orbiting the black hole and ii) the precession of the inner part of the accretion disk caused by misalignment of black hole and disk rotational axes. An alternative tempting hypothesis involves the presence of a jet: as PKS~0558–504 is a radio-loud object, beamed emission from a jet could partly contribute to the brightness and variability in X-rays. Radio-loud quasars are known to have flatter X-ray spectral indices than radio-quiet ones, and they are also known to be brighter X-ray sources (e.g. Yuan et al. (1998)). A widely accepted explanation for this fact is that the X-ray emission includes an additional harder component associated with the radio jet. In this scenario a naive interpretation of the temporal and the spectral variability associated with changes of the broad band flux, involves the presence of either a precessing jet or moving knots with helical trajectories in a relativistic magnetized jet (Camenzind & Krockenberger came (1992)). At this point it is worth noticing that also in one other of the few radio-loud NLS1, RX J0134.2–4258, the spectral variability (characterized, as in our case, by a hardening of the spectrum during the flux increase) was tentatively interpreted as due to the jet emission (Grupe et al. (2000)).

The exciting results from the temporal analysis of a few short exposures demonstrate the extraordinary capabilities of XMM, in particular in light of the long observations possible in virtue of the highly eccentric orbit.

###### Acknowledgements.

The XMM - Newton project is supported by the Bundesministerium für Bildung und Forschung / Deutsches Zentrum für Luft- und Raumfahrt (BMBF/DLR), the Max-Planck Society and the Heidenhain-Stiftung. EPIC was developed by the EPIC Consortium led by the Principal Investigator, Dr. M.J.L. Turner. The consortium comprises the following Institutes: University of Leicester, University of Birmingham, (UK); CEA/Saclay, IAS Orsay, CESR Toulouse, (France); IAAP Tuebingen, MPE Garching, (Germany); IFC Milan, ITESRE Bologna, IAUP Palermo, (Italy). EPIC is funded by: PPARC, CEA, CNES, DLR and ASI. MG acknowledges support from the European Commission under contract number ERBFMRX-CT98-0195 (TMR network “Accretion onto black holes, compact stars and protostars”).


text:
# Lattice polytopes with distinct pair-sums

Let $`𝒫`$ be a lattice polytope in $`^n`$, the convex hull of a finite set in $`^n`$, and let

$$\left(𝒫\right):=𝒫^n=\{v_1,\mathrm{},v_N\},$$

where $`N=N(𝒫):=|(𝒫)|`$. Suppose the $`N+\left(\genfrac{}{}{0pt}{}{N}{2}\right)`$ points in $`(𝒫)+(𝒫)`$,

$$2v_1,\mathrm{},2v_N;v_1+v_2,v_1+v_3,\mathrm{},v_{N1}+v_N$$

are distinct. In this case, we say that $`𝒫`$ is a distinct pair-sum or dps polytope. Our interest in dps polytopes comes from the study of the representation of polynomials as a sum of squares of polynomials.

The following lemma offers two other geometrical characterizations of dps polytopes.

###### Lemma 1.

Let $`𝒫`$ be a lattice polytope. Then the following are equivalent:

$`\left(1\right)`$ $`\left(𝒫\right)`$ is a dps polytope.
$`\left(2\right)`$ $`\left(𝒫\right)`$ does not contain the vertices of a (nondegenerate) parallelogram, and does not contain three collinear points.
$`\left(3\right)`$ Suppose $`vv^{}`$ and $`ww^{}`$ are in $`\left(𝒫\right)`$. Then $`v^{}v`$ and $`w^{}w`$ are parallel only if $`\{v,v^{}\}=\{w,w^{}\}`$.

###### Proof.

$`\left(1\right)\left(2\right)`$. Suppose $`v_1,v_2,v_3,v_4\left(𝒫\right)`$ are the vertices of a parallelogram. Then $`v_1v_2=v_3v_4`$ implies $`v_1+v_4=v_2+v_3`$, so that $`𝒫`$ is not dps. Now suppose $`v_1,v_2,v_3\left(𝒫\right)`$, and $`v_2`$ is interior to the line segment $`\overline{v_1v_3}`$. If $`v_2`$ is the midpoint of the segment, then $`v_2+v_2=v_1+v_3`$, so $`𝒫`$ is not dps. Otherwise, we may assume that $`v_2`$ is closer to $`v_1`$ than to $`v_3`$. Then $`v_4=v_2+\left(v_2v_1\right)`$ will also be a lattice point on the line segment $`\overline{v_1v_3}`$, and $`v_2`$ is the midpoint of $`\overline{v_1v_4}`$; again, $`𝒫`$ is not dps.

$`\left(2\right)\left(3\right)`$. For $`u^n`$, let $`g\left(u\right)=\text{gcd}(u_1,\mathrm{},u_n)`$. Suppose $`g\left(u^{}u\right)=d>1`$. Then $`u^{}u=du^{\prime \prime }`$ for $`u^{\prime \prime }^n`$, and the line segment $`\overline{uu^{}}`$ contains the lattice points $`u,u+u^{\prime \prime },\mathrm{},u+du^{\prime \prime }=u^{}`$. Thus, if (2) holds and $`u,u^{}\left(𝒫\right)`$, $`uu^{}`$, we have $`g\left(u^{}u\right)=1`$. Suppose $`w^{}w=\alpha \left(v^{}v\right)`$. Then $`\alpha =p/q`$ for nonzero integers $`p,q`$, and $`q\left(w^{}w\right)=p\left(v^{}v\right)`$. Hence $`\left|q\right|=g\left(q\left(w^{}w\right)\right)=g\left(p\left(v^{}v\right)\right)=\left|p\right|`$, so $`\alpha =\pm 1`$. Now the parallelogram condition in (2) implies that $`\{v,v^{}\}=\{w,w^{}\}`$.

$`\left(3\right)\left(1\right)`$. If (3) holds for $`𝒫`$, and $`v_i,v_j,v_k,v_{\mathrm{}}\left(𝒫\right)`$ with $`i\{k,\mathrm{}\}`$, then $`v_iv_kv_{\mathrm{}}v_j`$, and so $`v_i+v_jv_k+v_{\mathrm{}}`$. This proves (1). ∎

Our main results are these: if $`𝒫`$ in $`^n`$ is a dps polytope, then $`N(𝒫)2^n`$, and, for every $`n`$, we construct dps polytopes in $`^n`$ for which $`N(𝒫)=2^n`$.

Example 1. Let $`𝒫^2`$ be the triangle with vertices $`\{(0,1),(1,2),(2,0)\}`$. Then $`𝒫`$ is a dps polytope, because

$$\left(𝒫\right)=\{(0,1),(1,2),(2,0),(1,1)\},$$

and

$$\left(𝒫\right)+\left(𝒫\right)=\{(0,2),(1,2),(1,3),(2,1),(2,2),(2,3),(2,4),(3,1),(3,2),(4,0)\}.$$

We can view $`𝒫`$ as the projection onto the first two coordinates of the triangle with vertices $`\{(0,1,2),(1,2,0),(2,0,1)\}`$, which lies in the hyperplane $`x_1+x_2+x_3=3`$. (In this example, we could have just as well taken the triangle with vertices $`\{(0,0),(1,2),(2,1)\}`$; again, $`\left(𝒫\right)`$ will consist of the vertices of $`𝒫`$ and $`(1,1)`$.)

Example 2. Let

$$\begin{array}{c}\text{ }𝒜=\{(4,1,0,0),(0,4,1,0),(0,0,4,1),(1,0,0,4)\},\text{ }\hfill \\ \text{ }=\{(2,1,1,1),(1,2,1,1),(1,1,2,1),(1,1,1,2)\};\text{ }\hfill \end{array}$$

and let $`𝒫=cvx\left(𝒜\right)^4`$ be the convex hull of $`𝒜`$. By construction, $`𝒫`$ is cyclically symmetric with respect to its coordinates. It is not hard to show that $`\left(𝒫\right)=𝒜`$. Suppose $`w=(w^{\left(1\right)},w^{\left(2\right)},w^{\left(3\right)},w^{\left(4\right)})\left(𝒫\right)`$. Since $`w`$ is a convex combination of $`𝒜`$, we have $`w^{\left(i\right)}0`$ and $`_iw^{\left(i\right)}=5`$. If $`w^{\left(i\right)}1`$ for all $`i`$, then $`w`$ must be a permutation of $`(2,1,1,1)`$ and so lies in $``$. Otherwise, $`w^{\left(i\right)}=0`$ for some $`i`$, and by cycling the coordinates, we may assume that $`w^{\left(1\right)}=0`$. But then $`w`$ must be a convex combination of $`(0,4,1,0)`$ and $`(0,0,4,1)`$ and so $`w𝒜`$. A routine check, which we omit, shows that the $`8+\left(\genfrac{}{}{0pt}{}{8}{2}\right)=36`$ sums in $`\left(𝒫\right)+\left(𝒫\right)`$ are distinct. By projecting $`𝒫`$ onto its first three coordinates, we obtain a dps polytope in $`^3`$ with $`N\left(𝒫\right)=8`$.

###### Theorem 2.

Suppose $`𝒫`$ is a dps polytope in $`^n`$. Then $`N(𝒫)2^n`$.

###### Proof.

If $`N\left(𝒫\right)>2^n`$, then by the Pigeonhole Principle, there exist $`v_iv_j`$ so that $`v_i`$ and $`v_j`$ are component-wise congruent modulo 2. This means that $`v_k=\frac{1}{2}\left(v_i+v_j\right)=v_i+\frac{1}{2}\left(v_jv_i\right)`$ is also a lattice point, and it follows from Lemma 1 that $`𝒫`$ is not a dps polytope. ∎

This argument is essentially the same one used to solve Putnam Problem 1971-A1 (see ): “Let there be given nine lattice points (points with integral coordinates) in three dimensional Euclidean space. Show that there is a lattice point on the interior of one of the line segments joining two of these points.” The proof of Theorem 2 also applies to the less restrictive class of convex polytopes which do not contain three lattice points on a line. One such polytope is the $`n`$-cube $`𝒞_n=\{0,1\}^n`$, which has many lattice parallelograms.

We shall say that a dps polytope $`𝒫^n`$ for which $`N(𝒫)=2^n`$ is maximal. The proof of Theorem 2 implies that no two points in a dps polytope are component-wise congruent modulo 2; hence a maximal dps polytope contains one representative from every congruence class modulo 2 (and at most one representative from every congruence class modulo $`m`$, $`m3`$).

Suppose $`M`$ is an $`n\times n`$ unimodular matrix with integer entries. Then $`M`$ defines a linear mapping on $`^n`$ (viewed as column vectors) by matrix multiplication. Since linear mappings preserve inclusions and both $`M`$ and $`M^1`$ have integer entries, it is easy to see that $`(M(𝒫))=M((𝒫))`$ for any lattice polytope $`𝒫`$, and since linear mappings preserve sums, it is then clear that $`𝒫`$ is dps if and only if $`M(𝒫)`$ is dps.

###### Theorem 3.

There exist maximal dps polytopes in $`^n`$ for every $`n`$.

###### Proof.

For $`n=1`$, let $`𝒫=[0,1]`$; for $`n=2,3`$, consider Examples 1 and 2. Suppose now that $`𝒫`$ is a maximal dps polytope in $`^n`$, $`n3`$. Write $`=\left(𝒫\right)`$ and define the (finite) set of differences

$$𝒟=\left(\right)^{}:=\{vv^{}:v,v^{},vv^{}\}.$$

Let $`M`$ be a unimodular integer matrix such that if $`u𝒟`$, then $`M\left(u\right)𝒟`$. (We shall construct such an $`M`$ below.)

We define the polytope $`𝒫^{}`$ in $`^{n+1}`$ as follows. Let

$$𝒜=\{(v,0)^{n+1}:v\left(𝒫\right)\},=\{(M\left(v\right),1)^{n+1}:v\left(𝒫\right)\},$$

and let $`𝒫^{}=cvx\left(𝒜\right)`$. If $`w=(w^{\left(1\right)},\mathrm{},w^{\left(n+1\right)})\left(𝒫^{}\right)`$, then $`0w^{\left(n+1\right)}1`$, hence $`w^{\left(n+1\right)}`$ equals 0 or 1. Thus, $`w`$ lies either on the face determined by $`𝒜`$, in which case $`w=(v,0)`$, or on the face determined by $``$, in which case $`w=(M\left(v\right),1)`$. It follows that $`\left(𝒫^{}\right)=𝒜`$, so $`N\left(𝒫^{}\right)=2^{n+1}`$.

Now consider $`\left(𝒫^{}\right)+\left(𝒫^{}\right)`$; this consists of three disjoint sets of points:

$$\left\{(v_i,0)+(v_j,0)\right\},\left\{(v_i,0)+(M\left(v_j\right),1)\right\},\left\{(M\left(v_i\right),1)+(M\left(v_j\right),1)\right\},$$

where $`v_i,v_j\left(𝒫\right)`$. Since both $`𝒫`$ and $`M\left(𝒫\right)`$ are dps, the sums in the first and the third set are distinct. For the second set, we suppose that

$$(v_i,0)+(M\left(v_j\right),1)=(v_k,0)+(M\left(v_{\mathrm{}}\right),1),$$
(1)
or equivalently,

$$v_iv_k=M\left(v_{\mathrm{}}\right)M\left(v_j\right)=M\left(v_{\mathrm{}}v_j\right).$$

If $`j=\mathrm{}`$, then $`v_iv_k=0`$, so $`i=k`$, which is the only possible way for (1) to hold in a dps polytope. Otherwise, $`j\mathrm{}`$, so $`M\left(v_{\mathrm{}}v_j\right)=v_iv_k𝒟`$, a contradiction to the choice of $`M`$. Thus, $`𝒫^{}`$ is a maximal dps polytope in $`^{n+1}`$.

We now construct a matrix $`M`$ with the desired properties. First, let

$$R=\mathrm{max}\left\{\left|u_j^{\left(k\right)}\right|:u_j𝒟,\mathrm{\hspace{0.33em}1}kn\right\}.$$

and let $`M`$ be the $`n\times n`$ matrix given below:

$$M=\left(\begin{array}{ccccccc}1+\left(R+1\right)^2& R+1& 0& 0& \mathrm{}& 0& 0\\ R+1& 1& R+1& 0& \mathrm{}& 0& 0\\ 0& 0& 1& R+1& \mathrm{}& 0& 0\\ 0& 0& 0& 1& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& 0& \mathrm{}& 1& R+1\\ 0& 0& 0& 0& \mathrm{}& 0& 1\end{array}\right).$$

(In words, the only non-zero entries in $`M`$ are the diagonal, the superdiagonal, and the first entry in the second row.) It is easy to see that $`M`$ is unimodular.

We show now that for every $`u𝒟`$, at least one entry of $`w=M\left(u\right)`$ has absolute value greater than $`R`$. This implies that $`M\left(u\right)𝒟`$, and will complete the proof. Write $`u=(u^{\left(1\right)},\mathrm{},u^{\left(n\right)})`$ and suppose that $`k`$ is the smallest index such that $`u^{\left(k\right)}0`$. (Such an index exists because $`0𝒟`$.)

If $`k=1`$, then $`w^{\left(1\right)}=\left(1+\left(R+1\right)^2\right)u^{\left(1\right)}+\left(R+1\right)u^{\left(2\right)}`$, and hence

$$\left|w^{\left(1\right)}\right|\left|\left(1+\left(R+1\right)^2\right)u^{\left(1\right)}\right|\left(R+1\right)\left|u^{\left(2\right)}\right|1+\left(R+1\right)^2R\left(R+1\right)=R+2.$$

If $`k2`$, then $`u^{\left(1\right)}=\mathrm{}=u^{\left(k1\right)}=0`$, so $`w^{\left(k1\right)}=\left(R+1\right)u^{\left(k\right)}`$ and $`\left|w^{\left(k1\right)}\right|R+1`$. Finally, we remark that the same proof applies in the case $`n=2`$, if we take as our matrix the $`2\times 2`$ submatrix at the upper left of $`M`$. ∎

Example 3. We illustrate the last construction by applying it to the polytope in Example 1, for which

$$𝒟=\{\pm (0,1),\pm (1,2),\pm (1,1),\pm (1,0),\pm (1,1),\pm (2,1)\},$$

so $`R=2`$ and

$$M=\left(\begin{array}{cc}10& 3\\ 3& 1\end{array}\right).$$

Thus, $`cvx\left(𝒜\right)`$ is a maximal dps polytope in $`^3`$, where

$$\begin{array}{c}\text{ }𝒜=\{(0,1,0),(1,1,0),(1,2,0),(2,0,0)\},\text{ }\hfill \\ \text{ }=\{(3,1,1),(13,4,1),(16,5,1),(20,6,1)\}.\text{ }\hfill \end{array}$$

We could now apply the shear $`(x_1,x_2,x_3)(x_13x_25x_3+5,x_2x_3,x_3)`$, which maps $`𝒜`$ and $``$ to

$$𝒜^{}:=\{(2,1,0),(3,1,0),(0,2,0),(7,0,0)\}$$

and

$$^{}:=\{(0,0,1),(1,3,1),(1,4,1),(2,5,1)\},$$

respectively, in order to reduce the magnitude of the coordinates in the example.

Since any translate of a dps polytope is also dps, we may always assume, as we have done in the examples, that $`𝒫`$ lies in the non-negative orthant of $`^n`$. In this case, we define $`s(𝒫)`$, the size of $`𝒫`$:

$$s\left(𝒫\right)=\mathrm{max}\left\{v_j^{\left(1\right)}+\mathrm{}+v_j^{\left(n\right)}:v_j\left(𝒫\right)\right\}.$$

If $`s=s(𝒫)`$, then $`𝒫`$ can be viewed as a projection onto the first $`n`$ coordinates of a polytope in $`^{n+1}`$ which lies in the simplex

$$\mathrm{\Delta }_{n+1}\left(s\right):=\{u=(u^{\left(1\right)},\mathrm{},u^{\left(n+1\right)}):u^{\left(i\right)}0,\underset{i=1}{\overset{n+1}{}}u^{\left(i\right)}=s\}.$$

Let $`s_n`$ denote the minimum size of any maximal dps polytope in $`^n`$. Examples 1 and 2 show that $`s_23`$ and $`s_35`$. It is not difficult to show that these estimates are sharp. The first case can be done by hand: if $`𝒫`$ is a maximal dps polytope with size 2 in $`^2`$, then $`(𝒫)`$ must consist of four points chosen from

$$\{(0,0),(0,1),(0,2),(1,0),(1,1),(2,0)\}.$$

Since each congruence class is represented in $`(𝒫)`$, it must contain $`(0,1)`$, $`(1,0)`$ and $`(1,1)`$. These three points form a parallelogram with each of the points $`(0,0),(0,2)`$ and $`(2,0)`$. Hence no fourth point can exist in $`(𝒫)`$ while preserving the dps property. The second case is similar, but much more complicated. Computer-aided calculations can be used to conclude that no dps polytope in $`^3`$ has size 4 or less. (We thank Dr. Bruce Carpenter for doing the Mathematica coding.)

It can also be shown, using the style of argument of \[6, Ch. 3\], that every maximal dps polytope in $`^2`$ is the image of the triangle in Example 1 under an affine unimodular linear mapping, and consists of a triangle with area 3/2, and a single lattice point inside, which will always be the centroid of the triangle. The tetrahedron determined by $``$ in Example 2 lies within the tetrahedron determined by $`𝒜`$, whereas in Example 3, each point in $``$ is on the boundary of the polytope. Thus there are at least two distinct combinatorial types of maximal dps polytopes in $`^3`$.

We make no serious conjecture about the growth of $`s_n`$. On the one hand, any maximal dps polytope must contain a lattice point with odd coordinates, so $`s_nn`$. In the other direction, it is not difficult to use the proof of Theorem 3 to obtain a doubly-exponential bound for $`s_n`$. Since this bound is likely to be very crude, we do not present it explicitly. Another open question is to determine the minimum volume of a maximal dps polytope in $`^n`$ for $`n3`$. We also do not know the answer to the following question: is every dps polytope a subset of a maximal dps polytope?

We now discuss our original interest in this subject. Given $`u_+^n`$, define the monomial $`x^u[x_1,\mathrm{},x_n]`$ by

$$x^u=x_1^{u^{\left(1\right)}}\mathrm{}x_n^{u^{\left(n\right)}}.$$

Suppose $`𝒰_+^n`$ and consider the polynomial

$$p(x_1,\mathrm{},x_n)=\underset{u𝒰}{}b_ux^u.$$

In , the present authors developed an algorithm for determining whether $`p`$ can be written as a sum of squares of polynomials. A necessary condition is that $`p`$ is psd; that is, $`p(x_1,\mathrm{},x_n)0`$ for all $`x^n`$. Suppose $`p`$ is psd and let

$$𝒞\left(p\right)=cvx\{u:b_u0\}.$$

Then $`𝒞(p)`$ is a lattice polytope; in fact it can be shown that the vertices of $`𝒞(p)`$ lie in $`(2)^n`$, so that $`𝒫:=\frac{1}{2}𝒞(p)`$ is a lattice polytope. Let

$$\left(𝒫\right)=\{v_1,\mathrm{},v_N\},$$

and for $`u𝒞(p)`$, let $`D(u)=\{(i,j):v_i+v_j=u\}`$. It is proved in \[4, Thm. 2.4\] that $`p`$ can be written as a sum of at most $`r`$ squares of polynomials if and only if there is a real $`N\times N`$ symmetric psd matrix $`A=[a_{ij}]`$ of rank at most $`r`$, so that

$$\underset{(i,j)D\left(u\right)}{}a_{ij}=b_u\text{for all }u𝒞\left(p\right)\text{.}$$

If $`𝒫`$ is a dps polytope in $`^n`$, then either $`|D(u)|1`$ or $`D(u)=\{(i,j),(j,i)\}`$. In either case, $`a_{ij}`$ is completely determined by $`b_u`$. In particular, if

$$h_𝒫(x_1,\mathrm{},x_n):=\underset{i=1}{\overset{N}{}}\left(x^{v_i}\right)^2,$$
(2)
then $`A`$ must equal $`I_N`$, the $`N\times N`$ identity matrix, so that $`p`$ is a sum of $`N`$ squares, and no fewer.

Finally, we note that the homogenization of polynomials with $`n`$ variables into forms with $`n+1`$ variables is precisely analogous to the embedding of polytopes in $`_+^n`$ into the hyperplane $`\mathrm{\Delta }_{n+1}(s)`$.

Example 4. (See \[4, Ex. 3.9\])

We return to Example 1, in its homogeneous version. Let $`A=[a_{ij}]`$ be a real symmetric $`4\times 4`$ matrix and let

$$f(t_1,t_2,t_3,t_4)=\underset{i=1}{\overset{4}{}}\underset{j=1}{\overset{4}{}}a_{ij}t_it_j$$

be its associated quadratic form. We use the substitution suggested by $`(𝒫)`$ and define the ternary sextic form

$$p(x_1,x_2,x_3)=f(x_2x_3^2,x_1x_2^2,x_1^2x_3,x_1x_2x_3).$$

Then $`p`$ is a sum of squares of polynomials (cubic forms) if and only if $`f`$ is a psd quadratic form; that is,

$$f(t_1,t_2,t_3,t_4)0\text{for all }(t_1,t_2,t_3,t_4)^4.$$

Since $`t_4^3=t_1t_2t_3`$, the condition for $`p`$ to be a psd form is weaker:

$$f(t_1,t_2,t_3,\left(t_1t_2t_3\right)^{1/3})0\text{for all }(t_1,t_2,t_3)^3.$$

If $`f(t_1,t_2,t_3,t_4)=t_1^2+t_2^2+t_3^23t_4^2`$, then $`f`$ is not psd, but $`f(t_1,t_2,t_3,(t_1t_2t_3)^{1/3})0`$ by the arithmetic-geometric inequality. It follows that

$$p(x_1,x_2,x_3)=x_2^2x_3^4+x_1^2x_2^4+x_1^4x_3^23x_1^2x_2^2x_3^2$$

is a form which is psd, but not a sum of squares of polynomials. This particular example was discussed in . For a history and bibliography of this subject and its relation to Hilbert’s 17th Problem, see .

More generally, the Pythagoras number of a ring $`A`$, $`P(A)`$, is the smallest number $`n\mathrm{}`$ such that any sum of squares in $`A`$ can be expressed as a sum of at most $`n`$ squares in $`A`$. Pfister  proved in 1967 that $`P((x_1,\mathrm{},x_n))2^n`$. It is easy to see that $`P([x_1])=2`$. Since maximal dps polytopes exist in $`^n`$ for every $`n`$, a consideration of $`h_𝒫`$ (c.f. (2)) shows that $`P([x_1,\mathrm{},x_n])2^n`$. This is not the strongest result possible: in \[2, p.60\], using other methods, Dai and the present authors have shown that $`P([x_1,\mathrm{},x_n])=\mathrm{}`$ for $`n2`$.


text:
# STATUS OF RING-LIKE CORRELATIONS AND WAVELETS

## Abstract

The problem of large-scale correlations of particles produced in high-energy collisions is discussed. Among them are, e.g., those correlations which lead to ring-like and elliptic flow shapes of individual high-multiplicity events in the polar+azimuthal angles plane. The wavelet method of analysis which allows to disentangle various patterns is proposed and applied to central Pb-Pb collisions at 158AGeV.

Various correlations of particles produced in high-energy collisions are known. Especially well studied are the two-particle correlations due to the decay of two-particle resonances and due to Bose-Einstein effect for identical particles. The common correlation function technique is well suited for this purpose. It is more difficult to apply it if several particles are involved. Nevertheless, correlations leading to clusters (or mini-jets, non-reducible to resonances) and to jets have been studied as well, however mostly in $`e^+e^{}`$-collisions where jet and subjet structures are clearly visible. Otherwise one has to apply the averaging procedure to get knowledge of dynamical correlations from studies of the moments of various distributions and their behavior in different regions of the phase space (for recent reviews, see ). It has lead to understanding of some global features of QCD .

However, the genuine multiparticle correlations originating from some collective effects may appear in individual events. There exists a method of the wavelet analysis which allows to recognize patterns due to correlations at different scales and damp down the statistical noise even if it is very large in an initial sample.

Before delving into applications of this method, let us consider examples of possible correlations which one would seek for. Everybody knows about Cherenkov radiation of photons by a bunch of electrons traversing a medium whose refractivity index exceeds 1. These photons form a ring in the plane perpendicular to electrons motion, i.e. they are emitted at a definite polar angle. As a hadronic analogue, one may treat an impinging nucleus as a bunch of confined quarks each of which can emit gluons when traversing a target nucleus. I speculated  about possible Cherenkov gluons relying on experimental observation of the positive real part of the elastic forward scattering amplitude of all hadronic processes at high energies. This is a necessary condition for such process because in the commonly used formula for the refractivity index its excess over 1 is proportional to this real part. However, later  I noticed that for such thin targets as nuclei the similar effect can appear due to small confinement length thus giving us a new tool for its estimate. If several gluons are emitted and each of them generates a mini-jet, the ring-like substructure will be observed in the target diagram. If the density of mini-jets within the ring is so high that they overlap, then they form a ridge pattern (or a wall ). If the number of emitted gluons is not large, we’ll see several jets (tower structure ) correlated in their polar, but not in the azimuthal angle. Central collisions of nuclei are preferred for observation of such effects because of a large number of participating partons. If the number of correlated particles within the ring is large enough, it would result in spikes in the pseudorapidity distributions. However, the usual histogram method is not always good to verify these spikes because it may split a single spike into two bins thus diminishing its role. However, some hints to such structure can be found from these histograms. The wavelet analysis is well suited for this purpose because it clearly resolves the local properties of a pattern on the event-by-event basis by choosing the so-called Heisenberg windows (see, e.g., ).

Another example is provided by the elliptic flow i.e. the azimuthal asymmetry in individual events. It may be related to a collective classical sling-effect  of the rotation of colliding nuclei after peripheral collisions initiated by the directed pressure  at some impact parameter or to some jetty structures.It can tell us about the equation of state of the hadronic matter by the shape of the created squeezed states. By passing, let me say that we also observed elliptic flow patterns (”cucumber”) corresponding to the large value of the second Fourier coefficient, and the ”three leaves flower” pattern with large third coefficient but I do not discuss them here.

The event-by-event analysis of patterns in experimental and Monte Carlo events becomes especially important with the advent of RHIC and LHC. The 4$`\pi `$ acceptance of detectors will be crucial for the event-by-event analysis .

When individual events are imaged visually, the human eye has a tendency to observe different kinds of intricate patterns with dense clusters (spikes) and rarefied voids. However, the observed effects are often dominated by statistical fluctuations. The method of factorial moments was proposed  to remove the statistical background in a global analysis and it shows fractal properties even in the event-by-event approach.

Both in cosmic rays  and at accelerators  some ring-like structures were observed in individual events. Especially astonishing was the famous NA22 event of pion-nucleon interaction with the pseudorapidity spike 60 times exceeding the average level.

Some event-by-event attempts (see, e.g., ) to treat nuclei collisions by methods of the traditional inclusive and correlation measures revealed that there are some jetty structures in individual events which lead to spikes in the angular (pseudorapidity) distribution. Various Monte Carlo simulations of the process were compared to the data and failed to describe this jettyness in its full strength. More detailed analysis  of large statistics data on hadron-hadron interactions (unfortunately, however, for rather low multiplicity) with dense groups of particles well separated from other particles in an event showed some ”anomaly” in the angular distribution of these groups. It was found that the centers of these groups prefer to be positioned at a definite polar angle. This feature favors the above interpretation in terms of Cherenkov gluons.

The wavelet analysis reveals the local properties of any pattern in an individual event at various scales and, moreover, avoids smooth polynomial trends and underlines the fluctuation patterns. By choosing the strongest fluctuations, one hopes to get those dynamical ones which exceed the statistical component.

First attempts to use wavelet analysis in multiparticle production go back to P. Carruthers  who used wavelets for diagonalisation of covariance matrices of some simplified cascade models. The proposals of correlation studies with the help of wavelets were promoted , and used, in particular, for looking for the disoriented chiral condensate . The wavelet transform of the pseudorapidity spectra of JACEE events was done in .

The wavelet analysis of patterns in central Pb-Pb collisions at 158AGeV was first done in . Five events with the highest registered multiplicities from 1034 to 1221 charged particles were chosen from 150 events of the emulsion chamber experiment EMU15 at CERN by the group from Lebedev Physical Institute.

The target diagrams of secondary particles distributions for these events are shown in Fig.1, where the radial distance from the center measures the polar angle $`\theta `$, and the azimuthal angle $`\varphi `$ is counted around the center. One can sum over the azimuthal angle and plot the corresponding pseudorapidity ($`\eta =\mathrm{log}\mathrm{tan}\theta /2`$) distributions shown in Fig.2. The pronounced peaks ($`\eta `$-spikes) strongly exceeding expected statistical fluctuations are seen in individual events. This inhomogeneity in pseudorapidity can arise either due to a very strong jet i.e. a large group (tower) of particles close both in polar and azimuthal angles or due to a ring-like (ridge) structure when several jets with smaller number of particles in each of them have similar polar angle but differ in their azimuthal angles.

In  wavelets were first used to analyze two-dimensional patterns of fluctuations of the event 19 in Fig.1. The results of the one-dimensional analysis  of the two-dimensional target diagram are presented. To proceed in this way, the whole azimuthal region was divided into 24 sectors. The pseudorapidity distributions in each of them were separately analyzed. Neighboring sectors were connected afterwards. Both jet and ring-like structures are found from the values of squared wavelet coefficients as seen from Fig.3 . At small scales $`a`$, the wavelet analysis reveals individual particles. At larger scales, the clusters or jets of particles are resolved. Finally, at ever larger scale one notices the ring-like structure around the center of the target diagram which penetrates from one azimuthal sector to another at nearby values of the polar angle (pseudorapidity), thus forming an elliptic ridge. This structure approximately corresponds to the peak in the pseudorapidity distribution (for more detail, see ). Let us note, that it is not easy to notice in the target diagram of this event any increase of density within the ring just by eye because of the specific properties of the $`\theta \varphi `$ plot where the density of particles decreases fast toward the external region of large polar angles.

To reveal these patterns in more detail one should perform the two-dimensional local analysis. It is strongly desirable to get rid of such drawback of the histogram method as fixed positions of bins that gives rise to splitting of a jet into pieces contained in two or more bins. This chance is provided by the wavelet transform of particle densities on the two-dimensional plot. Wavelets choose automatically the size and shapes of bins depending on particle densities at a given position.

In principle, the wavelet coefficients $`W_{j_1,k_1,j_2,k_2}`$ of the two-dimensional function $`f(\theta ,\varphi )`$ are found from the formula

$$W_{j_1,k_1,j_2,k_2}=f(\theta ,\varphi )\psi (2^{j_1}\theta k_1;2^{j_2}\varphi k_2)𝑑\theta 𝑑\varphi .$$
(1)
Here $`\theta _i,\varphi _i`$ are the polar and azimuthal angles of particles produced, $`f(\theta ,\varphi )=_i\delta (\theta \theta _i)\delta (\varphi \varphi _i)`$ with a sum over all particles $`i`$ in a given event, $`(k_1,k_2)`$ denote the locations and $`(j_1,j_2)`$ the scales analyzed. The function $`\psi `$ is the analyzing wavelet. The higher the density fluctuations of particles in a given region, the larger are the corresponding wavelet coefficients.

In practice, the discrete wavelets obtained from the tensor product of two multiresolution analyses of standard one-dimensional Daubechies 8-tap wavelets were used. Then the corresponding $`ss,sd`$ and $`dd`$ coefficients in the two-dimensional matrix were calculated (see ). The common scale $`j_1=j_2=j`$ was used.

As stressed above, the ring-like structure should be a collective effect involving many particles and large scales. Therefore, to get rid of the low-scale background due to individual particles and analyze their clusterization properties, we have chosen the scales $`j>5`$ where both single jets and those clustered in ring-like structures can be revealed as is seen from Fig.3. Therefore all coefficients with $`j<6`$ are put equal to 0. The wavelet coefficients for any $`j`$ from the interval $`6j10`$ are now presented as functions of polar and azimuthal angles in a form of the two-dimensional landscape-like surface over this plane i.e. over the target diagram. Their inverse wavelet transform allows to get modified target diagrams of analyzed events with large-scale structure left only. Higher fluctuations of particle density inside large-scale formations and, consequently, larger wavelet coefficients correspond to darker regions on this modified target diagram shown in Fig.4. Here we demonstrate two events (numbered 3 and 6) from those five shown in Figs.1 and 2. They display both jet and ring-like structures which are different in different events. To discard the methodical cut-off at $`\eta 1.61.8`$ the region of $`\eta >1.8`$ was only considered.

Even though the statistics is very low, it was attempted to plot the pseudorapidity distribution of the maxima of wavelet coefficients with the hope to see if it reveals the peculiarities observed in high statistics but low multiplicity hadron-hadron experiments . In Fig. 5, the number of highest maxima of wavelet coefficients exceeding the threshold value $`W_{j,k}>210^3`$ is plotted as a function of their pseudorapidity for all five events considered. It is quite peculiar that positions of the maxima are discrete. They are positioned quite symmetrically about the value $`\eta 2.9`$ corresponding to $`90^0`$ in cms as it should be for two Pb nuclei colliding. Difference of heights is within the error bars. More interesting, they do not fill in this central region but are rather separated. Qualitatively, it coincides with findings in .

For comparison, there were generated 100 central Pb-Pb interactions with energy 158AGeV according to Fritiof model and the same number of events according to the random model describing the inclusive rapidity distribution shape. The fluctuations in these events are much smaller than in experimental ones and do not show any ring-like structure.

Thus I conclude that even on the qualitative level there is the noticeable difference between experimental and simulated events with larger and somewhat ordered fluctuations in the former ones. My aim here is to show the applicability and power of the two-dimensional wavelet analysis , the qualitative features and differences leaving aside quantitative characteristics till higher statistics of high multiplicity AA-events becomes available. In particular, the special automatic complex for emulsion processing with high space resolution in Lebedev Physical Institute (www.lebedev.ru/structure/pavicom/index.htm) is coming into operation, and it will enlarge the statistics of central Pb-Pb collisions quite soon. STAR Collaboration of RHIC should also give soon some data at ever higher energy. Wavelets provide a powerful tool for event-by-event analysis of fluctuation patterns in such collisions.

The target diagrams of five events of central Pb-Pb collisions at energy 158 GeV/nucleon obtained by EMU-15 collaboration.

The pseudorapidity distributions of particles in five events shown in Fig.1.

Fig. 3 The wavelet coefficients for the event 19 analyzed in . Dark regions correspond to large values of the coefficients. Four of 24 sectors are shown. Rapidities are along x-axis, the scales increase down the vertical axis.

|  |  |
| --- | --- |

Fig.4 The modified large-scale target diagrams of two events (3 and 6). Darker regions correspond to larger particle density fluctuations.

Fig.5 The pseudorapidity distribution of the maxima of wavelet coefficients.
The irregularity in the maxima positions, the empty voids between
them and absence of peaks at $`\eta 2.9`$ are noticed.