format: text

text:
# A Map-Making Algorithm for the Planck Surveyor

## 1 Introduction

Making CMB maps out of balloon (e.g. QMAP, de Oliveira-Costa et al. 1998; MAXIMA-1, Hanany et al. 2000; BOOMERanG, de Bernardis et al. 2000) or space borne (MAP111http://map.gsfc.nasa.gov/ and the Planck Surveyor222http://astro.estec.esa.nl/SA-general/Projects/Planck/) experiments is an important step of the standard data analysis pipeline. It allows for: a major lossless compression of the raw data with a minimum number of assumptions; checks on systematics; tests for and/or estimation of foreground contamination; etc. Map-making is also important when simulating a CMB experiment. It allows to: look for possible systematics; optimize the focal plane design and/or the experiment scanning strategy; etc. It is then highly desirable to develop tools able to attack the map-making problem also for forthcoming, high-resolution CMB space experiments.

Map-making a la COBE (see e.g. Lineweaver 1994) can be extended to the differential, high resolution MAP experiment (Wright, Hinshaw and Bennett, 1996). Moreover Wright (1996) has discussed how to perform map-making in the case of one-horned experiments significantly affected by $`1/f`$ noise, which introduces spurious correlations in the data time stream (see e.g. Janssen et al. 1996). To our knowledge this algorithm has not yet been implemented for the Planck Surveyor mission. To date, maps for Planck have only been produced by means of destriping algorithms (Delabrouille et al. 1998; Maino et al. 1999). We remind that Planck spins at 1 rpm, has a boresight angle of $`85^{}`$ and observes the same circle of the sky for 1 hour (every hour the spin axis is moved, say along the ecliptic, by $`2.5^{}`$). The destriping algorithms implemented for Planck assume that, after averaging data taken along a given circle (i.e. in one hour), the residual noise along the circle is well approximated by white noise plus an offset. These offsets are determined by minimizing the set of all possible differences between the antenna temperatures of the same sky pixel observed in two different scan circles. This method rests on a number of assumptions, the stronger being that the knee frequency of the $`1/f`$ noise is smaller than or at most comparable to the spin frequency. Wright’s method does not suffer this limitation because it assumes to know the statistical properties of the noise, directly derived from the data themselves (see e.g. Ferreira and Jaffe, 2000; Prunet et al. 2000).

The purpose of this paper is to present the first implementation of Wright’s method to the Planck Surveyor, for the moment assuming a symmetric beam profile. In fact, we want to show, to our knowledge for the first time, the analysis of the entire time stream (14 months) of Planck simulated data to produce Maximum Likelihood, minimum variance CMB maps (we stress that maps obtained from destriping algorithms are not necessarily minimum variance). The parallel, Message Passing Interface (hereafter MPI) implementation of the algorithm has been tested on a SGI Origin 2000 (hereafter O2K) and runs on time scales that might render Monte Carlo simulations feasible. This opens up the possibility of evaluating the CMB angular power spectrum via Monte Carlo techniques (Wandelt, Hivon & Górski 2000) rather than by maximizing a Likelihood (see e.g. Bond, Jaffe & Knox 1998) or by directly evaluating the map angular correlation function (Szapudi et al. 1999). We want to stress that our implementation does not assume any exact periodicity of sky pointing within single circle observation (i.e. no average on the circle is performed, as instead required by destriping algorithms). If, on the contrary, it turns out that pointing periodicity is a good assumption (i.e. average on circles is performed) then the same code runs very efficiently on medium sized workstations (e.g. Pentium based PCs), again in quite short time scales.

Although this paper might seem rather technical, we think that it can be of interest to a large community of CMB data analysts, involved either in the Planck collaboration or in the forthcoming, new generation balloon experiments.

The plan of this paper is as follows. In Sect. 2 and 3 we will briefly discuss the method and its implementation. In Sect. 4 we will show tests and benchmarks of our implemented software. In Sect. 5 we will apply our tools to Planck and BOOMERanG simulated data. Finally, in Sect. 6 we will briefly review our conclusions.

## 2 Method

### 2.1 Statement of the problem

Let us for completeness outline here the map-making algorithm and its assumptions. The primary output of a CMB experiment are the Time Ordered Data (TOD), $`𝐝`$, which consist of $`𝒩_d`$ sky observations made with a given scanning strategy and at a given sampling rate (three points per FWHM, say). A map, $`𝐦`$, can be thought as a vector containing $`𝒩_p`$ temperature values, associated with sky pixels of dimension $``$ FWHM/3. Following Wright (1996) we assume that the TOD depend linearly on the map:

$$𝐝=\mathrm{𝐏𝐦}+𝐧,$$
(1)
where $`𝐧`$ is a vector of random noise and $`𝐏`$ is some known matrix333Note that the elements of $`𝐦`$ need not to be sky temperature values. Other parameters may be fitted provided that they also depend linearly on the TOD. Sometimes these parameters are called “virtual” pixels.. The rectangular, $`𝒩_d\times 𝒩_p`$ matrix $`𝐏`$ is dubbed a pointing matrix. That is, applying $`𝐏`$ on a map “unrolls” the latter on a TOD according to a given scanning strategy. Conversely, applying $`𝐏^T`$ on the TOD “sums” them into a map444The value of a pixel of this map is the sum of all the observations of that pixel made at different times according to a given scanning strategy.. The structure of $`𝐏`$ depends on what we assume for $`𝐦`$. If $`𝐦`$ contains a pixelized but unsmeared image of the sky then $`𝐏`$ must account for beam smearing. This is the most general assumption and allows one to properly treat, for instance, an asymmetric beam profile. In this case, applying $`𝐏`$ to $`𝐦`$ implies both convolving the sky pattern with the detector beam response and unrolling $`𝐦`$ into a “signal-only” time stream. If, on the other hand, the beam is –at least approximately– symmetric then it is possible –and certainly more convenient– to consider $`𝐦`$ as the beam smeared pixelized sky. The structure of $`𝐏`$ for a one-horned experiment would then be very simple. Only one element per row would be different from zero, the one connecting the observation of $`j`$-th pixel to the $`i`$-th element of the TOD. Hereafter we will restrict ourselves to this case and we will address the implementation of an asymmetric beam profile in a forthcoming paper.

### 2.2 Least Squares Approach to Map-making

Many methods have been proposed to estimate $`𝐦`$ in Eq. (1) \[for a review see, e.g. Tegmark (1997)\]. Since the problem is linear in $`𝐦`$, the use of a Generalized Least Squares (GLS) method appears well suited. This involves the minimization of the quantity

$$\chi ^2=𝐧^T\mathrm{𝐕𝐧}=(𝐝^T𝐦^T𝐏^T)𝐕(𝐝\mathrm{𝐏𝐦})$$

for some nonsingular, symmetric matrix $`𝐕`$. Deriving with respect to $`𝐦`$ yields an estimator, say $`\stackrel{~}{𝐦}`$, for the map:

$$\stackrel{~}{𝐦}=(𝐏^T\mathrm{𝐕𝐏})^1𝐏^T\mathrm{𝐕𝐝}$$
(2)
The proof that this estimator is unbiased is straightforward. Just note that:

$$\stackrel{~}{𝐦}𝐦=(𝐏^T\mathrm{𝐕𝐏})^1𝐏^T\mathrm{𝐕𝐧}.$$

So, provided that $`𝐧=0`$, we have $`\stackrel{~}{𝐦}=𝐦`$ (the symbol $``$ indicates, as usual, an average over the ensemble). The map covariance matrix is, then:

$`𝚺^1`$ $``$ $`(𝐦\stackrel{~}{𝐦})(𝐦^T\stackrel{~}{𝐦}^T)=`$
$`=`$ $`(𝐏^T\mathrm{𝐕𝐏})^1𝐏^T𝐕\mathrm{𝐧𝐧}^T\mathrm{𝐕𝐏}(𝐏^T\mathrm{𝐕𝐏})^1`$
In order to have a “low noise” estimator, one has to find $`𝐕`$ that minimizes the variance of $`\stackrel{~}{𝐦}`$. This is attained if we take $`𝐕`$ to be the noise inverse covariance matrix, i.e. $`𝐕^1=𝐍\mathrm{𝐧𝐧}^T`$. Then $`\stackrel{~}{𝐦}`$ has the very nice property of being, amongst all linear and unbiased estimators, the one of minimum variance 555To see why this is the case, remember that any linear estimator can be written in the form $`\widehat{𝐦}=\mathrm{𝐀𝐝}=\mathrm{𝐀𝐏𝐦}+\mathrm{𝐀𝐧}`$. So, the condition $`\mathrm{𝐀𝐏}=𝐈`$ must hold for it to be unbiased. To prove that the variance of $`\widehat{𝐦}`$ can not be smaller than the variance of $`\stackrel{~}{𝐦}`$, it helps (see, e.g. Lupton 1993) to write $`\widehat{𝐦}=\stackrel{~}{𝐦}+(\widehat{𝐦}\stackrel{~}{𝐦})`$. It follows that
$$\mathrm{Var}(\widehat{𝐦})=\mathrm{Var}(\stackrel{~}{𝐦})+\mathrm{Var}(\widehat{𝐦}\stackrel{~}{𝐦})+2\mathrm{C}\mathrm{o}\mathrm{v}(\stackrel{~}{𝐦},\widehat{𝐦}\stackrel{~}{𝐦}).$$
By direct substitution, one finds that $`\mathrm{Cov}(\stackrel{~}{𝐦},\widehat{𝐦}\stackrel{~}{𝐦})=(\mathrm{𝐓𝐏}𝐈)(𝐏^T𝐍^1𝐏)^1=0`$. Thus, the variance of $`\widehat{𝐦}`$ is equal to the variance of $`\stackrel{~}{𝐦}`$ plus a nonnegative quantity, which is all we need to show. . The GLS solution to the map-making problem is, then:

$$\stackrel{~}{𝐦}=𝚺^1𝐏^T𝐍^1𝐝$$
(3)
where

$$𝚺=𝐏^T𝐍^1𝐏$$
(4)
### 2.3 ML Solutions

The statistical properties of detector noise are, as usual, described by a multivariate Gaussian distribution. This fact has two consequences. The first one is rather obvious: if the noise distribution is Gaussian, $`\stackrel{~}{𝐦}`$ is indeed the ML estimator. In fact, in the Gaussian case, the Likelihood of the data time stream given the (true) map is

$`(𝐝|𝐦)=(2\pi )^{𝒩_p/2}\times `$
$`\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}[(𝐝^T𝐦^T𝐏^T)𝐍^1(𝐝\mathrm{𝐏𝐦})+\mathrm{Tr}(\mathrm{ln}𝐍)]\right\}`$ (5)
Solving for the maximum would clearly give Eq. (3) and Eq. (4) again. The second remark has to do with the notion of a lossless map, that is, a map which contains all the relevant cosmological information contained in the TOD. A good way to prove that a method is lossless is to show that the TOD and the related map have the same Fisher information matrix. The estimator defined in Eq. (3) leads to a lossless map (Tegmark 1997).

## 3 Implementation

The method outlined in Sect. 2 has the advantage of being simple, linear and, hence, computationally appealing for Planck. In this Section we discuss our implementation of the map-making algorithm to the 30 GHz channel of the Planck Low Frequency Instrument (LFI), with a nominal resolution $`30^{}`$ FWHM. This implies $`𝒩_d10^9`$ (14 months of observation, i.e. about two sky coverages) and $`𝒩_p10^6`$. Note, however, that the implementation we present here is rather general since the algorithm does not take specific advantage of the details of Planck’s scanning strategy and instrumental performances. In fact, we want to stress that our implementation of Wright’s algorithm is, in principle, well suited for any one-horned, balloon- or space-borne CMB experiment.

### 3.1 The Noise Covariance Matrix

As already mentioned in the Introduction, the GLS method assumes to know the statistical properties of the noise. They are fully described by $`𝐍`$, the noise covariance matrix in the time domain. It is of course highly desirable to estimate $`𝐍`$ directly from the data since there is no *a priori* guarantee that the noise statistical properties match with what is measured during ground testing. The measured TOD is a combination of signal and noise. Thus, one has to estimate the signal, subtract it from the TOD and only then estimate $`𝐍`$. Ferreira and Jaffe (2000) and Prunet et al. (2000) discuss iterative methods to attack this issue.

The problem of evaluating the noise properties directly from the TOD is a bit far from the main point we want to make here: the possibility, given $`𝐍`$, to analyze the entire Planck simulated data set to produce a ML, minimum variance map. In any case, we will show in Sect. 5 that for the Planck case, to zero-th order, we can estimate the signal by projecting the TOD onto the sky. This implies applying $`𝐏^T`$ \[summing the TOD into a map\], $`(𝐏^T𝐏)^1`$ \[dividing the pixel values by the number of hits\] and $`𝐏`$ \[unrolling a map into a time stream\] to $`𝐝`$. The noise estimator, $`\stackrel{~}{𝐧}`$, can then be written as follows:

$$\stackrel{~}{𝐧}=𝐝𝐏(𝐏^T𝐏)^1𝐏^T𝐝=[𝐈𝐏(𝐏^T𝐏)^1𝐏^T]𝐧$$
(6)
and does not include any contribution from the signal. In a forthcoming paper (Natoli et al. 2001) we will address in more detail, and specifically for the Planck Surveyor, the problem of estimating the noise statistical properties directly from flight data. Here, in Sect. 5, we will only show that $`\stackrel{~}{𝐧}`$ allows a reasonably good estimate of the “in-flight” noise behavior.

### 3.2 Stationary noise

It is customary to assume that the statistical properties of detector noise do not change over the mission life time. The formal way to state this property is to write the elements of the noise covariance matrix as follows:

$$N_{ij}=\xi (|ij|)$$
(7)
or, equivalently, to say that $`𝐍`$ is a Toeplitz matrix. If $`\xi `$ vanishes for $`|ij|>𝒩_\xi 𝒩_d`$, then $`𝐍`$ is band diagonal and it can be well approximated by a circulant matrix:

$$N_{i+1,j+1}=N_{ij}\mathrm{mod}𝒩_\mathrm{d}.$$
(8)
Obviously, circularity does not hold exactly for $`𝐍`$: there is no physical reason for the noise correlation function to wrap around itself past the edges. However, the number of elements to modify in order to turn $`𝐍`$ into a circulant matrix is a mere $`𝒩_\xi (𝒩_\xi +1)𝒩_d^2`$. A circulant matrix is diagonal in Fourier space: $`𝐍=𝐅^{}𝚵𝐅`$, where $`𝐅`$ ($`𝐅^{}`$) is the discrete Fourier transform (antitransform) operator (DFT), and $`𝚵`$ is diagonal. It is trivial to show that the inverse of a circulant matrix is still circulant and, in particular, that $`𝐍^1=𝐅^{}𝚵^1𝐅`$. Circularity ensures that we will not have edge problems: this is exactly the boundary condition imposed by the DFT.

### 3.3 Non-stationary noise

The ML solution of the map-making problem \[see Eq.(3) and Eq. (4)\] can be easily generalized to the case of piecewise stationary noise. In this case, the whole TOD can be split into a set $`\{𝐝_{(\kappa )}\}`$ of smaller pieces $`\kappa =1,𝒩_c`$. Clearly,

$$𝐝_{(\kappa )}=𝐏_{(\kappa )}𝐦+𝐧_{(\kappa )},$$
(9)
where $`𝐏_{(\kappa )}`$ is a pointing matrix and $`𝐧_{(\kappa )}`$ is a vector of random, stationary noise with covariance matrix $`𝐍_{(\kappa )}`$. Neglecting cross-correlations between different $`𝐧_{(\kappa )}`$’s (a reasonable assumption if $`𝒩_\xi 𝒩_d`$) reduces the Likelihood given in Eq.(2.3) to the product of the Likelihood of the different pieces, which of course do not need to share the same noise covariance matrix. Thus, in this case, the ML map is given by

$$\stackrel{~}{𝐦}=\left[\underset{\kappa =1}{\overset{𝒩_c}{}}𝐏_{(\kappa )}^T𝐍_{(\kappa )}^1𝐏_{(\kappa )}\right]^1\underset{\kappa =i}{\overset{𝒩_c}{}}𝐏_{(\kappa )}^T𝐍_{(\kappa )}^1𝐝_{(\kappa )}.$$
(10)
Since neglecting the correlation between different $`𝐧_{(\kappa )}`$’s is *de facto* equivalent to treating different sections of the same TOD as independent experiment outcomes, we can obviously think of Eq. (10) as a simple recipe to produce maps from $`𝒩_c`$ different horns whose sky coverage is at least partially overlapped.

### 3.4 The Noise Inverse Covariance Matrix

In what follows we do not need to explicitly use $`𝐍^1`$, although we have in principle all the relevant information to do so. In fact, the first row of $`𝐍^1`$, usually called a noise filter in the literature, can be computed by taking a DFT of $`𝚵^1`$, the noise inverse spectral density. Then, by exploiting the circularity of $`𝐍^1`$ we could in principle reconstruct the whole matrix. Needless to say, handling a $`𝒩_d\times 𝒩_d`$ matrix is, for the Planck case, unconceivable even for the 30 GHz case.

Fortunately, under the assumption of circularity the application of the noise inverse covariance matrix to the TOD is simply a convolution,

$$𝐍^1𝐝=𝐅^{}𝚵^1\mathrm{𝐅𝐝},$$
(11)
and we do not need to store a huge matrix.

Asking how long a noise filter should be is a very important question. In fact, from the one hand, it is desirable to convolve the TOD with a filter properly sampled in $`𝒩_d`$ points (something already numerically feasible even for a Planck-like TOD). On the other hand, Eq. (11) works for real data in the limit $`𝒩_\xi 𝒩_d`$ (so that $`𝐍^1`$ can be considered circulant). Also note that an obvious bound on $`𝒩_\xi `$ comes from the instrument: it is difficult to imagine that the Planck detectors will remain “coherent” for the whole mission life-time. A further benefit in having $`𝒩_\xi 𝒩_d`$ is the possibility of performing a piecewise convolution in Eq. (11). In any case, in Sect. 5 we will show the sensitivity of our results to $`𝒩_\xi `$.

### 3.5 The Conjugate Gradient Solver

Evaluating $`\stackrel{~}{𝐦}`$ requires solving Eq. (3), that for convenience we rewrite as

$$𝐇^1𝐏^T𝐍^1(𝐏\stackrel{~}{𝐦})=𝐇^1𝐏^T𝐍^1𝐝$$
(12)
where $`𝐇^1`$ is a preconditioning matrix (or preconditioner). This linear system can be solved by means of a Conjugate Gradient (CG) minimization technique (see e.g. Axelsson & Barker, 1984). CG methods are known to converge efficiently \[i.e. in less than the $`𝒪(𝒩_p^3)`$ operations required by standard matrix inversion\] if the eigenvalues of $`𝚺`$ spawn a few orders of magnitudes. If this is not the case, the method should be preconditioned. This implies finding a matrix $`𝐇`$ such that its inverse approximates666To be quantitative, an optimal preconditioner would give $`𝐇^1𝚺=𝐈+𝐑`$ with all the eigenvalues of $`𝐑`$ less than 1 (Axelsson & Barker, 1984). It is then clear that multiplying $`𝚺`$ by the preconditioner changes its spectrum to give an equivalent but better conditioned system. $`𝚺^1`$. The method then solves the equivalent system given in Eq. (12). Other than being a good approximation to the original matrix, the preconditioner should also be fast to compute and invert. These often are conflicting requirements: finding a good preconditioner may be hard. In our case, we find more than satisfactory to precondition our system with the diagonal part of the symmetric and positive definite matrix $`𝚺`$. This choice is motivated by the fact that $`𝚺`$ is diagonally dominant. Furthermore, the diagonal part of $`𝚺`$ is very close to the number of hits per pixel $`𝐏^T𝐏`$ so in practice we are normalizing each row of the inverse covariance matrix to the “redundancy” of the corresponding pixel777One other problem should be mentioned: $`\mathrm{\Sigma }`$ is remarkably ill conditioned. In fact, any monopole term belongs to its kernel since applying $`𝐍^1`$ kills the corresponding TOD offset. Our preconditioner does not cure this problem which is instead fixed by specifically imposing that $`\stackrel{~}{𝐦}`$ has a given (usually zero) average. This prescription, critical if one attempts a matrix inversion, is not needed in our case: we can, of course, always subtract out the monopole..

A CG algorithm does not need to explicitly invert $`𝚺`$. We want to stress that earlier methods implemented to solve the map-making problem for balloon-borne experiments (Borrill, 1999) require the inversion of the same matrix, surely an unpleasant task when considering large maps. Such implementations usually rely on the fact that $`𝚺`$ is (analytically speaking) positive definite and, as such, a candidate to be Cholesky decomposed in $`𝒪(𝒩_p^3)`$ operations. This fact makes the procedure prohibitively expensive –even for present day supercomputers– when $`𝒩_p`$ becomes large ($`10^5`$). Also, $`𝚺`$ may not be exactly positive definite due to small inaccuracies when estimating the noise correlation properties, making Cholesky decomposition critical. Note that using an iterative solver could potentially reduce the above operation count by a factor $`𝒩_p/𝒩_{iter}`$. However, storage requirement for $`𝚺`$ would still be $`𝒪(𝒩_p^2)`$, prohibitive for Planck

To conclude this Section, let us stress that the map we obtain, $`\stackrel{~}{𝐦}`$, has the correct noise covariance matrix, $`𝚺`$, even if we never evaluate nor store it directly.

### 3.6 Parallelization

The scheme outlined above naturally lends itself to a multi-processor environment. In fact, it is straightforward to divide the entire TOD in a set of smaller, partially overlapped, time streams and force each of the Processor Elements (PE) to perform map-making with its own time stream. This is accomplished by explicitly assigning work loads to each PE by means of MPI calls embedded in the Fortran 95 code.

As stated above, the working assumption that validates Eq. (11) is that $`𝒩_\xi 𝒩_d`$. Since the number of PE, $`𝒩_{PE}`$, is at most $`𝒪(100)`$ the previous condition is even fulfilled for the single PE. Thus, each PE performs the convolution of Eq. (11) by executing its own serial FFT as discussed in Sect. 3.4. This is more efficient than spreading over the PE’s an intrinsically parallel FFT, the advantage being that we limit inter-processor communications. However, as a drawback, each PE has to handle a partial map of size $`𝒩_p`$ and cross talk among PE’s is necessary to merge $`𝒩_{PE}`$ of these maps at the end of each CG iteration. Note that increasing the number of processors $`𝒩_{PE}`$ while keeping $`𝒩_d`$ constant may result in the operation count being dominated by this merging (see however Sect. 4.2 below).

## 4 Performances

In this Section we benchmark an implementation of the map-making algorithm described above. The benchmarks are produced within the framework of the Planck’s simulation pipeline. However, since the algorithm does not take specific advantage of Planck’s scanning strategy and instrumental performance, they are indeed very general. To show that this is the case we will apply our code to a simulation of the BOOMERanG experiment (see Sect. 5.3 below).

### 4.1 Scaling with $`𝒩_\xi `$, $`𝒩_d`$ and $`𝒩_p`$

Part of the l.h.s. of Eq.(12) can be thought as the application of the $`𝐏^T𝐍^1`$ operator (the same acting on $`𝐝`$ on the r.h.s.) to $`𝐏\stackrel{~}{𝐦}^{(k)}`$, the data stream obtained by unrolling the temptative solution produced by the CG at the $`k`$-th step. While the application of $`𝐏^T𝐍^1`$ to $`𝐝`$ is done just once, the application of the same operator to $`𝐏\stackrel{~}{𝐦}^{(k)}`$ must be done for every iteration of the CG algorithm. Therefore, speed of execution becomes critical.

If we use a radix-2 Fast Fourier Transform (FFT) to perform a piecewise convolution of the TOD, the operation count is expected $`𝒩_d\mathrm{log}_2𝒩_\xi `$. The lion share of CPU time is taken by the (real) DFT transform. The use of an efficient FFT library greatly speeds up the calculation. We used the publicly available FFTW library (Frigo & Johnson 1998). This code has the nice property of tailoring itself to the architecture over which is executed, therefore greatly enhancing cross-platform efficiency. In Fig. 1 (upper panel) we plot CPU time per CG iteration, $`\tau `$, versus $`𝒩_\xi `$. The behavior shown in the figure is due to the non constant efficiency attained by the FFTW library when performing transforms of different lengths. FFT routines are usually more efficient when processing power-of-two transforms and FFTW is no exception. Thus, we fix $`𝒩_\xi `$ to a power of two.

The application of $`𝐏`$ ($`𝐏^T`$) to a map (TOD) is expected to scale linearly with $`𝒩_d`$. In Fig. 1 (middle panel) we give $`\tau `$ as a function of $`𝒩_d`$ for a given map size \[$`𝒩_p=\mathrm{786\hspace{0.17em}432}`$, as expected by choosing $`𝒩_{\mathrm{side}}=256`$ in the HEALPix pixelization (Górski et al. 1999)\] and the Planck baseline scanning strategy. The scaling with $`𝒩_d`$ is indeed linear: $`\tau =[(𝒩_d/1,400)+50]`$ms.

Since we never build up the $`𝒩_p\times 𝒩_p`$ inverse covariance matrix $`𝚺`$, but rather perform the piecewise multiplication $`𝐏^T𝐍^1(\mathrm{𝐏𝐦})`$ each time we need to do so, this operation is not expected to scale with the number of pixels. However, multiplying by the preconditioner, $`𝐇^1`$, does. When $`𝒩_p𝒩_d`$ (a condition usually fulfilled for a CMB experiment) the preconditioner should produce a negligible effect on the operation count. Conversely, we expect an increase in the operation count as, keeping $`𝒩_d`$ fixed, $`𝒩_p`$ increases. This is confirmed in Fig. 1 (lower panel). Note that, for the last two points in the plot, $`𝒩_p𝒩_d`$. Another factor that should in principle contribute to the scaling of $`\tau `$ with $`𝒩_p`$ is the performance overhead resulting from handling the arrays containing the maps. Note, however, that $`\tau `$ increases by only a factor $`2`$ when $`𝒩_p`$ is boosted by a factor $`>200`$.

All the scalings given in Fig.1 have been obtained from a single processor job using a 500 MHz Pentium III CPU.

### 4.2 Scaling with $`𝒩_{PE}`$

In order to test our parallel software we use the following simulated data sets (relative to the LFI 30 GHz channel), for which $`𝒩_d`$ changes by a couple of orders of magnitudes: $`𝒫`$60, the time stream obtained by averaging on circles the full Planck TOD ($`𝒩_d=20,196,000`$); $`𝒫`$3, a fictitious time stream obtained assuming that each circle is observed only 20 times ($`𝒩_d=403,920,000`$); $`𝒫`$, the full Planck TOD ($`𝒩_d=1,211,760,000`$). In the ideal case, the algorithm speed $`\tau ^1`$ should scale linearly with $`𝒩_{PE}`$ and should be inversely proportional to $`𝒩_d`$. So, we can define an efficiency index

$$=\frac{1}{\tau }\frac{𝒩_d}{𝒩_{PE}}$$
(13)
that we measure on an O2K. For $`𝒫`$60, increasing $`𝒩_{PE}`$ from $`1`$ to $`8`$ results in a loss of efficiency of $`50\%`$ due to inter-processor communications. If we keep $`𝒩_{PE}=8`$ but change $`𝒫`$60 to $`𝒫`$3 and $`𝒫`$, $``$ increases. We find that $``$ is basically the same when comparing $`𝒫`$60 and $`𝒩_{PE}=1`$ with $`𝒫`$ and $`𝒩_{PE}=8`$. Thus, for these cases, our code shows the ideal scaling given in Eq. 13.

### 4.3 Convergence criterion and precision

A natural criterion to stop the CG solver is to require that the fractional difference

$$ϵ^{(k)}=\frac{𝐏^T𝐍^1𝐏\stackrel{~}{𝐦}^{(k)}𝐏^T𝐍^1𝐝}{𝐏^T𝐍^1𝐝}$$
(14)
obtained at the $`k`$-th iteration step is less than a supplied tolerance. The norm $`||||`$ in Eq. (14) represents the usual Euclidean norm888An alternative convergence criterion is the more stringent $`L_{\mathrm{}}`$ norm, i.e. the absolute value of the largest element of a vector. We have tested that, for our applications, the CG solver returns almost identical results using either the Euclidean or $`L_{\mathrm{}}`$ norm. However, in the latter case, convergence is significantly slower..

To test the CG convergence efficiency we plot in Fig. 2 $`ϵ^{(k)}`$ for two datasets of different length: $`𝒩_d=20,196,000`$ (corresponding to $`𝒫`$60) and $`𝒩_d=1,211,760,000`$ (corresponding to $`𝒫`$). Both datasets have been benchmarked for the cases of “signal-only” and “noise-only” TOD. In the latter case, the noise properties are kept independent of the TOD length to facilitate comparisons. Note that, after a few steps, the precision $`ϵ^{(k)}`$ decreases exponentially and drops to the $`10^6`$ level in a few tens of iterations. Although always exponential, the rate of convergence of the CG solver depends on the TOD length. This is a consequence of having kept fixed $`𝒩_\xi `$ and $`𝒩_p`$ while increasing $`𝒩_d`$. In fact, the number of elements of $`𝚺`$ connected by the noise filter is lower for a longer TOD, a fact that makes this matrix more efficiently preconditioned by its diagonal part.

To quote a single number, the algorithm converges, to $`10^6`$ precision, in about 90 minutes ($`𝒩_{PE}=8`$) to produce a map out of the entire (*not* averaged over circles) 14 months time stream, expected from the Planck Surveyor @ 30 GHz.

Of course, solver precision and accuracy (that is, “distance” from the true solution) are not the same thing. In other words, being confident that the algorithm converges quickly is only half the story. We also want it to converge to the exact solution. This is a rather important issue because we would not like the algorithm to kill modes in the map. Discussion of this point is deferred until the next Section.

## 5 Numerical results

### 5.1 “Signal-only”

In the limiting case of a noiseless TOD (i.e. $`𝐝=\mathrm{𝐏𝐦}`$) Eq. (3) and Eq. (4) yield $`\stackrel{~}{𝐦}=𝐦`$. Thus, any good map-making code should return the “true”, pixelized sky when fed with a “signal-only” TOD. We tested our code with a pure CMB anisotropy pattern 999We use the standard Cold Dark Model with $`\mathrm{\Omega }_{CDM}=0.95`$, $`\mathrm{\Omega }_B=0.05`$ and $`h=0.5`$. No dipole was included. and the same CMB pattern contaminated by the foregrounds most relevant at 30 GHz 101010Specifically, synchrotron and free free emission were included as well as simulated emission from SZ clusters and millimetric bright sources.. The following scheme was employed to prepare the TOD. First, a “signal map” is built, pixelized at $``$ FWHM/3 using the HEALPix pixelization ($`𝒩_{\mathrm{side}}=256`$) and FWHM-smeared. The latter is “read” by simulating a 30 GHz Planck/LFI time stream. This TOD is the input to the map-making code.

In order for the code to run a noise filter, $`𝐍^1`$, is needed. First note that for this noiseless case the filter should, in principle, be irrelevant: as stated above, the GLS solution trivially gives back the “true” map independently of $`𝐍^1`$. We must keep in mind, however, that we only consider $`𝒩_\xi `$ elements of the filter. In Fourier space this translates into estimating the noise inverse spectral density precisely at $`𝒩_\xi `$ frequencies, the lowest, $`f_{\mathrm{min}}`$ say, being $`𝒩_\xi `$ times smaller than the sampling frequency. One obvious limitation on $`𝒩_\xi `$ is that it must be fixed accordingly, so that the Fourier representation of $`𝐍^1`$ well comprises the spin frequency. The amplitude of modes with $`f<f_{\mathrm{min}}`$, if any, is not recovered unless the corresponding frequencies are considered in the analysis. However, there are two reasons why one should not worry about this effect: first, and more importantly, it is small and becomes utterly negligible for $`\mathrm{}`$ a few. Furthermore, as previously stated, this contribution is likely to get lost in any case since it falls in a region of the time Fourier spectrum strongly dominated by the noise and/or, more realistically, suppressed by hardware-induced decoherence (e.g. instrumental high passing). We show (for the typical case $`𝒩_\xi =4096`$) the angular power spectra of the CMB only input and output maps together with their difference (lower panel of Fig. 3). It is clear that our map-making code recovers, as expected, the input power spectrum to an impressively high precision over the whole range of explored multipoles (although with a somewhat coarser accuarcy at low $`\mathrm{}`$’s). The same happens when considering the case of CMB plus foregrounds (see Fig. 4, upper panel): the level of the residuals between the input and output map is remarkably low, underlying the high accuracy of our code in recovering the input signal even in the presence of sharp features (see Fig. 4, lower panel).

As shown above, the code also recovers to very high precision even the broadest features in the map, i.e. the ones corresponding to low $`\mathrm{}`$’s. We want to stress that the CG solver converges to the corresponding modes in a very reasonable number of iterations ($`50`$, see Fig. 2). Other map-making codes, implementing different solving algorithms, seem not to share this property (Prunet et al. 2000) requiring more sophisticated methods to speed up the convergence of low multipoles and/or a change of variable to solve for the featureless noise-only map (see Doré et al. 2001).

A final comment before concluding this Section. As mentioned at the end of Sect. 4.3, precision and accuracy are not the same thing. We want to be sure that the CG algorithm converges to the “true” solution. So, from a given “signal-only” TOD we generate different maps by changing the solver precision. For each of these maps we evaluate the maximum difference in pixel value between them and the input map. The results are shown in Fig. 5 where we plot this maximum difference, normalized to the $`L_{\mathrm{}}`$ norm of the input map, as a function of the required solver precision $`ϵ`$. There are two regimes separated by $`ϵ10^8`$: the map accuracy decreases when $`ϵ`$ gets larger and saturates when $`ϵ`$ gets smaller. So, $`ϵ10^8`$ is just about the right number to have a good compromise between speed of convergence and accuracy. In any case, any residual map error smaller than $`1\mu \mathrm{K}`$ is acceptable in most applications, given the Planck sensitivity goal.

### 5.2 ”Noise-only”

In the opposite limiting case of a “noise-only” TOD (i.e. $`𝐦0`$), the map-making solution \[cf. Eq.(3) and Eq.(4)\] gives $`\stackrel{~}{𝐦}=(𝐏^𝐓𝐍^\mathrm{𝟏}𝐏)^\mathrm{𝟏}𝐏^𝐓𝐍^\mathrm{𝟏}𝐧`$. The efficiency of a map-making code can be tested, in this case, by examining the quality of the map obtained, a natural figure of merit being the variance of the reconstructed map or -better- its angular power spectrum.

We generate a noise time stream assuming the following form for the noise spectral density:

$$P(f)=A[1+(|f|/f_k)^\alpha ],$$
(15)
where $`f_k`$ is the knee frequency. We choose $`f_k=0.1`$ Hz (the Planck goal), $`\alpha =1`$ and an amplitude A corresponding to the expected white noise level of the $`30`$ GHz receivers. The minimum and maximum frequencies are set by the inverse of the mission life-time and the sampling frequency ($`33`$ Hz). We then generate a “noise-only” data stream that is the input to the map-making code.

Fig. 6 shows the map obtained by just averaging different observations of the same pixel (upper panel), the ML solution (middle panel) and the angular power spectra of the two (lower panel). It is clear from the middle panel of this figure that the map-making algorithm strongly suppresses the stripes due to $`1/f`$ noise, very visible in the upper panel. This is confirmed by the angular power spectrum of the ML map (see lower panel of Fig. 6) which is basically flat (as expected in the case of white and isotropic noise) for $`\mathrm{}100`$. The increasing power at $`\mathrm{}100`$ is due to two effects: residuals of spurious correlations in the noise and nonuniform sky coverage due to the Planck scanning strategy. In spite of that, the overall level of instrumental noise has been considerably lowered by the map-making algorithm.

The noise filter or, equivalently, the inverse of the noise spectral density, has been evaluated as explained in Sect. 3.1. In order to assess how good $`\stackrel{~}{𝐧}`$ is as a noise estimator, we use the same “noise-only” TOD to produce two maps: the first, by using the “true” theoretical noise spectral density \[cf. Eq.(15)\]; the second, by using the noise spectral density estimated with $`\stackrel{~}{𝐧}`$. In both cases $`𝒩_\xi =4096`$ as for the map in Fig. 6. In Fig. 7 we show, as a function of $`\mathrm{}`$, the percentage difference between the spectra derived from the two maps. So, although not optimized to the specific case of the Planck Surveyor, the noise estimator $`\stackrel{~}{𝐧}`$ reproduces to quite a good extent the in-flight noise properties. In any case, as also shown in Fig. 7, the uncertainties introduced in the estimate of the noise angular power spectrum are less than those due to cosmic variance.

As already discussed in the previous Sections, the length of the noise filter is an important issue to carefully address. One should expect that in the “noise-only” case the dependence on $`𝒩_\xi `$ could be stronger because of the $`1/f`$ tail in the noise power spectrum. To address this issue in detail and to show the sensitivity of our results to $`𝒩_\xi `$, the following scheme is employed: we generate the “true”, fully correlated noise time stream, accordingly to Eq. (15). Then, we do map-making using a noise filter with $`𝒩_\xi =𝒩_d/2`$, and evaluate the r.m.s. of this map, $`\sigma `$.111111 Note that we simulate the whole “noise-only”, $`𝒩_d`$ long, time stream by using Eq. (15) and standard FFT techniques. This time stream is obviously periodic. Thus, $`𝐍^1`$ is strictly circulant even for $`𝒩_d=𝒩_\xi /2`$. The same would not be true for real data.  Then, we generate a set of maps out of the same fully correlated noise time stream, but imposing a noise filter length, $`𝒩_\xi `$, shorter and shorter w.r.t. $`𝒩_d/2`$. In Fig. 8 we plot the fractional variation of the map r.m.s. w.r.t. $`\sigma `$ as a function of $`𝒩_\xi `$. It is clear from the figure that, for the Planck Surveyor, we can lower $`𝒩_\xi `$ by a factor of $``$ thousand w.r.t. $`𝒩_d/2`$ and change the final r.m.s. of the map by much less than one percent. This shows that increasing $`𝒩_\xi `$ above a given threshold does not induce any significant difference. This result holds even if we consider knee frequencies significantly higher than $`f_k=0.1`$ Hz: in Fig. 8 we also plot the effect of having $`f_k=1`$ Hz and $`f_k=10`$ Hz, surely two very pessimistic assumptions in the Planck context.

In principle, the noise properties of the Planck detectors, as described by Eq. (15), suggest the existence of spurious correlations in the TOD over the whole mission life-time. In practice, the error introduced by neglecting them is negligible, at least when the r.m.s. is used as a figure of merit. To stress this point even further, we evaluate the angular power spectrum from two maps generated with the same input (i.e. the same TOD), but using noise filters of different lengths. Again, we plot the percentual difference between the spectrum obtained with $`𝒩_\xi =𝒩_d/2`$ and the spectra obtained with noise filters shorter by a factor of $`64`$ and $`1024`$, respectively. Fig. 9 shows that this difference is below $`1\%`$ for $`\mathrm{}100`$, well under cosmic variance in the region of interest for recovering the cosmological parameters.

To conclude this section let us note that the algorithm discussed here is intrinsically linear \[cf. Eq.(1)\] since the map estimator can be written as the “true” map plus a noise term:

$$\stackrel{~}{𝐦}=𝐦+(𝐏^T𝐍^1𝐏)^1𝐏^T𝐍^1𝐧.$$

Our code shares this property. In fact, if we produce a map from a “signal-only” TOD and sum it to the map produced from a “noise-only” TOD the result is virtually indistinguishable (to machine precision) from the map obtained by summing the two (“signal-only” and “noise-only”) time streams.

### 5.3 Balloon Borne Experiments

As already mentioned, our code does not take any advantage from the specific Planck scanning strategy. It is then straightforwardly applicable to other one-horned experiments, such as BOOMERanG (de Bernardis et al. 2000). Here we do not want to enter in any data analysis issue. Rather, we want to show the flexibility and the efficiency of our code to process data coming from an experiment completely different from Planck.

As far as the signal is concerned, we simulate a theoretical input map 9 pixelized at 1/3 of the BOOMERanG FWHM ($`10^{}`$ @ 150 GHz) and properly smeared. Knowing its scanning strategy, we extract a typical BOOMERanG “signal-only” time stream. As far as the noise is concerned, we generate a “noise-only” time stream with the noise properties described in Eq.(15) and choosing $`A=150\mu \mathrm{K}\sqrt{\mathrm{s}}`$ and $`f_k=0.07`$ Hz. In what follows we consider a noise filter of length $`𝒩_\xi =2^{17}`$.

On the very same line of the previous subsections, we compare the theoretical input map with $`\stackrel{~}{𝐦}`$, which is reconstructed at a level of accuracy of $`10^5`$ when the CG solver precision is chosen to be $`10^6`$. The length of the noise filter does not affect, even in this case, the final results.

In the above tests we only consider the 1 degree per second (d.p.s.) section of the BOOMERanG scan. However, in order to facilitate comparisons with Planck, we truncated the BOOMERanG simulated TOD at the length of $`𝒫`$60. Not surprisingly, our code runs the BOOMERanG map-making in $`15`$ minutes on a 500 MHz Pentium III workstation, while the whole 1 d.p.s. BOOMERanG scan is processed on the same machine in about 20 minutes.

Given these time scales, a parallel machine does not seem necessary for the analysis of a single channel. However, a parallel environment becomes quite appealing if one performs a combined analysis of more than one receiver. Also, having a parallel implementation may be important to perform extensive Monte Carlo simulations of the experiment. Specific applications of this code to BOOMERanG will be discussed elsewhere.

## 6 Summary

The purpose of this work was to present a parallel implementation of a map-making algorithm for CMB experiments. In particular, we have shown for the first time Maximum Likelihood, minimum variance maps obtained by processing the entire data stream expected from the Planck Surveyor, i.e. a TOD covering the full mission life span (14 months). Here we restrict ourselves to the simple case of the Planck/LFI 30 GHz channel. However, the extension of our implemented software to other, higher frequency, channels is straightforward. At the moment our implementation is limited to a symmetric antenna beam profile.

The code was shown to scale linearly with the number of time samples $`𝒩_d`$, logarithmically with the noise filter length and very slowly with the number of pixels $`𝒩_p`$, the latter scaling being mostly influenced by the choice of the preconditioner. Furthermore, on a multiprocessor environment the code scales nearly optimally with the number of processor elements.

The code is extremely accurate and fast. We have shown that it recovers the true, input map to a remarkably high accuracy, over the whole range of considered multipoles. Moreover, we run the entire 30 GHz simulated TOD in $``$ 90 minutes on an 8 processor job on a O2K machine, and the BOOMERanG case in only 20 minutes on a 500 MHz Pentium III workstation. On the latter machine, the Planck case runs in about 15 minutes if we perform the average on circles.

A key problem for a successful map-making pipeline is to obtain a reliable estimate of the noise properties directly from flight data. It is straightforward to further iterate the GLS solution implemented here to get more accurate estimates of the underlying noise. Here we have explicitly chosen not to do so. Rather, we stop to the zero-th order solution, $`\stackrel{~}{𝐧}`$, showing that this is more than enough for the purposes discussed here.

The possibility of running map-making algorithms on reasonably short time scales opens up the possibility of evaluating the CMB angular power spectrum via Monte Carlo simulations. We will address this point in a forthcoming paper, together with the extension to non-symmetric antenna beam profile.

###### Acknowledgements.

We are indebted with C. Burigana and D. Maino for helping us with the generation of the TOD sets. We thank A. Balbi, P. Cabella, D. Marinucci and C. di Fiore for useful suggestions and comments. We also thank P. de Bernardis and the BOOMERanG collaboration for having provided us with the BOOMERanG scan and instrumental performances. PN acknowledges fruitful discussions with J. Staren and E.L. Wright. We acknowledge use of the HEALPix package and of the FFTW library. The supercomputing resources used for this work have been provided by Cineca.


text:
# Neutron Diffusion and Nucleosynthesis in an Inhomogeneous Big Bang Model

## 0.1 Introduction

In 1999 the author wrote a big bang nucleosynthesis code corresponding to a universe with an inhomogeneous baryon distribution. The code envisions the universe as a lattice of spheres and treats one sphere as characteristic of the rest of the universe. The sphere has an initial radius ( the distance scale $`r_i`$ ) of 25000 cm. at temperature $`T=`$ 100 GK. The sphere consists of a high density inner region and a low density outer region. The boundary radius $`r_b`$ between the regions is set to $`r_i/2`$ and the density contrast ratio $`R_\rho =`$ 800:1. The baryon-to-photon ratio $`\eta `$ in this model is 3 $`\times 10^{10}`$ ( $`\eta _{10}=`$ 3 ). The sphere is divided into a core and 63 spherical shells.

The code solve for the number density $`n(i,s)`$ of isotope species $`i`$ contained in spherical shell $`s`$. The evolution of $`n(i,s)`$ obeys the equation (?)

$`{\displaystyle \frac{dn(i,s)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{n_b(s)}}{\displaystyle \underset{j,k,l}{}}N_i\left({\displaystyle \frac{n^{N_i}(i,s)n^{N_j}(j,s)}{N_i!N_j!}}[ij]+{\displaystyle \frac{n^{N_k}(k,s)n^{N_l}(l,s)}{N_k!N_l!}}[kl]\right)`$
$`3\dot{\alpha }n(i,s)+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left(r^2D_n{\displaystyle \frac{n(i,s)}{r}}\right)`$
The first term accounts for nuclear and weak reactions occuring between isotopes within shell $`s`$. The second term corresponds to the expansion of the universe. And the last term accounts for diffusion of isotope $`i`$ between shells, though in the model only neutrons diffuse.

## 0.2 Results for $`r_i=`$ 25000 cm, $`\eta _{10}=`$ 3

Figure (1) shows the number densities of free neutrons ( $`n(\mathrm{n},s)`$ ), protons ( $`n(\mathrm{p},s)`$ ) and 4 $`\times ^4`$He ( $`4n(^4\mathrm{H}e,s)`$ ). The horizontal axes correspond to the temperature $`T`$ from 100 GK down to 0.01 GK and the distance scale $`r_i`$ normalized to 64 dimensionless units, and the vertical axis is number density $`n(i,s)`$ normalized to the initial baryon density ( $`n_{b0}`$ ). At $`T=`$ 100 GK number densities in the inner shells vs. the outer shells are large. Weak reactions convert neutrons to protons until they no longer become significant at around $`T=`$ 13 GK. Figure (2a) shows the neutron-proton conversion rates for shell 3, a high density inner shell ( normalized to $`n_{b0}`$ and to the expansion rate of the universe $`\dot{\alpha }`$ ). Figure (2a) also shows that around $`T=`$ 6 GK neutrons start to diffuse significantly into the lower density outer shells. But since the conversion rates are no longer effective at this temperature the protons in the inner shells remain there. At around $`T=`$ 2.5 GK the neutrons are homogeneously distributed throughout the whole sphere.

Nucleosynthesis occurs at about $`T=`$ 0.9 GK, when the nuclear reaction n + p $``$ d + $`\gamma `$ falls out of equilibrium. Other nuclear reactions then build heavier nuclei, but most baryons wind up as part of 4He nuclei. Because of higher density nucleosynthesis occurs earlier in the inner shells, depleting neutrons. So neutrons diffuse back into the inner regions as shown on the right of Figure (2a). Figure (2b) shows the same reaction rates for shell 62, a low density outer shell. The rate for n + p $``$ d + $`\gamma `$ is lower and flatter than in Figure (2a). Figure (1) shows that most neutrons back diffuse into the inner shells and then undergo nucleosynthesis, leading to 4He production overwhelmingly concentrated in the inner shells.

The final mass fraction $`X_{{}_{}{}^{4}\mathrm{H}e}`$ of the whole model is 0.243664. This is the result for $`r_i=`$ 25000 cm and $`\eta _{10}=3`$ Table (1) shown below lists $`X_{{}_{}{}^{4}\mathrm{H}e}`$ for $`r_i=`$ 100 cm up to 6.3 $`\times 10^6`$ cm. For $`r_i`$ less than 25000 cm neutron diffusion occurs early enough to coincide at least partly with the time when the conversion reactions are significant. In the inner regions the reactions convert protons into neutrons that then diffuse to the outer regions where they get converted back to protons. But as $`r_i`$ gets larger more protons remain in the inner regions and $`X_{{}_{}{}^{4}\mathrm{H}e}`$ increases. $`r_i=`$ 25000 cm is the turning point when all the protons in the innermost regions remain there. For larger $`r_i`$ the back diffusing neutrons can’t reach the innermost regions before nucleosynthesis is completed. So with fewer neutrons in the inner regions $`X_{{}_{}{}^{4}\mathrm{H}e}`$ decreases. At $`r_i=7.9\times 10^5`$ cm though diffusion becomes significant at $`T=`$ 0.9 GK, right during nucleosynthesis. For larger $`r_i`$ nucleosynthesis occurs before neutron diffusion and the inner regions, with high proton and neutron densities, produce enough 4He to raise $`X_{{}_{}{}^{4}\mathrm{H}e}`$ significantly.

| $`r_i`$ in cm. | $`X_{{}_{}{}^{4}\mathrm{H}e}`$ |
| --- | --- |
| $`10^2`$ | 0.239029657 |
| 1.58489 $`\times 10^3`$ | 0.239369667 |
| 2.51189 $`\times 10^4`$ | 0.243663794 |
| $`10^5`$ | 0.241847338 |
| 7.94328 $`\times 10^5`$ | 0.229781933 |
| 6.30957 $`\times 10^6`$ | 0.250711334 |

Table (1): 4He, where $`\eta =3.0\times 10^{10}`$.

## 0.3 Conclusions and Future Research

Figure (3a) shows $`X_{{}_{}{}^{4}\mathrm{H}e}`$ for a range of $`\eta `$ from $`10^{10}`$ to $`3\times 10^9`$ as well as the range of $`r_i`$ listed in Table (1) For all values of $`\eta `$ shown one can see the turning points at which $`X_{{}_{}{}^{4}\mathrm{H}e}`$ increases and then decreases and then increases again. The code can generate contour maps for 67 other isotope species as well. Comparing the contour maps of $`X_{{}_{}{}^{4}\mathrm{H}e}`$, $`Y(\mathrm{d})/Y(\mathrm{p})`$ ( $`Y(\mathrm{d})`$ is the abundance of deuterium ), and $`Y(^7\mathrm{L}i+^7\mathrm{B}e)/Y(\mathrm{p})`$ to the $`2\sigma `$ ranges of recent observations ( IT, 1998; (?) BT, 1998; (?) PWSN, 1998 (?) ).

$`X_{{}_{}{}^{4}\mathrm{H}e}`$ $`=`$ $`0.244\pm 0.002`$
$`Y(\mathrm{d})/Y(\mathrm{d})`$ $`=`$ $`(3.4\pm 0.3)\times 10^5`$
$`\mathrm{log}[Y(^7\mathrm{L}i)/Y(\mathrm{d})]`$ $`=`$ $`9.45\pm 0.20`$
the author determined an observational range for $`\eta `$ from 4.7 $`\times 10^{10}`$ to 6 $`\times 10^{10}`$.

A planned article will describe 4He production in greater detail, describing results for values of $`r_i`$ other than 25000 cm. $`R_\rho `$, $`r_b`$ and the geometry of the model can be varied in the code. In the future the author will look at results with different values for these parameters. Also, the code has written into it subroutines to calculate the neutrino heating effect. Electrons and positrons annihilate each other from $`T=`$ 5 GK to $`T=`$ 1 GK. Through annihilation and scattering they can transfer some energy to neutrinos, which in turn can alter the neutron-proton conversion reaction rates and ultimately 4He production. Figure (3b) shows the change $`\mathrm{\Delta }X_{{}_{}{}^{4}\mathrm{H}e}`$ of the mass fraction of 4He for the same ranges of $`\eta `$ and $`r_i`$ as in Figure (3a). $`\mathrm{\Delta }X_{{}_{}{}^{4}\mathrm{H}e}`$ imitates $`X_{{}_{}{}^{4}\mathrm{H}e}`$, but not exactly. The change also remains within an order of $`10^4`$. The author will determine if the change is equivalent to a constant shift of Figure (3a) to the left, and if the coding for the neutrino heating effect is correctly represented in the code.


text:
# Contents

## 1 Introduction

The connection of the orthogonal polynomials theory with the group theory as well as the connection the theory of quantum (deformed) groups and algebras with the theory of basic hypergeometric functions (deformed polynomials)  is well known . In particular, the Hermite polynomials system (after multiplication on $`\mathrm{exp}(x^2)`$) is a system of eigenfunctions for quantum-mechanical harmonic oscillator  energy operator. Besides, these polynomials are connected with representations of the algebra generated by Heisenberg commutation relations (and with representations of the appropriate group as well). Similarly, many known $`q`$-deformed Hermite polynomials  make up eigenfunctions systems of the energy operators for the related deformed oscillators . Moreover, these polynomials naturally arise from the analysis of irreducible representations of appropriate deformed algebras.

At the standard approach it is suggested that a representation of considered algebra (or group) is prescribed. Then we look for a system of orthogonal polynomials which make up the eigenfunctions system of a significant operator of this algebra (for example, Hamiltonian). In other words, we look for the basis in the representation space, in which this operator take the diagonal form.

In  we considered in some sense an inverse task . Let $`\mu `$ be a positive Borel measure on the real line $`R^1`$. We consider a set of real orthonormal polynomials in the space $`L^2(R^1;\mu (dx))`$. For the orthogonal polynomial system one wants to construct oscillator like algebra, such that these polynomials make up the eigenfunctions system of the energy operator corresponding to this oscillator. By the full solution of this problem it is suggested that a differential or a difference equation for the given system of polynomials is known. In work  (strictly) classical polynomials, (namely, the Hermite, Jacobi and Laguerre polynomials) were considered as examples. From these examples it follows that the central problem with obtaining the differential (or difference) equation for the given system of polynomials was in finding a realization of the annihilation (lowering) operator by a differential operator. Note that the operator of a finite order in work  occurres only for the standard quantum mechanical oscillator connecting with the orthonormal Hermite polynomial system. In all other cases considered in  the annihilation operator connected with the orthonormal classical polynomial system was realized by a differential operator of the infinite order. More exactly these operators were described as rather simple functions of the first order differential operator. Other examples of realization of the annihilation operator by the known functions of the first order differential operator are considered in works .

In the present work we shall obtain a expression for an annihilation operator of the generalized deformed oscillator algebra related to a complete system of the real orthonormal polynomials. The above-mentioned examples from  are concrete realizations of the general expressions obtained here.

Let us remark that the operators of the finite order arise for the Hermite-Chihara polynomial systems (see,   and ) and the Hermite-Hahn polynomial systems .

It is more interesting to apply the obtained results to investigating a generalized oscillator connected with symmetric basic number $`[a]={\displaystyle \frac{q^aq^a}{qq^1}}`$ of quantum groups theory. This case is still not clearly understood. The research of this case is in a stage of completion now.

For the convenience of the reader in the following section we recall some necessary information from the theory of Jacobi matrices and the classical moment problem.

## 2 Background information

Let a operator $`X`$ is given by action on the standard orthonormal basis $`\left\{e_n\right|nZ_+\}`$ in the Hilbert space $`H=l^2(Z_+)`$:

$$Xe_n=b_ne_{n+1}+a_ne_n+b_{n1}e_{n1},a_n0,b_n.$$
(1)
Then $`X`$ can be represent in $`H`$ by an infinite 3-diagonal matrix

$$X=\left(\begin{array}{cccccc}a_0& b_0& 0& 0& \mathrm{}& \mathrm{}\\ b_0& a_1& b_1& 0& \mathrm{}& \mathrm{}\\ 0& b_1& a_2& b_2& \mathrm{}& \mathrm{}\\ 0& 0& b_2& a_3& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right),$$
(2)
known as the Jacobi matrix.

With the Jacobi matrix it is related the set of polynomials $`P_n(x)`$ of degree $`N`$, satisfying the recurrent relation,

$$b_nP_{n+1}(x)+a_nP_n(x)+b_{n1}P_{n1}(x)=xP_n(x),$$
(3)
and subjecting to the following ”initial conditions”

$$P_0(x)=1,P_1(x)=0.$$
(4)
Further, we shall assume that $`a_n,b_n`$ and $`a_n=0,b_n>0`$. Under such conditions polynomials $`P_n(x)`$ have real coefficients and fulfill the following parity conditions $`P_n(x)=(1)^nP_n(x).`$ Notice that the recurrent relations (3) have two linearly independent solutions. The polynomials $`P_n(x)`$, which solve to (3) and subject to the initial conditions (4), are called polynomials of the first kind. Besides, there are independent set of solutions consisting from polynomials $`Q_n(x)`$, which satisfy with the following initial conditions

$$Q_0(x)=0,Q_1(x)=\frac{1}{b_0}.$$
(5)
For these polynomials the standard agreement $`\mathrm{deg}Q_n(x)=n1`$ is carried out. These polynomials are known as polynomials of the second kind for the Jacobi matrix $`X`$ (2).

The polynomials of the first and second kind are related by the following equation

$$P_{n1}(x)Q_n(x)P_n(x)Q_{n1}(x)=\frac{1}{b_{n1}},(n=1,2,3,\mathrm{}).$$
(6)
It is known that polynomials $`\left\{P_n(x)\right\}_{n=0}^{\mathrm{}}`$ are orthonormal

$$_{}P_n(x)P_m(x)𝑑\sigma (x)=\delta _{n,m}$$
(7)
with respect to a positive Borel measure $`\sigma `$ on $``$. In order to evaluate the measure $`\sigma `$ one must solve the following task.

By $`s_n`$ we denote the factor $`\alpha _0`$ at $`P_0(t)`$ in the expansion

$$t^n=\underset{k=0}{\overset{n}{}}\alpha _kP_k(t).$$

For the given numerical sequence $`\left\{s_n\right\}_{n=0}^{\mathrm{}}`$ it is required to find a positive measure $`\sigma `$ on the real line such that the following equalities:

$$s_n=_{\mathrm{}}^{\mathrm{}}t^n𝑑\sigma (t),n=0,1,2,\mathrm{}$$
(8)
are hold. This task is known as the degree Hamburger moment problem .

If the system (8) defines the measure uniquely, then the moment problem is called determined one. Otherwise there exists an infinite set of such measures, and the moment problem is called undetermined one.

The matrix (2) determines a symmetrical operator $`X`$ in $`H.`$ It is known that the deficiency indices for the operator $`X`$ are equal $`(0,0)`$ or $`(1,1).`$

###### Proposition 2.1.

Let $`\overline{X}`$ be a closure of an operator $`X`$ in $`H`$. The following conditions are equivalent:

1. the deficiency indices of an operator $`X`$ are equal $`(0,0)`$;
2. an operator $`\overline{X}`$ is selfadjoint in the Hilbert space $`H`$;
3. the moment problem (8) for the given number set $`\left\{s_n\right\}_{n=0}^{\mathrm{}}`$ such that

$$s_n=(e_0,X^ne_0)$$
(9)
is determined;
4. a series $`_{n=0}^{\mathrm{}}\left|P_n(z)\right|^2`$ is divergent for all $`z`$ such that $`\mathrm{Im}z0`$.

Let the operator $`X`$ has the deficiency indices $`(1,1)`$. Then the operator $`\overline{X}`$ not selfadjoint and there exists infinite number of its selfadjoint extensions. In this case series $`_{n=0}^{\mathrm{}}\left|P_n(z)\right|^2<\mathrm{}`$ is divergent at all $`z`$ such that $`\mathrm{Im}z0`$, and moment problem (8) for a set of numbers $`s_n`$ from (9) is indetermined.

The following proposition gives a sufficient condition under which the operator $`X`$ has the deficiency indices $`(1,1)`$.

###### Proposition 2.2.

Under the following conditions

1. there exists $`N`$ such that for all $`n>N`$ the inequality

$$b_{n1}b_{n+1}b_n^{\mathrm{\hspace{0.17em}\hspace{0.17em}2}};$$
(10)
2. $$\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{b_n}<\mathrm{};$$
(11)

the operator $`X`$ has the deficiency indices equal to $`(1,1).`$

Moreover, we have the following statement.

###### Proposition 2.3.

(see ). If the conditions (10) and (11) are fullfiled, there is an infinite number of selfadjoint extensions of the operator $`\overline{X}`$, and related moment problem (8) for a set of numbers $`\left\{s_n\right\}_{n=0}^{\mathrm{}}`$ from (9) is undetermined.

The polynomials $`P_n`$ and $`Q_n`$ can are represented in the form

$$P_n(x;q)=\underset{m=0}{\overset{ϵ(\frac{n}{2})}{}}\frac{(1)^m}{\sqrt{\left\{n\right\}!}}b_0^{2mn}\alpha _{2m1,n1}x^{n2m}$$
(12)
$$\alpha _{1;n1}1;\alpha _{2m1;n1}=\underset{k_1=2m1}{\overset{n1}{}}\left\{k_1\right\}\underset{k_2=2m3}{\overset{k_12}{}}\left\{k_2\right\}\mathrm{}\underset{k_m=1}{\overset{k_{m1}2}{}}\left\{k_m\right\},m1$$
(13)
$$Q_{n+1}(x;q)=\underset{m=0}{\overset{ϵ(\frac{n}{2})}{}}\frac{(1)^m}{\sqrt{\left\{n+1\right\}!}}b_0^{2mn}\beta _{2m,n}x^{n2m}$$
(14)
$$\beta _{0;n}1;\beta _{2m;n}=\underset{k_1=2m1}{\overset{n}{}}\left\{k_1\right\}\underset{k_2=2m2}{\overset{k_12}{}}\left\{k_2\right\}\mathrm{}\underset{k_m=2}{\overset{k_{m1}2}{}}\left\{k_m\right\},m1,$$
(15)
where $`\{s\}={\displaystyle \frac{b_{s1}^2}{b_0^2}}`$, and the integral part of $`x`$ is denoted by $`ϵ(x)=\mathrm{Ent}(x)`$. Because the representation (12) - (15) of the polynomials $`P_n`$ and $`Q_n`$ have not been adequately explored, we shall recall some properties of coefficients $`\alpha _{m;n}`$ and $`\beta _{m;n}.`$

From the recurrent relation (3) it follows that the coefficients $`\alpha _{m;n}`$ and $`\beta _{m;n}`$ satisfy the conditions:

$`\alpha _{2m1;n}`$ $`=\left\{n\right\}\alpha _{2m3;n2}+\alpha _{2m1;n1};`$ (16)
$`\beta _{2m;n}`$ $`=\left\{n\right\}\beta _{2m2;n2}+\beta _{2m;n1}.`$ (17)
From the definitions (13) and (15) of coefficients $`\alpha _{k,n}`$ and $`\beta _{k,n}`$ it follows that

$$\alpha _{2k1,2k1}=\{2k1\}!!,\beta _{2k2,2k2}=\{2k2\}!!.$$
(18)
Here

$$\left\{2n\right\}!!=\left\{2n\right\}\left\{2n2\right\}\mathrm{}\left\{2\right\},\left\{2n1\right\}!!=\left\{2n1\right\}\left\{2n3\right\}\mathrm{}\left\{1\right\}.$$

## 3 Statement of a problem

Let $`=L^2(;\mu (\mathrm{dx}))`$ be a Hilbert space, and $`\left\{\mathrm{\Psi }_n(x)\right\}_{n=0}^{\mathrm{}}`$ be a system of polynomials, which are orthonormal with respect to the measure $`\mu `$. The recurrent relations of these polynomials take the following form:

$$b_n\mathrm{\Psi }_{n+1}(x)+b_{n1}\mathrm{\Psi }_{n1}(x)=x\mathrm{\Psi }_n(x),b_n>0,$$
(19)
$$\mathrm{\Psi }_0(x)=1,\mathrm{\Psi }_1(x)=0.$$
(20)
In the work  it was shown, that one can construct the oscillator-like algebra $`𝒜_\mathrm{\Psi }`$ corresponding to this polynomial system. In addition, the space $``$ is a space of the Fock representation, and the polynomials $`\left\{\mathrm{\Psi }_n(x)\right\}_{n=0}^{\mathrm{}}`$ make up the Fock basis. The generators $`a_{\mu }^{}{}_{}{}^{+},a_{\mu }^{}{}_{}{}^{},N`$ of this algebra acts in the Fock representation as follows

$`a_{\mu }^{}{}_{}{}^{+}\mathrm{\Psi }_n(x)`$ $`=\sqrt{2}b_n\mathrm{\Psi }_{n+1}(x),`$ $`n0;`$
$`a_{\mu }^{}{}_{}{}^{}\mathrm{\Psi }_n(x)`$ $`=\sqrt{2}b_{n1}\mathrm{\Psi }_{n1}(x),`$ $`n1a_{\mu }^{}{}_{}{}^{}\mathrm{\Psi }_0(x)=0;`$ (21)
$`N\mathrm{\Psi }_n(x)`$ $`=n\mathrm{\Psi }_n(x),`$ $`n0;`$
From these relations it follows that

$$a_{\mu }^{}{}_{}{}^{}a_{\mu }^{}{}_{}{}^{+}\mathrm{\Psi }_n(x)=2b_{n}^{}{}_{}{}^{2}\mathrm{\Psi }_n(x),a_{\mu }^{}{}_{}{}^{+}a_{\mu }^{}{}_{}{}^{}\mathrm{\Psi }_n(x)=2b_{n1}{}_{}{}^{2}\mathrm{\Psi }_{n}^{}(x),n0,$$
(23)
and $`b_1=0.`$ Then we have the following commutation relations:

$$[a_{\mu }^{}{}_{}{}^{},a_{\mu }^{}{}_{}{}^{+}]=2\left(B(N+1)B(N)\right),[N,a_{\mu }^{}{}_{}{}^{\pm }]=\pm a_{\mu }^{}{}_{}{}^{\pm },$$
(24)
between the generators of the algebra $`𝒜_\mathrm{\Psi }.`$ The function $`B(N)`$ is defined by the formula

$$B(N)\mathrm{\Psi }_n(x)=b_{n1}{}_{}{}^{2}\mathrm{\Psi }_{n1}^{}(x).$$
(25)
If the sequence $`\left\{b_n\right\}_{n=0}^{\mathrm{}}`$ satisfies the recurrent relation

$$b_n{}_{}{}^{2}Qb_{n1}{}_{}{}^{2}=C(n),$$
(26)
then, as follows from (25), the relation

$$a_{\mu }^{}{}_{}{}^{}a_{\mu }^{}{}_{}{}^{+}Qa_{\mu }^{}{}_{}{}^{+}a_{\mu }^{}{}_{}{}^{}=2C(N).$$
(27)
is fulfilled also.

According to (12) and (13), we have

$$\mathrm{\Psi }_n(x)=\underset{m=0}{\overset{ϵ(n/2)}{}}\frac{(1)^m}{\sqrt{\{n\}!}}\alpha _{2m1,n1}b_0^{2mn}x^{n2m},n0.$$
(28)
For simplification of further formulas we define the operator $`A`$ by a relation $`a^{}=\sqrt{2}A.`$ Then

$$A\mathrm{\Psi }_n=b_{n1}\mathrm{\Psi }_{n1},n1,A\mathrm{\Psi }_0=0;$$
(29)
$$A\mathrm{\Psi }_0=0a_{n0}=0,n0$$
(30)
The main purpose of the present work is to find, using the formulas (29) and (30), the conditions on the coefficients $`a_{ks}`$, under which the operator $`A`$ take the form

$$A=\underset{k,s=0}{\overset{\mathrm{}}{}}a_{ks}x^k\frac{\mathrm{d}^\mathrm{s}}{\mathrm{d}x^s}.$$
(31)
Substituting (31) and (3) into (29), one obtains

$$\begin{array}{c}\underset{k,s=0}{\overset{\mathrm{}}{}}a_{ks}\underset{m=0}{\overset{ϵ(n/2)}{}}\frac{(1)^m}{\sqrt{\left\{n\right\}!}}\alpha _{2m1,n1}b_0^{\mathrm{\hspace{0.17em}2}mn}x^{n+k2ms}\frac{\left(n2m\right)!}{\left(n2ms\right)!}=\hfill \\ \hfill b_{n1}\underset{m=0}{\overset{ϵ(\frac{n1}{2})}{}}\frac{(1)^m}{\sqrt{\left\{n1\right\}!}}\alpha _{2m1,n2}b_0^{\mathrm{\hspace{0.17em}2}m+1n}x^{n12m}\end{array}$$
(32)
Equating the coefficients at $`x^t`$ in both sides of a relation (32), we receive

$$a_{n,0}=0,n0$$

$$\begin{array}{c}\underset{s=0}{\overset{\mathrm{}}{}}\underset{m=\delta (\frac{ns1}{2})}{\overset{ϵ(\frac{ns}{2})}{}}(1)^ma_{t+s+2mn,s}\alpha _{2m1,n1}b_0^{\mathrm{\hspace{0.17em}2}mn}\frac{\left(n2m\right)!}{\left(n2ms\right)!}=\hfill \\ \\ \hfill \{\begin{array}{ccc}0& if& n1t=2p+1,p0\\ \sqrt{\left\{n\right\}}b_{n1}(1)^p\alpha _{2p1,n2}b_0^t& if& n1t=2p,p0\end{array},\end{array}$$
(33)
where the least integer greater than $`a`$ or equal to $`a`$ is denoted by $`\delta (a)`$.

From (33) it follows that the problem of a realization of the annihilation operator by a differential operator (31) was rather cumbersome in general case. For this reason we slightly change the statement of a problem. Namely, we will seek for not an annihilation operator but a ”lowering operator” $`A`$ (denoted by the same symbol). In other words from (33) we should look for a operator $`A`$ such that:

$$A\mathrm{\Psi }_n=\gamma _n\mathrm{\Psi }_{n1},n1,A\mathrm{\Psi }_0=0,$$
(34)
under the supposition, that the operator (31) works in space $`H`$ as a shift operator. Thus we supposed that among coefficients $`\left\{a_{k,s}\right\}_{k,s=0}^{\mathrm{}}`$ in the matrix of this operator only elements on the first upper diagonal are different from zero

$$a_{n,k}=\{\begin{array}{ccc}0& if& kn+1\\ a_{n,k}0& if& k=n+1\end{array}$$
(35)
Let us remark that condition

$$a_{n,0}=0,n0$$
(36)
and relations (32) and (33) (with a replacement $`b_{n1}\gamma _n`$ in a right hand side) are hold.

In the following section we shall obtain some consequence of the relations (34) and (35).

## 4 Conditions under which the matrix of an operator $`A`$ has non zero elements only upper above diagonal

Taking into account that all elements $`a_{k,s}`$ except for $`a_{k,k+1}`$ are equal to zero, and using the relation 33 for $`t=0,1,2,\mathrm{}`$ and $`P=0,1,2,\mathrm{}`$, we obtain the following conditions on coefficients $`\left\{b_n\right\}_{n=0}^{\mathrm{}}`$ and $`\left\{\gamma _n\right\}_{n=1}^{\mathrm{}}:`$

$$\sqrt{\{2p+1\}}\gamma _{2p+1}=\frac{\alpha _{2p1,2p}}{\alpha _{2p1,2p1}}\gamma _1,p1$$
(37)
$$\sqrt{\left\{2p+2\right\}}\gamma _{2p+2}=\frac{\alpha _{2p1,2p+1}}{\alpha _{2p1,2p}}\sqrt{\left\{2\right\}}\gamma _2,p0$$
(38)
$$\frac{a_{0,1}}{k!}+\frac{a_{1,2}}{(k1)!}+\mathrm{}+\frac{a_{k1,k}}{1!}+a_{k,k+1}=\frac{b_0}{(k+1)!}\sqrt{\left\{k+1\right\}}\gamma _{k+1}$$
(39)
where $`k0.`$ Besides, we have the following consistency conditions (for $`p0,`$ $`k0`$)

$$\frac{\alpha _{2p+1,2p+k+1}}{\alpha _{2p+1,2p+k}}=\frac{\alpha _{2p1,2p+k+1}}{\alpha _{2p1,2p+k}}\frac{\alpha _{1,k+1}}{\alpha _{1,k}}$$
(40)
As a result we obtain the following theorem.

###### Theorem 4.1.

For the operator $`A`$ defined by (31) fullfil the conditions (34) and (35) it is necessary and sufficient to have the conditions (37) and (38). Here $`\alpha _{2m1,n1}`$ are the coefficients of the function $`\mathrm{\Psi }_n(x)`$ defined by (40). Under these conditions the elements $`a_{k,k+1}`$ of the matrix of the operator $`A`$ are determined by the relations (39).

###### Corollary 4.2.

Under the additional conditions

$$a_{0,1}0,0=a_{1,2}=a_{2,3}=\mathrm{}=a_{k,k+1}=\mathrm{},$$

we have

$$\sqrt{\left\{k\right\}}\gamma _k=k\gamma _1,k1.$$
(41)
###### Theorem 4.3.

The coefficients $`\alpha _{2m1,n1}`$ defined by (13) satisfy the conditions (40) iff there exist a positive sequence $`\left\{v_k\right\}_1^{\mathrm{}},`$ such that

$`\alpha _{2m1,n1}={\displaystyle \frac{\left\{2m1\right\}!!\left(v_1v_2\mathrm{}v_{n1}\right)}{\left(v_1v_2\mathrm{}v_{2m1}\right)\left(v_1v_2\mathrm{}v_{n2m1}\right)}},m1,n3`$
$`\alpha _{1,n1}=1,\alpha _{1,1}=\left\{1\right\}=1.`$ (42)
Under these conditions the sequence $`\left\{v_k\right\}_1^{\mathrm{}}`$ has the following properties

1. $`1=v_0v_1v_2\mathrm{}`$
2. $`v_{n2}v_{2p1}+v_{2p3}v_{n2p}=v_nv_{2p3}+v_{2p1}v_{n2p}`$ for all $`n2,`$ $`p1,`$ $`2pn`$ $`(v_10)`$.

Proof.

The sufficiency of the condition (42) is checked by the direct substitution (42) into (40).

Necessity. Let the coefficients $`\alpha _{2m1,n1}`$ determined by relations (13) are satisfied the conditions (40). Let us put $`v_0=a_{1,1}=1,v_11`$, and define the remaining terms of a sequence $`\left\{v_k\right\}_1^{\mathrm{}}`$ from the equalities

$$\alpha _{1,n}=\frac{v_{n1}v_n}{v_1},n1.$$
(43)
Using relations (40), we immediately check the validity of the formulas (42).

An example of a polynomial system of the above-mentioned type give us the generalized Hermite polynomials (entered in one of notes in the monography  and explicitly investigated in ). We consider these polynomials in the next section.

## 5 Generalized Hermite polynomials

The more complete exposition of the material considered in this section can be found in .

Let us consider the Hilbert space

$$𝙷_\gamma =L^2(R;|x|^\gamma (\mathrm{\Gamma }(\frac{1}{2}(\gamma +1)))^1\mathrm{exp}(x^2)dx),\gamma 1.$$
(44)
According to the method suggested in , we shall construct a canonical orthonormal polynomial system $`\left\{\psi _n(x)\right\}_{n=0}^{\mathrm{}}`$, which is complete in the space $`𝙷_\gamma `$. The polynomials $`\psi _n(x)`$ will satisfy the recurrent relations (19) and (20) with coefficients $`\left\{b_n\right\}_{n=0}^{\mathrm{}}`$, which are given by the following formulas:

$$b_{n1}=\frac{1}{2}\{\begin{array}{cc}\sqrt{n}& n=2m,\\ \sqrt{n+\gamma }& n=2m+1.\end{array}.$$
(45)
The polynomials $`\psi _n(x)`$ are defined by the relations (42), (18) and the formulas (12), where $`\{s\}={\displaystyle \frac{b_{s1}^2}{b_0^2}}`$. Besides, the sequence $`\left\{v_n\right\}_{n=0}^{\mathrm{}}`$ is given by the formulas

$$v_n=\{\begin{array}{cc}\frac{\gamma +n+1}{\gamma +1}& n=2m,\\ \frac{n+1}{\gamma +1}& n=2m+1.\end{array}.$$
(46)
It is clear that $`v_0=1,v_1={\displaystyle \frac{2}{\gamma +1}}=b_0^2.`$

The appropriate ”lowering” operator $`A`$ is defined by relations (31) and (35), where

$$a_{0,1}=1,a_{m1,m}=\frac{(2)^{m1}}{m!}\frac{\gamma }{\gamma +1},m2.$$
(47)
The sequence $`\left\{\gamma _n\right\}_{n=0}^{\mathrm{}}`$ included in (34) is given by the expressions:

$$\gamma _n=\frac{\sqrt{2}}{\gamma +1}\{\begin{array}{cc}\sqrt{n}& n=2m,\\ \sqrt{n+\gamma }& n=2m+1.\end{array}.$$
(48)
The annihilation operator $`a_\mu ^{}`$ of the generalized oscillator algebra is equal to

$$a_\mu ^{}=\frac{\gamma +1}{\sqrt{2}}A.$$
(49)
It is valid the following commutation relations for the position operator $`X_\mu `$ and the number operator $`N`$:

$$X_\mu \frac{d}{dx}N=(a_\mu ^{})^2.$$
(50)
Using the eigenvalue equation for a Hamiltonian $`H_\mu =a_\mu ^{}a_\mu ^++a_\mu ^+a_\mu ^{}`$, and taking into account the equality $`a_\mu ^+a_\mu ^{}=2B(N)`$, one can obtain the following relation :

$$a_\mu ^{}a_\mu ^+=2B(N+I).$$
(51)
From (50) – (51), it follows the differential equation

$$x\psi _n^{\prime \prime }+(\gamma 2x^2)\psi _n^{}+(2nx\frac{\theta _n}{x})\psi _n=0,n0,$$
(52)
where $`\theta _n={\displaystyle \frac{\gamma (1(1)^n)}{2}}.`$ This equation coincides with the known differential equation for generalized Hermite polynomials ( see also ).

###### Remark 5.1.

For the generators $`a_\mu ^+,a_\mu ^{}`$ of the Heisenberg algebras $`A_\mu `$, corresponding to the system of the generalized Hermite polynomials, we have (see ):

$$[a_\mu ^{},a_\mu ^+]=(\gamma +1)I2(2B(N)N).$$
(53)
”The energy levels” for Hamiltonian of the related oscillator are equal to:

$$\lambda _0=\gamma +1,\lambda _n=2n+\gamma +1,n1.$$
(54)
Finally, we have for the momentum operator the following formula:

$$P_\mu =ı(\frac{d}{dx}+X_\mu ^1\mathrm{\Theta }_NX_\mu ),$$
(55)
where

$$\mathrm{\Theta }_N=2B(N)N.$$
(56)
## 6 Conditions under which the annihilation operator $`A`$ can be realized by a differential operator

Now we consider the relations (32) and (33) without the supposition (35)relative to the coefficients $`a_{k,s}`$.

Let us introduce

$$C_{k+2m2w,k+2p+12w}=\frac{\left(k+2p+12m\right)!}{b_0^{k+2p+12m}}a_{k+2m2w,k+2p+12w}$$
(57)
Then the following assertion is hold.

###### Theorem 6.1.

The operator $`A`$ defined by the formula (31) satisfy the relation (29) iff the elements $`a_{k,s}`$ of this operator matrix satisfy the following conditions:

* For all $`t=2l,n=2p+2l+1,p0,l0`$ the following relation

$$\begin{array}{c}_{w=0}^{l1}_{m=0}^w\left(1\right)^m\frac{\alpha _{2m1,2p+2l}}{\left(2w2m\right)!}B_1+_{w=l}^{p+l}_{m=wl}^w\left(1\right)^m\frac{\alpha _{2m1,2p+2l}}{\left(2w2m\right)!}B_2=\\ \\ =(1)^p\frac{\left\{2l+2p+1\right\}!!}{b_0^{2l}},\end{array}$$
(58)
where

$$B_1=\left(C_{2l+2m2w,2p+2l+12w}+\frac{\left(2p+2l+12m\right)}{b_0}C_{2l+2m2w1,2p+2l2w}\right),$$

$$B_2=\left(C_{2l+2m2w,2p+2l+12w}+\frac{\left(2p+2l+12m\right)}{b_0}C_{2l+2m2w1,2p+2l2w}\right)$$

must be hold.
* For all $`t=2l+1,n=2p+2l+2,p0,l0`$ the following relation

$$\begin{array}{c}_{w=0}^{l1}_{m=0}^w\left(1\right)^m\frac{\alpha _{2m1,2p+2l+1}}{\left(2w2m\right)!}B_3+_{w=l}^{p+l}_{m=wl}^w\left(1\right)^mB_4\frac{\alpha _{2m1,2p+2l+1}}{\left(2w2m\right)!}=\\ \\ (1)^p\left\{2l+2p+2\right\}\frac{\alpha _{2p1,2p+2l}}{b_0^{2l+1}}\end{array}$$
(59)
should be fulfilled, where

$$B_3=\left(C_{2l+2m2w+1,2p+2l+22w}+\frac{\left(2p+2l+22m\right)}{\left(2w+12m\right)b_0}C_{2l+2m2w,2p+2l+12w}\right),$$

$$B_4=\left(C_{2l+2m2w+1,2p+2l+22w}+\frac{\left(2p+2l+22m\right)}{\left(2w+12m\right)b_0}C_{2l+2m2w,2p+2l+12w}\right),$$

and the coefficients $`\alpha _{2mp1,n1}`$ are given by (13).

## 7 Conclusion

From the results given above it follows that in the general case a orthonormal polynomial system satisfies differential (or difference) equation of a finite order. Therefore it is interesting to describe such orthonormal polynomial system, for which one can obtain a differential (or difference) equation of a finite order. It is desirable that we have to deal only with differential equation of the second order. The elementary example of such class of polynomials selected by the condition (35), is given by Hermite - Chihara polynomials  considered in detail in .

Unfortunately, the relations (5759) are rather complicated in the general case. Therefore generally local formulas representing an annihilation operator (and also for a momentum operator and Hamiltonian) by a differential operator are too cumbersome.

Finally, we note the that the formalism put forward in the present work holds true also after replacement usual derivation by some of its $`q`$-analog. The last circumstance provides a useful guide to investigating of deformed polynomials.