text:
# Time resolved tracking of a sound scatterer in a complex flow: non-stationary signal analysis and applications

## I Introduction

In several areas of fluid dynamics research, it is desirable to study the motion individual fluid particles in a flow, i.e. the Lagrangian dynamics of the flow. The properties of this motion governs the physics of mixing, the behavior of binary flows and the Eulerian complexity of chaotic and turbulent flows. Lagrangian studies are possible in numerical experiments where chaotic Fountain00  and turbulent Yeung89 ; Squires91 ; Yeung94 ; Yeung97  flows have been studied. For turbulence, the numerical studies are limited to small Reynolds number flows whose evolution is only followed during a few large-eddy turnover times. In addition only the small scales properties of homogeneous turbulence are captured; the influence of inhomogeneities (such as large scale coherent structures) are not taken into account. It does not seem possible at the moment to extend high resolution turbulent DNS computations to long periods of time or to high Reynolds number flows. Experimental studies are thus needed. They differ from the numerical studies because one cannot tag and follow individual fluid particles; most techniques aim at recording the motion of solid particles carried by the flow motion. The degree of fidelity with which solid particles can act as Lagrangian tracers is an open problem; it depends on the size and density of the particle. While the interaction between the particle and its wake can be important for large particles or particles with a large density difference with the surrounding fluid Maxey83 ; Sridhar95 ; Mordant00 , it is generally admitted that density matched particles with a size smaller that the Kolmogorov length follow the fluid. Measurements of small particle motion have been made, using optical techniques that follow individual particle motion over short times/distances Virant97 ; Voth98 . We propose here an acoustic technique that can resolve an individual particle motion over long periods of time (compared with the characteristic time of flow forcing).

The principle of the technique is to monitor the Doppler shift of the sound scattered by a particle which is continuously insonified. This is an extension of the pulsed Doppler method that has been developed to measure velocity profiles and that has many applications in fluid mechanics and medicine Omer99 . The main advantage of the continuous insonification is to improve the time resolution of the measurement, although it is limited to the tracking of a very small number of particles (the tests reported here are made with only one particle in the flow). The measurement relies on the ability to track a Doppler frequency and its variation in time. For this signal processing problem two classes of approaches have been developed: (i) time-frequency analysis and (ii) high resolution parametric spectral analysis. Time frequency methods mainly rely upon the quadratic Wigner-Ville transform, or smoothed versions of it. Numerous studies and papers have recently been published, in which the theoretical issues are presented (see e.g. the textbooks by Flandrin Flandrin98  or Cohen Cohen95 ). These non parametric techniques are convenient and well-suited for weakly non-stationary signals with a good signal-to-noise ratio (SNR). However, time frequency representations present numerous drawbacks when it comes to extract trajectory information. Their quadratic nature give rise to numerous spurious interference terms that require post processing. For signals with a faster frequency modulation and a low SNR, we show here that an optimized parametric approach is a better choice. Parametric high resolution spectral analysis methods take advantage of an a priori knowledge of the spectral content of the recorded signal, namely the emitted signal frequency plus one or many doppler-shifted echoes. Furthermore, a time-recursive frame for the estimation of the Doppler shift is proposed here, where the evolution of the frequency is taken into account in the algorithm.

The two methods are tested in two experiments, in which the acoustic signals have different time scales and noise levels. The first experiment is a study of the transient acceleration of a heavy sphere settling under gravity in a fluid at rest. In this case the characteristic time scale of velocity variations is slow ($`\tau 50`$ ms) and the signal to noise ratio is fair (about 20 dB); we show that a technique of reassignment of the spectrogram gives good results. The second experiments deals with the motion of a neutrally buoyant sphere embedded in a turbulent flow. In this case, velocity variations occur over times of about 1 ms and the signal to noise ratio is low (less than 6 dB). We show that the AML parametric method yields very good results in that situation.

The paper is organized as follows: in section II we present the acoustic technique and measurement procedure. In section III we describe the signal processing techniques, with a particular emphasis on the AML method which has been developed and optimized to this particle tracking problem. Examples of applications to measurements in real flows are given in section IV.

## II Acoustical set-up

### II.1 Principle of the measurement

In the experimental technique proposed here, a particle is continuously insonified. It scatters a sound wave whose frequency is shifted from the incoming sound frequency due to the Doppler effect. This Doppler shift is directly related to the particle velocity $`𝐯_p`$:

$$\mathrm{\Delta }\omega =𝐪𝐯_p,$$
(1)
where $`𝐪`$ is the scattering wavevector (the difference between the incident and scattered wavevectors $`𝐪=𝐤_{\mathrm{scat}}𝐤_{\mathrm{inc}}`$) and $`\omega `$ is the wave pulsation.

We choose a backscattering geometry (see figure 1) so that $`𝐪=2𝐤_{\mathrm{inc}}`$ and the frequency shift becomes

$$\mathrm{\Delta }\omega (t)=2\frac{v(t)}{c}\omega _0,$$
(2)
where $`c`$ is the speed of sound, $`\omega _0`$ is the incident pulsation, and $`v(t)`$ is the component of the velocity on the incident direction at time $`t`$. We continuously insonify the moving particle and record the scattered sound. If need be, the particle position can be obtained by numerical integration of the velocity signal.

### II.2 Transducers characteristics and acquisition

We use a Vermon array of ultrasonic transducers made of individual elements of size 2$`\times `$2 mm each, separated by 100 $`\mu `$m. Their resonant frequency is about 3.2 MHz and their bandwidth at -3 dB is 1.5 MHz. Sound emission is set at 3 MHz or 3.5 MHz; experiments are performed in water so that the wavelength is $`\lambda `$=0.50 mm or 0.43 mm. The corresponding emission cone for each $`d=2`$ mm square element is $`29^{}`$ at 3 MHZ and $`24^{}`$ at 3.5 MHz. In our measurements, the particle to transducer distance lies between 5 cm and 40 cm, so that measurements are made in the far field ($`d^2/\lambda >10`$ mm). Given maximum flow and particle velocities of the order of 1.5 m.s-1, we expect a maximum sound frequency shift of the order of 5 kHz or 6 kHz, depending on whether the emission is at 3 MHz or 3.5 MHz. This yields a frequency modulation rate of at most 0.25%. One element of the transducer array is used for continuous sound emission and another for scattered sound detection. As the operation is continuous (as opposed to pulsed) and the elements are located close to one-another, we observe a coupling between the emitter and the receiver of the order of 60 dB (this is due both to electromagnetic and acoustic surface waves cross-talk).

The sound scattered by the moving particle is detected by a piezoelectric transducer. Upon connection to a 50 $`\mathrm{\Omega }`$ impedance, it yields an electrical signal of about 2 to 30 $`\mu `$V. In comparison, the noise is 1 $`\mu `$V and the electromagnetic coupling with the emitter is 8 mV. Hence the signal to noise ratio is between 0 dB and 30 dB.

The transducer output is sampled at 10 MHz over a 21 bit dynamical range (input range 31.25 mV) and numerically heterodyned at the emitting frequency. Then it is decimated at the final sampling frequency of 19531 Hz. The acquisition device is a HP-e1430A VXI digitizer.

### II.3 Scattering by an elastic sphere

The study of sound scattered by a fixed solid sphere is a classic but continuing area of study and difficulties arise in the interpretation of observed phenomena especially when trying to deal with elasticity and absorption Gaunard83 ; Thorne94 ; Hefner00 . Complex behaviour is observed linked with resonances of Rayleigh waves at the surface of the sphere. As a consequence the scattered pressure distribution varies both in directivity and amplitude. A generic expression for the far field pressure is the following :

$$p_{scat}(r,\theta )=p_{inc}\frac{af(ka,\theta )}{2r}e^{ikr},$$
(3)
where $`r`$ is the distance from the center of the sphere, $`a`$ its radius, $`p_{inc}`$ the incident pressure on the sphere, $`k`$ the incident wavenumber in the fluid, $`\theta `$ the scattering angle and $`f`$ is a form function which depends on the physical properties of the solid medium. Under very general assumptions, $`f`$ can be developed as a series of partial waves:

$$f(ka,\theta )=\frac{2}{ika}\underset{n=0}{\overset{\mathrm{}}{}}(2n+1)\frac{B_n(ka)}{D_n(ka)}P_n(\mathrm{cos}\theta ),$$
(4)
where $`P_n`$ is the a Legendre polynomial, $`B_n`$ and $`D_n`$ are determinants of matrices composed of spherical Bessel and Hankel function and their derivatives Gaunard83 . Physically, $`f`$ represents the sum of the specular echo and of interferences due to the radiation by Rayleigh waves Thorne94 ; Hefner00 . As a result, $`f`$ is a strongly varying function, particularly for high values of $`ka`$. In our experiments we used spheres of different material (polypropylene PP, steel, tungsten carbide, glass) with corresponding $`ka`$ between 7 and 15. The flow acts on the sphere motion, thus causing its acceleration and, eventually, its rotation. These effects may change the radiation diagram: first there is Doppler shift for the sound received by the sphere, and, perhaps more importantly, the sphere rotation may change the Rayleigh emission. For these reasons, the evolution of the amplitude of the scattered sound during the particle motion is quite complex. However, the observed amplitude modulation (see figures 2 and 8) varies slowly enough to allow a correct estimate of the frequency modulation of the scattered sound.

## III Signal processing

Numerous spectral estimation techniques are based on the ideas behind Fourier analysis of linear time invariant (LTI) differential equations. These techniques may be divided into (i) non-parametric techniques where the basis functions are implicitly the harmonically related complex exponentials of Fourier analysis and (ii) parametric techniques whose task is the estimation of the parameters of a (sub)set of complex exponentials. The spectrogram and the reassigned spectrogram belong the former category, whereas the maximum likelihood and its approximate form belong to the latter.

### III.1 Time-Frequency analysis

The most common time frequency distribution (TFD), the spectrogram, involves a moving time window. This window attempts to capture a portion of the signal which is sufficiently restricted in time so that stationarity and LTI assumptions are approximately met. To overcome the inherently poor localization in the time-frequency plane, a method has been proposed by Gendrin et al. Gendrin79 , and extended more recently by Auger and Flandrin Auger95 ; Flandrin98 . The idea is to locally reassign the energy distribution to the local center of gravity of the Fourier transform. Despite its ability to exhibit clear and well localized trajectories in the time-frequency plane, this technique requires an additional image processing step to extract the TF trajectory. For rapidly fluctuating frequency modulations and/or low SNR spurious clusters appear which makes this extraction difficult. The parametric method presented below is more robust.

### III.2 AML spectral estimation

This approach is largely based upon maximum likelihood spectral estimation (see e.g. KayKay88 ). The fundamentals are briefly recalled, as they serve as a basis for the approximate likelihood scheme, originally developed by Clergeot and Tressens Clergeot90 . This work is extended here within a recursive estimation frame, thus allowing to track the variations of the Doppler frequency shift induced by fast velocity changes of a scattering sphere imbedded in a turbulent flow. Michel and Clergeot have developed a similar approach for non stationary spectral analysis in an array processing frameMichel91a ; Michel91b .

#### III.2.1 Introduction

In this section, we address the problem of estimating the frequencies $`f_1,\mathrm{},f_M`$ of $`M`$ harmonic signals embedded in noise, from a small number of samples

$$x(t)=\underset{m=1}{\overset{M}{}}a_m(t)\mathrm{exp}(j(2\pi f_mt+\varphi _m)+n(t).$$
(5)
As the number of sampling points that are supposed to be available is low, classical Fourier based approaches fail to provide good results. We focus our attention on parametric approach, where an a priori knowledge about the structure of the signal is taken into account to improve the analysis.

The following assumptions are made:

* The time series is regularly sampled with time period $`T_s`$, so as to insure $`\frac{1}{T_s}>\frac{f_{max}}{2}`$ where $`f_{max}`$ stands for the bandwidth of the anti-aliasing filter used in the recording process. For convenience, $`T_s`$ will be set to $`T_s=1`$ and the term frequency will refer to normalized frequency (i.e. the actual frequency, divided by $`F_s=\frac{1}{T_s}`$.
* The amplitudes $`a_i(t)`$ are deterministic but unknown.
* The noise is a complex white gaussian circular and iid (independent increment identically distributed) with (unknown) variance $`\sigma ^2`$; the distribution function of a $`k`$-dimensional vector $`\stackrel{}{N}`$ defined by

$$\stackrel{}{N}(t)=[n(t),n(t+1),n(t+2),\mathrm{},n(t+(K1))]^T$$
(6)
reads

$$p(\stackrel{}{N})=\frac{1}{(\sqrt{2\pi }\sigma )^K}\mathrm{exp}\left(\frac{|\stackrel{}{N}|^2}{2\sigma ^2}\right).$$
(7)
Furthermore, the noise and the signal are independent.
* The term observation refers to a set of $`Q`$ $`K`$-dimensional vectors constructed from the sampled time series, according to

$$\stackrel{}{X}(t_j)=[x(t_j),x(t_j+1),x(t_j+2),\mathrm{},x(t_j+(K1))]^T\text{ }j=1,\mathrm{},Q$$
(8)
* The frequency $`f`$ is supposed to remain constant during an observation.

Under the assumption that the noise process is iid, and that the observed vectors are corrupted by independent realisations of the noise process, the likelihood of the observation is given by the product of the likelihood of each vector. Let $`𝒫`$ be the set of searched parameters ($`𝒫`$ contains $`\sigma ^2`$, $`F=\{f_1,\mathrm{},f_M\}`$ and the $`\stackrel{}{A}=[a_1\mathrm{exp}(j\varphi _1),\mathrm{},a_M\mathrm{exp}(j\varphi _M]^T`$), the loglikelihood of an observation is simply given by

$$(𝒫)=KQ\mathrm{log}(2\pi \sigma ^2)\frac{1}{2\sigma ^2}\underset{q=1}{\overset{Q}{}}|\stackrel{}{X}(q)\text{S}(\stackrel{}{F})\stackrel{}{A}(q)|^2,$$
(9)
where

$$\text{S}(\stackrel{}{F})=[\stackrel{}{S}_1,\mathrm{},\stackrel{}{S}_M]=\left(\begin{array}{cccc}1& \mathrm{exp}(2\pi f_1)& \mathrm{}& \mathrm{exp}(2\pi (K1)f_1\\ \mathrm{}& & & \mathrm{}\\ 1& \mathrm{exp}(2\pi f_M& \mathrm{}& \mathrm{exp}(2\pi (K1)f_M\end{array}\right)^T.$$
(10)
According to the maximum likelihood principle, the set $`𝒫`$ of parameters must be chosen in order to maximize expression (9).

#### III.2.2 Reduced expression

Minimizing (9) jointly for all the parameters is usually untractable. Most authors propose a separate maximisation for each of the parameters. For our application, the spectral components (i.e. $`\stackrel{}{A}`$ and $`\stackrel{}{F}`$) are the relevant variables. We first maximize with respect to $`\stackrel{}{A}`$ and derive an expression for the optimal $`\stackrel{}{F}`$; $`\sigma ^2`$ is estimated independently.

The value of vector $`\stackrel{}{A}`$ which minimizes the norm $`|\stackrel{}{X}(q)\text{S}(\stackrel{}{F})\stackrel{}{A}(q)|^2`$ is easily obtained :

$$\stackrel{}{A}(q)=(\text{S}^+\text{S})^1\text{S}^+(\stackrel{}{F})\stackrel{}{X}(q).$$
(11)
Note that the “signal only” vector $`\stackrel{}{Y}=\stackrel{}{X}\stackrel{}{N}`$ appears to be the orthogonal projection of $`\stackrel{}{X}`$ on the signal subspace spanned by the row vectors of S:

$$\stackrel{}{Y}=\text{S}(\stackrel{}{F}).\stackrel{}{A}(q)=\text{S}((\text{S}^+\text{S})^1\text{S}^+(\stackrel{}{F})\stackrel{}{X}=\mathrm{\Pi }_s(\stackrel{}{F})\stackrel{}{X},$$
(12)
where $`\mathrm{\Pi }_s(\stackrel{}{F})`$ stands for the parametric projector on the signal subspace111It is straightforward to establish that $`\mathrm{\Pi }_s\stackrel{}{Y}=\stackrel{}{Y}`$ for any vector $`\stackrel{}{Y}`$ lying in the signal subspace. Here ‘parametric projector’ must be understood as the projector calculated for the vector of frequencies $`\stackrel{}{F}`$ . Let $`\mathrm{\Pi }_n(\stackrel{}{F})=\text{I}\mathrm{\Pi }_s(\stackrel{}{F})`$ be the noise subspace, I is the identity matrix. By substituting (12) and using the definition of $`\mathrm{\Pi }_n(\stackrel{}{F})`$ in the expression of the log-likelihood (9), one gets the following simplified expression to minimise

$$L(\stackrel{}{F})=\frac{1}{\sigma ^2}\underset{q=1}{\overset{Q}{}}|\mathrm{\Pi }_n(\stackrel{}{F})\stackrel{}{X}|^2.$$
(13)
Using the properties of the trace operator (hereafter denoted $`\text{T}r`$) and those of the projection matrix $`\mathrm{\Pi }_n(\stackrel{}{F})`$, the maximum likelihood estimation of $`\stackrel{}{F}`$ takes the common form : minimize

$$L(\stackrel{}{F})=\frac{Q}{\sigma ^2}\text{T}r\left[\mathrm{\Pi }_n(\stackrel{}{F})\widehat{\text{R}}_x\right],$$
(14)
where $`\widehat{\text{R}}_x`$ is an estimate of the correlation matrix $`\text{R}_x`$ of the vector process $`\stackrel{}{X}(q)`$. Minimizing $`L(\stackrel{}{F})`$ in (14) leads to the exact value $`\stackrel{}{F}_{ML}`$ which has the maximum likelihood. It is important here to emphasize the following : if the vectors $`\stackrel{}{X}`$ are obtained by time-shift over the recorded time series, an observation runs over $`K+Q1`$ samples, i.e. the actual duration of one observation is $`T_{obs}=(Q+K2)T_s`$. In this case, the observed vectors may not be considered as being corrupted by independent realisations of the noise process, as some ’time integration’ is performed in the estimation of $`\text{R}_x`$. The consequences and interest of such smoothing have been studied by Clergeot and TressensClergeot90 , and OuamriOuamri86 , in the frame of array processing (in this context, ’time integration’ becomes ’spatial smoothing’). In the remainder of this paper, the development are based on equation (14), no matter how $`\text{R}_x`$ is estimated; see appendix for the practical implementation.

#### III.2.3 Approximate Max-likelihood

Equation (14) is still too complicated to be solved analytically in a simple way. A minimization can be easily performed if $`L(\stackrel{}{F})`$ has a quadratic dependance in SClergeot90 . Let $`\text{R}_y`$ be the correlation matrix of the signal vectors $`\stackrel{}{Y}(q)`$, the assumption that signal and noise are independent allow to establish the following equalities

$`\widehat{\text{R}}_x=\text{R}_y+\widehat{\sigma }^2\text{I},`$ (15)
$`\text{R}_y=\text{S}\text{P}\text{S}^+,`$ (16)
$`\text{P}=[\stackrel{}{A}\stackrel{}{A}^+],`$ (17)
where $``$ stands for the mathematical expectation. Substituting in equation (14) leads to:

$$L(\stackrel{}{F})=\frac{Q}{\widehat{\sigma }^2}\text{T}r\left[\mathrm{\Pi }_n(\stackrel{}{F})\text{S}\widehat{\text{P}}\text{S}^+\right].$$
(18)
Clergeot and TressensClergeot90  propose a second order approximation of $`L(\stackrel{}{F})`$:

$$L_{AML}(\stackrel{}{F})=\frac{Q}{\widehat{\sigma }^2}\text{T}r\left[\widehat{\mathrm{\Pi }}_n\text{S}(\stackrel{}{F})\widehat{\text{P}}\text{S}^+(\stackrel{}{F})\right],$$
(19)
in which $`\widehat{\mathrm{\Pi }}_n`$ is estimated by computing the projector spanned by the $`(KN)`$ smallest eigenvalues of the estimated covariance matrix $`\widehat{\text{R}}_x`$. They prove that this approach leads to more reliable estimates of $`\stackrel{}{F}`$ at low signal to noise ratio (SNR), and that the minimization of $`L_{AML}`$ is asymptotically efficient. In practice, the following set of equations is used

$`\widehat{\sigma }^2={\displaystyle \frac{1}{KM}}\text{T}r(\widehat{\mathrm{\Pi }}_n\widehat{\text{R}}_x),`$ (20)
$`\mathrm{\Pi }_s(\stackrel{}{F})=\text{S}(\stackrel{}{F})(\text{S}^+(\stackrel{}{F}).\text{S}(\stackrel{}{F}))^1.\text{S}^+(\stackrel{}{F}),`$ (21)
$`\text{S}(\stackrel{}{F}).\widehat{\text{P}}.\text{S}^+(\stackrel{}{F})=\mathrm{\Pi }_s(\stackrel{}{F})(\widehat{\text{R}}_x\sigma ^2\text{I}).\mathrm{\Pi }_s(\stackrel{}{F}).`$ (22)
The approximately quadratic dependence of $`L_{AML}`$ in $`\text{S}(\stackrel{}{F})`$, allows a fast convergence of the minimization algorithm by using a simple Newton-Gauss algorithm:

$$\stackrel{}{F}(k+1)=\stackrel{}{F}(k)\text{H}^1.\stackrel{}{\mathrm{grad}}(L_{AML})|_{\stackrel{}{F}=\stackrel{}{F}(k)},$$
(23)
where $`k`$ stands for the iteration step in the minimization process, $`\stackrel{}{\mathrm{grad}}`$ and H are the gradient and hessian respectively (see expressions in the appendix).

#### III.2.4 Combining new measurements and estimates

In this section, it is assumed that new measurements do not allow by itself the derivation of a good estimate. The variance of such an estimate varies as $`\frac{1}{T_{obs}}(K+Q1)^1`$, whereas integrating new measurements to this estimate allow to derive a better estimation.
Let $`\widehat{\stackrel{}{F}}(t)`$ be an estimate of $`\stackrel{}{F}`$ at time $`t`$, and $`𝒩(\widehat{\stackrel{}{F}}(t),\mathrm{\Gamma }(t))`$ its density, assumed to be normal with variance $`\mathrm{\Gamma }(t)`$222This is generally not the case, but this assumption remains valid as long as the loglikelihood is well approximated by its second order expansion around $`\widehat{\stackrel{}{F}}`$.. If a linear evolution model is known for $`\stackrel{}{F}(t)`$, one has

$`\stackrel{}{F}(t+1)=\text{M}\stackrel{}{F}(t)+\epsilon (t),`$ (24)
$`p_{t+1|t}(\stackrel{}{F})=𝒩(\text{M}\widehat{\stackrel{}{F}}(t),\text{M}\mathrm{\Gamma }(t)\text{M}^++\text{R}_\epsilon ),`$ (25)
where M is the evolution matrix; $`\epsilon `$ is a perturbation term, which is statistically independent from $`\stackrel{}{F}`$, and , $`\text{R}_\epsilon `$ is its covariance matrix. $`p_{t+1|t}`$ is the probability density function that can be derived for time $`t+1`$, if the observations are made until time $`t`$ only. As such an evolution equation is usually unknown, M will be set to the identity matrix in the rest of the paper (see MichelMichel91b ) for a detailed discussion). Applying the Bayes rule over conditional probabilities gives :

$$p_{t+1|t+1}(\stackrel{}{F})=\frac{p_{t+1|t}(\stackrel{}{F}).p_{t+1}(\stackrel{}{X}|\stackrel{}{F})}{p_{t+1}(\stackrel{}{X})}.$$
(26)
Noting that $`\mathrm{log}(p_{t+1}(\stackrel{}{X}|\stackrel{}{F}))`$ is the loglikelihood function for which a reduced expression has been derived in the previous section, one gets after all reductions and identifications the simple following expressions

$`\widehat{\stackrel{}{F}}(t+1|t)=\widehat{\stackrel{}{F}}(t),`$ (27)
$`\mathrm{\Gamma }(t+1|t)=\mathrm{\Gamma }(t)+\text{R}_\epsilon ,`$ (28)
$`\mathrm{\Gamma }(t+1)^1=\text{H}+\mathrm{\Gamma }(t+1|t)^1,`$ (29)
$`\widehat{\stackrel{}{F}}(t+1|t+1)=\widehat{\stackrel{}{F}}(t+1)=\widehat{\stackrel{}{F}}(t+1|t)\mathrm{\Gamma }(t+1)^1.\stackrel{}{\mathrm{grad}},`$ (30)
where it can be shown that the gradient function has the same expression as in the previous section. $`\text{R}_\epsilon `$ is an unknown matrix which will be practically set to $`v^2\text{I}`$, where $`v^2`$ will be tuned in order to allow the algorithm to take slight changes in $`\stackrel{}{F}`$ into account. Furthermore, it is interesting that the set of expression above expresses a generalized Kalman filter for estimating $`\stackrel{}{F}`$ (in the sense that it relies upon second order expansion of the loglikelihood functions). The statistical convergence properties and numerical efficiency of these approaches are described in the work of Michel & ClergeotMichel91a  and MichelMichel91b .

## IV Experimental results

We first describe the simple case of a particle settling in a fluid at rest. It is well adapted to the reassigned spectrogram method because the acoustic signal has a good SNR and a slow frequency modulation. We show that it allows to extract the subtle interaction between the falling particle and its wake. We then study the more complicated case of the motion of a particle embedded in a turbulent flow, where the dynamics of motion is much faster and the SNR is poor. We show that the AML method is well suited.

### IV.1 The settling sphere

#### IV.1.1 Motivation and experimental setup

When a particle is released in a fluid at rest, its developing motion creates a wake. The particle velocity is then set by the balance between buoyancy forces and drag, and additional subtle effects: first, ‘added mass’ corrections because the particles ‘pushes’ the fluid, and second, a ‘history’ force because the wake reacts back on the particle. Formally, one can write the equation of motion as Maxey83 ; Lawrence95 ; Mordant00 :

$$(m_p+\frac{1}{2}m_f)\frac{\mathrm{d}𝐯_p}{\mathrm{d}t}=(m_pm_f)𝐠\frac{1}{2}\pi a^2\rho _f𝐯_p𝐯_pc_D(Re)+𝐅_{\mathrm{history}}.$$
(31)
where $`m_p`$ is the particle mass, $`m_f`$ is the mass of a fluid particle of the same size, $`v_p`$ is the particle velocity, $`𝐠`$ is the acceleration of gravity, $`a`$ is the sphere radius, $`\rho _f`$ is the fluid density, $`c_D`$ is the static empiric drag coefficient, $`Re`$ is the Reynolds number $`Re=\frac{2av_p}{\nu }`$ and $`𝐅_{\mathrm{history}}`$ is the so-called history force. In this expression, the drag coefficient is usually obtained from measurement of the forces acting on a body at rest in an hydrodynamic tunnel. The history term, however, is largely unknown. Analytic expression can only be derived in the limit of small Reynolds numbers (less than 10) and cannot be applied for real flow configurations (e.g. multiphase flows) where $`Re1`$.

We perform measurements of the motion of a settling sphere, with the aim of evaluating the influence of the history forces. We use a water tank of size 1.1 m$`\times `$0.75 m and depth 0.65 m, filled with water at rest (figure 3). The bead is held by a pair of tweezers, five centimeters below the transducers. It is released a time t=0 without initial velocity and its trajectory is about 50 cm long. The data acquisition is started before the bead is released in order to capture the onset of motion.

#### IV.1.2 Results

Let us use as a first example, the fall of steel bead, 0.8 mm in diameter. The Doppler shift during the bead motion is detected using the spectrogram representation and a subsequent reassignment scheme. The simple spectrogram and reassigned version are shown in figure 4.

The reassignment technique drastically improves the localization of the energy in the time-frequency plane. In this case, the image processing step computes $`v_p(t)`$ as the line of maxima. The precision of the overall measurement depends on two factors : first on the intrinsic precision of the reassignment method and second on the dispersion of the measurements (the reproducibility of the bead motion over several experiments). The intrinsic precision of the reassignment method has been empirically studied using synthetic signals modelling the particle dynamics plus a noise that mimics the experimental data. We observed that for our choice of parameters (a time-frequency picture with 256$`\times `$256 pixels) the $`rms`$ precision is about one half pixel both in time and frequency directions. The method thus allows a precise analysis of the dynamics of the fall; we describe below two sets of experiment that illustrate the potential of the reassignment technique.

First, we show in figure 5 the velocity of a 1 mm steel bead (average over ten falls) together with two numerical simulations based on equation 31, first without the memory force and second with the expression of the memory force derived at low Reynolds numbers (called the Stokes memory term, as in Maxey & RileyMaxey83 ). The precision of the detection technique is sufficient for the measured profile to be compared to the simulated curves and to draw physical conclusions about the hydrodynamical forces. At early times, the trajectory is close to the simulation with memory force. This is due to the diffusion away from the bead surface of the vorticity generated at the boundary Maxey83 ; Lawrence95 ; Mordant00 . However, as the instantaneous Reynolds number increases, the curve deviates from this simple regime: vorticity is advected into the wake. Memory is progressively lost and the sphere reaches a terminal velocity in a finite time as does the simulation without memory.

The measurement and signal processing techniques are then tested on a more non-stationary motion, as in the case of a bead whose density is close to that of the fluid. In this situation a stronger interaction is expected between the particle motion and the development of its wake. Formally, this traces back to differences in the effective inertial mass and buoyancy mass of the particle – see equation (31). In figure 6, we show the velocity variation for a light glass sphere (density 2.48) compared to a tungsten bead (density 14.8). We observe that the velocity of the glass oscillates before reaching a constant terminal value whereas the other particle has a regular acceleration. In the case of light beads the hydrodynamic forces may be large enough to overcome the gravity and change the sign of the acceleration . This is linked with the non-stationarity of the wake, as vortex shedding is known to occur for Reynolds number above critical ($`Re_c250`$).

### IV.2 Turbulent flow : Lagrangian velocity measurement

#### IV.2.1 Experimental set-up

The turbulent flow is generated in a von Kármán geometry : the water is set into motion by two coaxial counter rotating disk in a cylindrical tank (figure 7). The Reynolds number $`Re=\frac{2\pi R^2f}{\nu }`$ (where $`\nu =0.8910^6`$ m2s-1 is the water kinematic viscosity) is equal to $`10^6`$. To prevent cavitation in the flow, we boil the water before filling the tank by lowering the pressure with a vacuum pump and during the experiment the pressure is increased to two bars. For the acoustic measurement, we use the same array of transducers as in the previous experiments, at emitting frequency 3 MHz. The cylinder and the surface of the disks are covered by 3 cm of Ciba Ureol 5073A and 6414B. Its density is 1.1 and the sound velocity is 1460 m.s-1 so that its acoustic impedance is close to that of the water, reducing drastically the reflections at the interface water/ureol compared to water/steel. The attenuation at 2.5 MHz is about 6 dB per cm. With a 3 cm layer and after the reflection on steel the total absorption is about 36 dB. The total reflection at the interfaces is reduced by a factor 60. The particle is a polypropylene (PP) sphere of radius 1 mm and density $`0.9`$.

#### IV.2.2 Results

We show in figure 8 the time series when one particle is in the ultrasonic beam and the corresponding spectrogram and reassigned spectrogram. The signal to noise ratio is very poor, typically less than 6dB (to give an idea, in figure 8(a) the bead enters the ultrasonic beam at $`t20`$ ms). One can also see some events localized in the time frequency plane that may be considered as noise and that may have several origin (noise of the motors, external electromagnetic noise, cavitation in the flow…). Altogether, the time-frequency pictures show the trajectory of the particle but the low SNR prevents it from being easily extracted. In particular, the trajectory in the reassigned picture becomes quite lacunar and extracting it would require sophisticated (and CPU greedy) image processing techniques.

The result of the AML algorithm is plotted in figure 9. The extracted frequency modulation is of course within the estimation in the spectrogram as in figure 8(b), but one observes that fine variations in the velocity of the bead are now detected. The algorithm provides also an estimate of the amplitude of the source (figure 9(c)). It can be seen that there is a strong amplitude modulation and that the SNR is at most 6 dB and may become less than 0 dB.

As the hessian is related to the Fisher information matrix Scharf , its inverse square root is linked with the variance of the estimation: a large value of the hessian indicates an accurate estimation of the modulation frequency and, hence, of the bead velocity. The inverse square root of the hessian is plotted in figure 9(b): very large values are calculated in the absence of a bead in the measurement volume at the beginning and end of the time series (as a signature of the mismatch between the model which is composed of one source at least and the reality: no source). Local lower values (typically less than 0.1) are observed when the variance on the estimation is small. Spurious effects are generated when the frequency modulation approaches zero as the hessian also becomes very small because of the filtering operation made in order to get rid of the coupling part of the signal. Finally, one observes that the hessian decreases as the signal to noise ratio increases (see at time 0.55 s).

## V Concluding remarks

As can be seen in the previous section, both methods, time-frequency analysis and parametric spectral analysis are suited for extracting the time-varying frequency modulation due to a Doppler effect. The domain of application of each method depends on the degree of non-stationarity and on the SNR.

For high SNR and weakly non-stationary signals, the time-frequency approach yields very good results. One drawback is the need of a second processing stage to extract the trajectory from the time-frequency picture. This stage may become increasingly difficult if there is more than one spectral component or if the SNR degrades. In both cases the quadratic nature of the algorithm produces interference patterns in the image: spurious clusters and a lacunar trajectory result. Another, more fundamental, limitation is that the length of the time window must be long enough to preserve an acceptable frequency resolution, even with the reassigned spectrogram. This limits the methods to weakly non-stationary signals.

For signals with a rapid frequency modulation, the AML spectral estimation is well suited, as long as the noise is near iid. The size of the time window can be decreased because of the parametric nature of the method, since a priori knowledge has been taken into account. The performance is further increased by the use of a Kalman-like filter. The drawback is the necessity to find a good dynamical model for the evolution of the spectral components. We have chosen here the simplest model which works well for our experiments but the approach can be refined by increasing the number of parameters in order to consider more precisely the variation of the frequency. The AML algorithm also provides a quantitative estimation of the quality of the demodulation and the instantaneous power of the spectral component. Finally, the AML method has the advantage to provide directly frequency modulation as a function of time, in one stage.

###### Acknowledgements.

We are indebted to Pascal Metz for the development of the signal conditioning electronics. We thank Marc Moulin for his help in the design of the von Kármán setup, VERMON for continuous assistance in the development of the transducer array. This work is partially supported by ACI grant No. 2226.

## AML algorithm

* First step : calculate $`\widehat{\text{R}}_x`$ using the following expressionMichel91b  :

$$\widehat{\text{R}}_x=\frac{1}{2Q}\underset{i=t+1}{\overset{t+Q}{}}\left(\stackrel{}{X}(i)\stackrel{}{X}(i)^T+\stackrel{~}{\stackrel{}{X}}(i)\stackrel{~}{\stackrel{}{X}}(i)^T\right),$$
(32)
with

$$\stackrel{~}{\stackrel{}{X}}(i)=[x(i+K1),x(i+K2),\mathrm{},x(i)]_{}^{}{}_{}{}^{T},$$
(33)
where  stands for complex conjugate. $`\stackrel{~}{\stackrel{}{X}}`$ is the complex conjugate of the time reversed version of $`\stackrel{}{X}`$.
* Second step : diagonalize $`\widehat{\text{R}}_x`$; one obtains the eigenvectors $`(\stackrel{}{V}_i)_{i=1..K}`$ and eigenvalues $`(\lambda _i)_{i=1..K}`$ sorted in decreasing order.
* Third step : Compute $`\widehat{\mathrm{\Pi }}_n`$ and $`\widehat{\sigma }^2`$, using the set of equations

$`\widehat{\mathrm{\Pi }}_n={\displaystyle \underset{i=M+1}{\overset{K}{}}}\stackrel{}{V}_i\stackrel{}{V}_i^T,`$ (34)
$`\widehat{\sigma }^2={\displaystyle \frac{1}{KM}}\text{T}r(\widehat{\mathrm{\Pi }}_n\widehat{\text{R}}_x)={\displaystyle \frac{1}{KM}}{\displaystyle \underset{i=M+1}{\overset{K}{}}}\lambda _i.`$ (35)
* Forth step : choose $`\stackrel{}{F}=\widehat{\stackrel{}{F}}(t)`$ as candidate value. Compute $`\stackrel{}{\mathrm{grad}}`$ and H usingMichel91a

$`\stackrel{}{\mathrm{grad}}={\displaystyle \frac{2Q}{\sigma ^2}}\mathrm{R}_\mathrm{e}\left\{\mathrm{Diag}(\text{S}^+(\stackrel{}{F}).\mathrm{\Pi }_n(\stackrel{}{F}).\widehat{\mathrm{\Pi }}_n.\text{S}(\stackrel{}{F}).\widehat{\text{P}})\right\},`$ (36)
$`\text{H}={\displaystyle \frac{2Q}{\sigma ^2}}\mathrm{R}_\mathrm{e}\{\mathrm{Diag}\left((\text{S}^+(\stackrel{}{F}).\mathrm{\Pi }_n(\stackrel{}{F}).\widehat{\mathrm{\Pi }}_n.\mathrm{\Pi }_n(\stackrel{}{F}).\text{S}^{})\right)\widehat{\text{P}}^{}\},`$ (37)
where the operator $``$ stand for the term to term matrix multiplication, and $`\text{P}^{}`$ is the conjugate of P, and

$$\text{S}^{}=[\frac{d\stackrel{}{S}_1}{df_1},\mathrm{},\frac{d\stackrel{}{S}_M}{df_M}]^T.$$
(38)
* Fifth step : using equations (27) to (30), compute $`\widehat{\stackrel{}{F}}(t+1)`$ and $`\mathrm{\Gamma }(t+1)`$.

The initialization of the algorithm is done by either (i) setting an initial value of $`\stackrel{}{F}(1)`$ or (ii) estimating this value using the maxima of the amplitude of the FFT of a small window of signal (of length 64 or 128 samples) and using the iterative algorithm described in section III.2.3 to converge towards $`\stackrel{}{F}(1)`$.

For example, the extracted velocity of figure 9 is obtained by starting at the maximum of energy of the signal and applying the algorithm forward and backward in time. The algorithm is stopped as the mean of the inverse square root of the hessian over a window of size 400 samples exceeds 0.5 for more than 400 samples.


format: text

text:
# 𝑆𝑈(1,1) Random Polynomials

## 1 Introduction

In this paper we are interested in statistical properties of zeros of random polynomials and random analytic functions associated with the pseudoeuclidean group of symmetries $`SU(1,1)`$. The motivation for the study of zeros of random polynomials and random analytic functions comes from different applications, most importantly from the theory of quantum chaos (see papers of Bogomolny, Bohigas, and Leboeuf \[BBL\], Leboeuf and Shukla \[LS\], Hannay \[Han\], Korsch, Miller, and Wiescher \[KMW\], Nonnenmacher and Voros \[NV\], Forrester and Honner \[FH\], Leboeuf \[L\], Shiffman and Zelditch \[SZ1\], and others). There are different ensembles of random polynomials associated with different groups of symmetries. In this respect we can make a parallel to the theory of random matrices. In the latter, the basic ensembles include orthogonal, unitary, symplectic, circular, and some others (see \[M\]). Making the parallel, we may think of the $`O(n+1)`$ ensemble of real polynomials as of an analogue of the orthogonal ensemble of random matrices. The $`O(n+1)`$ ensemble of random polynomials consists of multivariate homogeneous real polynomials of $`(n+1)`$ variable $`z=(z_0,z_1,\mathrm{},z_n)`$ of the form

$$\psi (z)=\underset{|m|=N}{}\sqrt{C_N^m}a_mz^m\text{,}$$
(1)
where $`m=(m_0,m_1,\mathrm{},m_n)`$ is a multiindex, $`|m|=m_0+m_1+\mathrm{}+m_n`$, $`z^m=z_0^{m_0}z_1^{m_1}\mathrm{}z_n^{m_n}`$, $`C_N^m`$ is the multinomial coefficient,

$$C_N^m\frac{N!}{m_0!m_1!\mathrm{}m_n!}\text{,}$$
(2)
and $`a_m`$ are independent standard, $`N(0,1)`$, real Gaussian random variables. Consider a set $`Z`$ of joint real zeros of $`k`$, $`kn`$, independent copies of the $`O(n+1)`$ random polynomial in the projective space $`^n`$. This is a random real algebraic variety of dimension $`nk`$. It is nondegenerate almost surely and it possesses a natural volume element induced by the standard metric in $`^n`$. The joint distribution functions of the zeros are invariant with respect to the action of the group $`O(n+1)`$, and, in particular, the density function of the zeros is constant (cf. \[ShSm\], \[BSZ2\], \[BSZ3\]).

Similarly, the $`SU(n+1)`$ ensemble of random polynomials consists of multivariate homogeneous complex polynomials of $`(n+1)`$ complex variable $`z=(z_0,z_1,\mathrm{},z_n)`$ of form (1) where $`a_m`$ are independent standard complex Gaussian random variables (cf. \[BBL\], \[Han\], \[BSZ1\]-\[BSZ4\]). It corresponds to the unitary ensemble of random matrices. The distribution of joint zeros of $`k`$ independent copies of the $`SU(n+1)`$ random polynomial is invariant with respect to the action of the group $`SU(n+1)`$ on the complex projective space $`^n`$. As the degree $`N`$ of the $`SU(n+1)`$ random polynomial goes to infinity, the scaled correlation functions of zeros approach a limit, which is represented by the correlation functions of the $`W_n`$ ensemble of random analytic functions (see \[Han\], \[BSZ1\]-\[BSZ4\], \[L\]). The latter consists of random functions of the form

$$\psi (z)=\underset{m}{}\sqrt{\frac{1}{m!}}a_mz^m\text{,}$$
(3)
where the sum runs over multiindices $`m=(m_1,\mathrm{},m_n)`$ with $`m_j0`$, $`m!=m_1!\mathrm{},m_n!`$, and $`a_m`$ are independent standard complex Gaussian random variables. The existence of the limit for the scaled correlation functions of zeros is valid in a very general framework of random sections of powers of line bundles over compact manifolds and the limit is universal (see \[BSZ1\]-\[BSZ4\]). The explicit combinatorial formulae for the limit of two-point correlation functions are obtained in \[BSZ2\]-\[BSZ4\]. They demonstrate a quadratic repulsion if $`k=n=1`$, neutrality if $`k=n=2`$, and attraction if $`k=n>2`$.

In this paper we are interested in the pseudoeuclidean, $`SU(1,1)`$ ensemble of random analytic functions. The general SU$`(n,1)`$ ensemble consists of multivariate analytic functions of the form

$$\psi (z)=\underset{m}{}\sqrt{C_{|m|+L1}^m}a_mz^m\text{,}z=(z_1,\mathrm{},z_n)\text{,}$$
(4)
where $`L1`$ is a fixed integer, a parameter of the ensemble, the sum runs over multiindices $`m=(m_1,\mathrm{},m_n)`$ with $`m_j0`$, $`z^m=z_1^{m_1}\mathrm{}z_n^{m_n}`$,

$$C_{|m|+L1}^m\frac{(|m|+L1)!}{(L1)!m_1!\mathrm{}m_n!}\text{,}$$
(5)
and $`a_m`$ are independent standard complex Gaussian random variables. We will restrict our study to the case $`n=1`$. The basic calculations are extended to the case of any $`n`$, and we are going to return to this extension in subsequent publications. Another possible extension concerns the universality result like in \[BSZ2\], \[SZ2\], and \[BSZ3\].

The plan of the paper is the following. In Section 2 we consider the ensemble of $`SU(1,1)`$ random polynomials, which is obtained by restricting $`m`$ in (4) to be bounded by $`N`$ (with $`n=1`$). We calculate the density function for the distribution of zeros of the $`SU(1,1)`$ random polynomial and we show that the zeros are concentrated in a narrow annulus around the unit circle, of the width of the order of $`1/N`$. We find the scaled profile of the density function along the radius in the limit $`N\mathrm{}`$. In Subsection 2.2 we rederive the result of Leboeuf \[L\] for the density of the $`SU(1,1)`$ random analytic function. In Section 3 we derive different formulae for the correlation functions between zeros of the $`SU(1,1)`$ random analytic function. First, we apply the general result of \[BSZ2\] to get the Kac-Rice type expression for the $`n`$-point correlation function. We then specify it for the case of a Gaussian random field (cf. \[Han\], \[BD\], \[BSZ2\]). This expresses the $`n`$-point correlation function in terms of some Gaussian averages, which is further represented as a sum over Feynman diagrams, or in a different approach developed in \[BSZ4\], as a supersymmetric integral. In Subsection 3.3 we carry out concrete calculations for the two point correlation function in the $`SU(1,1)`$ ensemble. They demonstrate a quadratic repulsion at small distances and fast decay at infinity. In the limit when the main parameter $`L`$ of the $`SU(1,1)`$ polynomial goes to infinity we recover the $`W_1`$ correlation function of Hannay. There is an appendix at the end of the paper, in which we consider correlations between outer zeros $`z_j`$, $`|z_j|>1`$, of the $`SU(1,1)`$ random polynomial. We show that in the limit $`N\mathrm{}`$, the outer zeros are independent of the inner zeros $`|z_j|<1`$, and after changing the variable $`z`$ to $`1/z`$, the correlations between outer zeros coincide with the correlations between zeros in the $`SU(1,1)`$ ensemble with $`L=1`$.

## 2 Basic Statistical Properties

### 2.1 Scaled Density Function

Consider the following random polynomial, associated with the pseudoeuclidean group $`SU(1,1)`$:

$$\psi (z)=\underset{m=0}{\overset{N}{}}\sqrt{C_{m+L1}^m}a_mz^m\text{,}C_{m+L1}^m\frac{(m+L1)!}{(L1)!m!}\text{,}$$
(6)
where $`L`$ is a fixed positive integer, a parameter of the problem, $`C_{m+L1}^m`$ are Newton’s binomial coefficients, and $`a_m`$ are independent standard complex Gaussian random variables, so that

$$Ea_m=0\text{,}Ea_m\overline{a_n}=\delta _{mn}\text{,}Ea_ma_n=0\text{.}$$
(7)
When $`L=1`$, (6) reduces to the classical form

$$\psi (z)=\underset{m=0}{\overset{N}{}}a_mz^m\text{.}$$
(8)
We wish to investigate the density of zeros $`p_N(z)`$ for the $`SU(1,1)`$ polynomial (6) as $`N\mathrm{}`$. The density $`p_N(z)`$ is determined by the condition that for any test function $`\phi (z)`$, which is infinitely differentiable and compactly supported,

$$E\left(\underset{j=1}{\overset{N}{}}\phi (z_j)\right)=_{}p_N(z)\phi (z)𝑑z\text{,}dzdxdy\text{,}$$
(9)
where $`z_j`$ are zeros of the random polynomial $`\psi (z)`$. Since the total number of zeros is equal to N, we have that

$$_{}p_N(z)𝑑z=N\text{.}$$
(10)
Observe that the polynomial

$$\psi (e^{i\theta }z)=\underset{m=0}{\overset{N}{}}\sqrt{C_{m+L1}^m}a_me^{im\theta }z^m$$
(11)
has the same probability distribution as $`\psi (z)`$, hence

$$p_N(re^{i\theta })=p_N(r)\text{.}$$
(12)
The general formula for $`p_N(z)`$ is given by the Poincaré-Lelong type expression

$$p_N(z)=\frac{1}{\pi }\frac{^2}{z\overline{z}}[\mathrm{ln}E(\psi (z)\overline{\psi (z)})]$$
(13)
(see \[BSZ1\], \[L\]). From (6),

$$E(\psi (z)\overline{\psi (z^{})})=\underset{m=0}{\overset{N}{}}C_{m+L1}^m(z\overline{z^{}})^m=\frac{1}{(L1)!}\underset{m=0}{\overset{N}{}}(m+1)\mathrm{}(m+L1)(z\overline{z^{}})^m\text{.}$$
(14)
In particular,

$$E(\psi (z)\overline{\psi (z)})=_{N,L}(x)\underset{m=0}{\overset{N}{}}C_{m+L1}^mx^m=\frac{1}{(L1)!}\underset{m=0}{\overset{N}{}}(m+1)\mathrm{}(m+L1)x^m\text{.}$$
(15)
where here and in what follows we use the notation

$$x=z\overline{z}=|z|^2\text{ .}$$
(16)
It is obvious that

$$_{N,L}(x)>0,x0\text{ .}$$
(17)
Substituting (15) into (13) gives that

$$p_N(z)=\frac{1}{\pi }\left[x\left(\frac{_{N,L}^{}(x)}{_{N,L}(x)}\right)^{}+\frac{_{N,L}^{}(x)}{_{N,L}(x)}\right]\text{ .}$$
(18)
We further define

$$𝒢_{N,L}(x)=\underset{m=0}{\overset{N}{}}x^{m+L1}=x^{L1}\frac{x^{N+1}1}{x1}\text{ .}$$
(19)
Then from (15),

$$_{N,L}(x)=\frac{1}{(L1)!}\frac{d^{L1}𝒢_{N,L}(x)}{dx^{L1}}\text{ .}$$
(20)
Numerical simulations show that most of the roots of $`SU(1,1)`$ polynomials (6) are concentrated in a small annulus (of the width of the order of $`1/N`$) near the unit circle. This property is illustrated in Fig. 1, which contains a point-plot of the zeros of 150 $`SU(1,1)`$ polynomials of degree $`N=200`$, with $`L=30`$.

To get the asymptotics of the density $`p_N(z)`$ near the unit circle, we introduce a scaling of the variable $`x`$ in the form

$$x=1+\frac{s}{N}\text{ .}$$
(21)
Then,

$$\frac{d}{dx}=N\frac{d}{ds}\text{ .}$$
(22)
In subsequent calculations, we will assume that

$$AsA\text{ ,}$$
(23)
for some arbitrary fixed $`A>0`$. The notation $`R(s)=O(N^1)`$ used below means that there exists some $`C(A)>0`$ so that $`|R(s)|C(A)N^1`$ for all $`s[A,A]`$. The main result of this section is the following theorem.

###### Theorem 2.1.

As $`N\mathrm{}`$,

$$N^2p_N\left(\left(1+\frac{s}{N}\right)^{1/2}e^{i\theta }\right)=\frac{1}{\pi }\left[\frac{g^{(L)}(s)}{g^{(L1)}(s)}\right]^{}+O(N^1)\text{ .}$$
(24)
where

$$g^{(j)}(s)=\frac{d^j}{ds^j}\left(\frac{e^s1}{s}\right)\text{ .}$$
(25)
Remark. Observe that one $`N`$ in the normalization of $`p_N(z)`$ on the left is due to the rescaling (21) and one is due to the integral condition (10).

Proof. By substituting (21) into (19) we obtain that as $`N\mathrm{}`$,

$$𝒢_{N,L}(x)=\left(1+\frac{s}{N}\right)^{L1}\frac{(1+\frac{s}{N})^{N+1}1}{\frac{s}{N}}=N\frac{e^s1}{s}\left(1+O(N^1)\right)=Ng(s)\left(1+O(N^1)\right)\text{ ,}$$
(26)
where we define

$$g(s)=\frac{e^s1}{s}=\underset{j=0}{\overset{\mathrm{}}{}}\frac{s^j}{(j+1)!}\text{ .}$$
(27)
Observe that for $`j=0,1,2,\mathrm{}`$,

$$g^{(j)}(s)\frac{d^jg(s)}{ds^j}>0\text{ ,}\mathrm{}<s<\mathrm{}\text{ .}$$
(28)
It is obvious from (27) for $`s0`$. For negative $`s`$, use the identity

$$g^{(j)}(s)=\frac{j!}{s^{j+1}}\left[1e^s\left(1+\frac{s}{1!}+\mathrm{}+\frac{s^j}{j!}\right)\right]\text{ .}$$
(29)
Formula (26) holds obviously in a small complex neighborhood of the segment $`AsA`$. Therefore, we can differentiate it in $`s`$, so that

$$\frac{d^j𝒢_{N,L}(x)}{dx^j}=N^{j+1}g^{(j)}(s)\left(1+O(N^1)\right)\text{ ,}$$
(30)
Thus, from (20), we obtain that

$$_{N,L}(x)=\frac{N^Lg^{(L1)}(s)}{(L1)!}\left(1+O(N^1)\right)\text{,}g^{(L1)}(s)\frac{d^{L1}g(s)}{ds^{L1}}\text{ ,}$$
(31)
and

$$_{N,L}^{}(x)=\frac{N^{L+1}g^{(L)}(s)}{(L1)!}\left(1+O(N^1)\right)\text{ ,}$$
(32)
hence

$$\frac{_{N,L}^{}(x)}{_{N,L}(x)}=\frac{Ng^{(L)}(s)}{g^{(L1)}(s)}\left(1+O(N^1)\right)\text{ .}$$
(33)
Let us go back to formula (18). In this formula, $`x=1+O(N^1)`$, and the first term in the brackets has an extra derivative so this is the leading term, and we can neglect the second term. This gives that

$$N^2p_N(z)=\frac{1}{\pi }\left[\frac{g^{(L)}(s)}{g^{(L1)}(s)}\right]^{}+O(N^1)\text{,}z\overline{z}=1+\frac{s}{N}\text{ ,}$$
(34)
which was stated. Theorem 2.1 is proved.

Equation (24) can also be expressed in the following form:

$$N^2p_N(z)=\frac{1}{\pi }\left[\mathrm{log}g^{(L1)}(s)\right]^{\prime \prime }+O(N^1)\text{.}$$
(35)
From (24),

$$\underset{N\mathrm{}}{lim}N^2p_N\left(\left(1+\frac{s}{2N}\right)e^{i\theta }\right)=p(s)\frac{1}{\pi }\left[\frac{g^{(L)}(s)}{g^{(L1)}(s)}\right]^{}=\frac{1}{\pi }\left[\mathrm{log}g^{(L1)}(s)\right]^{\prime \prime }\text{.}$$
(36)
Consider now the distribution function of zeros. Define

$$P_N(x)=N^1E\mathrm{\#}\{j:|z_j|^2x\}\text{,}$$
(37)
which gives the expected value of the fraction of zeros in the disk of radius $`\sqrt{x}`$. Then

$$P_N^{}(x)=N^1\pi p_N(z)\text{,}x=|z|^2\text{.}$$
(38)
Hence (24) implies that

$$P_N\left(1+\frac{s}{N}\right)=\frac{g^{(L)}(s)}{g^{(L1)}(s)}+O(N^1)\text{.}$$
(39)
As follows from (29), the limiting distribution function,

$$\underset{N\mathrm{}}{lim}P_N\left(1+\frac{s}{N}\right)=P(s)\frac{g^{(L)}(s)}{g^{(L1)}(s)}\text{,}$$
(40)
has the following asymptotics at $`\mathrm{}`$:

$$P(s)=\frac{L}{s}+O\left((s)^{L1}e^s\right)\text{,}s\mathrm{}\text{.}$$
(41)
At $`+\mathrm{}`$ we use the formula

$$g^{(j)}(s)=\frac{e^s}{s}\left[1\frac{j}{s}+\frac{j(j1)}{s^2}\mathrm{}+\frac{(1)^jj!}{s^j}\right]+\frac{(1)^{j+1}j!}{s^{j+1}}\text{,}$$
(42)
which gives that

$$P(s)=\frac{1}{s}\frac{L1}{s^2}+O\left(s^3\right)\text{,}s\mathrm{}\text{.}$$
(43)
For $`L=1`$, (40) reduces to the well-known result (see \[BS\])

$$P(s)=\frac{se^se^s+1}{s(e^s1)}\text{.}$$
(44)
For $`L=1`$, the limiting distribution is symmetric, so that

$$P(1s)=1P(s)$$
(45)
(it is related to the symmetry of the zeros of $`\psi (z)`$ in (8) with respect to the transformation $`z1/z`$). Since

$$g(s)=1+\frac{s}{2!}+\frac{s^2}{3!}+\frac{s^3}{4!}+\mathrm{}\text{,}$$
(46)
we have that

$$g^{(j)}(s)=\frac{1}{j+1}+\frac{s}{1!(j+2)}+\frac{s^2}{2!(j+3)}+\frac{s^3}{3!(j+4)}+\mathrm{}\text{.}$$
(47)
Thus,

$$P_N\left(1+\frac{s}{N}\right)=\frac{{\displaystyle \frac{1}{L+1}}+{\displaystyle \frac{s}{1!(L+2)}}+{\displaystyle \frac{s^2}{2!(L+3)}}+\mathrm{}}{{\displaystyle \frac{1}{L}}+{\displaystyle \frac{s}{1!(L+1)}}+{\displaystyle \frac{s^2}{2!(L+2)}}+\mathrm{}}+O(N^1)\text{.}$$
(48)
In particular,

$$P_N(1)=\frac{L}{L+1}+O(N^1)\text{,}$$
(49)
which means that the expected value of the fraction of zeros inside the unit disk is asymptotically equal to $`L/(L+1)`$.

In our numerical simulations we generated a large number of $`SU(1,1)`$ random polynomials of degree $`N`$, calculated the zeros using standard techniques, and counted the number of zeros in annuli of fixed width, which are concentrically spread around the origin, per one polynomial. The scaled density was obtained by dividing each of these numbers by the area of the corresponding annulus, as well as by adjusting it with respect to $`N`$ (according to the previously established result). Figure 2 shows the results of this procedure for $`L=4`$ and $`N=150`$, in comparison with the theoretical value as given by Eq. (24).

### 2.2 $`SU(1,1)`$ Ensemble of Random Analytic Functions

In this section we consider the $`SU(1,1)`$ random analytic function,

$$\psi (z)=\underset{m=0}{\overset{\mathrm{}}{}}\sqrt{C_{m+L1}^m}a_mz^m\text{,}$$
(50)
which is obtained from (6) by setting $`N=\mathrm{}`$. Here again $`C_{m+L1}^m`$ are Newton’s binomial coefficients, and $`a_m`$ are independent standard complex Gaussian random variables.

###### Proposition 2.2.

Random series (50) converges almost surely for all $`|z|<1`$. For all mutually distinct $`z_1,\mathrm{},z_n,`$ in the disk $`|z|<1`$ the covariance matrix

$$A_n=\left(E\left(\psi (z_p)\overline{\psi (z_p^{})}\right)\right)_{p,p^{}=1,\mathrm{},n}$$
(51)
is positive definite.

Proof. Almost sure convergence. Observe that $`(m+j)/mj+1`$ hence

$$\frac{(m+1)\mathrm{}(m+L1)}{m^{L1}}L!$$
(52)
and

$$C_{m+L1}^m=\frac{(m+L1)\mathrm{}(m+1)}{(L1)!}Lm^{L1}\text{.}$$
(53)
Consider the set

$$\mathrm{\Lambda }_N=\{a=(a_0,a_1,\mathrm{}):|a_m|mmN\}\text{.}$$
(54)
Then for all $`a\mathrm{\Lambda }_N`$ series (50) converges and

$$\underset{N\mathrm{}}{lim}\mathrm{Prob}\mathrm{\Lambda }_N=1\text{,}$$
(55)
which implies the almost sure convergence.

Positive definiteness. Consider the quadratic form

$$A_n(\mu )=\underset{p,p^{}=1}{\overset{n}{}}E\left(\psi (z_p)\overline{\psi (z_p^{})}\right)\mu _p\overline{\mu _p^{}}=E\left|\underset{p=1}{\overset{n}{}}\mu _p\psi (z_p)\right|^2\text{.}$$
(56)
Assume that $`A_n`$ is not positive definite. Then, for some nonzero vector $`\mu =(\mu _1,\mathrm{},\mu _n)`$,

$$\underset{p=1}{\overset{n}{}}\mu _p\psi (z_p)=0$$
(57)
for almost all $`\psi `$. Hence, (57) holds for almost $`\psi `$ with coefficients vector $`a\mathrm{\Lambda }_N`$. Assume that $`N>n`$. There exists a polynomial

$$\psi _0(z)=\underset{m=0}{\overset{n1}{}}\sqrt{C_{m+L1}^m}a_m^0z^m\text{,}$$
(58)
such that

$$\psi _0(z_p)=\overline{\mu _p}\text{,}$$
(59)
so that

$$\underset{p=1}{\overset{n}{}}\mu _p\psi _0(z_p)=\underset{p=1}{\overset{n}{}}|\mu _p|^20\text{.}$$
(60)
Let $`a^0=(a_0^0,a_1^0,\mathrm{},a_{n1}^0,0,0,\mathrm{})`$. From (54), it is obvious that if $`a\mathrm{\Lambda }_N`$ then $`(a+a^0)\mathrm{\Lambda }_N`$. The shift $`aa+a^0`$ is a measurable one-to-one transformation in $`\mathrm{\Lambda }_N`$, hence for almost all vectors $`a\mathrm{\Lambda }_N`$ we have that (57) holds along with

$$\underset{p=1}{\overset{n}{}}\mu _p[\psi (z_p)+\psi _0(z_p)]=0\text{.}$$
(61)
But this implies that

$$\underset{p=1}{\overset{n}{}}\mu _p\psi _0(z_p)=0\text{,}$$
(62)
which is in contradiction with (60). This contradiction proves the positive definiteness of $`A_n`$. Proposition 2.2 is proved.

Remark. The above proof gives also the positive definiteness of the covariance matrix $`A_n`$ corresponding to the random polynomial (6), provided $`N>n`$.

In this section, we will be interested in the density function $`\rho (z)`$ of the distribution of zeros of (50) in the disk $`|z|<1`$. We have from (8) that

$$E(\psi (z)\overline{\psi (z^{})})=\underset{m=0}{\overset{\mathrm{}}{}}C_{m+L1}^m(z\overline{z^{}})^m=\frac{1}{(1z\overline{z^{}})^L}\text{ .}$$
(63)
In particular,

$$E(\psi (z)\overline{\psi (z)})=f(x)\frac{1}{(1x)^L}\text{,}xz\overline{z}=|z|^2\text{ .}$$
(64)
In terms of $`f`$ the density function $`\rho (z)`$ of the zeros of (50) is given by:

$$\rho (z)=\frac{1}{\pi }\left[x\left(\frac{f^{}(x)}{f(x)}\right)^{}+\frac{f^{}(x)}{f(x)}\right]\text{ .}$$
(65)
Substituting (64) into this formula gives that

$$\rho (z)=\frac{1}{\pi }\left[x\left(\frac{L}{1x}\right)^{}+\frac{L}{1x}\right]=\frac{L}{\pi (1x)^2}\text{ ,}$$
(66)
or going back to the complex variable $`z`$, we obtain the result of Leboeuf \[L\],

$$\rho (z)=\frac{L}{\pi \left(1|z|^2\right)^2}\text{ .}$$
(67)
Numerical simulations were similar to those of the previous section. The density was simply obtained by dividing the number of counted zeros by the area of the corresponding annulus. Figure 3 shows the results of this procedure for the degree $`N`$ equal to 50, 75, 100, and 150, in comparison with the theoretical value as given by Eq. (67).

## 3 Correlations between Zeros

### 3.1 The $`SU(1,1)`$ Symmetry

We have the following important property:

###### Theorem 3.1.

The distribution of zeros of the random analytic function (50) is invariant with respect to the action of the group $`SU(1,1)`$,

$$z\frac{az+b}{cz+d},\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)SU(1,1)\text{ .}$$
(68)
Remark. This implies that all joint distribution functions of zeros of $`\psi (z)`$ are $`SU(1,1)`$ invariant.

Proof. Consider the following homogeneous analytic function of two variables:

$$\mathrm{\Psi }(z_0,z_1)=z_0^L\underset{m=0}{\overset{\mathrm{}}{}}\sqrt{C_{m+L1}^m}a_m\left(\frac{z_1}{z_0}\right)^m\text{,}|z_1|<|z_0|\text{.}$$
(69)
Then

$$\mathrm{\Psi }(1,z_1)=\psi (z_1)\text{.}$$
(70)
Let us find the covariance function of $`\mathrm{\Psi }(z_0,z_1)`$:

$`E\left(\mathrm{\Psi }(z_0,z_1)\overline{\mathrm{\Psi }(z_0^{},z_1^{})}\right)`$ $`=`$ $`\left(z_0\overline{z_0^{}}\right)^L{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}C_{m+L1}^m\left({\displaystyle \frac{z_1\overline{z_1^{}}}{z_0\overline{z_0^{}}}}\right)^m`$ (71)
$`=`$ $`\left(z_0\overline{z_0^{}}\right)^L\left(1{\displaystyle \frac{z_1\overline{z_1^{}}}{z_0\overline{z_0^{}}}}\right)^L=\left(z_0\overline{z_0^{}}z_1\overline{z_1^{}}\right)^L\text{.}`$ (72)
The action of a matrix $`ASU(1,1)`$ preserves the (1,1) indefinite scalar product $`z_0\overline{z_0^{}}z_1\overline{z_1^{}}`$. Thus, the covariance function of the Gaussian random analytic function $`\mathrm{\Psi }(z_0,z_1)`$ is $`SU(1,1)`$ invariant. In addition,

$$E\mathrm{\Psi }(z_0,z_1)=0\text{,}E\mathrm{\Psi }(z_0,z_1)\mathrm{\Psi }(z_0^{},z_1^{})=0\text{.}$$
(73)
For the Gaussian random function, the first two moments determine it uniquely. This implies that the distribution of zeros of $`\mathrm{\Psi }(z_0,z_1)`$ is $`SU(1,1)`$ invariant. Restricting this to $`z_0=1`$ we get that the distribution of zeros of $`\psi (z)`$ is $`SU(1,1)`$ invariant as well. Theorem 3.1 is proved.

Recall that a general formula for a matrix $`ASU(1,1)`$ is the following:

$$A=\left(\begin{array}{cc}a& b\\ \overline{b}& \overline{a}\end{array}\right)\text{,}|b|^2|a|^2=1\text{,}$$
(74)
so that it depends on three real parameters. An $`SU(1,1)`$-invariant metric is

$$ds^2=\frac{4(dx^2+dy^2)}{(1x^2y^2)^2}\text{,}$$
(75)
and the corresponding distance $`\tau (z_1,z_2)`$ on the disk $`\{|z|<1\}`$ is determined by the equation

$$\mathrm{tanh}(\tau /2)=\frac{|z_1z_2|}{|1z_1\overline{z_2}|}\text{.}$$
(76)
The corresponding $`SU(1,1)`$-invariant volume element is

$$\frac{4dxdy}{(1x^2y^2)^2}\text{.}$$
(77)
Theorem 3.1 implies the $`SU(1,1)`$ invariance for the normalized correlation functions.

### 3.2 Correlation Functions - Preliminaries

The $`n`$-point correlation function $`K_n(z_1,\mathrm{},z_n)`$ is determined by the condition that for any test functions $`\phi _1(z),\mathrm{},\phi _n(z)`$, which are infinitely differentiable and compactly supported in the disk $`\{|z|<1\}`$, such that their supports do not intersect,

$$E\underset{p=1}{\overset{n}{}}\left(\underset{j=1}{\overset{\mathrm{}}{}}\phi _p(\zeta _j)\right)=_^nK_n(z_1,\mathrm{},z_n)\underset{p=1}{\overset{n}{}}[\phi (z_p)dz_p]\text{,}$$
(78)
where $`\zeta _j`$ are zeros of the $`SU(1,1)`$ random analytical function $`\psi (z)`$ (we denote them by $`\zeta _j`$ to distinguish them from the variables $`z_j`$). The sum over $`j`$ is, in fact, finite because $`\phi _p(z)`$ has a compact support. The general formula for $`K_n(z_1,\mathrm{},z_n)`$ is given by the Kac-Rice expression

$$K_n(z)=𝑑\xi D_n(0,\xi ;z)\underset{p=1}{\overset{n}{}}(\xi _p\xi _p^{})\text{,}z=(z_1,\mathrm{},z_n)\text{,}\xi =(\xi _1,\mathrm{},\xi _n)\text{,}$$
(79)
where $`D_n(x,\xi ;z)`$, $`x,\xi C^n`$, is the distribution density of the two random vectors

$$X=(\psi (z_1),\mathrm{},\psi (z_n))\text{,}\mathrm{\Xi }=(\psi ^{}(z_1),\mathrm{},\psi ^{}(z_n))\text{.}$$
(80)
Formula (79) is derived in \[BSZ2\] in a much more general situation of sections of powers of a line bundle over a complex manifold, as a generalization of the original formula by Kac and Rice \[Ka\], \[Ri\] (see also \[Han\] and \[BD\]).

For the Gaussian random functions, formula (79) can be specialized as follows (see \[Han\], \[BD\], and \[BSZ2\]):

$$K_n(z)=\frac{1}{\pi ^ndetA_n}\underset{p=1}{\overset{n}{}}(\xi _p\xi _p^{})_{\mathrm{\Lambda }_n}\text{ ,}$$
(81)
where

$$A_n=\left(E\psi (z_p)\overline{\psi (z_p^{})}\right)_{p,p^{}=1,\mathrm{},n}\text{,}$$
(82)
and $`_{\mathrm{\Lambda }_n}`$ stands for averaging with respect to Gaussian complex random variables $`\xi _1,\mathrm{},\xi _n`$, such that

$$\left(E\xi _p\overline{\xi _p^{}}\right)_{p,p^{}=1,\mathrm{},n}=\mathrm{\Lambda }_n\text{,}E\xi _p=0\text{,}E\xi _p\xi _p^{}=0\text{,}p,p^{}=1,\mathrm{},n\text{,}$$
(83)
where

$$\mathrm{\Lambda }_n=C_nB_n^{}A_n^1B_n$$
(84)
and

$$B_n=\left(E\psi (z_p)\overline{\psi ^{}(z_p^{})}\right)_{p,p^{}=1,\mathrm{},n}\text{,}C_n=\left(E\psi ^{}(z_p)\overline{\psi ^{}(z_p^{})}\right)_{p,p^{}=1,\mathrm{},n}\text{.}$$
(85)
For $`n=1`$, formula (81) reduces to

$$K_1(z)=\frac{1}{\pi }\frac{A_n(z)C_n(z)B_n(z)\overline{B_n(z)}}{A_n^2(z)}\text{,}$$
(86)
where

$$A_n(z)=E\psi (z)\overline{\psi (z)}\text{,}B_n(z)=\frac{A_n(z)}{z}\text{,}C_n(z)=\frac{^2A_n(z)}{z\overline{z}}\text{.}$$
(87)
The 1-point correlation function $`K_1(z)`$ is nothing else than the density function of the distribution of zeros and formula (86) is equivalent to the Poincaré-Lelong type formula (13).

The normalized correlation function $`k_n(z_1,\mathrm{},z_n)`$ is defined as

$$k_n(z_1,\mathrm{},z_n)=\frac{K_n(z_1,\mathrm{},z_n)}{K_1(z_1)\mathrm{}K_1(z_n)}\text{.}$$
(88)
It satisfies the following theorem.

###### Theorem 3.2.

The function $`k_n(z_1,\mathrm{},z_n)`$ is invariant with respect to the action of the group $`SU(1,1)`$,

$$k_n(Az_1,\mathrm{},Az_n)=k_n(z_1,\mathrm{},z_n)\text{ ,}ASU(1,1)\text{,}$$
(89)
where

$$Az=\frac{az+b}{cz+d}\text{.}$$
(90)
Remark. This implies that $`k_n(z_1,\mathrm{},z_n)`$ is a function of pairwise distances $`\tau (z_p,z_q)`$ defined by (75).

Proof. By Theorem 3.1, the distribution $`K_n(z_1,\mathrm{},z_n)dz_1\mathrm{}dz_n`$ is $`SU(1,1)`$ invariant. Also, the distribution $`K_1(z_1)\mathrm{}K_1(z_n)dz_1\mathrm{}dz_n`$ is $`SU(1,1)`$ invariant. Hence their quotient (the Radon - Nikodim derivative),

$$\frac{K_n(z_1,\mathrm{},z_n)dz_1\mathrm{}dz_n}{K_1(z_1)\mathrm{}K_1(z_n)dz_1\mathrm{}dz_n}=k_n(z_1,\mathrm{},z_n)$$
(91)
is $`SU(1,1)`$ invariant, which was stated. Theorem 3.2 is proved.

The normalized correlation function can be expressed as a supersymmetric (Berezin) integral,

$$k_n(z_1,\mathrm{},z_n)=\frac{1}{detA_n}\frac{1}{det[I+\mathrm{\Lambda }_n(z)\mathrm{\Omega }]}𝑑\eta \text{,}$$
(92)
where $`\mathrm{\Omega }`$ is $`n\times n`$ matrix,

$$\mathrm{\Omega }=\left(\delta _{pp^{}}\eta _p^{}\overline{\eta }_p\right)_{p,p^{}=1,\mathrm{},n}\text{,}$$
(93)
and the $`\eta _p,\overline{\eta }_p`$ are anti-commuting (fermionic) variables, with $`d\eta =_pd\overline{\eta }_pd\eta _p`$. The integral in (92) is a Berezin integral, which is evaluated by simply taking the coefficient of the top degree form of the integrand $`\frac{1}{det[I+\mathrm{\Lambda }_n(z)\mathrm{\Omega }]}`$. For a derivation of formula (92), see \[BSZ4\].

### 3.3 The Two-Point Correlation Function

From (50), (67), (82), and (85), we can directly obtain the following expressions for $`A_n`$, $`B_n`$, and $`C_n`$ in the case of a two-point correlation function $`k_2(z_1,z_2)`$:

$`A_2`$ $`=`$ $`\left(\begin{array}{cc}(1|z_1|^2)^L& (1z_1\overline{z_2})^L\\ & \\ (1\overline{z_1}z_2)^L& (1|z_2|^2)^L\end{array}\right)`$ (96)
$`B_2`$ $`=`$ $`\left(\begin{array}{cc}Lz_1(1|z_1|^2)^{L1}& Lz_1(1z_1\overline{z_2})^{L1}\\ & \\ Lz_2(1\overline{z_1}z_2)^{L1}& Lz_2(1|z_2|^2)^{L1}\end{array}\right)`$ (99)
$`C_2`$ $`=`$ $`\left(\begin{array}{cc}L(1+L|z_1|^2)(1|z_1|^2)^{L2}& L(1+Lz_1\overline{z_2})(1z_1\overline{z_2})^{L2}\\ & \\ L(1+L\overline{z_1}z_2)(1\overline{z_1}z_2)^{L2}& L(1+L|z_2|^2)(1|z_2|^2)^{L2}\end{array}\right)\text{ .}`$ (102)
By Theorem 3.2, $`k_2(z_1,z_2)`$ is $`SU(1,1)`$ invariant, that is

$$k_2(Az_1,Az_2)=k_2(z_1,z_2)\text{,}ASU(1,1)\text{.}$$
(103)
The action of $`SU(1,1)`$ is transitive, hence we can move $`z_1`$ to the origin, and by rotation, we can then move $`z_2`$ to the positive half-axis. Therefore, we will assume that $`z_1=0`$ and $`z_2=r>0`$. In that case, equations (99) become

$`A_2`$ $`=`$ $`\left(\begin{array}{cc}1& 1\\ & \\ 1& (1r^2)^L\end{array}\right)`$ (106)
$`B_2`$ $`=`$ $`\left(\begin{array}{cc}0& 0\\ & \\ Lr& Lr(1r^2)^{L1}\end{array}\right)`$ (109)
$`C_2`$ $`=`$ $`\left(\begin{array}{cc}L& L\\ & \\ L& L(1+Lr^2)(1r^2)^{L2}\end{array}\right)\text{ .}`$ (112)
From (109) and (84), we further obtain

$$\mathrm{\Lambda }_2=\left(\begin{array}{cc}LL^2r^2\frac{1}{(1r^2)^L1}& LL^2r^2\frac{(1r^2)^{L1}}{(1r^2)^L1}\\ \multicolumn{2}{c}{}\\ LL^2r^2\frac{(1r^2)^{L1}}{(1r^2)^L1}& L(1+Lr^2)(1r^2)^{L2}L^2r^2\frac{(1r^2)^{2L2}}{(1r^2)^L1}\end{array}\right)\text{ .}$$
(113)
For the two-point case, formula (81) reduces to the following expression:

$$K_2(z_1,z_2)=\frac{1}{\pi ^2detA_2}\left(\mathrm{\Lambda }_{2,11}\mathrm{\Lambda }_{2,22}+\mathrm{\Lambda }_{2,12}\mathrm{\Lambda }_{2,21}\right)\text{ ,}$$
(114)
which, when combined with (88), (86), and (67), results in the following formula for the normalized two-point correlation function:

$`k_2(z_1,z_2)`$ $`=`$ $`[(1r^2)^{3L+2}+((L^22L2)r^4+(4L+4)r^21)(1r^2)^{2L}`$ (115)
$`+((L+1)^2r^4(4L+2)r^21)(1r^2)^L+1]/(1(1r^2)^L)^3\text{ ,}`$
where according to (76),

$$r=\mathrm{tanh}\left(\frac{\tau }{2}\right)=\frac{|z_1z_2|}{\left|1z_1\overline{z_2}\right|}\text{.}$$
(116)
When $`L=1`$, formula (115) simplifies to

$$k_2(z_1,z_2)=r^2(2r^2)\text{,}L=1\text{.}$$
(117)
Plots of $`k_2\left(\mathrm{tanh}\left(\frac{\tau }{2}\right)\right)`$ for $`L=1,5`$ and 50 are shown in Fig. 4.

As it can be seen from the plot, the two-point correlation function goes to 0 as $`\tau 0`$ (or, in other words, as $`r0`$). The limiting behavior can be obtained through a series expansion of (115). The following expression was obtained using Maple<sup>TM</sup>:

$$k_2(r)=\frac{1}{2}\frac{(L+1)^2}{L}r^2\frac{1}{4}\frac{(L+1)^2}{L}r^4\frac{1}{36}\frac{(L^21)^2}{L}r^6\frac{1}{72}\frac{(L^21)^2}{L}r^8+O(r^{10})\text{ ,}$$
(118)
which shows dominating quadratic behavior in the neighborhood of $`r=0`$. Thus, there is a quadratic repulsion between zeros.

We are also interested in the asymptotic behavior of the correlation function as $`L\mathrm{}`$. For this purpose, we introduce the scaling

$$r=\frac{u}{\sqrt{L}}\text{ .}$$
(119)
Substituting this into (115) and taking the limit $`L\mathrm{}`$, we obtain the following expression:

$$k_2(u)=\frac{e^{3u^2}+(u^44u^21)e^{2u^2}+(u^4+4u^21)e^{u^2}+1}{\left(e^{u^2}1\right)^3}\text{ ,}$$
(120)
which can also be written as

$$k_2(u)=\frac{\left(\mathrm{sinh}^2t+t^2\right)\mathrm{cosh}t2t\mathrm{sinh}t}{\mathrm{sinh}^3t}\text{ ,}t=\frac{u^2}{2}\text{ ,}$$
(121)
and agrees with the result obtained by Hannay \[Han\]. A plot of (121) is shown in Fig. 5.

In fact, under the scaling (119), the $`SU(1,1)`$ random analytic function converges, as $`L\mathrm{}`$, to the $`W_1`$ random analytic function (cf. \[L\]). Indeed,

$$\psi \left(\frac{u}{\sqrt{L}}\right)=\underset{m=0}{\overset{\mathrm{}}{}}\sqrt{\frac{L(L+1)\mathrm{}(L+m1)}{L^mm!}}a_mu^m\text{.}$$
(122)
As $`L\mathrm{}`$, the expression under the radical approaches $`1/m!`$, so that $`\psi (u/\sqrt{L})`$ approaches

$$\psi _0(u)=\underset{m=0}{\overset{\mathrm{}}{}}\sqrt{\frac{1}{m!}}a_mu^m\text{,}$$
(123)
which is the $`W_1`$ random analytic function.

## Appendix A  Appendix. Correlations between Inner and Outer Zeros

In this appendix, we consider the limit of correlations between zeros of the random polynomial (6) as $`N\mathrm{}`$. In the open disk $`\{|z|<1\}`$, the random polynomial (6) approaches the $`SU(1,1)`$ random analytic function (50), and the correlations between zeros of the random polynomial approach the ones of the random analytic function. However, according to formula (49), there is a $`1/(L+1)`$ fraction of the zeros outside of the unit disk. The limiting correlations between those, outer zeros, and between inner and outer zeros are described by the following theorem. Let $`k_n^L(z_1,\mathrm{},z_n)`$ denote the normalized $`n`$-point correlation function (88) corresponding to the parameter $`L`$. Indication of $`L`$ is important for the theorem. Let furthermore $`k_{nN}^L(z_1,\mathrm{},z_n)`$ denote the normalized $`n`$-point correlation function for the zeros of polynomial (6).

###### Theorem A.1.

Assume that $`|z_1|,\mathrm{},|z_m|<1`$ and $`|z_{m+1}|,\mathrm{},|z_n|>1`$. Then

$$\underset{N\mathrm{}}{lim}k_{nN}^L(z_1,\mathrm{},z_n)=k_m^L(z_1,\mathrm{},z_m)k_{nm}^1(z_{m+1}^1,\mathrm{},z_n^1)\text{,}$$
(124)
so that in the limit $`N\mathrm{}`$, the inner and outer zeros become independent, and the limiting correlations between the outer zeros coincide, after the change of variable $`z1/z`$ in the argument, with the correlations between zeros of the $`SU(1,1)`$ random analytic function with the parameter $`L=1`$.

Proof. We will first find correlations between the outer zeros and then prove the independence of the inner and outer zeros.

Correlations between outer zeros. Consider the random polynomial (6) and define another random polynomial,

$$\phi (z)=\frac{1}{\sqrt{C_{N+L1}^N}}z^N\psi (z^1)\text{.}$$
(125)
Observe that if $`z_j`$, $`|z_j|>1`$, is a zero of $`\psi (z)`$ then $`z_j^1`$ is a zero of $`\phi (z)`$ and $`|z_j^1|<1`$. Consider $`\phi (z)`$ in the disk $`|z|<1`$. From (6),

$$\phi (z)=a_N+\sqrt{\frac{N\mathrm{}(N+L2)}{(N+1)\mathrm{}(N+L1)}}a_{N1}z+\sqrt{\frac{(N1)\mathrm{}(N+L3)}{(N+1)\mathrm{}(N+L1)}}a_{N2}z^2+\mathrm{}\text{.}$$
(126)
As $`N\mathrm{}`$, the expressions under the radical approaches 1 from below. In addition, we can replace $`a_{Nm}`$ by $`a_m`$, because they are the same standard random variables. Thus, as $`N\mathrm{}`$, $`\phi (z)`$ approaches the random function

$$\phi (z)=a_0+a_1z+a_2z^2+\mathrm{}\text{,}|z|<1\text{,}$$
(127)
which is the $`SU(1,1)`$ random analytic function with the parameter $`L=1`$. Observe that $`\phi (z)`$ is a Gaussian random polynomial and for its correlations we have formula (81). From this formula, we obtain the convergence of correlations between zeros of $`\phi (z)`$ as $`N\mathrm{}`$ to the ones of the $`SU(1,1)`$ random analytic function with $`L=1`$.

Independence. We introduce the random function $`\eta (z)`$ such that $`\eta (z)=\psi (z)`$ for $`|z|1`$ and

$$\eta (z)=\frac{1}{\sqrt{C_{N+L1}^N}}z^N\psi (z)\text{,}|z|>1\text{.}$$
(128)
Then the zeros of $`\eta (z)`$ and $`\psi (z)`$ coincide. In addition, $`\eta (z)`$ is a Gaussian random field and its covariance function $`E(\eta (z)\overline{\eta (z^{})})`$ coincides with $`E(\psi (z)\overline{\psi (z^{})})`$ when $`|z|,|z^{}|<1`$, and hence, as $`N\mathrm{}`$, it approaches the correlation function (63) of the $`SU(1,1)`$ random analytic function with the parameter $`L`$. Similarly, as we saw above in this appendix, when $`|z|,|z^{}|>1`$, $`E(\eta (z)\overline{\eta (z^{})})`$ approaches the $`SU(1,1)`$ covariance function with $`L=1`$, if we replace $`z,z^{}`$ by $`z^1,z_{}^{}{}_{}{}^{1}`$, respectively. Consider now the correlation function $`E(\eta (z)\overline{\eta (z^{})})`$ when $`|z|>1`$ and $`|z^{}|<1`$. From (128),

$$E(\eta (z)\overline{\eta (z^{})})=\frac{1}{z^N\sqrt{C_{N+L1}^N}}\underset{m=0}{\overset{N}{}}C_{L+m1}^m(z\overline{z^{}})^m\text{ .}$$
(129)
Assume that $`|zz^{}|1`$. Then we can estimate the correlation function as follows:

$$\left|E(\eta (z)\overline{\eta (z^{})})\right|\frac{1}{|z|^N\sqrt{C_{N+L1}^N}}\underset{m=0}{\overset{N}{}}C_{L+m1}^m\frac{N^L}{|z|^N}\text{ .}$$
(130)
If $`|zz^{}|>1`$, then we similarly obtain that

$$\left|E(\eta (z)\overline{\eta (z^{})})\right|=\left|\frac{z_{}^{}{}_{}{}^{N}}{\sqrt{C_{N+L1}^N}}\underset{m=0}{\overset{N}{}}C_{L+m1}^m(z\overline{z^{}})^{mN}\right|N^L|z^{}|^N\text{ .}$$
(131)
Combining the two cases we can write

$$\left|E(\eta (z)\overline{\eta (z^{})})\right|N^L\left(\mathrm{max}\{|z|^1,|z^{}|\}\right)^N\text{ .}$$
(132)
This shows that the values of $`\eta (z)`$ inside and outside of the unit disk become independent as $`N\mathrm{}`$. Hence their zeros become independent. Explicit estimates for the correlations between the inner and outer zeros follow from formula (81) applied to $`\eta (z)`$.

Acknowledgement. Research of P.B. is partially supported by NSF grant #DMS-9970625 and this support is gratefully acknowledged.


text:
# HIGH ENERGY NEUTRINO FLUXES FROM COSMIC ACCELERATORS

## 1 Introduction

As several cubic-kilometer neutrino telescopes are being planned, $`^\mathrm{?}`$ it is instructive to ask about the maximum event rate one should expect for these telescopes. Evidently any answer to this questions depends on some crucial assumptions, especially concerning the nature of possible neutrino sources.

Following the discussion in Mannheim *et al*., $`^\mathrm{?}`$ we derive an upper flux limit under the assumptions that (i) protons do not escape from the acceleration region (whereas neutral particles do), (ii) the power index of the target flux density spectrum (i.e. the synchrotron photons which accelerated electrons and protons can interact with photohadronically) is given by $`\alpha =1`$, and (iii) the cosmic ray proton proton flux has the value

$$N_{p,\mathrm{obs}}(E_p)=0.8\times (E_p/1\mathrm{GeV})^{2.75}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{GeV}^1.$$
(1)
Eq. 1 is an upper limit consistent with all observations for proton energies $`E_p`$ between $`3\times 10^6`$ and $`10^{12}\mathrm{GeV}`$. Assumption (ii) implies that the bolometric gamma-ray luminosity $`L_\gamma `$ of the sources considered equals twice the neutrino luminosity $`L_\nu `$$`^\mathrm{?}`$

$$L_\gamma =2L_\nu .$$

Using particle physics, one may obtain the muon neutrino production rate from the corresponding neutron production rate: $`^\mathrm{?}`$

$$Q_{\nu _\mu }(E)83.3Q_n(25E).$$
(2)
Here and in the following particles and antiparticles are not distinguished. The cosmic ray proton flux (which is due to the decay of neutrons leaving the source) is given by

$$Q_{cr}(E_p)Q_n(E_p)\times P_{\mathrm{esc},n}(E_p),$$
(3)
where $`P_{\mathrm{esc},n}`$ denotes the neutron escape probability.

## 2 Constructing neutrino flux bounds

The recipe for constructing a generic upper bound on the neutrino flux can now be stated as follows: Given some neutron production rate one may compute the corresponding production rates of muon neutrinos and cosmic ray protons by means of Eqs. 2 and 3. In order to compute the neutrino and proton fluxes $`I`$ at Earth resulting from these production rates, one has to integrate over all sources in the universe, bearing in mind that interactions during the flight (e.g. Bethe-Heitler processes and pion production) change the particle spectra:

$$I(E)\frac{1}{4\pi }_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}M(E,z)\frac{(1+z)^2}{4\pi d_L^2}\frac{dV_\mathrm{c}}{dz}\frac{dP_{\mathrm{source}}}{dV_\mathrm{c}}Q((1+z)E,z)𝑑z.$$
(4)
Here $`z`$ denotes the redshift factor, $`V_\mathrm{c}`$ the comoving volume, $`d_L`$ the luminosity distance, and $`dP_{\mathrm{source}}/dV_\mathrm{c}`$ the redshift distribution. $`M(E,z)`$ takes into account the modification of the spectrum due to interactions during the flight. For neutrinos $`M(E,z)=1`$ holds valid.

Finally, the expressions for $`I_{\nu _\mu }`$ and $`I_{\mathrm{cr}}`$ thus obtained are normalized so that $`I_{\mathrm{cr}}`$ is tangent to the observational cosmic ray proton flux limit (Eq. 1). If gamma-rays are overproduced for this normalization (i.e. if $`L_\gamma =2L_\nu `$ is greater than the bolometric diffuse gamma-ray luminosity $`^\mathrm{?}`$, the proton and neutrino fluxes are reduced appropriately.

Fig. 1 shows several bounds obtained by means of the formalism just described:

* Generic bound for optically thin sources with $`Q_n(E_n)E_n^1\mathrm{exp}(E_n/E_{\mathrm{max}})`$, the redshift distribution of which equals that of AGN and galaxies. $`^\mathrm{?}`$ The generic bound (denoted by $`\tau _{n\gamma }<1`$) is the envelope of the flux limits computed for each $`E_{\mathrm{max}}`$.
* Generic bound for optically thick sources. Here the bound (denoted by $`\tau _{n\gamma }1`$) is derived by demanding consistency with the diffuse gamma ray background only.
* Generic bound for blazars. The construction parallels that of the generic bound for optically thin sources, but we assume a spectral break between $`10^7`$ GeV and $`10^{11}`$ GeV in the escaping cosmic ray flux due to opacity effects (still allowing for $`>\mathrm{GeV}`$ photon emission).
* Bounds for EGRET blazars and BL Lac objects. In computing the flux limit from Eq. 4 $`Q_{\mathrm{cr}}`$ and $`Q_\nu `$ are averaged over the luminosity function and redshift distribution of EGRET-detected blazars and BL Lacs, respectively.

Details of the calculation of these bounds may be found in Mannheim *et al.* $`^\mathrm{?}`$

## 3 Event rates in neutrino telescopes

So far, we have considered muon neutrinos only. However, as neutrinos are produced predominantly via pion decays, one should expect the ratio $`Q_{\nu _e}:Q_{\nu _\mu }:Q_{\nu _\tau }1:2:0`$ for the various neutrino production rates. In addition, flavor oscillations during the flight to Earth lead to a flux ratio of $`I_{\nu _e}:I_{\nu _\mu }:I_{\nu _\tau }=1:1:1`$$`^\mathrm{?}`$ Hence in the following, we will assume equal numbers of neutrinos for the three flavors and change the upper bounds accordingly.

Furthermore, the neutrino fluxes are altered due to neutrino-nucleon interactions during the crossing of the inner Earth: $`^\mathrm{?}`$

$$\frac{dI_{\nu _i}(E)}{dt}=\sigma _{i,\mathrm{tot}}(E)I_{\nu _i}(E)+\underset{k}{}_E^{\mathrm{}}𝑑E^{}\frac{d\sigma _{ki}}{dE}(E^{},E)I_{\nu _k}(E^{})+\text{decay of }\tau .$$

In this equation, $`i`$ and $`k`$ denote neutrino flavors, $`t`$ the column number density. $`d\sigma _{ki}(E^{},E)/dE`$ constitutes the differential cross section for turning a neutrino of flavor $`k`$ and energy $`E^{}`$ into one of flavor $`i`$ and energy $`E`$. $`\sigma _{i,\mathrm{tot}}`$ is the total (charged and neutral current) cross section of $`\nu _i`$.

Limiting ourselves to muon neutrinos, we may write the differential muon event rate $`d\dot{N}/dE`$ in a neutrino telescope with an effective area $`A_{\mathrm{eff}}`$ as a sum of two terms covering the muons due to neutrino-nucleon interactions inside and outside the detector:

$$\frac{d\dot{N}}{dE}=\frac{\rho }{m_p}A_{\mathrm{eff}}\left(L_E^{\mathrm{}}𝑑E^{}\frac{d\sigma _{\mathrm{CC}}}{dE}(E^{},E)I_{\nu _\mu }(E^{})+_0^{\mathrm{}}𝑑x_E^{\mathrm{}}𝑑E^{}\frac{d\sigma _{\mathrm{CC}}}{dE}(E^{},E_0(E,x))I_{\nu _\mu }(E^{})\right).$$

Here $`d\sigma _{\mathrm{CC}}/dE`$ constitutes the differential cross section for charged current interactions, $`\rho `$ the density of the detector medium (i.e. $`\rho =1\mathrm{g}/\mathrm{cm}^3`$ for water), $`L`$ the detector length, and $`m_p`$ the proton mass. $`E_0(E,x)`$ denotes the energy a muon must have so that after crossing a distance $`x`$ of the detector medium (and thus suffering radiative losses) it retains the energy $`E`$.

The event rates corresponding to various neutrino flux bounds are shown in Fig. 2 for horizontal and vertical incidence, assuming a cubic kilometer telescope.

## 4 Discussion

Fig. 2 clearly shows that whereas EGRET-blazar-like neutrino spectra might imply an event rate of about ten events per year and steradian, generic blazar models could allow neutrino event rates of several hundred events per year and steradian (most sources would then not be gamma ray emitters above GeV). This exceeds the limit obtained by Waxman and Bahcall $`^\mathrm{?}`$ for optically thin sources (which is comparable to our BL Lac case) by one to two orders of magnitude, reflecting the different assumptions.

Any derivation of an upper neutrino flux limit faces two important unknowns: $`^\mathrm{?}`$ Firstly, intergalactic magnetic fields (which affect protons but have no influence on neutrinos) might considerably change the number of neutrinos compatible with cosmic ray observations. Secondly, the same applies to exotic processes such as the decay of superheavy particles.

The planned cubic-kilometer neutrino telescopes will reach the sensitivity necessary to probe the assumptions that enter the flux-limit calculations. If photohadronic processes play an important role in blazars, we have shown that they will very likely be detected by such telescopes.

## Acknowledgments

This work was supported by the Studienstiftung des deutschen Volkes. We gratefully acknowledge the support of the ESF Neutrino Network, and we would like to thank the Moriond organizers.

## References