{'fraction_non_alphanumeric': 0.029869798315037017, 'fraction_numerical': 0.018964951311134615, 'mean_word_length': 4.347113884555382, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The trispectrum of the Cosmic Microwave Background on sub-degree angular scales: an analysis of the BOOMERanG data

The trispectrum of the cosmic microwave background can be used to assess the level of non-Gaussianity on cosmological scales. It probes the fourth order moment, as a function of angular scale, of the probability distribution function of fluctuations and has been shown to be sensitive to primordial non-gaussianity, secondary anisotropies (such as the Ostriker-Vishniac effect) and systematic effects (such as astrophysical foregrounds). In this paper we develop a formalism for estimating the trispectrum from high resolution sky maps which incorporates the impact of finite sky coverage. This leads to a series of operations applied to the data set to minimize the effects of contamination due to the Gaussian component and correlations between estimates at different scales. To illustrate the effect of the estimation process, we apply our procedure to the BOOMERanG data set and show that it is consistent with Gaussianity. This work presents the first estimation of the CMB trispectrum on sub-degree scales.

from homogeneity which have been detected by over two dozen different experiments can be directly related to the primordial origin of perturbations in the early universe and therefore to fundamental physics at very high energies. A new threshold was crossed in the experimental forum with the high resolution, high sensitivity mapping of significant fractions of the CMB sky by the BOOMERanG (de Bernardis et al. 2000) and MAXIMA (Hanany et al. 2000) experiments. A careful analysis of the variance of fluctuations in these maps has led to accurate estimates of the angular power spectrum, far surpassing previous experimental analyses on small angular scales. More recent results from BOOMERanG , Ruhl et al. 2002, MAXIMA (Lee et al. 2001), and from other experiments like DASI (Halverson et al. 2002), CBI (Pearson et al. 2002), VSA (Grainge et al. 2002), ACBAR (Kuo et al. 2002) and ARCHEOPS (Benoit et al. 2002) have posed our knowledge of the CMB angular power spectrum on even more solid ground. However, there is more information in the CMB fluctuations than what is provided by its power spectrum alone. The standard way to extract this information is to analyse high signal to noise maps of the CMB field like those already produced by BOOMERanG or the ones that the MAP † satellite is expected to provide shortly.
There is, therefore, strong motivation to attempt a more detailed study of the CMB sky; in principle one would like a complete characterisation of the probability distribution function of the fluctuations in the CMB with the hope that it might probe more fundamental features of the origin of structure in the universe. For example one relatively stringent test of whether the origin of fluctuations is due to a standard, single-field inflationary model is if the CMB is a realization of a nearly Gaussian distribution (Coles and Barrow 1987), while other models like the curvaton might lead to a much larger non-Gaussian contribution (Lyth & Wands 2002;Bartolo & Liddle 2002;Bernardeau & Uzan 2002;Bernardeau & Uzan 2002). Many secondary anisotropies like the Ostriker-Vishniac (OV) effect (Ostriker & Vishniac 1986) and the Sunyaev-Zeldovich effect (Sunyaev & Zeldovich 1980) can introduce measurable non-Gaussianity while foregrounds and systematics may contribute as well.
A standard method of parametrising an arbitrary probability distribution function is in terms of all its higher order moments. In the case of statistically homogeneous and isotropic fields it is more convenient to consider moments of the Fourier transform of the field; these symmetries will impose a set of selection rules which pick out the true degrees of freedom. The past few years have seen initial attempts at constraining these moments by measuring them with the available data. A series of analyses of the bispectrum (the cubic moment) of the COBE data have revealed the presence of a non-Gaussian systematic (Ferreira, Magueijo, & Gorski 1998;Banday, Zaroubi & Gorski 2000;Komatsu et al. 2001). This has been confirmed with an analysis of the trispectrum, the quartic moment (Kunz et al 2001;Komatsu 2002). But no primordial non-Gaussianity was detected. On smaller angular scales analyses of the QMAP and QMASK (Park et al. 2001) and MAXIMA (Wu et al 2001;Santos et al 2002) data have shown that their observations are consistent with the assumption that the CMB anisotropies are the result of an isotropic Gaussian random process. The analyses of these data sets have also revealed that statistical methods may be sensitive to the data processing pipeline. Moreover, different technical issues must be confronted if one is considering finite sky coverage as opposed to full sky coverage.
The trispectrum, which we consider in this study, probes a rather different kind of non-Gaussianity than the bispectrum; (Aghanim et al 2003) found it to be very sensitive to point sources and (Castro 2003) showed it to be a powerful probe of the OV effect. It can also be used to detect weak lensing in the CMB (Cooray & Kesden 2002). In particular, as described in (Bernardeau 1997), weak lensing does not introduce a three-point correlation function, meaning that its expected bispectrum is zero. The first non trivial higher order correlation function to detect weak lensing effects is the trispectrum. Using the trispectrum to test for non-Gaussianity in high-resolution CMB maps complements therefore other analyses using the bispectrum, like (Santos et al 2002). The purpose of this paper is to improve, extend and apply the method for estimating the trispectrum first proposed in Kunz et al (2001) where it was applied to the full sky map produced by the COBE satellite.
Using techniques developed for performing the operation in pixel space (Ferreira, Magueijo, & Silk 1997;Spergel & Goldberg 1999;Hu 2001) we present a formalism which can easily be extended to high ℓs and discuss a robust method for identifying an orthogonal set of estimators in the full-sky case. We then discuss the various numerical and statistical problems one faces when analysing finite sky coverage. We finally use the BOOMERanG data as a test case to extract the first estimate of the trispectrum on sub-degree angular scales.

FORMALISM AND NOTATION
In this section we present the notation that will be used throughout this work. We start with a temperature anisotropy field defined on the sphere, ∆T T (n); it may be expanded in terms of spherical harmonic functions, Y ℓm (n): For any theory of structure formation, the a ℓm coefficients are a set of random variables; we shall restrict ourselves to theories which are statistically homogeneous and isotropic. The power spectrum C ℓ of the temperature anisotropies is then defined by If we consider the 3-point function of the temperature field, we obtain the bispectrum, defined as The term (· · ·) is a Wigner 3J symbol, which arises due to the "selection rules" of the moments. Following the same steps, we can construct the 4-point function and the associated trispectrum. We represent the rotationally invariant solution for the trispectrum as in (Hu 2001): Using the orthogonality properties of the Wigner 3J symbols and the relation Q = T + G, we can invert the equation (4) to obtain the estimator T ℓ 1 ℓ 2 ℓ 3 ℓ 4 (L) = (2L + 1) a ℓ 1 m 1 a ℓ 2 m 2 a ℓ 3 m 3 a ℓ 4 m 4 −Ĝ l 1 l 2 l 3 l 4 (L) The term G l 1 l 2 l 3 l 4 (L) represents the unconnected Gaussian contribution and it is given in (Hu 2001) as: G l 1 l 2 l 3 l 4 (L) = (−1) (l 1 +l 3 ) (2l 1 + 1)(2l 3 + 1)C l 1 C l 3 δ l 1 l 2 δ l 3 l 4 δ L0 + (2L + 1)C l 1 C l 2 × [(−1) l 2 +l 3 +L δ l 1 l 3 δ l 2 l 4 + δ l 1 l 4 δ l 2 l 3 ] The termT ℓ 1 ℓ 2 ℓ 3 ℓ 4 (L) is the connected part of the angular trispectrum and its expectation value is exactly zero for a Gaussian field. This means that the connected part is sensitive to the presence of non Gaussianities. The unconnected term is non-zero only for L = 0 or l 1 = l 2 = l 3 = l 4 , but only with full sky coverage. In the case of incomplete sky coverage the unconnected terms can contaminate all other modes. We will discuss this situation in section (3).
For the purpose of this work we have not computed the possible trispectrum components. We concentrated only on the simpler case ℓ 1 = ℓ 2 = ℓ 3 = ℓ 4 = ℓ, i.e. the diagonal component. Recent papers have investigated in detail the power of the diagonal trispectrum in the presence of some non-Gaussian signals mentioned in the introduction. In particular, (Aghanim et al 2003) has shown that for simulated point-source maps the diagonal trispectrum is much more powerful than the nearly diagonal estimator (ℓ, ℓ + 1, ℓ + 2, ℓ + 3), even though the latter does not contain a Gaussian contribution. Also, in (Castro 2003) it is discussed how the OV effect generates a signature on the diagonal tispectrum which could easily be detected on the arcminute scales probed by the Planck ‡ mission.
Finally it should be noticed that computing all components of the trispectrum is a serious computational challenge. Many modes are also correlated due to the limited sky coverage. For these reasons we have decided to restrict this analysis to the case of ℓ 1 = ℓ 2 = ℓ 3 = ℓ 4 = ℓ.
We start with the method described in (Kunz et al 2001). We define a ℓm 1 a ℓm 2 a ℓm 3 a ℓm 4 = n a=0 T ℓ;a τ a;ℓ m 1 m 2 m 3 m 4 where the T ℓ,a are the components of the trispectrum that we wish to estimate, and τ is a tensor which we have to determine in order to construct an estimator for T ℓ,a . The geometrical considerations stated above, together with the required symmetries with respect to the interchange of {ℓ i , m i } pairs suggest Although theτ s satisfy all the correct symmetries, they define an over-complete basis. To correct for this deficiency we define an orthonormalised set of tensors where the matrix L aα ℓ is derived from the required property that the τ be orthogonal with respect to the product given in equation (5) of Kunz et al (2001). The estimator of the trispectrum is then given bŷ Note that there are only int(ℓ/3) independent estimators due to the symmetry properties of the a ℓm . In this paper we will consider the "normalised" trispectrum used in Kunz et al (2001) where we divide each estimate of the trispectrum by (Ĉ ℓ ) 2 , whereĈ ℓ = 1 2ℓ+1 |a ℓm | 2 . Its statistical properties are equivalent to the ones of the unnormalised estimator, and it is more robust with respect to fluctuations in the power spectrum (Aghanim et al 2003).

APPLICATION TO HIGH RESOLUTION MAPS WITH INCOMPLETE SKY COVERAGE
In this paper we will be focusing on a high-resolution map with incomplete sky coverage, in particular the BOOMERanG data set. This leads to a set of algorithmic problems which did not have to be addressed in Kunz et al (2001). The three problems we wish to highlight are: Speed: The numerical evaluation of Wigner 3J coefficients for large values of ℓ becomes time consuming and practically unfeasible. Indeed beyond the COBE resolution of a maximum ℓ of approximately 25 it is not possible to estimate theT α;ℓ sufficiently rapidly for a robust Monte Carlo assessment of the statistics. Gaussian Contamination: The finite sky coverage will induce correlations between the estimators with different values of ℓ and a (or α). As a consequence, all estimators may be heavily contaminated by the Gaussian (or disconnected) contributions to the maps.
Correlations: The correlations between modes in the cut sky mean that theT α;ℓ will be even more correlated than in the full sky case.
We shall now focus on the solutions we propose to these three problems

Speed
We have opted to use the method described in (Hu 2001;Spergel & Goldberg 1999) for calculatingT α;ℓ : we define a new set of sky maps weighted in rings centred around a point q: To implement this method we start with the relation (11) and we use the relation (1) to express the temperature T as a function of spherical harmonics and the relation: to express also the Legendre polynomials as a function of spherical harmonics. Combining them with equation (11) we obtain: The e ℓ calculation is quite fast because we can use the fast Fourier transform on rings of equal latitude (Muciaccia et al 1997).
We can then rewrite equation (5) in terms of this new set of sky maps (Komatsu 2002): where If we expand the Wigner 3J symbols in terms of spherical harmonics and use the addition theorem we obtain: where We have thus computed a set ofT α;ℓ which we can orthonormalise to get the estimator for the trispectrumT a;ℓ . This method is very fast, especially when estimating the trispectrum at high values of ℓ. Note that direct evaluation of the Wigner 3J coefficients, e.g. by recurrence relations, would result in an O(ℓ 5 ) problem, requiring ∼ 10 16 operations for ℓ ∼ 1000.

Gaussian Contamination
Kunz et al (2001) found that the purely Gaussian contribution to the trispectrum (the disconnected part) corresponds to the α = 0 term. By orthogonalising all other estimators with respect to this tensor it is possible to remove the Gaussian contribution exactly on a map by map basis. The resulting estimators are only sensitive to non-Gaussian contributions, i.e. in the case of Gaussian skies they would have a zero expectation value. Additionally, the variance of this estimator is shown to be minimal (Kunz et al 2001). In the cut sky case, there are strong cross correlations between components of the trispectrum with different values of ℓ and α. In this case, the orthogonalisation method fails. Given that the Gaussian contribution to the trispectrum may be much larger than the non-Gaussian contribution, it is essential that we remove it as completely as possible nonetheless.
To overcome these problems we have chosen to employ the following Monte Carlo scheme: we generate an ensemble of maps with the same angular power spectrum, sky coverage and noise as the data maps we want to analyse. We estimate theT α;ℓ from each map and calculate the mean of these quantities over the whole ensemble. Let us denote this mean byT α;ℓ G . We then use this quantity to correct for the Gaussian contamination in the estimate of the trispectrum from the data by definingT α;ℓ GC =T α;ℓ −T α;ℓ G . Note that, by so doing, we are removing the Gaussian contamination before performing the full-sky orthogonalisation, i.e. before multiplying by L.

Correlations
The fact that there is only limited sky coverage also implies that there will be correlations between values of the trispectrum at different ℓs. This is shown in section 5 below. To strongly suppress the correlations and to end up with a simple covariance matrix, we consider band averaged values of the trispectrum. Since there is no a priori given band width, we study the cases of ∆ℓ = 40, 50 and 60 and ∆a = 10 and 15. The choice ∆ℓ ≃ 50 is consistent with previous analyses of the BOOMERanG power spectrum (see, e.g., (Ruhl et al. 2002)). In this way we can check the sensitivity of the results to the chosen band size.

NUMERICAL IMPLEMENTATION AND CONSISTENCY TESTS
The process we use is basically the same as in a number of previous analyses (Kogut et al 1996;Ferreira, Magueijo, & Gorski 1998;Kunz et al 2001;Santos et al 2002). We generate an ensemble of Gaussian maps with the same angular power spectrum and noise property as estimated from the BOOMERanG data and the same sky coverage. We then apply our estimators to the set of maps to obtain a distribution for each estimator in the Gaussian case. In particular we characterise the full distribution in terms of the mean values of the estimators and the covariance matrix between them. These quantities are used to define a standard multivariate χ 2 as a goodness of fit. The estimators of the trispectrum are then evaluated from the BOOMERanG data. In section 5 we will discuss their behaviour. The goodness of fit of these estimators is compared against its distribution for a new ensemble of Gaussian maps. From this comparison we can quantify the confidence with which the data can be said to be Gaussian from the point of view of our estimator. It is clear that the numerical details of this process must be well understood if we are to believe in our results. We focus on the particularities of the analysis in this paper which differ from previous analyses.
It is important to compare the results using this hybrid pixel/harmonic analysis with the standard methods which have been used before. We do so by looking at the two lower order statistics, i.e. we have calculated the power spectrum C ℓ and the bispectrum B ℓℓℓ using this new approach as well as summing up the 3J symbols. The relevant expressions are: and (Spergel et al. 1999) ℓ ℓ ℓ 0 0 0 B ℓℓℓ = dqe ℓ (q)e ℓ (q)e ℓ (q) We have compared the C ℓ and B ℓℓℓ using these expressions with the standard results obtained using a ℓm and the Wigner 3J symbols, for a maximum ℓ value of 1000. Using a set of CMB Gaussian maps with the best fit power spectrum measured by Boomerang ) and a pixel resolution of 7 ′ we have found that the bispectrum obtained with e ℓ is affected by a pixelisation effect for high values of ℓ (while the power spectrum shows no difference). To check for this, we have done the same analysis with higher resolution (≃ 3 ′ ) and have found that the pixelisation effect vanishes. Given that we are restricted to the pixelisation level of the data, we can use the comparison of the two estimates of the bispectrum to define a maximum ℓ out to which we can trust the new estimate of the trispectrum. Note that it is computationally intractable to perform such a comparison in terms of the trispectrum, although this would be preferable. From Fig. 1 one can see that discrepancies arise for ℓ > 800 and we chose not to estimate the trispectrum beyond ℓ = 700, leaving a conservative margin as we did not test the trispectrum itself. Furthermore, we chose not to consider any ℓ below 100 because the BOOMERanG data are not very sensitive to these modes, due to limited sky coverage and data filtering.
Another novelty in our analysis (as compared to that of Kunz et al. 2001) is the method Figure 1. In this figure we represent the normalised bispectrum I 3 ℓ = B ℓ /C 3/2 ℓ calculated for a full sky Gaussian map with a pixel size of 7 ′ . The red crosses show I 3 ℓ obtained with the e ℓ method, the black boxes I 3 ℓ with a ℓm and the 3J symbols. We can see that at ℓ ≃ 800 there is a quite evident difference between two plots, due to a pixelisation effect. We limit therefore our analysis to ℓ ≤ 700. for constructing L a;α ℓ . There, a Gram-Schmidt (GS) procedure was used to calculate the orthonormal transformation matrix L a;α ℓ . Due to the inherent instability of the GS procedure, it is not applicable to large matrices, i.e. for large ℓ. We have opted therefore to use an alternative orthonormalisation method. We subtract the a = 0 part and use a Jacobi routine to obtain a spectral decomposition (SD) of the remaining matrix. The eigenvectors of the non-vanishing eigenvalues form then the transformation matrices. This method is robust and, moreover, gives us an unambiguous procedure for ordering the estimators through the different eigenvalues. As a strong consistency test we have applied our trispectrum code to the COBE data and compared the results with those of (Kunz et al 2001). The results of the SD method lie, up to a possible sign change, very close to the original ones (see figure 2). In any case, the statistical significance of the results (and the conclusions one can extract) are the same as in Kunz et al (2001). We advocate the use of the SD method from now on, even in the case of analyses limited to low ℓs.
For our analysis we have used the best four of the six 150 GHz channels of the BOOMERanG 1998 flight; we naively coadd the data taken at the scan speed of 1 degree per second (1 dps). We simulate three sets of 1000 Gaussian maps each. In fact, we need three statistically independent ensembles of simulated maps for our analysis: one to estimate the Gaussian contribution described above, another to estimate our estimator's covariance matrix and a third one to perform the actual comparison with real data. To produce the maps we follow the very same steps used during the real Boomerang data reduction. To generate a map we use timestream simulations created with the actual flight pointing and transient flagging. The signal component of these time streams is generated from simulated gaussian CMB maps, while the noise component is from gaussian realizations of the measured detector noise power spectrum. The noise spectrum is determined using the iterative procedure described in Ferreira & Jaffe (2000) and Netterfield et al. (2002).
To reduce the effects of 1/f noise on this naively binned map, a brick-wall highpass Fourier filter is first applied to the timestream at a frequency of 0.1 Hz. A notch filter is also applied between 8 and 9 Hz to eliminate a non-stationary spectral line in the data. This only affects angular scales well above ℓ = 1000 and is therefore irrelevant for our analysis.
The coaddition of several channels is achieved by averaging the maps (both from the data, and from the Monte-Carlos of each channel). Each channel has a slightly different Figure 2. In this figure we reproduce the COBE results for the trispectrum with the GS and SD orthonormalisation methods (see text) and compare them with (Kunz et al. 2001).The black crosses are the (Kunz et al. 2001) results, green squares are GS method results and blue triangles are SD method results.
beam size, which is taken into account in the generation of the simulated maps. We select the most central region of the scan by applying an elliptical mask as in . This corresponds to ≃ 1.8% of the full sky. The mask selects a region of approximately uniform coverage and high signal to noise and comprises ∼ 57000 pixels of size ∼ 7 ′ each in the Healpix pixelization scheme (Górski et al. 1998). We refer the reader to Ruhl et al. (2002) for a thorough description of a simulation pipeline (based on the MASTER/FASTER algorithms described in (Hivon et al. 2000)) very similar to the one used here.

RESULTS AND APPLICATION TO THE BOOMERANG DATA
To show in detail the method proposed in Sections 2 and 3 we are going to discuss the results obtained at each step from both the data and the Monte-Carlo simulations. We start with the estimate of theT α;ℓ without Gaussian corrections. In the top panel of Figure 3 we plot T α;ℓ as a function of ℓ for selected values of α. We can highlight two features. Firstly thē T α;ℓ are highly correlated for adjacent values of ℓ due to the finite sky coverage, as expected. Secondly, and because of the finite sky coverage, there is a strong contamination from the disconnected component of the trispectrum. This is evident in the fact that the values of T α;ℓ scatter about the (C ℓ ) 2 and that the 95% confidence limits are not centred about zero. As one would expect the lower the value of α, the more contaminated the estimate is by the disconnected part. As advocated in Section 3, we correct for the contamination due to the disconnected component by using a Monte Carlo ensemble (of 1000 realizations) to generate a correction. This can be seen as a bias which must be subtracted off all estimates ofT α;ℓ with α > 0. In the bottom panel of Figure 3 we plot the "Gaussian corrected" estimate ofT α;ℓ with corresponding 95% confidence limits. As expected the estimates now scatter around zero, while the confidence limits, although not necessarily symmetric around the ℓ axis are much more centred. The remaining asymmetry is merely a manifestation that for low αs the distribution of theT α;ℓ is slightly skewed.
Given that we will be working with normalised estimates of the trispectrum (as in Kunz  et al 2001) it is illustrative to plot the ℓ dependence ofT α;ℓ /(Ĉ ℓ ) 2 for a few values of α. We do this in Figure 4. One can see the dependence on α of the 95% confidence intervals, i.e. the ℓ value of minimum scatter depends on α.
Let us now proceed to the orthornomalised estimators,T a;ℓ ; a selection of estimators are plotted for a choice of as in Figure 5. As noted above, the a are limited to a ≤ int(l/3), and we see a clear suppression of the high a values for each ℓ (or, in the case of the figure, of the low ℓ values for fixed a), as the maps with limited sky coverage contain less information than full sky maps.   Once we have calculated the trispectrum estimators both for the BOOMERanG data and for the Monte-Carlo simulations, we can proceed to get the χ 2 distribution for the simulated Gaussian maps and compare it to the data. We used two different approaches: one taking as estimatorT a;ℓ (the orthogonalised trispectrum corrected for Gaussian contamination and normalised to C 2 l ) and the other one taking its absolute value |T a;ℓ |. We construct a standard multivariate χ 2 as: χ 2 = ℓ,ℓ ′ ,a,a ′ ( T l;a G −T l;a )C −1 ℓ,a,ℓ ′ ,a ′ ( T l ′ ;a deriving the expectation values T l;a G and the covariance matrix C ℓ,a,ℓ ′ ,a ′ = T l;aT l ′ ;a ′ G − T l;a G T l ′ ;a we don't sum over all ℓ and a, but bin both a and ℓ, varying the size of the bins. Finally we calculate the χ 2 distribution from the last ensemble of Gaussian maps. In figures 6 and 7 we show the χ 2 distribution derived from 1000 Gaussian realizations compared to the BOOMERanG data for both estimators and for different bin widths. The probability P (χ 2 > χ 2 B ) that a Gaussian map has a larger χ 2 than the Boomerang map is given in table 1. Although the values vary considerable with the choice of bin-widths, none of them is below 5% or 2σ. We conclude that the trispectrum does not detect any non-Gaussianity in the coadded Boomerang 150 GHz maps.

CONCLUSIONS
We have applied an improved version of the method of (Kunz et al 2001) for measuring the trispectrum to the four best 150 GHz Boomerang maps. To this end, we used maps containing only one multipole each to avoid computing the Wigner 3J symbols and subtracted the average Gaussian contribution using an ensemble of simulated maps. We then orthogonalised these maps and normalised them to C ℓ . We then binned the resulting trispectrum values with a variety of different bin sizes, and computed the χ 2 value, using a full covariance matrix estimated from a second ensemble of simulated Gaussian maps. When comparing the data χ 2 value to the Gaussian realizations (obtained from a third ensemble of simulated Gaussian maps) we concluded that the trispectrum does not detect any deviations from Gaussianity.
This work complements the pixel-space analysis (Polenta et al. 2002) and the bispectrum analysis (Contaldi et al. in preparation) of the BOOMERanG data.
In this paper we have studied for the first time the trispectrum of real CMB data with sub-degree resolution. The main problem encountered was the limited sky coverage, which introduces strong correlations, and prevents the use of orthogonalisation to remove the Gaussian (connected) part of the trispectrum. We expect therefore that the MAP and Planck satellites will be able to improve on this analysis considerably, but this study provides a proof of feasibility for measuring the trispectrum of full sky high resolution maps as well as first results on small angular scales.