text:
# An iterative destriping technique for diffuse background polarization data

## 1 introduction

Low frequency noise is known to affect all radiometers and to induce correlations among successive samples of the measured signal, leading to typical striping effects when producing sky maps. It is characterised by a power-law spectrum $`S(f)(1/f)^\beta `$, with $`\beta `$ roughly in the range $`1÷2.5`$ depending on the noise source. The total instrumental noise power spectrum is usually specified in terms of the knee frequency, $`f_\mathrm{k}`$, at which the white and low frequency components of the noise are equal:

$$S(f)=\sigma _0^2\left[1+\left(\frac{f_\mathrm{k}}{f}\right)^\beta \right]$$
(1)
where $`\sigma _0^2`$ represents the white noise variance per second. As shown by Janssen et al. (1996), when data are taken from spinning spacecrafts, where the spin acts as an ideal switch or chopper, most of the low frequency noise can be removed by software, provided $`f_\mathrm{k}`$ is lower than the satellite spin frequency, $`f_\mathrm{s}`$.
For present and future space experiments primarily dedicated to Cosmic Microwave Background (CMB) anysotropy measurements, like MAP111MAP homepage: http://map.gsfc.nasa.gov and PLANCK222PLANCK homepage: http://astro.estec.esa.nl/Planck, destriping techniques have already been set up to clean the Time Ordered Data (TOD) from low frequency noise. The impact of different scanning strategies on the quality of final maps has been studied as well (Delabrouille 1998; Maino et al. 1999).
The problem of removing $`1/f`$ noise from the data stream is even more important for experiments attempting measurements of CMB polarization (CMBP), where even a tiny amount of residual correlated noise, at a level negligible for investigations of CMB temperature anysotropy ($`\mathrm{\Delta }T/T10^5`$), might distort the signal characteristics: in the most favourable case the expected polarization level is only in the order of 10% of the temperature fluctuations. The first algorithm specifically designed for destriping polarization data has been proposed by Revenu et al. (2000), who extended to polarization data ideas previously studied to destripe PLANCK anisotropy data (Delabrouille 1998).

The SPOrt333SPOrt homepage: http://sport.bo.iasf.cnr.it experiment onboard the ISS (Cortiglioni et al. 2002; Carretti et al. 2003) is equipped with 4 correlation polarimeters directly measuring the $`Q`$ and $`U`$ Stokes parameters in the 22-90 GHz frequency range, with $`\mathrm{FWHM}=7^{}`$. Besides providing polarization maps of Galactic synchrotron, SPOrt will try a first detection of CMBP on large angular scales ($`\mathrm{}25`$), where the CMB E-mode power spectrum is particularly sensitive to the value of the optical depth of the reionized medium (Seljak 1997; Zaldarriaga et al. 1997). Such information cannot be extracted from CMB temperature anysotropy data, nor from CMBP data taken by balloon and ground-based experiments looking at small sky patches. This makes CMBP measurements on large angular scale of utmost importance.

The algorithm presented here is particularly suited to SPOrt, acting directely on $`Q`$ and $`U`$ Stokes parameter data.

The SPOrt radiometers are based on InP HEMT amplifiers, for which the low frequency noise is known to have a $`1/f`$-like spectrum ($`\beta 1`$) dominated by transistor gain fluctuations (Wollack 1995), with knee frequencies in the range 100-1000 Hz. However, when using correlation techniques, the knee frequency is reduced by a factor $`(T_{\mathrm{offset}}/T_{\mathrm{sys}})^2`$, with $`T_{\mathrm{offset}}`$ and $`T_{\mathrm{sys}}`$ the instrumental offset and the system noise temperature, respectively (Wollack & Pospiesalski 1998; Carretti et al. 2001). The present state of the instrument already guarantees a value $`f_\mathrm{k}<2.5\times 10^5`$ Hz (Carretti et al. 2003), close to the goal knee frequency of one tenth of the ISS orbit frequency, i.e. $`f_\mathrm{k}^{\mathrm{goal}}=1.8\times 10^5`$ Hz.
The SPOrt scanning strategy is bound to the ISS motion. The ISS orbit is tilted by $`51.6^{}`$ with respect to the Celestial equator and is characterized by a frequency $`f_\mathrm{o}=1.8\times 10^4`$ Hz. The SPOrt antennae, looking at the zenit, will cover about 80% of the sky under precessing circles intersecting each other. The nominal sampling rate is 1 Hz, and a sky map will be provided every 70 days in a mission lasting at least 18 months. Further details about the SPOrt sky coverage are shown in Fig. 1.
The iterative algorithm described in this paper has been specifically designed to destripe $`Q`$ and $`U`$ Stokes parameter data. It is general enough to be valid for any experiment where data are taken by scanning the sky at a frequency higher than the receiver knee frequency, and where different scans intersect each other in a sufficient number of points. Measurements of scalar quantities are easily handled.
The description of the algorithm is provided in Sect. 2, whereas Sect. 3 presents various tests of its performances and Sect. 4 summarises our conclusions.

## 2 Destriping method

The only assumptions our algorithm relies upon is that the radiometer is stable during the signal modulation period (the time needed to complete one orbit in case of SPOrt), so that the instrumental offset does not change significantly during this time, and there is enough overlap between different orbits. In such a case the noise can be split in two parts: for time scales shorter than the orbit period (high frequencies) it is essentially white, whereas for longer timescales (low frequencies) it also contains the $`1/f`$ component, which we model as a different constant offset for each orbit. A simple iterative procedure is then applied to remove these offsets from the TOD before map-making. The map-making itself consists in averaging the cleaned TOD corresponding to the same sky pixel after rotating them into a fixed reference frame. Parallel transport to the center of each pixel is also implemented, as described by Bruscoli et al. (2002), before averaging, and can optionally be switched on, though it is not strictly necessary for the portion of sky covered by SPOrt. It would be mandatory if polar caps were included in the accessible region. The time needed to simulate 1 year of SPOrt data taking with the nominal sampling rate of 1 Hz, destripe the TOD and produce $`Q`$ and $`U`$ clean maps is about $`10`$ minutes on a 1.3 GHz Pentium III equipped with a 4 GB RAM when using 10 destriping cycles, and increases by 2.5 minutes for each additional bunch of 10 iteration cycles.

### 2.1 Formalism

Each datum consists of a ($`\stackrel{~}{Q}`$,$`\stackrel{~}{U}`$) pair as measured in the polarimeter reference frame: the same sky point, when seen from different orbits, shows different ($`\stackrel{~}{Q}`$,$`\stackrel{~}{U}`$) values even in absence of noise. The standard longitude-latitude fixed reference frame (Berkhuijsen 1975) is used when considering the real sky emission, which is projected into an $`N`$-pixel map. Our $`M`$ TOD, collected during $`R`$ orbits about the Earth, are used to build an $`M`$-vector $`Y`$ of measured ($`\stackrel{~}{Q}`$,$`\stackrel{~}{U}`$) pairs:

$$Y=AX+BOFF+N.$$
(2)
$`X`$ is an $`N`$-vector whose elements are the ($`Q`$,$`U`$) pairs of the true sky emission; $`A`$ is an $`M\times N`$ pointing matrix connecting sky pixels to observations (Wright 1996; Tegmark 1997a): its elements are $`2\times 2`$ matrices rotating the Stokes parameters from the fixed to the polarimeter reference frame; $`OFF`$ is an $`R`$-vector containing a pair of radiometric offsets for each orbit; $`B`$ is an $`M\times R`$ pointing matrix connecting orbits to observations, and $`N`$ is an $`M`$-vector containing pairs of white noise. In more details, the elements of the above mentioned vectors and matrices are:

$`Y_t=\left(\begin{array}{c}\stackrel{~}{Q}_t\\ \stackrel{~}{U}_t\end{array}\right)`$ (5)
$$A_{tp}=\{\begin{array}{cc}(\alpha _t)=\left(\begin{array}{cc}\mathrm{cos}(2\alpha _t)& \mathrm{sin}(2\alpha _t)\\ \mathrm{sin}(2\alpha _t)& \mathrm{cos}(2\alpha _t)\end{array}\right)\hfill & \begin{array}{c}\text{if pixel }p\hfill \\ \text{is observed}\hfill \\ \text{at time }t\hfill \end{array}\hfill \\ & \\ \left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right)\hfill & \mathrm{otherwise}\hfill \end{array}$$
(6)
where $`\alpha _t`$ is the angle between the polarimeter and the fixed reference frames at time $`t`$.

$`X_p=\left(\begin{array}{c}Q_p\\ U_p\end{array}\right)`$ (9)
$$B_{to}=\{\begin{array}{cc}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\hfill & \text{if }o\text{ is the orbit run at time }t\hfill \\ & \\ \left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right)\hfill & \mathrm{otherwise}\hfill \end{array}$$
(10)
$`Off_o=\left(\begin{array}{c}Off_o^{\stackrel{~}{Q}}\\ Off_o^{\stackrel{~}{U}}\end{array}\right)`$ (13)
$`N_t=\left(\begin{array}{c}N_t^{\stackrel{~}{Q}}\\ N_t^{\stackrel{~}{U}}\end{array}\right)`$ (16)
### 2.2 Iteration cycle

If we knew the sky emission vector $`X`$ we could use it to obtain the offset vector $`OFF`$ by inverting Eq. (2). Of course, estimating the noise offsets would not be worthwhile if we already knew the true sky emission. However, to start we can guess a 0-order sky emission $`𝐗^{(\mathrm{𝟎})}`$, the most neutral choice being a null vector; then obtain a 0-order offset estimate, and subtract it from our measurements to get a 0-order cleaned TOD; from this compute $`𝐗^{(\mathrm{𝟏})}`$, the 1st-order estimate of the real sky emission, and iterate the cycle.
Starting from Eq. (2) and following Tegmark (1997b), the best estimators for the offset vector and the signal matrix can be written:

$`OFF`$ $`=`$ $`W_{\mathrm{obs}}(YAX)`$ (17)
$`X`$ $`=`$ $`W_{\mathrm{pix}}(YBOFF)`$ (18)
where, in case of uniform white-noise variance, as expected for the SPOrt receivers,

$`W_{\mathrm{obs}}`$ $`=`$ $`N_{\mathrm{obs}}^1B^\mathrm{T}`$ (19)
$`W_{\mathrm{pix}}`$ $`=`$ $`N_{\mathrm{pix}}^1A^\mathrm{T}`$ (20)
with $`N_{\mathrm{obs}}=B^\mathrm{T}B`$ and $`N_{\mathrm{pix}}=A^\mathrm{T}A`$ diagonal matrices containing the number of observations in each orbit and in each pixel, respectively. In such a case Eq (10) simply corresponds to rotating into the fixed reference frame the offset-subtracted observations and averaging those corresponding to the same sky pixel. When starting from a null guess map, i.e. $`X^{(0)}=0`$, the first offset estimate is a simple average of the measurements taken over each orbit.

In the more general case of non-uniform white noise, i.e. when the noise covariance matrix $`C=NN^\mathrm{T}`$ has unequal diagonal elements, the matrices $`W_{\mathrm{obs}}`$ and $`W_{\mathrm{pix}}`$ become (Tegmark 1997b)

$`W_{\mathrm{obs}}`$ $`=`$ $`(B^\mathrm{T}C^1B)^1B^\mathrm{T}C^1`$ (21)
$`W_{\mathrm{pix}}`$ $`=`$ $`(A^\mathrm{T}C^1A)^1A^\mathrm{T}C^1`$ (22)
meaning that the simple averages mentioned above should be replaced with noise-weighted averages.

In any case the first offset estimate writes:

$`OFF^{(0)}`$ $`=`$ $`W_{\mathrm{obs}}^1Y`$ (23)
The 0-order cleaned TOD are:

$$Y^{(0)}=YBOFF^{(0)}$$
(24)
and the 1st-order estimate of the true sky emission is

$$X^{(1)}=W_{\mathrm{pix}}Y^{(0)}$$
(25)
The $`i^{th}`$ cycle writes:

$`OFF^{(i)}`$ $`=`$ $`W_{\mathrm{obs}}(YAX^{(i)})`$ (26)
$`X^{(i+1)}`$ $`=`$ $`W_{\mathrm{pix}}(YBOFF^{(i)})`$ (27)
and must be iterated until the difference between the sky emission estimates of order $`i`$ and $`i+1`$ is negligible when compared to the white noise level.

The proof that the procedure is effective has been obtained numerically, as described in detail in the next section. Destriping measurements of total intensity (or, in general, any other scalar quantity) can be simply performed by replacing pairs with scalars in the previously defined vectors, and rotation matrices with scalar 1. The average sky signal is lost when destriping maps of scalar quantities. However, for maps of polarization data like $`Q`$ and $`U`$ it can be kept provided the polarimeter reference frame rotates about the fixed frame while running along each orbit, as in the case of SPOrt. This is a nice feature, especially when measuring foreground contributions.

## 3 Simulations and tests

We test the method with numerical simulations of the data stream expected from about one year of SPOrt data taking.
As real sky emission we consider both CMB and synchrotron radiation, both convolved with a $`7^{}`$ Gaussian beam. Other foregrounds are expected to be less important at the frequencies and angular scales covered by SPOrt (Tegmark et al. 2000; Lazarian & Prunet 2002) and are not included.
CMB emission is generated by the CMBFAST444http://www.physics.nyu.edu/matiasz/CMBFAST/cmbfast.html package according to a $`\mathrm{\Lambda }`$CDM cosmological model with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\tau =0.2`$, $`H_0=65\mathrm{Kms}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }_b=0.05`$ and no contributions from gravitational waves. For the synchrotron radiation we use $`Q`$ and $`U`$ map templates at 22 GHz developed by Bernardi et al. (2003) and featuring a $`P_{\mathrm{rms}}=\sqrt{Q^2+U^2}35\mu `$K. The polarized intensity peak emission is about 130 $`\mu `$K. The main ingredients of the template are total intensity low frequency data (Haslam et al. 1982; Reich 1982), from which the syncrothron polarized intensity is extracted, and optical starlight data (Heiles 2000), from which polarization angles are derived. The template frequency has been chosen to correspond to the SPOrt channel where synchrotron emission is expected to be most important.

Noise streams to be added to the sky signal are built in the frequency space, where their Fourier coefficients are generated according to the spectrum defined in Eq. (1), with $`\beta =1`$, and transported to the time domain via FFT. Our code allows us to calculate the angle between our polarimeters and the longitude-latitude fixed reference frame for each data sample, and to set the sampling frequency from the nominal rate (1 Hz) to any other rate while keeping the SPOrt nominal istantaneous sensitivity of $`1\mathrm{mK}\mathrm{s}^{1/2}`$ (Carretti et al. 2003). In the tests described here the sampling rate was set to 1 Hz, and the SPOrt observing time was about 1 year. The final output of our simulations are TOD and noise arrays, clean $`Q`$ and $`U`$ signal maps to investigate the effects of SPOrt partial sky survey, $`Q`$ and $`U`$ maps containing both signal and noise before and after destriping as well as a hit map to evaluate the time spent over each map pixel. We use the HEALPix555http://www.eso.org/science/healpix pixelisation scheme (Gorski et al. 1998) with $`Nside=64`$ to ensure the effects of pixelisation are negligible when compared to the antenna beam (the pixel size of the output maps is $`1^{}`$).
We test the destriping efficiency on simulated maps by evaluating the residual low frequency noise in terms of added rms noise with respect to the white noise case, angular correlation functions $`C^X(\theta )`$ ($`X=Q,U`$) and power spectra $`C_{\mathrm{}}^Y`$ ($`Y=E,B`$).

### 3.1 Correlation functions and power spectra

To avoid edge effects due to partial sky coverage and to efficiently take into account the presence of noise, even in case of non uniform sky coverage, we primarly test our destriping technique by studying the two-point correlation functions $`C^X(\theta )`$ ($`X=Q,U`$) measured from our simulated maps. We sample the correlation functions with about one degree of angular resolution according to the expression:

$$\stackrel{~}{C}^X(\theta )=\frac{_{ij}w_{ij}\mathrm{\Delta }_i^X\mathrm{\Delta }_j^X}{_{ij}w_{ij}}.$$
(28)
Here $`\mathrm{\Delta }_i^X`$ is the content of pixel $`i`$ of map $`X`$, and $`w_{ij}`$ is the weight of the $`ij`$ pixel pair

$$w_{ij}=\frac{2w_iw_j}{w_i+w_j}$$
(29)
where each pixel weight $`w_i`$ is proportional to the pixel observing time and $`_{i=1}^{N_{pix}}w_i=N_{pix}`$.

In the computation the Stokes parameters are properly rotated so that $`Q`$ is defined along the great circle connecting the points in each pair (Zaldarriaga 1998). The exact value of the angle $`\theta `$ at which the correlation functions are sampled is calculated as the weighted average of the angles between pixels contributing to the same measured point. As an example, the measured two-point correlation functions of $`Q`$,$`U`$ maps containing only white noise are on average zero everywhere but at 0-distance, where they are checked to be equivalent to the average pixel variance.
The $`C_{\mathrm{}}^{E,B}`$ power spectra can be recovered from the measured correlation functions by integration (Kamionkowski et al. 1997; Zaldarriaga 1998):

$$C_{\mathrm{}}^E=_0^\pi [C^Q(\theta )F_{\mathrm{},2}^1(\theta )+C^U(\theta )F_{\mathrm{},2}^2(\theta )]\mathrm{sin}\theta d\theta $$
(30)
$$C_{\mathrm{}}^B=_0^\pi [C^U(\theta )F_{\mathrm{},2}^1(\theta )+C^Q(\theta )F_{\mathrm{},2}^2(\theta )]\mathrm{sin}\theta d\theta $$
(31)
where the functions $`F_{\mathrm{},2}^i(\theta )`$ have been defined by Zaldarriaga (1998). We perform the integration by Gauss quadrature, where the correlation function values at the needed points $`\theta `$ are obtained by interpolation, with a second order polynomial, of the three closest measured points. This method of measuring power spectra is similar to that used by Szapudi et al. (2001), though here it has been extended to polarization data. It has been tested on full-sky maps of simulated CMB emission, convolved with a Gaussian filter with FWHM=7, by comparing the measured power spectra to those obtained by using the ANAFAST package: we get the same results, within numerical precision, up to $`\mathrm{}=25`$, i.e. in the full multipole range covered by SPOrt. On partial sky maps this method is found to work better than ANAFAST, thanks to the lack of edge effects: our average of power spectra extracted from CMB maps without noise and featuring the SPOrt sky coverage matches the input theoretical spectrum very well, whereas that measured with ANAFAST is flatter, as shown in Fig. 2. Furthermore, when trying to extract the real sky signal from noisy maps this method allows for pixel weighting in case of non-uniform sky coverage as well as noise subtraction, provided the noise statistical properties (correlation functions) are known. It has the additional advantages of being simple and rather fast: the required computing time practically corresponds to that needed to calculate the two-point correlation functions, i.e. about 35 minutes on our platform. The number of operations is $`O(N_{\mathrm{pix}}^2)`$, while optimal methods to measure power spectra require matrix inversions tipically implying a number of operations $`O(N_{\mathrm{pix}}^3)`$ (Wright 1996). The implementation of a new algorithm for computing N-point correlation functions is expected to further reduce the computing time (Moore et al. 2000).

### 3.2 Destriping signal-only maps

To check for the presence of spurious correlations introduced by our algorithm and depending on the underlying signal, as a first test we run the destriping code over CMB and Synchrotron maps containing no noise. The test confirms our algorithm subtracts the average signal from the input map when dealing with scalar quantities (CMB temperature in our case), whereas nothing is lost in our $`Q`$ and $`U`$ maps. We then analyse the difference between the destriped and the original maps as a function of the number of iterations of the destriping cycle. Should no residual effects be present, these differences would asimptotically become null for $`Q`$ and $`U`$ maps, and a constant number corresponding to the average of the input signal for scalar maps. Indeed, the maximum pixel-to-pixel difference, which depends on the level of the underlying signal, decreases with increasing number of iterations. In particular, the spurious noise introduced when destriping maps of synchrotron emission at 22 GHz, i.e. the SPOrt frequency where the Galactic contribution is expected to be higher, becomes negligible also for CMBP after about 20 loops, reaching the tiny level of few nK after 30 loops. More details are shown in Table 1, both in terms of peak-to-peak amplitude and rms of spurious noise. For comparison, the peak-to-peak amplitude and rms of the underlying signal maps are shown as well.

The same conclusions can be derived by inspecting the angular correlation functions $`C^X(\theta )`$ ($`X=Q,U`$) of the spurious noise: they are not flat, but deviations from a flat null function decrease with increasing number of iteration cycles and can be made negligible after about 10 (20) loops when destriping CMB (synchrotron) $`Q`$, $`U`$ maps.

### 3.3 Destriping noise maps

The effectiveness of our destriping technique is then tested on simulated maps containing only noise, in the following cases:

* white noise
* $`1/f`$ noise
* both white and $`1/f`$ noise

The good performances of the technique are made evident in Fig. 3 where we show a simulated noise map (both $`1/f`$ and white) before and after destriping.

One possible way to quantify the residual correlated noise after destriping is measuring the fractional excess pixel noise with respect to the case of purely white noise. Results are shown in Table 2 for two different values of the knee frequency, corresponding to the SPOrt goal knee frequency, $`f_\mathrm{k}=1.8\times 10^5`$ Hz, and the SPOrt orbit frequency, $`f_\mathrm{o}=1.8\times 10^4`$ Hz, the latter representing a very conservative case.

Another way to quantify the residual correlated noise is measuring and inspecting the two-point correlation functions C$`{}_{}{}^{X}(\theta )`$ of simulated $`Q`$ and $`U`$ noise maps. Averages and 1$`\sigma `$ bands of 500 correlation functions C$`{}_{}{}^{Q}(\theta )`$ measured from maps containing both white and $`1/f`$ noise, before and after destriping, are compared to the purely white noise case in Fig. 4 for the same knee frequencies as in the previous test. The correlation functions C$`{}_{}{}^{U}(\theta )`$ are similar and are not shown. As expected, the correlated noise is strongly reduced by the destriping procedure, the residuals falling within the statistical error of the white noise case for $`f_\mathrm{k}=f_\mathrm{k}^{\mathrm{goal}}`$. The correlation functions for the case of purely $`1/f`$ noise are similar to the white plus $`1/f`$ case everywhere but at 0 distance, and are not shown.

If the noise statistical properties are known, its two-point correlation function after destriping can be calculated with Monte Carlo techniques and subtracted from the measured values before integrating, the recovering of signal power spectra presenting no huge problems. However, the residual correlated noise also increases the expected error on measured quantities. In our case the rms of correlation functions measured from many noise maps is increased by roughly 4% (11%), for $`f_\mathrm{k}=1.8\times 10^5`$ Hz ($`f_\mathrm{k}=1.8\times 10^4`$ Hz), with respect to the case of purely white noise.

For completeness, and to ease the comparison with other methods, average noise power spectra, obtained from the correlation functions used to make Fig. 4, are shown in Fig. 5, though they carry the same information as the correlation functions. As expected, even for the largest knee frequency the $`C_{\mathrm{}}^E`$ power spectrum of the residual noise is close to that of purely white noise. The region of low multipoles is the most sensitive to low frequency residuals, some contributions being always found here also after the application of other destriping techniques (Maino et al. 1999). The excess noise, in percentage of the measured white noise level, is shown in the bottom panel of Fig. 5 as a function of the multipole. The theoretical white noise level can be calculated as follows:

$$C_{\mathrm{}}^{wn}=\frac{4\pi }{N_{pix}}\sigma ^2$$
(32)
where $`\sigma ^2`$ represents the average noise variance per pixel. Our measured value is equivalent to its expectation of $`0.37\mu `$K2, again confirming that our map-making algorithm is correct.

Finally, Fig. 6 shows the increment, due to the presence of residual correlated noise, in the rms of measured noise power spectra, after destriping, in percentage of the rms of purely white noise power spectra, again as a function of the multipole.

## 4 Conclusions

We presented a new iterative destriping technique working directely on $`Q`$ and $`U`$ Stokes parameter data as well as on scalar quantities. As for other methods, in order for the technique to work properly the receiver knee frequency must be lower than the experiment modulation frequency, and there must be sufficient overlap between different orbits. On the other hand, no requirements on the statistical properties of the noise are necessary here. Moreover, our observing strategy preserves average $`Q`$ and $`U`$ values in the destriped maps, though the map average is still lost when dealing with scalar quantities.

The performances of the technique, studied on simulated SPOrt data by analysing both the measured $`C^{Q,U}(\theta )`$ two-point correlation functions and the measured $`C_{\mathrm{}}^{E,B}`$ power spectra, are comparable to those of other methods based on $`\chi ^2`$ minimisation. Computer time is not an issue, no matrix inversions being involved.

Power spectra are recovered from the measured two-point correlation functions by integration, this method being implemented here for the first time for polarization data. It has the advantages of avoiding edge problems arising when using the ANAFAST code and to allow for pixel weighting in case of non uniform sky coverage.

###### Acknowledgements.

This work has been carried out in the context of the SPOrt program, which is funded by the Italian Space Agency. We thank B. Audone and F. Amisano for useful discussions and the anonimous referee for good suggestions. We acknowledge the use of CMBFAST and HEALPix packages.


format: text

text:
# Static potential in baryon

(Institute of Theoretical and Experimental Physics,
117218, B.Cheremushkinskaya 25, Moscow, Russia)

## Abstract

The baryon static potential is calculated in the framework of field correlator method and is shown to match the recent lattice results. The effects of the nonzero value of the gluon correlation length are emphasized.

1. The static potential is a key quantity in the calculations of spectra and wave functions of baryons in the potential models , as well as in the QCD string approach . On the other hand the static potential is generated by the gluon forces of QCD and therefore sheds the light on the profound properties of strong interactions. Accurate numerical calculations of static potential in baryon performed in lattice QCD last years ,  induced new significant growth of interest to the structure of the gluon forces in baryon.

The results of calculations of the baryon potential in the framework of the method of the field correlators (MFC)  are presented in this report. We set the correspondence of this potential to lattice one and underline fundamental effects of strong interactions — the confinement of color charges and definite value of the correlation length of the gluon fields, which reveal themselves in the behavior of baryon potential.

2. The static potential is expressed through the Wilson loop $`𝒲_B`$ according to the relation

$$V_B=\underset{T\mathrm{}}{lim}\frac{1}{T}\mathrm{ln}𝒲_B,$$
(1)
and that is why we are proceeding now to the calculations of Wilson loop. The baryon Wilson loop is shown in Fig. 1 and is defined as follows,

$$𝒲_B=\frac{1}{6}ϵ_{\alpha \beta \gamma }ϵ^{\alpha ^{}\beta ^{}\gamma ^{}}\mathrm{\Phi }_\alpha ^{}^\alpha (C_1)\mathrm{\Phi }_\beta ^{}^\beta (C_2)\mathrm{\Phi }_\gamma ^{}^\gamma (C_3),$$
(2)
where the parallel transporter or Schwinger line $`\mathrm{\Phi }`$ is given as

$$\mathrm{\Phi }_\beta ^\alpha (x,y,C)=(P\mathrm{exp}ig_CA_\mu 𝑑z_\mu )_\beta ^\alpha .$$
(3)
The average of the Wilson loop over the vacuum fields using the nonabelian Stokes theorem and bilocal (two-point) approximation is performed in , and the result reads as

$$𝒲_B=\mathrm{exp}\{\underset{i=1}{\overset{3}{}}\frac{1}{2}_{S_i}_{S_i}d\sigma _{\mu _1\nu _1}(x)d\sigma _{\mu _2\nu _2}(x^{})𝒟_{\mu _1\nu _1,\mu _2\nu _2}(xx^{})+$$

$$\underset{i<j}{}\frac{1}{2}_{S_i}_{S_j}d\sigma _{\mu _1\nu _1}(x)d\sigma _{\mu _2\nu _2}(x^{})𝒟_{\mu _1\nu _1,\mu _2\nu _2}(xx^{})\},$$
(4)
where $`𝒟`$ designates the bilocal field strength correlators,

$$𝒟_{\mu _1\nu _1,\mu _2\nu _2}(xx^{})\frac{g^2}{N_c}\mathrm{tr}F_{\mu _1\nu _1}(x)\mathrm{\Phi }(x,x^{})F_{\mu _2\nu _2}(x^{})\mathrm{\Phi }(x^{},x),$$
(5)
and it is assumed that the straight-line trajectories for parallel transporters are chosen. Integrations in (4) are taken over the surfaces $`S_i`$ of the Wilson loop, which are formed by the trajectories of the corresponding quark, $`C_i`$, and that of the string junction, see Fig. 1. The trajectory of the string junction is shown in the figure by the dotted line.

Relying on the phenomena of the Casimir scaling , which is confirmed in lattice simulations, one can expect that the bilocal approximation (4) is valid within the accuracy of 1%.

Performing the surface integration in (4) and substituting the result into (1) we arrive at the expression for the baryon potential in bilocal approximation,

$$V_B(R_1,R_2,R_3)=\left(\underset{a=b}{}\underset{a<b}{}\right)n_i^{(a)}n_j^{(b)}_0^{R_a}_0^{R_b}𝑑l𝑑l^{}_0^{\mathrm{}}𝑑t𝒟_{i4,j4}(z_{ab}),$$
(6)
where $`R_a`$ is separation of the string junction and corresponding ($`a`$-th) quark, $`𝐧^{(a)}`$ the unity vector directed from the string junction to this quark, and $`z_{ab}=(l𝐧^{(a)}l^{}𝐧^{(b)},t)`$.

According to MFC , the bilocal correlators are written in the general form containing two scalar formfactors $`D(z)`$ and $`D_1(z)`$,

$$𝒟_{\mu _1\nu _1,\mu _2\nu _2}(z)=(\delta _{\mu _1\mu _2}\delta _{\nu _1\nu _2}\delta _{\mu _1\nu _2}\delta _{\mu _2\nu _1})D(z)+$$

$$\frac{1}{2}\left(\frac{}{z_{\mu _1}}(z_{\mu _2}\delta _{\nu _1\nu _2}z_{\nu _2}\delta _{\nu _1\mu _2})+\frac{}{z_{\nu _1}}(z_{\nu _2}\delta _{\mu _1\mu _2}z_{\mu _2}\delta _{\mu _1\nu _2})\right)D_1(z).$$
(7)
The formfactor $`D(z)`$ decreases exponentially,

$$D(z)=D(0)\mathrm{exp}\left(\frac{|z|}{T_g}\right),$$
(8)
which reflects the stochastic properties of the nonperturbative background gluon fields and is justified by the lattice computations . This behavior leads to the asymptotic area law for the Wilson loop, where the string tension $`\sigma `$ is expressed through $`D(z)`$ as follows,

$$\sigma =\frac{\pi }{2}_0^{\mathrm{}}𝑑z^2D(z)=\pi D(0)T_g^2=0.18\text{ GeV}^2.$$
(9)
The string tension is the main nonperturbative parameter of MFC. Its value is defined phenomenologically by the slope of the meson Regge trajectories  and is directly related to the radius of confinement.

There is another parameter in (8), the correlation length of the background gluon field $`T_g`$. However, it is not an independent parameter. Its value is extracted from the energy of the gluon excitation of the hadron spectra  and may be calculated in the QCD string approach ,  using the only parameter $`\sigma `$. The energy of the gluon excitation is large, $`1.5`$ GeV , and therefore the gluon correlation scale is significantly less than the confinement scale.

The formfactor $`D_1(z)`$ is dominated by the contribution of one-gluon-exchange of the perturbation theory, which gives rise to the static color-Coulomb potential. In the case of the baryon with the quarks forming equilateral triangle with the side $`r`$ the perturbative potential reads as

$$V_{\mathrm{pert}}(r)=\frac{3}{2}\frac{C_F\alpha _s}{r},$$
(10)
where $`C_F=4/3`$.

3. We now proceed to the comparison of the potential (6)-(10) with the lattice results , . In Fig. 2 the results of the lattice simulations  of the baryon potential with the perturbative part subtracted is shown by points in dependence on the length of the string, $`L=R_1+R_2+R_3`$. Solid line in the figure shows the behavior of the nonperturbative potential calculated in MFC for the configurations of equilateral triangle. When $`L\begin{array}{c}>\hfill \\ \hfill \end{array}1`$ fm, the potential grows linearly having the slope $`\sigma `$. When the length of the string gets smaller, the slope of the potential diminishes. Dotted line presents a tangent to the MFC potential at $`L=0.7`$ fm. The slope of the tangent is $`0.9\sigma `$. Linear phenomenological potential with the same slope was used in the constituent model for the description of the spectra of baryons long time ago . One can see from the figure that our curve goes well through all lattice points. The configurations of quarks forming triangles with the angles in the region the region $`[\pi /20,\pi /2]`$, which were used in lattice work , do not allow to establish any dependence of the potential on the configuration with the accuracy given. The study of the MFC potential on quark configurations will be performed below (see Fig. 4).

In Fig. 3 the results of lattice computations of the potential in equilateral triangle  are presented (points) vs. the quark separations $`r`$. The MFC potential with the perturbative part included is shown by solid line. One can see that all the lattice data are completely described by the MFC potential.

Computations in ,  are performed in quenched approximation. In the MFC calculations sea quarks are not considered too. The studies of the effects of light dynamical quarks on the $`Q\overline{Q}`$ static potential in lattice were recently performed in . No clear evidence for a flattening of the potential was found up to distances 2.5 fm. One can therefore expect that the sea quarks would not change the baryon potential significantly at typical hadron sizes.

It is also interesting to study the dependence of the potential in baryon on the quark locations at fixed length of the string $`L`$. Let us consider isoceles triangles with the full length of the string $`L`$ and vertex $`\gamma `$ and denote $`V^L(\gamma )`$ the nonperturbative potential in these triangles. At large enough sizes $`L\begin{array}{c}>\hfill \\ \hfill \end{array}1`$ fm$`T_g`$ asymptotic relations follow from (6)-(9). When $`\gamma =0`$ and locations of two quarks coinside, these quarks combine in antitriplet. The string consists of one line and

$$V_1V^L(\gamma =0)=\sigma L\frac{4}{\pi }\sigma T_g.$$
(11)
When $`0<\gamma <2\pi /3`$, the string consists of three lines, and potential

$$V_2V^L(0<\gamma <2\pi /3)=\sigma L+\left(\frac{12}{\pi }+\frac{2}{\sqrt{3}}\right)\sigma T_g.$$
(12)
When $`\gamma 2\pi /3`$, the string consists of two lines. In this case

$$V^L(\gamma 2\pi /3)=\sigma L+\left(\frac{8}{\pi }+\frac{2}{\pi }(\pi \gamma )\mathrm{ctg}(\pi \gamma )\right)\sigma T_g.$$
(13)
It is not difficult to calculate difference between the potentials in different configurations (we use $`\sigma `$=0.18 GeV2, $`T_g`$=0.12 fm),

$$\mathrm{\Delta }V_1V_1V_2=\left(\frac{8}{\pi }\frac{2}{\sqrt{3}}\right)\sigma T_g150\text{ MeV},$$
(14)
$$\mathrm{\Delta }V_2V^L\left(\frac{2\pi }{3}\right)V_2=\left(\frac{4}{\pi }\frac{4}{3\sqrt{3}}\right)\sigma T_g55\text{ MeV},$$
(15)
$$\mathrm{\Delta }V_3V^L(\gamma \pi )V^L\left(\frac{2\pi }{3}\right)=\left(\frac{2}{\pi }\frac{2}{3\sqrt{3}}\right)\sigma T_g30\text{ MeV}.$$
(16)
In Fig. 4 the dependence of the baryon potential in isoceles triangle on the angle $`\gamma `$ at $`L=1.8`$ fm is shown. One can verify that for the curve shown in the figure the relations (13)-(16) are justified (taking $`\gamma =\pi /3`$ for $`V_2`$).

The nonperturbative potential in the configurations under consideration has a pike near $`\gamma `$=0. However, in this region the perturbative part of the total potential dominates. Indeed, the effect reveals itself when the quark separation $`r_{qq}`$ becomes less than the background correlation length. But at separation $`r_{qq}0.1`$ fm the color-Coulomb interaction of quarks attains 300 MeV and grows rapidly when $`r_{qq}`$ diminishes. That is why to answer the question about the relevance of quark-diquark configurations we are to consider the motion of two quarks in their perturbative potential. In particular, the radius of diquark may be estimated using the Bohr formula, $`r_{qq}=2/(m_qC_F\alpha _s)`$, which yields $`r_{bb}0.3`$ fm and $`r_{cc}1.0`$ fm. More accurate values were obtained in the relativistic quark model , where relativistic corrections to perturbative potential as well as linear nonperturbative potential were accounted for. The radii of diquarks calculated in  for the ground states are $`r_{bb}=0.37`$ fm, $`r_{cc}=0.56`$ fm. We are to state that the radii of diquarks are comparable with the size of baryon, which perhaps signals that the formation of quark-diquark configuration is unprobable. Nevertheless, the quark-diquark approximation turns out usefull for the computations of the baryon spectra (see e.g. ). Account of the nonperturbative effects considered in (11)-(16) and Fig. 4 would spoil this approximation only a little.

4. To summarize, we have calculated nonperturbatively the static potential in baryon and demonstrated that it completely describes recent lattice results. The latter means in particular that the problem of the form of the potential is resolved. The potential in our approach has apparent $`Y`$-type structure, i.e. it depends only on the distances from the quarks to the string junction, the location of the latter being determined by the condition of the minimal total length of the string. The formal answer of the problem of the potential law is that the slope of the nonperturbative potential in dependence on the length of the string $`L`$ grows from zero to $`\sigma `$ when $`L`$ changes from zero to $``$1.5 fm. Therefore an effective slope at typical hadronic sizes may be chosen $`0.9\sigma `$.

The behavior of the baryon potential considered is a consequence of two profound properties of the strong interaction, the confinement of quarks and the definite value of the correlation length of nonperturbative gluon fields. The latter is directly related to the energy scale of the gluon excitations of the hadron spectra. It is the correlation length that induces the change of the slope of the potential in baryon. Note that for mesons situation differs. The slope of the nonperturbative static potential in mesons almost does not change at small distances due to the interference with the perturbative fields .

Apart from the change of the slope of the potential in baryon, another related effect was studied, namely the behavior of the potential when the length of the string is fixed. It was demonstrated in particular that the difference of the potentials in configurations with the string consisting of different (one, two or three) numbers of lines turns out to be proportional to $`\sigma T_g`$. The combination of the parameters directs immediately that the effect is induced by both the scale of confinement and correlations of gluon fields. The influence of the effect on the creation of the quark-diquark configuration is shown to be small. Last but not least, it would be interesting to study this effect in devoted lattice calculations.

The author is grateful to Yu.A.Simonov for many fruitful discussions. This work has been supported by RFBR grants 00-02-17836, 00-15-96786, and INTAS 00-00110, 00-00366.


text:
# Hamiltonians separable in cartesian coordinates and third-order integrals of motion

## I Introduction

In classical mechanics, an n-dimensional Hamiltonian system is called Liouville integrable if it allows $`n`$ functionally independent integrals of motion in involution (including the Hamiltonian), that is

$$\begin{array}{cc}\hfill \{H,X_i\}& =0,\hfill \\ \hfill \{X_i,X_j\}& =0,i,j.\hfill \end{array}$$

The Hamiltonian $`H=H(x_1,\mathrm{}x_n,p_1,\mathrm{},p_n)`$ and the integrals of motion $`X_i=X_i(x_1,\mathrm{}x_n,p_1,\mathrm{},p_n)`$ must be well defined functions on phase space (). The system is superintegrable if it allows more than $`n`$ functionally independent integrals, $`n`$ of them in involution. It is called maximally superintegrable if it allows $`2n1`$ integrals of motion. The best known superintegrable systems in $`n`$ dimensions are the harmonic oscillator $`V=\omega r^2`$ and the Coulomb potential $`V=\frac{\alpha }{r}`$, and they are indeed maximally superintegrable. This may be closely related to Bertrand’s theorem () which states that these are the only rotationally invariant systems for which all finite trajectories are closed.

In quantum mechanics, a Hamiltonian system is said to be integrable if there exists a set $`\{X_i\}`$ of $`n`$ well defined, algebraically independent operators (including the Hamiltonian) that commute pairwise. It is superintegrable if it possesses further independent operators, $`\{Y_j\}`$ that commute with the Hamiltonian. The $`Y_j`$ do not necessarily commute with each other, nor with the $`X_i`$.

The independence of operators in quantum mechanics remains to be defined rigorously . Since we are dealing here only with polynomial differential operators, we can proceed by analogy with the classical case, keeping in mind that a rigorous definition will be needed as soon as we will want to make some more general statements. This choice of a definition will be used only for discussion purposes, since we will find all potentials that admit third-order integrals and all their integrals. The results obtained will therefore hold for any definition of the independence of operators.

Integrable and superintegrable systems, both in quantum and in classical mechanics, attracted considerable interest in the last years. Extensive literature exists about systems with second-order integrals of motion, either in euclidian space , or in spaces with nonzero constant  or nonconstant curvature . As long as there was no magnetic field in the Hamiltonian, the quantum and classical integrals of motion obeyed the same determining equations, and therefore quantum and classical integrability were very similar. Both properties were related to separation of the Hamilton-Jacobi or Schrödinger equations, and also to exact solvability  and generalized symmetries .

Systems with higher-order integrals have been studied and classified as early as 1935 in a well-known paper by Drach. This paper considered classical Hamiltonians in complex Euclidian space. Efforts were made recently to understand and classify more completely systems with higher-order integrals in classical  and quantum mechanics . In spite of these efforts still relatively few such systems are known. This paper is the logical sequel of a systematic search for superintegrable systems with higher order integrals started in . Here we consider two-dimensional real Euclidian space with a one-particle Hamiltonian;

$$H=\frac{1}{2}\left(p_x^2+p_y^2\right)+V(x,y).$$

We request the existence of two additional integrals of motion, one of second order in the momenta and the other of third order.

The condition of existence for second-order integral implies, both in classical and quantum mechanics, that the Hamiltonian be separable in cartesian, polar, parabolic or elliptic coordinates. In this paper we consider potentials that are separable in cartesian coordinates;

$$H=\frac{1}{2}\left(p_x^2+p_y^2\right)+V_1(x)+V_2(y).$$

We found all such systems which admit third-order integrals. Quantum and classical mechanics will be treated simultaneously, for the conditions of existence of integrals of motion, even though not equivalent, are quite similar in both cases.

## II Existence of a third-order integral

In quantum and classical mechanics, the general third-order commuting operator

$$X=\underset{i+j=0}{\overset{3}{}}P_{ij}(x,y)p_x^ip_y^j$$

can be reduced to a much simpler form,

$$\begin{array}{cc}\hfill X=& \underset{\begin{array}{c}i,j,k\\ i+j+k=3\end{array}}{}A_{ijk}\{L_3^i,p_x^jp_y^k\}+\{g_1(x,y),p_x\}+\{g_2(x,y),p_y\}\hfill \\ \hfill L_3=& xp_yyp_x.\hfill \end{array}$$
(1)
where the $`A_{ijk}`$ are real constants, and the $`g_i`$ real functions. The bracket is the anticommutator. It is not needed in classical mechanics, but in quantum mechanics it allows us to get rid of terms with even powers of the $`p_i`$ and to make sure the operator is self-adjoint. Furthermore, its use allows us to see clearly the relations between the quantum and the classical case. Indeed, it was found in  that the requirement that the operator commutes (or Poisson-commutes) with the Hamiltonian

$$H=\frac{1}{2}\left(p_x^2+p_y^2\right)+V(x,y).$$

implies equations that behave well in the classical limit. Namely, commutativity implies

$`\begin{array}{cc}\hfill 0& =g_1V_x+g_2V_y{\displaystyle \frac{\mathrm{}^2}{4}}(f_1V_{xxx}+f_2V_{xxy}+f_3V_{xyy}+f_4V_{yyy}\hfill \\ & +8A_{300}(xV_yyV_x)+2(A_{210}V_x+A_{201}V_y)),\hfill \end{array}`$ (2)
$`(g_1)_x`$ $`=3f_1(y)V_x+f_2(x,y)V_y,`$ (3)
$`(g_2)_y`$ $`=f_3(x,y)V_x+3f_4(x)V_y,`$ (4)
$`(g_1)_y+(g_2)_x`$ $`=2\left(f_2(x,y)V_x+f_3(x,y)V_y\right),`$ (5)
in quantum mechanics, where

$`f_1(y)`$ $`=A_{300}y^3+A_{210}y^2A_{120}y+A_{030},`$
$`f_2(x,y)`$ $`=3A_{300}xy^22A_{210}xy+A_{201}y^2+A_{120}xA_{111}y+A_{021},`$
$`f_3(x,y)`$ $`=3A_{300}x^2y+A_{210}x^22A_{201}xy+A_{111}xA_{102}y+A_{012},`$
$`f_4(x)`$ $`=A_{300}x^3+A_{201}x^2+A_{102}x+A_{003}.`$
The equations in classical mechanics are obtained by setting $`\mathrm{}=0`$. We may notice from the quantum equations, or directly from the condition $`[H,X]=0`$, that we can express all quantum integrable potentials as $`\mathrm{}^2\stackrel{~}{V}(x,y)`$ where $`\stackrel{~}{V}`$ does not depend on $`\mathrm{}`$. It is often more natural and interesting though to choose arbitrary parameters contained in $`\stackrel{~}{V}`$ to be depending on $`\mathrm{}`$, so that the potential $`V`$ does not vanish in the classical limit. One may always verify though that through an appropriate transformation of the arbitrary parameters one can write the potential as being proportional to $`\mathrm{}^2`$

The three last equations, identical in the quantum and classical cases, yield a linear compatibility condition for $`V`$, which reads

$$\begin{array}{cc}\hfill 0=& f_3V_{xxx}+\left(2f_23f_4\right)V_{xxy}+\left(3f_1+2f_3\right)V_{xyy}f_2V_{yyy}\hfill \\ & +2\left(f_{2y}f_{3x}\right)V_{xx}+2\left(3f_{1y}+f_{2x}+f_{3y}3f_{4x}\right)V_{xy}+2\left(f_{2y}+f_{3x}\right)V_{yy}\hfill \\ & +\left(3f_{1yy}+2f_{2xy}f_{3xx}\right)V_x+\left(f_{2yy}+2f_{3xy}3f_{4xx}\right)V_y.\hfill \end{array}$$
(6)
Further nonlinear compatibility conditions can be obtained from (2) to (5) for the potential, and these are listed in . They are quite complicated though and were not used for the results stated in this paper.

## III Potentials separable in cartesian coordinates

If we set $`V=V_1(x)+V_2(y)`$ in equations (2) to (5), we find

$`\begin{array}{cc}\hfill 0& =g_1V_{1x}+g_2V_{2y}{\displaystyle \frac{\mathrm{}^2}{4}}(f_1V_{1xxx}+f_4V_{2yyy}\hfill \\ & +8A_{300}(xV_{2y}yV_{1x})+2(A_{210}V_{1x}+A_{201}V_{2y})),\hfill \end{array}`$ (7)
$`(g_1)_x`$ $`=3f_1(y)V_{1x}+f_2(x,y)V_{2y},`$ (8)
$`(g_2)_y`$ $`=f_3(x,y)V_{1x}+3f_4(x)V_{2y},`$ (9)
$`(g_1)_y+(g_2)_x`$ $`=2\left(f_2(x,y)V_{1x}+f_3(x,y)V_{2y}\right),`$ (10)
with $`\mathrm{}=0`$ in the classical case. Equations (8) and (9) are readily integrated, so in the cartesian case two equations remain to be solved.

The compatibility condition (6) allows us to find ODEs for $`V_1`$ and $`V_2`$. If we set alternatively $`y=0`$ and $`x=0`$, we find

$$(A_{210}x^2+A_{111}x+A_{012})V_1^{(3)}(x)+4(2A_{210}x+A_{111})V_1^{\prime \prime }(x)+12A_{210}V_1^{}(x)=ax+b$$
(11)
$$(A_{201}y^2A_{111}y+A_{021})V_2^{(3)}(y)+4(2A_{201}yA_{111})V_2^{\prime \prime }(y)+12A_{201}V_2^{}(y)=cy+d$$
(12)
The solutions to the homogeneous part of these equations are easily found and brought to a simple form by translations in $`x`$ and $`y`$. If we take the first one for definiteness, we have four different types of solution. When $`A_{210}0`$, we have two possible solutions,

$$\begin{array}{cc}\hfill V_{1hom}& =\frac{c_1}{(x+\alpha )^2}+\frac{c_2}{(x\alpha )^2}\hfill \\ \hfill V_{1hom}& =\frac{c_1}{x^2}+\frac{c_2}{x^3}.\hfill \end{array}$$

If $`A_{210}=0`$ and $`A_{111}0`$ we get

$$V_{1hom}=\frac{c_1}{x^2}+c_2x.$$

Finally if only $`A_{012}0`$, the solution may be brought to the form

$$V_{1hom}=c_2x^2+c_1x.$$

Special solutions are also simple. If $`A_{210}0`$, we have $`V_{1part}=\alpha x^2+\beta x`$. Otherwise, when $`A_{111}0`$, $`V_{1part}=\alpha x^3+\beta x^2`$, finally, if only $`A_{012}0`$, $`V_{1part}=\alpha x^4+\beta x^3`$. Provided that (11) or (12) do not vanish trivially, we can choose $`V_1`$ or $`V_2`$, respectively, amongst the following functions:

1. $`f_1=\frac{c_1}{(x+\alpha )^2}+\frac{c_2}{(x\alpha )^2}+c_3x^2+c_4x`$
2. $`f_2=\frac{c_1}{x^2}+\frac{c_2}{x^3}+c_3x^2+c_4x`$
3. $`f_3=\frac{c_1}{x^2}+c_2x^3+c_3x^2+c_4x`$
4. $`f_4=c_1x^4+c_2x^2+c_3x`$
5. $`f_5=c_1x^3+c_2x`$
6. $`f_6=c_1x^2`$
7. $`f_7=c_1x`$,

and then solve (7) to (10). These long but rather straightforward calculations yield the $`15`$ superintegrable potentials included in Table 1. Some of them are obviously particular cases of others, but we listed them separately to account for their additional integrals. Only the third-order integrals are listed, some of them being trivial consequences of lower-order ones. With the exception of the harmonic oscillator, potentials that have first-order integrals are not listed here for they were already presented in  with all their third-order integrals. The complete integrals of motion can be found in Appendix I.

Many of these potentials were not known. The only classical potentials among these are indicated with a $``$. These are well-known superintegrable potentials (see e.g. ), and all of them, except $`V_i`$, are in fact quadratically superintegrable.

All the potentials are superintegrable in the quantum case. We therefore notice that classical nontrivial potentials can have many different quantum equivalents. The classical harmonic oscillator $`V_a`$ can be seen as a limiting case of the quantum potentials $`V_a`$, $`V_c`$ and $`V_d`$, and also, if we set $`\alpha =\sqrt{\mathrm{}}/\omega `$, of $`V_e`$, $`V_f`$ and $`V_g`$, not to mention the similar potentials that can be obtained by permutations of $`x`$ and $`y`$. The anisotropic harmonic oscillator with ratio $`1:3`$ also admits many quantum deformations but, interestingly, the anisotropic oscillator with ratio $`1:2`$ does not admit such deformations. Notice also that if we want to deal with real potentials only, $`\alpha `$ must be either real or purely imaginary in potentials $`V_e`$, $`V_f`$, $`V_g`$ and $`V_k`$. Therefore these have as a classical limit harmonic oscillators with $`a>0`$.

All quantum superintegrable potentials reduce to classical ones when the classical limit is considered, sometimes in more than one way. For example, potentials $`V_e`$, $`V_f`$ and $`V_g`$ give the free motion potential instead of the harmonic oscillator if $`\alpha `$ remains constant as $`\mathrm{}0`$.

These potentials all satisfy the linear equations (11) and (12), and can be expressed as sums of simple superintegrable potentials.

Let us now set $`A_{210}=A_{111}=A_{012}=0`$ so that (11) vanishes trivially. We may also assume that $`V_1`$ does not take one of the forms (A.1) to (A.7), for we have already worked these cases out. This is quite useful, for if we set $`y=1`$ in (6), we obtain for $`V_1`$ an equation of the same form as (11) with different coefficients. These coefficients must therefore vanish, so $`A_{300}=A_{201}=A_{102}=0`$. This is a significant simplification that allows us to restrict our attention, when considering equation (12), to the following three cases:

1. $`V_2=ay^2;`$
2. $`V_2=ay;`$
3. $`A_{120}=A_{021}=0`$

Before we consider each case separately, it is worth noticing that potentials of the form $`V=V_1(x)`$ that admit third-order integrals independent of $`y`$ and $`p_y`$ should appear as solutions here, for the integral remains if we add a function $`V_2(y)`$ to these potentials. These potentials were found in  and  to satisfy equation

$$\mathrm{}^2V_1^2=4V_1^3g_2V_1g_3,$$
(13)
and can therefore be written as

$$V_1=\mathrm{}^2𝒫(x),$$
(14)
where $`𝒫(x)`$ is the Weierstrass elliptic function. Since the $`y`$ variable plays no role here and these potentials admit integrals with leading-order terms proportional to $`p_x^3`$, these solutions will appear in cases (i) to (iii).

### A. Case i: $`V=V_1(x)+ay^2`$

When $`V=V_1(x)+ay^2`$ and $`a0`$, we find that $`A_{300}=0`$, and the following two equations must be satisfied:

$`0`$ $`=A_{030}\left(\mathrm{}^2V_1^{(3)}6(V_1^2)^{}\right)+\gamma _1V_1^{}`$ (15)
$`\begin{array}{cc}\hfill 0& =A_{120}\left(\mathrm{}^2V_1^{(4)}24a(xV_1)^{}+6(V_1^2)^{\prime \prime }4ax^2V_1^{\prime \prime }+8a^2x^2\right)\hfill \\ & +8aA_{021}\left(2ax(xV_1^{})^{}2V_1^{}\right)+4\eta (2aV_1^{\prime \prime }).\hfill \end{array}`$ (16)
where $`\gamma _1`$ and $`\eta `$ are arbitrary constants. When $`A_{030}0`$, equation (15) is equivalent to (13) (up to a translation of $`V_1`$ to get rid of $`\gamma _1`$), hence its solutions are of the form (14). These potentials cannot satisfy simultaneously equation (16) for nontrivial parameters. This can be observed by expanding (13) in series around $`x=0`$ and substituting the result in (16). Therefore solutions given by $`A_{030}0`$ are of no special interest here.

Let us now set $`A_{030}=\gamma _1=0`$. Equation (16) can be greatly simplified. We assume that $`A_{120}0`$, for otherwise equation (16) can be solved to give potential (1). Then by an appropriate translation of $`x`$ and $`V`$, we can get rid of the terms involving $`A_{021}`$ and $`\eta `$ and finally divide by $`A_{120}`$;

$$0=\mathrm{}^2V_1^{(4)}24a(xV_1)^{}+6(V_1^2)^{\prime \prime }4ax^2V_1^{\prime \prime }+8a^2x^2.$$
(17)
This equation admits a first integral, namely

$$k=\mathrm{}^2\left(xV_1^{\prime \prime \prime }V_1^{\prime \prime }\right)+4x\left(ax^23V_1\right)V_1^{}+6V_1^2+12ax^2V_12a^2x^4.$$
(18)
Both (17) and (18) can be simplified by setting $`V_1=W(x)+ax^2/3`$. Then

$$\mathrm{}^2W^{(4)}=12WW^{\prime \prime }+12(W^{})^2+bxW^{}+2bW\frac{1}{6}b^2x^2$$
(19)
with $`b=8a0`$ for (17), and

$$k_2=3\mathrm{}^2\left(xW^{\prime \prime \prime }W^{\prime \prime }\right)18x(W^2)^{}+2(2ax^2+3W)^2.$$
(20)
for the first integral. Equation (19) is well known. It is equivalent to equations (3.16) in  and (2.17) in , which were obtained by nonclassical reduction of the Boussinesq Equation. It was also shown in  to be a nonclassical reduction of the Kadomtsev-Petviashvili equation. It has the Painlevé property, and, when $`b0`$, its solution, given in  (equation 2.88) may be written in terms of the fourth transcendant function of Painlevé, namely

$$W=\frac{\mathrm{}}{2}b_1P_4^{}(x,\frac{b}{\mathrm{}^2})\frac{1}{2}bP_4^2(x,\frac{b}{\mathrm{}^2})\frac{1}{2}bxP_4(x,\frac{b}{\mathrm{}^2})\frac{1}{6}(\frac{b}{2}x^2+\mathrm{}^2K_1\mathrm{}b_1)$$
(21)
where $`b_1\pm \sqrt{b}=\pm \sqrt{8a}`$ and $`P_4(x,\frac{b}{\mathrm{}^2})=P_4(x,\frac{b}{\mathrm{}^2},K_1,K_2)`$ is the fourth transcendant function of Painlevé, and therefore satisfies equation

$$\begin{array}{cc}\hfill P_4^{\prime \prime }(x,\alpha )=& \frac{(P_4^{}(x,\alpha ))^2}{2P_4(x,\alpha )}\frac{3\alpha }{2}P_4(x,\alpha )^32\alpha xP_4(x,\alpha )^2\hfill \\ & (\frac{\alpha }{2}x^2+K_1)P_4(x,\alpha )+\frac{K_2}{P_4(x,\alpha )}.\hfill \end{array}$$
(22)
$`K_1`$ and $`K_2`$ are integration constants. The potential therefore reads

$$\begin{array}{cc}\hfill V(x,y)=& a(x^2+y^2)+\frac{\mathrm{}}{2}b_1P_4^{}(x,\frac{8a}{\mathrm{}^2})+4aP_4^2(x,\frac{8a}{\mathrm{}^2})\hfill \\ & +4axP_4(x,\frac{8a}{\mathrm{}^2})+\frac{1}{6}(\mathrm{}^2K_1+\mathrm{}b_1).\hfill \end{array}$$
(23)
This potential admits as special cases two anisotropic harmonic oscillators, $`V=a(x^2+y^2)`$ when $`K_2=0`$ (and $`P_4=0`$), and $`V=a(x^2/9+y^2)`$ when $`K_1=0`$ and $`K_2=1/18`$ (and $`P_4=x/3`$), as well as all their quantum deformations that have the form $`V=a(p^2x^2+y^2)+f(x)`$, that is potentials $`V_d`$, $`V_e`$, $`V_j`$ and $`V_k`$ (up to a permutation of $`x`$ and $`y`$).

The constant term in the potential, $`(\mathrm{}^2K_1+\mathrm{}b_1)`$, can be set to zero, but we will keep it in order to be able to write the quantum and classical integrals in a unified way.

In classical mechanics, the equation (18) with $`\mathrm{}=0`$ admits a first integral, which reads

$$c=\frac{(9V_1ax^2)(V_1ax^2)^3+\frac{k^2}{4}k(V_1ax^2)(3V_1+ax^2)}{x^2}.$$

We may therefore write the solution for $`V_1`$ implicitly as

$$cx^2d^2+2d(V_1ax^2)(3V_1+ax^2)=(9V_1ax^2)(V_1ax^2)^3,$$
(24)
where $`c`$ and $`d`$ are arbitrary constants. If $`c=d=0`$, we find either the familiar anisotropic harmonic oscillators, or a potential obtained by joining at $`x=0`$ two halves of anisotropic harmonic oscillators with different ratios. Even though this potential does not have a continuous second derivative, it can be obtained as a limiting case of the family of smooth superintegrable potentials (26).

In the general case the potential may be expressed as the root of a fourth-order polynomial with three arbitrary parameters ($`a`$, $`c`$ and $`d`$), although we can set $`a=1`$ and one of the other two coefficients to $`\pm 1`$ by scaling $`x`$ and $`V_1`$.

Since equation (24) may describe new two-dimensional potentials with bounded motion, it is worth studying for interesting special cases. Indeed, if we assume $`a>0`$ and $`d0`$, we can consider the case $`c=128a^4\stackrel{~}{d}^3/27^2`$ and $`d=4a^2\stackrel{~}{d}^2/27`$. Four solutions exist in that case, two of which are translated harmonic oscillators. The remaining two potentials are

$$V=ay^2+\frac{a}{9}\left(2\stackrel{~}{d}+5x^2\pm 4x\sqrt{(\stackrel{~}{d}+x^2)}\right).$$
(25)
Up to an additive constant, potentials (25) are equal to

$$V=ay^2+\frac{a}{9}\left(x\pm 2\sqrt{\stackrel{~}{d}+x^2}\right)^2.$$
(26)
Those potentials are smooth interpolations between anisotropic harmonic oscillators with ratio $`1:1`$ and $`1:3`$. When $`\stackrel{~}{d}=0`$ (and therefore $`c=d=0`$) they are the junctions of the two halves of harmonic oscillators mentioned above.

The integral of motion is similar in quantum and classical mechanics, and reads

$$X=\{L,p_x^2\}+\{ax^2y3yV_1(x),p_x\}\frac{1}{2a}\{\frac{\mathrm{}^2}{4}V_{1xxx}+(ax^23V_1)V_{1x},p_y\},$$
(27)
where $`V_1`$ is a solution to equation (24). Therefore the integral for (26) must be slightly modified to take into account the constant term we removed.

### B. Case ii: $`V=ay+V_1(x)`$

Here we find again two equations:

$`0`$ $`=\mathrm{}^2A_{120}V_1^{(3)}6A_{120}(V_1^2)^{}+12a^2A_{003}+\gamma _1V_1^{},`$ (28)
$`0`$ $`={\displaystyle \frac{A_{030}}{4}}\left(6(V_1^2)^{}+\mathrm{}^2V_1^{(3)}\right)^{}{\displaystyle \frac{a}{2}}A_{120}\left(6(xV_1)^{}+x^2V_1^{\prime \prime }\right)`$
$`aA_{021}\left((xV_1^{})^{}+2V_1^{}\right)\gamma _2V_1^{\prime \prime }+a\gamma _1.`$
This time we cannot treat the equations separately, but we can, when $`A_{120}0`$, translate $`V`$ to annihilate $`\gamma _1`$. Then we can substitute the first equation in the second to get rid of $`A_{030}`$, and finally translate $`x`$ to get rid of $`A_{021}`$. We can then solve the remaining system, first by solving the second, linear equation, and then by substituting the result in the first one. The only solution remaining is then potential $`V_n=ay+\frac{\mathrm{}^2}{x^2}`$. Hence we can set $`A_{120}=0`$, and therefore also $`A_{003}=\gamma _1=0`$ (as $`V_1bx`$ and $`a0`$). In the second equation we can set $`\gamma _2=0`$. If then $`A_{030}=0`$, we find $`V_1=a/x^2`$ which was already classified. If $`A_{030}0`$ and $`A_{021}=0`$, though, the potential $`V_1(x)`$ is solution to

$$\mathrm{}^2V_1^{\prime \prime }=6V_1^2+\lambda x+k$$
(29)
If $`\lambda 0`$ we can set $`k=0`$, and we find that $`V_1`$ can be expressed in terms of the first Painlevé transcendent:

Comme $`\lambda `$ est une constante arbitraire, il est plus pratique d’ crire $`\omega ^5=\frac{\lambda }{\mathrm{}^4},`$

$$V=ay+\mathrm{}^2\omega ^2P_1\left(\omega x\right),$$
(30)
with $`\omega ^5=\frac{\lambda }{\mathrm{}^4}.`$

The integral of motion is

$$X=2p_x^3+3\{V_2(x),p_x\}+\{\frac{\mathrm{}^4\omega ^5}{4a},p_y\}.$$

Notice that in order to consider the limiting case $`a=0`$, we can multiply the integral by $`a`$, and thus we find a potential that depends only on $`x`$ with a trivial $`p_y`$ integral. We can also set $`\omega ^5=a\omega _{}^{5}{}_{}{}^{}`$ first, in which case we find the potentials (14). If we look for the classical limit, we find

$$V=ay+b\sqrt{x}.$$
(31)
with the integral

$$X=2p_x^3+3b\{\sqrt{x},p_x\}\{\frac{3b^2}{2a},p_y\}.$$

Returning to equation (29) and assuming $`\lambda =0`$, we find again that $`V_1`$ has the form (14), and the integral depends only on $`x`$ and $`p_x`$.

Finally, if we set $`A_{030}A_{021}0`$, we have to solve the equation

$$0=\left(\frac{3}{2}(V_1^2)^{}+\frac{\mathrm{}^2}{4}V_1^{(3)}\right)^{}+b\left((xV_1^{})^{}+2V_1^{}\right),b=\frac{aA_{021}}{A_{030}}0.$$
(32)
It can be integrated once to give

$$C_1=\frac{3}{2}(V_1^2)^{}+\frac{\mathrm{}^2}{4}V_1^{(3)}+b\left((xV_1^{})+2V_1\right).$$

We can set $`C_1=0`$ by translations in $`x`$ and $`V`$. The equation admits a first integral which reads

$$2b\mathrm{}^2\left(V_1(x)bx\right)V_1^{\prime \prime }(x)+b\mathrm{}^2\left(2bV_1^{}(x)\right)V_1^{}(x)8bV_1(x)\left(V_1(x)bx\right)^2=k_1.$$
(33)
In the classical case, $`\mathrm{}=0`$, this is enough to solve. The solution, which can be written implicitly in a more compact form, is given by

$$d=V_1(V_1bx)^2,$$
(34)
where $`d`$ is an arbitrary constant. When $`d=0`$ we find the familiar case $`V=bx+ay`$ which admits a first-order integral. We may notice that the implicit form of the solution is somewhat similar to (24).

In the quantum case, we notice that the transformation $`W(x)=V_1(x)bx`$, that preserves the Painlevé property, simplifies equation (33) that becomes

$$\mathrm{}^2\left(2WW^{\prime \prime }W^2\right)8\left(W+bx\right)W^2=k_2.$$
(35)
or

$$W^{\prime \prime }=\frac{W^2}{2W}+\frac{4W^2}{\mathrm{}^2}+\frac{4bxW}{\mathrm{}^2}+\frac{k_2}{2\mathrm{}^2W}$$
(36)
We can also substitute $`Y(x)=\sqrt{W(x)}`$ to find

$$\mathrm{}^2Y^{\prime \prime }=2\left(Y^2+bx\right)Y+\frac{k_2}{4Y^3}.$$

If $`k_2=0`$, we can set $`\mathrm{}=1`$ and $`b=1/2`$ by the change of variables

$$Y=(2\mathrm{}b)^{\frac{1}{3}}Z,x=\left(\frac{\mathrm{}^2}{2b}\right)^{\frac{1}{3}}\xi .$$

The solution for $`Z(\xi )`$ is then a special case of the second Painlevé transcendant, defined by the equation

$$P_2^{\prime \prime }(x,\alpha )=2P_2(x,\alpha )^3+xP_2(x,\alpha )+\alpha .$$

The solution for $`V`$ is

$$V(x,y)=bx+ay+(2\mathrm{}b)^{\frac{2}{3}}P_2^2(\left(\frac{2b}{\mathrm{}^2}\right)^{\frac{1}{3}}x,0),$$
(37)
with the particular case $`V=bx+ay`$.

Now if $`k_20`$, we use (36), that can be normalized by a change of variables,

$$x=\left(\frac{\mathrm{}^2}{(4b)}\right)^{\frac{1}{3}}\xi ,W=\frac{\sqrt{k_2}}{(4\mathrm{}b)^{\frac{1}{3}}}Y.$$

We then find

$$Y^{\prime \prime }=\frac{Y^2}{2Y}+4\beta Y^2\xi Y\frac{1}{2Y},$$
(38)
where $`\beta =\sqrt{k_2}/(4\mathrm{}b)0`$ is an arbitrary constant. Equation (38) corresponds to case XXXIV, p. 340, in . The solution for $`Y(\xi )`$ reads

$$2\beta Y=P_2^{}(\xi ,2\beta \frac{1}{2})+P_2(\xi ,2\beta \frac{1}{2})^2+\frac{\xi }{2}$$
(39)
where $`P_2`$ is once again the second Painlevé transcendant. Since $`\beta `$ is an arbitrary constant, we can set $`\kappa =2\beta \frac{1}{2}`$.

Back to the original variables, we get

$$V(x,y)=ay+\left(2\mathrm{}^2b^2\right)^{\frac{1}{3}}\left(P_2^{}((4b/\mathrm{}^2)^{\frac{1}{3}}x,\kappa )+P_2^2((4b/\mathrm{}^2)^{\frac{1}{3}}x,\kappa )\right)$$
(40)
This potential admits $`V=ay`$ and $`V=ay+\mathrm{}^2/x^2`$ as particular cases.

In all cases with $`V=ay+V_1(x)`$ and $`A_{012}A_{030}0`$ the integral of motion is

$$2ap_x^32bp_x^2p_y+a\{3V_1(x)bx,p_x\}2b\{V_1(x),p_y\}.$$
(41)
Limiting values for $`a`$, $`b`$ and $`\mathrm{}`$ give either known potentials or trivial integrals.

### C. Case iii: $`A_{120}=A_{021}=0`$

Since we set here $`A_{120}=A_{021}=0`$, we assume $`A_{030}A_{003}0`$ (otherwise we would find potentials of the form (14)). The conditions for the existence of a third-order operator then read

$$\begin{array}{c}\hfill \mathrm{}^2V_1^{\prime \prime }(x)=6V_1^2(x)+A_{003}\sigma x\\ \hfill \mathrm{}^2V_2^{\prime \prime }(y)=6V_2^2(y)A_{030}\sigma y.\end{array}$$
(42)
If $`\sigma `$ vanishes, we find the potential

$$V=\mathrm{}^2(𝒫(x)+𝒫(y)),$$

with two integrals that each depend on only one variable, as could have been predicted from the results of equation (13). If $`\sigma `$ does not vanish, let us set $`b_1=A_{003}\sigma `$ and $`b_2=A_{030}\sigma `$.

If $`\mathrm{}=0`$, we find that

$$V=\pm \sqrt{\beta _1x}\pm \sqrt{\beta _2y}$$

is superintegrable, where the $`\beta _i=b_i/6`$ are arbitrary constants.

This defines similar real potentials on each quadrant of the plane. We can again patchwork the pieces to find real continuous potentials defined everywhere, e.g. $`V=c_1\sqrt{|x|}+c_2\sqrt{|y|}`$, although in that case neither the Hamiltonian nor the integral are differentiable at the equilibrium point.

Finally, if $`\mathrm{}0`$, we find that $`V_1`$ and $`V_2`$ can be both written using the first Painlevé transcendent.

$$V=\mathrm{}^2\omega _1^2P_1\left(\omega _1x\right)+\mathrm{}^2\omega _2^2P_1\left(\omega _2x\right)$$
(43)
where the $`\omega _i=(b_i/\mathrm{}^4)^{1/5}`$ are arbitrary constants.

The integral of motion is

$$X=2b_2p_x^32b_1p_y^3+3b_2\{V_1(x,b_1),p_x\}3b_1\{V_2(y,b_2),p_y\},$$
(44)
both in quantum and classical mechanics.

## IV Conclusion

We have found all systems in two-dimensional euclidian space that admit separation of variables in cartesian coordinates and at least one third-order integral, both in quantum and classical mechanics. Many new superintegrable potentials were found, and, interestingly, all the quantum superintegrable potential are found as solutions of equations having the Painlevé property, and this is probably not accidental. Many of the quantum potentials can in fact be written in terms of different transcendent functions of Painlevé, and many are related to group invariant solutions of the Korteweg-de Vries equation, and to reductions of the Boussinesq and KP equations. All classical integrable potentials were found to be limiting cases of quantum ones. They do not obey to equations having the Painlevé property, though, for many classical integrable potentials have movable branch points of the form $`\sqrt{xb}`$. Thus, in that respect, quantum integrable potentials behave more regularly than classical ones. A natural question is what does the Painlevé property, in quantum mechanics, tell us about classical integrable potentials. Since this paper deals mostly with the classification of superintegrable systems, the consideration of such questions regarding their properties and their solutions will be postponed to a future article.

Our investigation provided interesting new examples of the differences between quantum and classical integrability. A systematic search for systems with higher-order integrals is therefore a useful task since we still know little about such systems, and they are likely to share interesting properties.

Note added in proof: The future article mentioned earlier () indeed shows that superintegrability of separable systems tells a lot about the physical properties of the separated systems. Even though the systems considered here can all be separated in two independent one-dimensional systems, two-dimensional superintegrability provides information on them that is far from obvious from a one-dimensional perspective.

Equations (2) to (5) are the necessary and sufficient conditions for the existence of third-order integrals. The general solution to these equation is highly nontrivial, and it is not likely that a direct approach will lead to such a solution for integrals of order higher than four, without the use of new methods. In order to develop such methods, it would be useful to have a rigorous definition of quantum integrability.

## Acknowledgements

The author thanks gratefully P. Winternitz for his very helpful comments and suggestions, and also J. Bérubé for useful discussions. The author benefits from a NSERC ES-A fellowship.

## V Annex I

Here is the complete list of two-dimensional Hamiltonians separable in cartesian coordinates that admit at least one third-order integral. We also gave all their third-order integrals. Potentials depending on only one cartesian coordinates are not listed here, for they were classified in . Many of the potentials listed below were already known to be superintegrable, but we listed them here for completeness and in some cases to take into account additional third-order integrals. Potentials in boxes are those who do not admit enough first- or second-order integrals to make them superintegrable. Most of them were not known before.

Some potentials (namely $`V_p`$ and $`V_q`$) are not superintegrable according to the usual definition, for their integrals have simple though nontrivial relations with the Hamiltonian. For definitions of parameters the reader is referred to the body of the article. In order to write the integrals of motion in a compact form we often use the notation $`V(x,y)=V_1(x)+V_2(y)=V_1+V_2`$.

### A. Quantum potentials

1. $`V=a(x^2+y^2)X_1=L^3X_2=\{L,p_xp_y\}+a\{2x^2y,p_y\}a\{2xy^2,p_x\}X_3=\{L,p_y^2\}+2a\left(\{xy^2,p_y\}\{y^3,p_x\}\right)X_4=\{L,p_x^2\}2a\left(\{x^2y,p_x\}\{x^3,p_y\}\right)`$
2. $`V=a(x^2+y^2)+\frac{b}{x^2}+\frac{c}{y^2}X_1=\{L,p_xp_y\}+\{xy\left(\frac{2b}{x^3}+2ax\right),p_y\}\{xy\left(\frac{2c}{y^3}+2ay\right),p_x\}`$
3. $`V=a(x^2+y^2)+\frac{\mathrm{}^2}{x^2}+\frac{\mathrm{}^2}{y^2}X_1=2L^3\mathrm{}^2\{\frac{3x^2}{y}+2y+\frac{3y^3}{x^2},p_x\}+\mathrm{}^2\{\frac{3y^2}{x}+2x+\frac{3x^3}{y^2},p_y\}X_2=\{L,p_xp_y\}+\{xy\left(\frac{2\mathrm{}^2}{x^3}+2ax\right),p_y\}\{xy\left(\frac{2\mathrm{}^2}{y^3}+2ay\right),p_x\}`$
4. $`V=a(x^2+y^2)+\frac{\mathrm{}^2}{y^2}X_1=2L^3\mathrm{}^2\{2y+\frac{3x^2}{y},p_x\}+\mathrm{}^2\{\frac{3x^3}{y^2}+2x,p_y\}X_2=\{L,p_xp_y\}2\{axy^2\frac{\mathrm{}^2x}{y^2},p_x\}+2a\{x^2y,p_y\}X_3=\{L,p_y^2\}\{2ay^3+\mathrm{}^2\frac{1}{y},p_x\}+\{\mathrm{}^2\frac{3x}{y^2}+2axy^2,p_y\}`$
5. $`\overline{)V=\mathrm{}^2\left(\frac{1}{8\alpha ^4}(x^2+y^2)+\frac{1}{(x\alpha )^2}+\frac{1}{(x+\alpha )^2}\right)}X_1=2L^33\alpha ^2\{L,p_y^2\}+\frac{\mathrm{}^2}{4}\{y(8+\frac{3y^2}{\alpha ^2}\frac{24y^2(x^2+\alpha ^2)}{(x\alpha )^2(x+\alpha )^2}),p_x\}+\frac{\mathrm{}^2}{4}\{x(8\frac{3y^2(x^410\alpha ^2x^224\alpha ^4)}{\alpha ^2(x\alpha )^2(x+\alpha )^2}),p_y\}X_2=\{L,p_x^2\}+\mathrm{}^2\{y(\frac{4\alpha ^2x^2}{4\alpha ^4}\frac{6(x^2+\alpha ^2))}{(x^2\alpha ^2)^2}),p_x\}+\mathrm{}^2\{\frac{x(x^24\alpha ^2)}{4\alpha ^4}\frac{2x}{x^2\alpha ^2}+\frac{4x(x^2+\alpha ^2)}{(x\alpha )^2(x+\alpha )^2},p_y\}`$
6. $`\overline{)V=\mathrm{}^2\left(\frac{1}{8\alpha ^4}(x^2+y^2)+\frac{1}{y^2}+\frac{1}{(x+\alpha )^2}+\frac{1}{(x\alpha )^2}\right)}X_1=2L^33\alpha ^2\{L,p_y^2\}+\mathrm{}^2\{\frac{3y^3}{4\alpha ^2}+\frac{6y^3(x^2+\alpha ^2)}{(x\alpha )^2(x+\alpha )^2}\frac{3(x^2\alpha ^2)}{y}2y,p_x\}+3\mathrm{}^2\{x(\frac{x^23\alpha ^2}{y^2}\frac{3y^28\alpha ^2}{12\alpha ^2}\frac{2y^2}{x^2\alpha ^2}+\frac{4y^2(x^2+\alpha ^2)}{(x\alpha )^2(x+\alpha )^2}),p_y\}`$
7. $`\overline{)V=\mathrm{}^2\left(\frac{1}{8\alpha ^4}(x^2+y^2)+\frac{1}{(y\alpha )^2}+\frac{1}{(x\alpha )^2}+\frac{1}{(y+\alpha )^2}+\frac{1}{(x+\alpha )^2}\right)}X_1=2L^33\alpha ^2(\{L,p_x^2\}+\{L,p_y^2\})+\frac{\mathrm{}^2}{4}\{y(124+\frac{3(x^2+y^2)}{\alpha ^2}+\frac{24(x^25y^2)}{y^2\alpha ^2}\frac{144x^2}{x^2\alpha ^2}+\frac{24(3x^2y^2)(x^2+\alpha ^2)}{(x\alpha )^2(x+\alpha )^2}+\frac{48(y^2x^2)(y^2+\alpha ^2)}{(y\alpha )^2(y+\alpha )^2}),p_x\}\frac{\mathrm{}^2}{4}\{x(124+\frac{3(x^2+y^2)}{\alpha ^2}\frac{24(5x^2y^2)}{x^2\alpha ^2}\frac{144x^2}{y^2\alpha ^2}\frac{24(x^23y^2)(y^2+\alpha ^2)}{(y\alpha )^2(y+\alpha )^2}+\frac{48(x^2y^2)(x^2+\alpha ^2)}{(x\alpha )^2(x+\alpha )^2}),p_y\}`$
8. $`V=a(4x^2+y^2)+\frac{b}{y^2}+cxX_1=2p_xp_y^2+\{2ay^2+\frac{2b}{y^2},p_x\}+\{8axy+cy,p_y\}`$
9. $`\overline{)V=a(9x^2+y^2)}X_1=\{L,p_y^2\}+\frac{2}{3}a\{y^3,p_x\}6a\{xy^2,p_y\}`$
10. $`\overline{)V=a(9x^2+y^2)+\frac{\mathrm{}^2}{y^2}}X_1=\{L,p_y^2\}+\{\frac{2ay^3}{3}\frac{\mathrm{}^2}{y},p_x\}+\{3x\left(2ay^2+\frac{\mathrm{}^2}{y^2}\right),p_y\}`$
11. $`\overline{)V=\mathrm{}^2\left(\frac{1}{8\alpha ^4}(9x^2+y^2)+\frac{1}{(y+\alpha )^2}+\frac{1}{(y\alpha )^2}\right)}X_1=\{L,p_y^2\}+\mathrm{}^2\{y\left(\frac{y^2}{12\alpha ^4}\frac{8\alpha ^2}{(y^2\alpha ^2)^2}\frac{2}{y^2\alpha ^2}\right),p_x\}+\frac{3\mathrm{}^2}{4}\{x(\frac{8(y^2+\alpha ^2)}{(y^2\alpha ^2)^2}\frac{y^2}{\alpha ^4}),p_y\}`$
12. $`\overline{)V=\frac{\mathrm{}^2}{x^2}+\frac{a}{y^2}}X_1=\{L^2,p_x\}2\mathrm{}^2\{\frac{y^2}{x^2},p_y\}+\{3\mathrm{}^2\frac{y^2}{x^2}+2a\frac{x^2}{y^2}+\frac{\mathrm{}^2}{2},p_x\}X_2=\{L,p_xp_y\}2\mathrm{}^2\{\frac{y}{x^2},p_y\}+2a\{\frac{x}{y^2},p_x\}X_3=2p_x^3+\{\frac{3\mathrm{}^2}{x^2},p_x\}`$
13. $`\overline{)V=\frac{\mathrm{}^2}{x^2}+\frac{\mathrm{}^2}{y^2}}X_1=2L^3\mathrm{}^2\{2y+\frac{3x^2}{y}+\frac{3y^3}{x^2},p_x\}+\mathrm{}^2\{2x+\frac{3y^2}{x}+\frac{3x^3}{y^2},p_y\}X_2=\{L^2,p_x\}2\mathrm{}^2\{\frac{y}{x},p_y\}+\{3\mathrm{}^2\frac{y^2}{x^2}+2\mathrm{}^2\frac{x^2}{y^2}+\frac{\mathrm{}^2}{2},p_x\}X_3=\{L^2,p_y\}2\mathrm{}^2\{\frac{x}{y},p_x\}+\{3\mathrm{}^2\frac{x^2}{y^2}+2\mathrm{}^2\frac{y^2}{x^2}+\frac{\mathrm{}^2}{2},p_x\}X_4=\{L,p_xp_y\}2\mathrm{}^2\{\frac{y}{x^2},p_y\}+2\mathrm{}^2\{\frac{x}{y^2},p_x\}X_5=2p_x^3+\{\frac{3\mathrm{}^2}{x^2},p_x\}X_6=2p_y^3+\{\frac{3\mathrm{}^2}{y^2},p_y\}`$
14. $`V=ax+\frac{\mathrm{}^2}{y^2}X_1=\{L,p_y^2\}\{\frac{\mathrm{}^2}{y},p_x\}+\{3\frac{\mathrm{}^2x}{y^2}a\frac{y^2}{2},p_y\}X_2=2p_y^3+\{3\frac{\mathrm{}^2}{y^2},p_y\}X_3=2p_xp_y^2+2\mathrm{}^2\{\frac{1}{y^2},p_x\}+a\{y,p_y\}`$
15. $`\overline{)V=\mathrm{}^2𝒫(y)+V(x)}X_1=2p_y^3+\{3\mathrm{}^2𝒫(y),p_y\}`$
16. $`\overline{)V=\mathrm{}^2(𝒫(x)+𝒫(y))}X_1=2p_x^3+\{3\mathrm{}^2𝒫(x),p_x\}X_2=2p_y^3+\{3\mathrm{}^2𝒫(y),p_y\}`$
17. $`\overline{)V=\mathrm{}^2\omega _1^2P_1\left(\omega _1x\right)+\mathrm{}^2\omega _2^2P_1\left(\omega _2x\right)}X_1=2\omega _2^5p_x^32\omega _1^5p_y^3+3\omega _2^5\{V_1(x,\omega _1),p_x\}3\omega _1^5\{V_2(y,\omega _2),p_y\}`$
18. $`\overline{)V=a(x^2+y^2)+\frac{\mathrm{}}{2}b_1P_4^{}(x,\frac{8a}{\mathrm{}^2})+4aP_4^2(x,\frac{8a}{\mathrm{}^2})+4axP_4(x,\frac{8a}{\mathrm{}^2})+\frac{1}{6}(\mathrm{}^2K_1+\mathrm{}b_1)}X_1=\{L,p_x^2\}+\{ax^2y3yV_1,p_x\}\frac{1}{2a}\{\frac{\mathrm{}^2}{4}V_{1xxx}+(ax^23V_1)V_{1x},p_y\}`$
19. $`\overline{)V=ay+\mathrm{}^2\omega ^2P_1(\omega x)}X_1=2p_x^3+3\mathrm{}^2\omega ^2\{P_1(\omega x),p_x\}+\{\frac{\omega ^5\mathrm{}^4}{4a},p_y\}`$
20. $`\overline{)V=bx+ay+(2\mathrm{}b)^{\frac{2}{3}}P_2^2(\left(\frac{2b}{\mathrm{}^2}\right)^{\frac{1}{3}}x,0)}X_1=2ap_x^32bp_x^2p_y+a\{3V_1(x)bx,p_x\}2b\{V_1(x),p_y\}`$
21. $`\overline{)V=ay+\left(2\mathrm{}^2b^2\right)^{\frac{1}{3}}\left(P_2^{}((4b/\mathrm{}^2)^{\frac{1}{3}}x,\kappa )+P_2^2((4b/\mathrm{}^2)^{\frac{1}{3}}x,\kappa )\right)}X_1=2ap_x^32bp_x^2p_y+a\{3V_1(x)bx,p_x\}2b\{V_1(x),p_y\}`$

### B. Classical potentials

1. $`V=a(x^2+y^2)X_1=L^3X_2=\{L,p_xp_y\}+a\{2x^2y,p_y\}a\{2xy^2,p_x\}X_3=\{L,p_y^2\}+2a\left(\{xy^2,p_y\}\{y^3,p_x\}\right)X_4=\{L,p_x^2\}2a\left(\{x^2y,p_x\}\{x^3,p_y\}\right)`$
2. $`V=a(x^2+y^2)+\frac{b}{x^2}+\frac{c}{y^2}X_1=\{L,p_xp_y\}+\{xy\left(\frac{2b}{x^3}+2ax\right),p_y\}\{xy\left(\frac{2c}{y^3}+2ay\right),p_x\}`$
3. $`V=a(4x^2+y^2)+\frac{b}{y^2}+cxX_1=2p_xp_y^2+\{2ay^2+\frac{2b}{y^2},p_x\}+\{8axy+cy,p_y\}`$
4. $`\overline{)V=a(9x^2+y^2)}X_1=\{L,p_y^2\}+\frac{2}{3}a\{y^3,p_x\}6a\{xy^2,p_y\}`$
5. $`\overline{)V=\pm \sqrt{\beta _1x}\pm \sqrt{\beta _2y}}X_1=2\beta _2p_x^32\beta _1p_y^3+3\beta _2\{V_1(x,6\beta _1),p_x\}3\beta _1\{V_2(y,6b_2),p_y\}`$
6. $`\overline{)V=ay^2+V_1}`$ where $`V_1`$ satisfies equation (24)
$`X_1=\{L,p_x^2\}+\{ax^2y3yV_1,p_x\}\frac{1}{2a}\{(ax^23V_1)V_{1x},p_y\}`$
7. $`\overline{)V=ay+b\sqrt{x}}X_1=2p_x^3+3b\{\sqrt{x},p_x\}\{\frac{3b^2}{2a},p_y\}`$
8. $`\overline{)V=ay+V_1(x)}`$, where $`(V_1bx)^2V_1=d`$
$`X_1=2ap_x^32bp_x^2p_y+a\{3V_1(x)bx,p_x\}2b\{V_1(x),p_y\}`$