A Labatie

J L Starck

M Lachièze-Rey

DETECTING BARYON ACOUSTIC OSCILLATIONS

Keywords:

Detecting baryon acoustic oscillations

INTRODUCTION

Large-scale structures in the Universe provide crucial information and can be used to test different cosmological models. This study is complementary to other observations such as the Cosmic Microwave Background (CMB) or Type Ia supernovae. Combining those different observations enables to cross-test models, to break degeneracies, and better constrain cosmological parameters. The good agreement found recently between them and the now-standard Lambda-Cold Dark Matter (ΛCDM) model gives hope for this model to be a lasting foundation.

CDM models with baryons predict the existence of acoustic waves traveling in the hot plasma before recombination, when baryons and photons were coupled together. Those spherical waves originate from the competition between gravitation making over-densities collapse and the photon pressure. About 380 000 years after the Big Bang, baryons and photons decoupled and those spherical waves became frozen at the sound horizon scale r s . Because of their large size (≈ 153 Mpc) their imprints in the matter density field have mainly undertook linear evolution and they should be clearly seen antoine.labatie@cea.fr in current large-scale structures. Those acoustic waves which translate into an excess of correlation at the sound horizon scale are known as Baryon Acoustic Oscillations (BAOs, [START_REF] Peebles | [END_REF]; [START_REF] Bassett | Dark Energy: Observational and Theoretical Approach[END_REF]).

There are two different uses of BAOs that should be distinguished. First they can be used as a very distinct feature to confirm cosmological models. Indeed their detection at the expected scale in large-scale structures gives a strong support for CDM models, with a linear gravitational evolution from z ≈ 1000 and the existence of dark energy. Concretely, the detection is made through hypothesis testing, by finding that models with BAOs are strongly preferred to models without BAOs.

The first convincing detection of BAOs in large-scale structures was reported in the correlation function analysis (Eisenstein et al. 2005) of the Sloan Digital Sky Survey (SDSS, York et al. (2000)) Luminous Red Galaxies (LRG, Eisenstein et al. (2001)) survey. It gave a 3.4σ significance in Data Release 3 (DR3). It was followed by other detections, as in the 2-degree Field Galaxy Redshift Survey (2dFGRS, Colless et al. (2001)), using power spectrum analysis (Cole et al. 2005) with a 2.5σ significance. The power spectrum analysis is also applied in Hütsi (2006) with a 3.3σ detection in the SDSS-LRG DR4. In Percival et al. (2007) and Percival et al. (2010) the combined power spectrum of LRG and Main [START_REF] Strauss | [END_REF] samples of SDSS (with respectively DR5 and DR7), together with the 2dFGRS survey is used to obtain respective significances 3σ and 3.6σ. Very recently BAOs were also detected in the power spectrum [START_REF] Blake | [END_REF] of the WiggleZ Dark Energy Survey (Drinkwater et al. 2010) at a higher redshift z = 0.6 with a significance of 3.2σ.

We must keep in mind that the different detection levels cannot usually be compared. In Eisenstein et al. (2005) and Cole et al. (2005) it is calculated with respect to zero-baryon models (pure CDM models). In Hütsi (2006) and [START_REF] Blake | [END_REF], the significance is calculated with respect to the "no-wiggles" fits of Eisenstein & Hu (1998), which remove the baryon oscillations signature but keep the intermediate suppression of power due to baryons. Finally in Percival et al. (2007) and Percival et al. (2010), it is calculated with respect to power spectrum models where oscillations are smoothed out using splines.

The second use of BAOs consists in constraining cosmological parameters. Again BAOs are very useful, because they provide a statistical standard ruler (Seo & Eisenstein 2003) with an absolute size given with small uncertainty by CMB measurements (Komatsu et al. 2009). So they directly constrain the redshift-distance relation in redshift surveys. Besides, BAOs appear to have the lowest systematic uncertainties among current methods for studying cosmic expansion (Albrecht et al. 2006). They have been used in combination with other cosmological probes to constrain more efficiently cosmological parameters (Eisenstein et al. 2005;Tegmark et al. 2006;Percival et al. 2007;Sánchez et al. 2009;Percival et al. 2010;Reid et al. 2010;Kazin et al. 2010;[START_REF] Blake | [END_REF].

Note that in most studies, the aforementioned constraints not only come from BAOs, but from the whole information in the estimated correlation function or power spectrum. Note also that BAOs do not need to be detected before they can be used for parameter constraints (Cabré & Gaztañaga 2011). One could think that the BAO peak must be proved to be significant and not a random fluctuation, before it is used as a standard ruler. However all sources of uncertainty are normally taken into account when obtaining the constraints (e.g. in the covariance matrix of ξ when using the correlation function), so this argument is not valid. The real question is whether the accepted cosmological models are correct, and BAO detection in a given sample is just a contribution to support these models.

In this paper we focus on the first use of BAOs, i.e. on the BAO detection. We restrict the analysis to the correlation function, but most reasonings could also be applied to the power spectrum.

The plan of this paper is as follows: we start in section 2 by discussing the correlation function estimation and modeling, as well as the general procedure for BAO detection. In section 3 we present wavelet methods for BAO detection, which are mildly model-dependent. In the rest of the paper we focus on fully model-dependent methods. In section 4 we present the classical method used for BAO detection, based on the χ 2 statistic. We find that this method does not provide the correct sig-nificance. So we propose in section 5 a new method that we call the ∆l method, based on two modifications to the χ 2 method. In section 6 we explain the other use of BAOs, i.e. how parameter constraints can be obtained. Finally we illustrate the different methods and results using simulations in section 7.

BAO DETECTION IN THE CORRELATION FUNCTION

Correlation function

The correlation function ξ is a second order statistic that measures the clustering of a continuous field or point process. More precisely for the distribution of galaxies, it quantifies the excess of probability to find a pair of galaxies in volumes dV 1 and dV 2 separated by r, compared to a random unclustered distribution

dP 12 = n [1 + ξ(r)] dV 1 dV 2 (1)
with n the mean density of the point distribution. With the isotropy hypothesis, ξ only depends on the norm r = r of the separation vector r. However in redshift space the field is not rigorously isotropic. In this case, one is usually interested in the monopole ξ(r) of ξ(r). In the following we will make the abuse of language of referring to the monopole when speaking about the correlation function.

Given a galaxy survey, the correlation function can be estimated by comparing the number of pairs at distance r with a random catalogue in the same volume. Different estimators based on this method have been proposed and empirically compared [START_REF] Labatie | Statistical Methodology[END_REF]Kerscher et al. 2000;Pons-Bordería et al. 1999)). While Pons-Bordería et al. (1999) did not recommend one estimator for all cases, Kerscher et al. (2000) recommend to use the Landy-Szalay estimator. The recommendation is the same in the more recent study [START_REF] Labatie | Statistical Methodology[END_REF], where Landy-Szalay is found to be nearly unbiased for current galaxy surveys. It is given by

ξ(r) = 1 + N RR N DD DD(r) RR(r) -2 N RR N DR DR(r) RR(r) (2)
with DD(r), RR(r), and DR(r) the number of pairs at a distance in [r±dr/2] of respectively data-data, randomrandom and data-random points, and N DD , N RR and N DR the total number of corresponding pairs in the catalogues.

Modeling the galaxy correlation function

In the early universe prior to recombination, baryons are tightly coupled with photons. This results in acoustic waves traveling in the plasma (BAOs), but also in the suppression of power on small and intermediate scales compared to a CDM model without baryons. After the time of decoupling (also called drag epoch), the matter density field becomes pressureless, allowing the perturbations to grow by gravitational instability. This evolution can be analytically solved in the linear regime where fluctuations are small, and only the overall amplitude of the power spectrum is changing. Eisenstein & Hu (1998) provides fitting formulae for the linear power spectrum, with a dependence on cosmological parameters. From the linear power spectrum, the linear correlation function is simply obtained by Fourier transform. The effect of BAOs is clearly identified as series of wiggles in the power spectrum, and as a localized bump in the correlation function at the sound horizon scale r s (see figure 2).

To fully model the matter correlation function, one also has to take into account non-linear evolution. This can be done using N -body simulations, which empirically provide a correction from the linear to the non-linear correlation functions. Smith et al. (2003) provide corrections for scale-free power spectrum. So one also has to correct for non-linear degradation of the acoustic peak. Eisenstein et al. (2007) found that this is well approximated by a Gaussian smoothing of the acoustic feature in real and redshift spaces.

A last step to model the galaxy correlation function is to take into account redshift distortions and galaxy bias with respect to matter. Again this can be done using Nbody simulations, where dark matter halos are populated using a halo model.

Models of galaxy correlation with BAOs are constructed using the linear matter correlation function with non-zero baryon fraction Ω b > 0, and further applying the different corrections. On the other hand, models of galaxy correlation without BAOs can be obtained by setting the baryon fraction to zero Ω b = 0. One can also construct no-BAO models with Ω b > 0, by using only the non-oscillatory part of the power spectrum to remove the effect of BAOs (e.g. the no-wiggles fit of Eisenstein & Hu (1998)). In this case, BAOs are artificially erased and the models are non-physical. Yet they enable to test the existence of BAOs independently of other baryonic effects.

BAO detection by hypothesis testing

Let us define the two different hypotheses H 0 : no-BAO hypothesis H 1 : BAO hypothesis BAO detection is equivalent to this problem of hypothesis testing. The common procedure is to design a test statistic to assess the truth of the null hypothesis H 0 . From the test statistic obtained with the measurement, one computes a p-value and a significance. If the measurement if found to be more unlikely than a given threshold under H 0 , one rejects H 0 and accepts H 1 (see section 4.3).

We focus on the case where the data measurement is the correlation function estimated in different bins ξ = ( ξi ) 1≤i≤n . Such a binning is always present for the measurement, and thus for the model correlations in the hypotheses H 0 and H 1 . As a slight abuse of language, we will use the terms estimated correlation function and model correlation functions for designating their binned versions.

The hypotheses H 0 and H 1 are based on BAO and no-BAO models of correlation function, as the ones presented section 2.2. The hypotheses also include the noise of the measurement, i.e. the covariance matrix C = (C i,j ) 1≤i,j≤n of ξ.

We will see in section 3 that wavelet methods are mainly sensitive to the BAO feature in the correlation function, and not on the global shape of the model.

On the other hand, usual detection methods (e.g. the χ 2 method) are based on full modeling of the correlation function. In this case, BAO and no-BAO models of correlation function ξ BAO,θ and ξ noBAO,θ are parameterized by θ to account for variations of cosmological parameters.

The hypotheses are

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C noBAO,θ ) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C BAO,θ )
Most methods work with these hypotheses, where ξ is Gaussian. Ideally the hypotheses are sampled by Nbody simulations for each model, which does not force ξ to be Gaussian. Yet we will see in section 7.3 that the Gaussian approximation works very well for our lognormal simulations.

The classical χ 2 method used for BAO detection that we present in section 4 simplifies the hypotheses by assuming constant covariance matrices ( i.e. independent of the model). The reason is that it can be hard to evaluate the covariance matrix for all tested models. In a modified version of the χ 2 method that we propose in section 5 this simplification is not imposed.

The distribution of ξ is entirely specified for each hypothesis when fixing θ, which is not the case when allowing θ to vary. In the first case the hypotheses are said to be simple, whereas in the second case they are said to be composite1 .

BAO detection makes more sense when allowing variations of θ, since it takes into account uncertainties in cosmological parameters. However we will see in section 4 with the classical χ 2 method, that one has to be careful when testing composite hypotheses.

WAVELET FILTERING METHODS

As explained in section 2.2, BAOs manifest as a characteristic peak in the correlation function at the acoustic scale r s . For detecting this feature new methods have recently emerged, based on wavelet analysis (Xu et al. 2010;Tian et al. 2011;Arnalte-Mur et al. 2011). Wavelet transforms are widely used in many areas, especially in image analysis [START_REF] Mallat | A Wavelet Tour of Signal Processing, The Sparse Way[END_REF][START_REF] Starck | Sparse Image and Signal Processing[END_REF]. They are specially suited for the analysis of data at different scales, and identification of characteristic patterns or structures.

Here the characteristic structure is the BAO feature in the correlation, with uncertainty in its scale and shape. Uncertainty in the scale comes from a wrong fiducial cosmology to convert redshifts into distances, and a weak dependence of r s on cosmological parameters. Uncertainty in the shape is due to non-linear evolution, redshift distortions which are subject to modeling errors (see section 2.2).

Our focus here is BAO detection, so we present two different methods that have been developed for this purpose. In these methods, a wavelet w = (w(r i )) 1≤i≤n acts as a peak finder to detect excess in the measured correlation ( ξ(r i )) 1≤i≤n . The wavelet w(R, s) is parametrized by two parameters R and s, linked respectively to the scale and width of the peak. A simple way to do it is to consider an original peak-finding wavelet w 0 and rescale it as w i (R, s) = w 0 ri-R s (see figure 1). One obtains a filtered signal S w (R, s) for every pair (R, s)

S w (R, s) = w(R, s), ξ = n i=1 w i (R, s) ξ(r i ) (3)
The next step is to divide S w (R, s) by its noise σ w (R, s) under H 0 to obtain Z-scores Z w (R, s). The parameters (R max , s max ) giving the maximum Z-score can be used to estimate the BAO scale and width.

Z w (R, s) = S w (R, s) σ w (R, s) (4)
Z max w = Z w (R max , s max ) (5)
The two hypotheses are roughly that the maximum response Z max w is negligible under H 0 (no peak is found), and that there is a non-negligible signal under H 1 (a peak is found). To reject H 0 one performs simulations without BAOs, and computes how rarely a value of Z max w as high as in the data is observed. This gives a p-value and thus the significance of the detection.

The major advantage of using wavelets is that they are mainly sensitive to the BAO feature, and not to smooth changes in the correlation function. Furthermore, scale variations of the analyzing wavelets allows dilations of the correlation function (to probe different cosmologies), and variations of its width allows for variations in the shape of the BAO peak. The consequence is that the wavelet response is mainly related to the existence of a BAO peak, and mildly dependent on the whole modeling of the correlation function. In other words, wavelet methods are robust to small modeling errors in the correlation function.

The price to pay is that they are outperformed by some model-dependent methods, when there are no modeling errors in H 0 and H 1 . As we will see in section 5, the ∆χ 2 (respectively ∆l) statistic can be seen as a generalized likelihood ratio between H 0 and H 1 in the case of constant (respectively varying) covariance matrix. The interest of the likelihood ratio is that it is optimal in the Neyman-Pearson sense for simple hypotheses (see section 4.3 and appendix C). For composite hypotheses there is no such notion of optimality. But we also find in section 7 that the generalized likelihood ratio ∆l gives better results than ∆χ 2 in the case of varying covariance matrix. Since wavelet statistics are not directly linked to the likelihood ratio, we expect less good results than with a generalized likelihood ratio.

A wavelet detection method was used in Tian et al. (2011) on the SDSS Main sample. An angular average of the 3D anisotropic correlation function is computed, by applying a flat weighting in all directions, instead of the angular average on the sphere. The justification is that the BAO feature is sharpened in the line-of-sight direction, and this type of weighting is giving more importance to it. The correlation is further analyzed with a mexican hat wavelet (see figure 1). Using Gaussian simulations with a no-wiggles power spectrum, they find that only 0.2% of the simulations have the statistic Z max w at the same level as the data, which means a 3.1σ detection. In Arnalte-Mur et al. (2011) two different galaxy samples of the SDSS DR7 are used, to compute the mean density profile of the Main sample around LRGs. This signal is analyzed with a 3D wavelet called BAOlet (see figure 1). Because the wavelet is isotropic, it is equivalent to applying a 1D wavelet transform on the cross-correlation LRG-Main. Simulations of H 0 are made by replacing LRGs by random centers to show that LRGs are located at special positions. Using this hypothesis, a signal as high as Z max w is found with a probability p = 4.10 -5 in simulations, corresponding to a 4.1σ detection. Again this cannot be compared to other existing methods, because the tested hypothesis H 0 is very different.

χ 2 METHOD

The χ 2 method is the classical method used for BAO detection and can deal with the general case of varying cosmological parameters. Unlike wavelets methods presented in section 3 it is fully model-dependent, so it is mainly useful when all effects in the correlation function are well understood.

The χ 2 method is also designed for hypotheses where the measurement ξ is Gaussian. In the rest of the paper we will only consider such hypotheses. In section 7.3, we will see using simulations that the Gaussian approximation is well justified.

χ 2 = ξ -ξ m , C -1 ( ξ -ξ m ) (6) = 1≤i,j≤n ξ(r i ) -ξ m (r i ) C -1 i,j ξ(r j ) -ξ m (r j ) (7)
Supposing that the model is correct (i.e. ξ ∼ N (ξ m , C)), the χ 2 statistic follows a χ 2 n distribution, i.e. a chi-square distribution with n degrees of freedom. The χ 2 n distribution can be interpreted as the one followed by the sum of the squares of n independent standard normal variables. If X 1 , . . . , X n are n such i.i.d. Gaussian variables then n i=1 X 2 i follows a χ 2 n distribution.

4.2. χ 2 method for BAO detection Let us show how the BAO detection is usually performed. We note θ = (θ 1 , . . . , θ k ) ∈ Θ the dependence parameters for the model correlation functions with and without BAOs, ξ BAO,θ and ξ noBAO,θ . The hypotheses are given by

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C)
As mentioned in section 2.3, the χ 2 method tests hypotheses with a constant covariance matrix C. The parameters θ are not directly cosmological parameters but are linked to them. For example in Eisenstein et al. (2005), they are given by a dilation parameter α (to account for a wrong fiducial cosmology to convert redshifts into distances), an amplitude parameter b 2 (to account for galaxy bias, redshift distortions, and power spectrum normalization σ 8 ), and the parameter Ω m h 2 (determining the horizon scale at matter-radiation equality, the amplitude of the BAO peak, and more moderately the position of the peak). Other parameters also have an impact on the expected correlation function (Ω b h 2 and the spectral index n) but they are well constrained by CMB data and can be fixed as a good approximation.

The χ 2 has a dependence on the parameters θ

χ 2 BAO,θ = ξ -ξ BAO,θ , C -1 ( ξ -ξ BAO,θ ) (8) χ 2 noBAO,θ = ξ -ξ noBAO,θ , C -1 ( ξ -ξ noBAO,θ ) (9)
For each class of models one can look at the χ 2 bestfits, min θ χ 2 noBAO,θ and min θ χ 2 BAO,θ . It is a widely used result that the best-fit χ 2 value follows a χ 2 n-k distribution, assuming that the true model is inside the fitting class. So the number of degrees of freedom decreases by the number of parameters in the fit. We stress that this result (see appendix A) is only rigorously valid when the space of model correlations is affine. We recall that the measurement vector ξ and the model correlations ξ m (θ) are n-dimensional binned versions of their continuous counterparts. So the set of all model correlations (ξ m (θ)) θ∈Θ constitutes a subspace of R n , which needs to be a k-dimensional affine space for the previous result to hold (see appendix A).

This result can be used on min θ χ 2 BAO,θ to verify that data are compatible with H 1 . More precisely it can be tested whether min θ χ 2 BAO,θ is compatible with its distribution when the true model is in

H 1 min θ χ 2 BAO,θ ∼ χ 2 n-k (10)
For the rejection of H 0 , the usual procedure is more complex. We add an artificial parameter in the fit, which accounts for the presence of BAOs in the model correlation. For example this parameter β can be a weighting of ξ BAO,θ and ξ noBAO,θ in the model correlation function.

ξ β,θ = β ξ BAO,θ + (1 -β) ξ noBAO,θ (11)
Under H 0 , the expected correlation function is of the form ξ noBAO,θ for the true parameters θ = θ 0 , but it is also of the form ξ β,θ with β = 0 and θ = θ 0 . Thus the best-fit χ 2 value in the no-BAO class follows a χ 2 n-k distribution, and the best-fit χ 2 value in the extended class follows a χ 2 n-(k+1) distribution (since β is an additional parameter).

We are in the case of two nested classes of models, which both contain the true model. In this case, the difference of the best-fit values follows a chi-square distribution with number of freedom equal to the difference of parameters between the two classes. Again this result is not rigorously true in the general case, but only when the two spaces of model correlations are affine (see appendix B).

Here there is only one additional parameter, β, in the extended class, thus the difference of the best-fit values follows a χ 2 1 distribution.

∆χ 2 global = min θ χ 2 noBAO,θ -min β,θ χ 2 β,θ ∼ χ 2 1 (12)
This accounts for the fact that fitting an additional parameter, which is not required by the true model, only moderately decreases the best-fit value. To reject H 0 , one can look at this difference ∆χ 2 global and compute how unlikely it is under H 0 (i.e for a χ 2 1 distribution). In practice, the best-fit in the whole extended class is replaced by the best-fit in the BAO class (i.e. restricting to β = 1)

∆χ 2 = min θ χ 2 noBAO,θ -min θ χ 2 BAO,θ (13)
This difference ∆χ 2 is necessarily less than ∆χ 2 global of equation ( 12), which follows a χ 2 1 distribution under H 0 . Thus for a realization with an arbitrary value ∆χ 2 = x, we have

P (∆χ 2 > x | H0) ≤ P (∆χ 2 global > x | H0) (14)
A χ 2 1 distribution is simply the distribution followed by the square of a standard normal variable. Noting Φ the cumulative distribution function of a standard Gaussian variable, we get for x ≥ 0

P (∆χ 2 > x | H0) ≤ P (∆χ 2 global > x | H0) = 2Φ(- √ x) (15)
which corresponds to a number of σ equal to √ x for a normal distribution. Thus, when ∆χ 2 ≥ 0, one can evaluate the significance of the BAO detection as ∆χ 2 .σ. A significance is originally given in terms of a p-value, i.e. the probability of obtaining the measurement value under H 0 . When given as a number of σ, it is simply the corresponding number of standard deviation to the mean for a Gaussian variable.

In Eisenstein et al. (2005) the difference of chi-square equals ∆χ 2 = 11.7, corresponding to a 3.4σ detection using this method. In Percival et al. (2010) it is equal to ∆χ 2 = 13.1, corresponding to a 3.6σ detection.

Because of the inequality in equation ( 14) the method may seem conservative. However because the assumptions of the method are not verified, we will see in section 7.4.1 that the method actually overestimates the significance.

Limitations of the χ 2 method

The Neyman-Pearson lemma states that, when performing a hypothesis test between H 0 and H 1 , the most powerful tests are based on the likelihood ratio

Λ( ξ) = L H0 ( ξ) L H1 ( ξ) (16)
More precisely, the most powerful test of significance α is

• if Λ( ξ) ≤ η then accept H 1 (i.e. reject H 0 ) • if Λ( ξ) > η then accept H 0 (i.e. reject H 1 )
with α the probability of rejecting H 0 if it is true (type I error)

α = P Λ( ξ) ≤ η | H 0 (17)
The power of the test is defined as the probability of accepting H 1 if it is true. It is equal to 1 -β where β is the probability of type II error, i.e. the probability of not accepting H 1 if it is true.

β = P Λ( ξ) > η | H 1 (18)
In practice, such tests with given thresholds are not really used, and it is more common to cite the significance level for the realization value. For a realization with an arbitrary value Λ( ξ) = x, the significance (given as a p-value) is

α(x) = P Λ( ξ) ≤ x | H 0 (19)
As we show in appendix C, the Neyman-Pearson lemma implies that the expected significance under H 1 obtained with Λ( ξ) is better than with any other statistic. More precisely, under H 1 , the expected p-value of equation ( 19) is smaller, and the expected number of σ is larger for Λ( ξ) than for any other statistic.

Note that the statistic Λ( ξ) is optimal in this sense, but we need to know its distribution under H 0 (to compute the significance α(x) corresponding to a realization value x). Moreover in the case of composite hypotheses, the distribution of Λ( ξ) is not well-defined under H 0 (see section 2.3) so the significance is also not well-defined. The advantage of ∆χ 2 global is that its distribution is identical for every model in H 0 (a χ 2 1 distribution). In this case one is able to give a significance even with composite hypotheses. Yet this result is subject to a regularity assumption, that spaces of model correlation functions are affine. Because it is not verified, we will see with simulations in section 7.4.1 that the distribution of ∆χ 2 global is quite different from a χ 2 1 and that the χ 2 method only gives a rough estimate of the significance.

The estimate can be even more wrong when considering more realistic hypotheses, where covariance matrices depend on the model

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C noBAO,θ ) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C BAO,θ )
Let us consider a common case where the parameters θ are given by θ = (Ω m h 2 , α, b) as introduced in section 4.2. For illustrative purposes, we only take into account the dependence in the bias b of the covariance matrix. It gives a multiplicative factor b 2 in the expected correlation function, and at first approximation (up to shot noise) a multiplicative factor b 4 in the covariance matrix. So the covariance matrices are given by

C BAO,θ = C noBAO,θ ∝ b 4 C (20)
In the classical χ 2 method, χ 2 BAO,θ and χ 2 noBAO,θ are computed with a constant covariance matrix C. So for a realization A ξ with A ≥ 0, we get the BAO best-fit

min θ χ 2 BAO,θ (A ξ) = min θ A ξ -ξ BAO,θ , C -1 (A ξ -ξ BAO,θ ) = A 2 min θ ξ - 1 A ξ BAO,θ , C -1 ξ - 1 A ξ BAO,θ = A 2 min θ χ 2 BAO,θ ( ξ) (21)
The last equality comes from the role of b, which enables any positive multiplication of the model. The same reasoning can be applied to min θ χ 2 noBAO,θ which gets multiplied by A 2 , and thus the statistic ∆χ 2 also gets multiplied by A 2 .

Given the hypotheses with varying covariance matrix, H 0 realizations with θ = (Ω m h 2 , α, b 1 ) can be obtained from H 0 realizations with θ = (Ω m h 2 , α, b 0 ) by a multiplicative factor (b 1 /b 0 ) 2 . As a result, the distribution of ∆χ 2 gets dilated by a factor (b 1 /b 0 ) 4 . This creates very large differences between the different models in H 0 , so the classical χ 2 method provides very bad estimates of the significance. The conclusion is that the classical χ 2 method cannot be used in the case of varying covariance matrix.

MODIFIED χ 2 METHOD

In this section we propose two modifications to the χ 2 methods to overcome its limitations. A first modification enables to obtain the correct significance in all cases. So unlike the classical χ 2 method, our modified method can be applied to hypotheses with varying covariance matrices

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C noBAO,θ ) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C BAO,θ )
The way to rigorously compute the p-value for a measurement ∆χ 2 = x is to consider the "worst-case"

H 0 model p(x) = max θ∈Θ P (∆χ 2 ≥ x | H 0 , θ) (22)
So for every model in H 0 , the p-value of the measurement ∆χ 2 = x is less than p(x). When considering the significance s(x)σ corresponding to this p-value, we get that every model in H 0 is at least rejected at s(x)σ. Note that this is the best significance we can get to reject all H 0 models simultaneously. This way of obtaining the significance does not rely on any assumption, unlike in the classical χ 2 method. However it requires more work to compute the distribution of ∆χ 2 under every H 0 model. We will see how precisely this can be achieved in section 7.4.2, when applying the procedure on simulations.

The second modification we propose is on the statistic itself. As we saw in section 4.3 the optimal statistic to test simple hypotheses H 0 and H 1 is the likelihood ratio. However when working with composite hypotheses, likelihoods are not well defined. Indeed they can be defined for any model in H 0 or H 1 but not for H 0 and H 1 themselves. The idea of the ∆χ 2 statistic is that it can be thought as a generalized likelihood ratio between composite hypotheses. Indeed in the special case of a constant covariance matrix we have

∆χ 2 = min θ χ 2 noBAO,θ -min θ χ 2 BAO,θ = -2 max θ ln (L noBAO,θ ) -max θ ln (L BAO,θ ) = -2 ln max θ L noBAO,θ max θ L BAO,θ (23)
where we used the relation between χ 2 and the likelihood of equation ( 32). This is only valid for a constant covariance matrix. To extend this idea in the case of varying covariance matrices, we simply consider the statistic

∆l = -2 max θ ln (L noBAO,θ ) -max θ ln (L BAO,θ ) (24)
We use the notation ∆l because it refers to a difference of log-likelihoods. It is a slight abuse of notation since it is not strictly speaking a difference of log-likelihoods. Unlike ∆χ 2 , ∆l is still a generalized likelihood ratio for varying covariance matrices. So it should give better results in this case as we verify in section 7.5.

In the case of varying covariance matrices the likelihoods write

L BAO,θ ∝ |C BAO,θ | -1/2 e - χ 2 BAO,θ 2 (25) L noBAO,θ ∝ |C noBAO,θ | -1/2 e - χ 2 noBAO,θ 2 (26)
where χ 2 BAO,θ and χ 2 noBAO,θ are computed by taking into account variations of the covariance matrix

χ 2 BAO,θ = ξ -ξ BAO,θ , C -1 BAO,θ ( ξ -ξ BAO,θ ) (27) χ 2 noBAO,θ = ξ -ξ noBAO,θ , C -1 noBAO,θ ( ξ -ξ noBAO,θ ) (28)
Let us write l noBAO,θ and l BAO,θ for -2 ln (L noBAO,θ ) and -2 ln (L BAO,θ ) we get ∆l = min

θ l noBAO,θ -min θ l BAO,θ (29)
l BAO,θ = χ 2 BAO,θ + ln |C BAO,θ | + cst (30) l noBAO,θ = χ 2 noBAO,θ + ln |C noBAO,θ | + cst (31)
with the same additive constant for l BAO,θ and l noBAO,θ , which can be taken as 0.

Note that ∆l is not equivalent to ∆χ 2 even if χ 2 BAO,θ and χ 2 noBAO,θ are computed using equations ( 27) and (28). Indeed in the case of varying covariance matrices, one has to take into account variations of the matrix determinant as in equations ( 30) and ( 31).

We will refer to the ∆l method for the χ 2 method modified as we suggested: replacing ∆χ 2 by ∆l and modifying the procedure to obtain the correct significance.

COSMOLOGICAL PARAMETERS CONSTRAINTS

Let us describe the second use of BAOs, which can help constrain cosmological parameters. Here the true cosmological model is supposed to be in H 1 .

Constraints with constant covariance matrix

In the case of a constant covariance matrix the hypotheses H 1 is given by

H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C)
For a given model in H 1 , the measurement ξ is Gaussian and χ 2 BAO,θ is equivalent to the log-likelihood

χ 2 BAO,θ = -n ln (2π|C|) -2 ln L BAO,θ ( ξ) (32)
To define the posterior probability p(θ | ξ) one needs a prior p(θ) on θ

p(θ | ξ) ∝ p(θ) L BAO,θ ( ξ) (33) ∝ p(θ) exp - 1 2 χ 2 BAO,θ (34)
To obtain constraints only coming from ξ one can assume a constant prior p(θ)

p(θ | ξ) ∝ L BAO,θ ( ξ) (35)
Note that this choice is still arbitrary because it is linked to a given parameterization. Indeed a constant prior p(θ) can lead to a non-constant prior for a different parameterization.

To combine constraints from ξ with the ones from other independent experiments, one has to modify the prior. For example with CMB data the posterior is given by

p(θ | CMB, ξ) ∝ p(θ, CMB, ξ) ∝ p(θ, CMB) p( ξ | θ, CMB) ∝ p(θ | CMB) L BAO,θ ( ξ) ( 36
)
where we used the independence of ξ and CMB measurement. Adding the CMB measurement is thus equivalent to using a prior p(θ) = p(θ | CMB). Again we consider the parameters θ given by θ = (Ω m h 2 , α, b). The parameter α accounts for dilation of the correlation function and b corresponds to a multiplicative factor b 2 . The correlation function models are thus given by

ξ BAO,θ (r) = b 2 ξ BAO,Ωmh 2 (α r) (37)
α is linked to the dilation scale D V (z) at the mean redshift of the sample z by the relation α = D V (z)/D V,f id (z), with D V,f id (z) the dilation scale for the fiducial cosmology used to construct the 3D data catalogue. The dilation scale expresses how distances dilate when modifying the fiducial cosmology (Eisenstein et al. 2005). It depends on the Hubble parameter (lineof-sight dilation) and the transverse comoving distance D M (z) (transverse dilation)

D V (z) = D M (z) 2 cz H(z) 1/3 (38)
One is interested in constraining Ω m h 2 and D V (z). We consider a constant prior p(θ) to obtain constraints only from ξ. The posterior p(Ω m h 2 , α | ξ) is obtained by marginalizing over the multiplicative factor

B = b 2 p(Ωmh 2 , α | ξ) ∝ L BAO,(Ωmh 2 ,α,B) dB (39) ∝ exp - 1 2 χ 2 BAO,(Ωmh 2 ,α,B) dB (40)
The posterior of Ω m h 2 is obtained by marginalizing p(Ω m h 2 , α | ξ) over α, and the posterior of α by marginalizing over Ω m h 2 . For each parameter, the maximum in the posterior gives the parameter estimate and the standard deviation can give a 1σ interval.

Constraints with varying covariance matrix

In the case of varying covariance matrix the hypothesis H 1 is given by

H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C BAO,θ )
In this case one has to take into account the dependence of the covariance matrix on the model

χ 2 BAO,θ = ξ -ξ BAO,θ , C -1 BAO,θ ( ξ -ξ BAO,θ ) (41) L BAO,θ ∝ |C BAO,θ | -1/2 e - χ 2 BAO,θ 2 (42)
Let us consider the simple dependence C BAO,θ ∝ b 4 C as an illustration again. In this case, the marginalization over b 2 gives a different result compared to the result with constant covariance. So the obtained posteriors of Ω m h 2 and α are also different. We will see with simulations in section 7.6 that this changes indeed the constraints.

TESTS ON SIMULATIONS

Simulations

We use the same procedure for generating lognormal simulations of the SDSS DR7 LRG sample as in [START_REF] Labatie | Statistical Methodology[END_REF] with only small differences in the input power spectrum. Because we consider a volume-limited LRG sample with only the northern contiguous region, the volume is approximately half of the the full LRG sample and the number of galaxies a third. In particular the expected detection significance is lower than for the full LRG sample. So our focus is not on the expected detection significance for current surveys, but rather on comparing the different methods.

We use a ΛCDM power spectrum given by the iCosmo software [START_REF] Refregier | [END_REF], with parameters h = 0.7, Ω b = 0.045, Ω m = 0.27, Ω Λ = 0.73, n s = 1.0, σ 8 = 0.8, and taken at redshift z = 0.3. The transfer function is the form with wiggles of Eisenstein & Hu (1998) and the non-linear correction to the power spectrum is obtained using the fitting formula of Smith et al. (2003). We model the non-linear degradation of the acoustic peak by multiplying the part of the power spectrum containing the oscillations by the function exp(-a 2 k 2 /2) with a = 7h -1 Mpc (i.e. smoothing the oscillations in the correlation with a Gaussian of width a). This is found to be a good approximation in Eisenstein et al. (2007). It consists more precisely in constructing the power spectrum using the forms with and without wiggles

P (k) = Pno wig (k)+exp(-a 2 k 2 /2)[Pwig(k)-Pno wig (k)] (43)
We further apply a constant bias b 2 with b = 2.5 to this power spectrum so that the corresponding correlation function matches the one estimated on the SDSS LRG sample. Our simulations do not take into account the scale-dependence of the galaxy bias. However this is not a problem here since we only use simulations, and do not analyze real data.

We estimated the correlation function on 2000 independent simulations using each time the Landy-Szalay estimator and 100.000 random points. With this procedure Landy-Szalay has been shown to be the estimator with minimum variance, and to be nearly unbiased (see [START_REF] Labatie | Statistical Methodology[END_REF]).

We estimated the covariance matrix by the empirical covariance matrix of the measured correlation function in the simulations. We use bins of size 10 h -1 Mpc and perform the analysis in the range 20 to 200 h -1 Mpc. In this way we obtain n = 18 bins, corresponding to 171 free parameters in the covariance matrix. This is much smaller than the number of simulations so that the empirical covariance matrix gives a good estimate of the true covariance matrix (see e.g. Pope & Szapudi (2008)). Another reason for not using a strong binning, is that most methods use the inverse of the covariance matrix. Since the bins are very correlated at small separation, the inverse matrix is very oscillating for too small binning, and the result becomes non robust to small modeling errors.

We find a small statistical bias in the correlation function estimation, due to the limited resolution of the lognormal simulations or the integral constraint (see [START_REF] Labatie | Statistical Methodology[END_REF]). This bias is negligible in absolute value (≈ 5 × 10 -4 ), but when multiplied by C -1 as in Z w , χ 2 BAO,θ or χ 2 noBAO,θ it can slightly affect the result. Here the lognormal simulations will only be used for computing the covariance matrix and verifying the Gaussianity of ξ. However when working with real data, one should verify that the integral constraint is not biasing results.

Models

Unless otherwise stated, we use for the BAO and no-BAO models in H 1 and H 0 the transfer function of Eisenstein & Hu (1998) with respectively the form with wiggles and without wiggles. At some points we also quote the results obtained with no-BAO models constructed with zero baryon.

The procedure for generating the power spectrums is the same as for the lognormal simulations. Obviously the non-linear degradation of the BAO peak is only per-formed for BAO models, by smoothing with a kernel of size a in comoving coordinates.

The respective correlation functions ξ BAO,θ and ξ noBAO,θ are obtained from the power spectrums by inverse Fourier transform in 3 dimensions, i.e. by a Hankel transform in the isotropic case. The power spectrum has bins with exponential sizes in k (i.e. the ln(k i )'s are spaced linearly) since it is smooth in that space. For doing the Hankel transform with this spacing we use the FFTLog code 2 . The correlation is finally binned equivalently as when it is estimated by pair counting, i.e. for a bin [r i -dr/2, r i + dr/2] ξ(r i ) = ri+dr/2 ri-dr/2 ξ(r) r 2 dr ri+dr/2 ri-dr/2 r 2 dr (44)

We fix the parameters Ω b h 2 = 0.0315, n s = 1.0 and σ 8 = 0.8 as in the lognormal simulations, and choose the same parametrization as before θ = (Ω m h 2 , α, b) for the ΛCDM correlation functions

ξ BAO,θ (r) = b 2 ξ BAO,Ωmh 2 (α r) ξ noBAO,θ (r) = b 2 ξ noBAO,Ωmh 2 (α r)
To obtain these functions we vary Ω m h 2 , adjust the value Ω Λ = 1 -Ω m , and perform the dilation in α of the correlation function obtained in comoving coordinates. There are different choices for varying Ω m h 2 , and we choose to keep h = 0.7 constant and vary Ω m . Another choice would be to keep Ω m = 0.27 and vary h. These choices lead to different amplitudes of the correlation due to different growth factors at redshift z = 0.3. We verified that our results are only slightly affected by this choice.

Finally the functions are rescaled by a factor b 2 for modeling the galaxy bias. We have to take into account that no-BAO models with zero baryon have different amplitude than the other models for the range of scales considered (for the same σ 8 ). Therefore we rescale them by an amplitude factor of 1.29, which is found to minimize the distance between the no-BAO and BAO models for the parameter values of the lognormal simulations

ξ BAO -ξ noBAO , C -1 (ξ BAO -ξ noBAO ) (45)
The lognormal simulations correspond to parameters Ω m h 2 = 0.1323, α = 1, and b = 2.5. In figure 2 we plot the correlation function of the simulations, for the two corresponding no-BAO models and for different values of Ω m h 2 . The no-wiggles model has only the BAO peak smoothed out, whereas the zero-baryon model (with Ω b = 0) has a different global shape. Among BAO models, Ω m h 2 controls the proportion of baryons in the total matter. When Ω m h 2 increases, the proportion of baryons decreases, and baryonic effects such as the BAO peak are reduced. On the contrary when Ω m h 2 decreases, baryonic effects are amplified.

Verification of the Gaussianity of ξ

In this section we verify the Gaussian hypothesis ξ ∼ N (ξ BAO , C) on our lognormal simulations. Even if we 2 http://casa.colorado.edu/∼ajsh/FFTLog/ Fig. 2.-Correlation function ξ(r) as observed in the fiducial cosmology. In inset we plot r 2 ξ(r) for a better visualization. For the lognormal simulation parameters (Ωmh 2 , α, b) = (0.1323, 1, 2.5), we plot the BAO model (black), the no-wiggles model (blue), and zero-baryon model (purple). The no-wiggles model has just the BAO peak smoothed out, whereas the zero-baryon model has a different global shape. We also plot the lognormal simulations mean and error bars (red) which shows simulations are very precise. Finally we plot two other BAO models by changing Ωmh 2 = 0.1423 (orange) and Ωmh 2 = 0.1223 (green). Increasing Ωmh 2 reduces the baryonic effects such as the BAO peak, whereas decreasing Ωmh 2 amplifies these effects. do not expect large differences, the next step would be to verify it on N -body simulations, which are more realistic.

First we look at the χ 2 BAO statistic

χ 2 BAO = ξ -ξBAO, C -1 ( ξ -ξBAO) = 1≤i,j≤n ξ(ri) -ξBAO(ri) C -1 i,j ξ(rj) -ξBAO(rj)
With the Gaussian hypothesis ξ ∼ N (ξ BAO , C), χ 2 BAO follows a chi-square distribution with n degrees of freedom

χ 2 BAO ∼ χ 2 n (46)
We compare the histogram of χ 2 BAO on our lognormal simulations to the probability density function (pdf) of a χ 2 n variable where n = 18. We show figure 3 the very good agreement between the two distributions.

For the next tests, we look at wavelet methods of section 3. In these methods we obtain a wavelet response S w (R, s) and a Z-score Z w (R, s) for every parameter (R, s)

S w (R, s) = w(R, s), ξ = n i=1 w i (R, s) ξ(r i ) Z w (R, s) = S w (R, s) σ w (R, s)
We only consider the mexican hat filter with parameters R = 113.6 h -1 Mpc, s = 20 h -1 Mpc and the BAOlet filter with parameters R = 116 h -1 Mpc, s = 36 h -1 Mpc. These parameters maximize the Z-score obtained on the data measurement, in the respective studies Tian et al. (2011) and Arnalte-Mur et al. (2011). In order to obtain Z w on our simulations, we compute the noise σ w of S w using the covariance matrix of the simulations

σ w = w, Cw (47)
With the Gaussian hypothesis ξ ∼ N (ξ BAO , C), Z w is Gaussian with mean E[Z w ] and standard deviation equal to 1. We plot in figures 4 and 5 the histogram of Z w on our lognormal simulations, respectively for the mexican hat and for the BAOlet filter. As we already mentioned, there is a small bias between Z w on simulations and E[Z w ]. Here we are only interested in the Gaussianity of Z w and not in this small bias, so we compare the histogram to the pdf of a Gaussian N Z w , 1 . Again we find a very good agreement between the simulations and the Gaussian prediction. The tested hypotheses are given by

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C)
The Gaussian hypothesis on ξ is well justified (at least for our lognormal simulations) as we have seen in section 7.3. The chi-square quantities are functions of θ

χ 2 BAO,θ = ξ -ξ BAO,θ , C -1 ( ξ -ξ BAO,θ ) χ 2 noBAO,θ = ξ -ξ noBAO,θ , C -1 ( ξ -ξ noBAO,θ )
An extended model of correlation function is implicitly defined, which mixes the BAO and the no-BAO models, e.g.

ξ β,θ = β ξ BAO,θ + (1 -β) ξ noBAO,θ
Let us write as before

∆χ 2 = min θ χ 2 noBAO,θ -min θ χ 2 BAO,θ ∆χ 2 global = min θ χ 2 noBAO,θ -min β,θ χ 2 β,θ
The basic assumption of the method is that ∆χ 2 global follows a χ 2 1 distribution under H 0 . Since we have ∆χ 2 ≤ ∆χ 2 global by construction, a conservative estimate of the significance is given by ∆χ 2 .σ when ∆χ 2 ≥ 0.

We recall that ∆χ 2 global ∼ χ 2 1 is subject to the assumption that the spaces of model correlation functions (ξ noBAO,θ ) θ∈Θ and (ξ β,θ ) β∈R,θ∈Θ are affine. Since it is not easy to verify we want to test with our simulations that we have indeed ∆χ 2 global ∼ χ 2 1 under H 0 . For a model θ in H 0 , we generate realizations as

C 1/2 g + ξ noBAO,θ (48)
where g is a standard multivariate Gaussian. For each realization, we find the best-fit model by testing all the Ω m h 2 and α values on a grid. The remaining parameters, which are the bias b 2 and the parameter β, are found analytically for the best-fit.

For the grid (Ω m h 2 , α) we take Ω m h 2 ∈ [0.0423, 0.2923] with grid step 0.005 and α ∈ [0.5, 1.5] with grid step 0.01. We also allow any b 2 ≥ 0 and β. We test two models in H 0 with values Ω m h 2 = 0.1323, α = 1.0 as in the lognormal simulations, and with other values Ω m h 2 = 0.0823, α = 0.9. Each time we generate 10 000 realizations using equation ( 48) to estimate the distribution of ∆χ 2 global and ∆χ 2 . We show in tables 1 and 2 for different thresholds t, the p-value and corresponding significance for a χ 2 1 variable, for ∆χ 2 global , and for ∆χ 2 . Our results show that the assumption ∆χ 2 global ∼ χ 2 1 is clearly wrong. In particular, the mean of ∆χ 2 global for the two different models are respectively 2.85 and 2.23, whereas the mean of a χ 2 1 variable is equal to 1.

In both cases one grossly overestimates the significance when identifying ∆χ 2 global with a χ 2 1 distribution. The fact that the classical χ 2 method is conservative because it uses the value of ∆χ 2 instead of ∆χ 2 global can compensate this overestimation. For one H 0 model, identifying ∆χ 2 with a χ 2 1 distribution still gives a small overestimation of the significance. For the other H 0 model, it gives an underestimation of the significance.

Let us stress that these significances are not strictly speaking the significances of the BAO detection, but only the rejection of particular H 0 models. Indeed the BAO detection consists in the rejection of all H 0 models simultaneously.

TABLE 1 χ 2 1 ∆χ 2 global ∆χ 2 P (X ≥ 1.0) 0.32 (1σ) 0.81 (0.25σ) 0.39 (0.85σ) P (X ≥ 2.25) 0.13 (1.5σ) 0.51 (0.65σ) 0.18 (1.35σ) P (X ≥ 4.0) 4.5×10 -2 (2σ) 0.23 (1.2σ) 6.8×10 -2 (1.8σ) P (X ≥ 6.25) 1.2×10 -2 (2.5σ) 8.3×10 -2 (1.75σ) 2.1×10 -2 (2.3σ) P (X ≥ 9.0) 2.7×10 -3 (3σ) 1.9×10 -2 (2.35σ) 4.3×10 -3 (2.85σ) NOTES.
-p-values and corresponding significances for different distributions and for the rejection of the particular H0 model with Ωmh 2 = 0.1323, α = 1.0. We show the χ 2 1 distribution and the ∆χ 2 global , ∆χ 2 distributions. The assumption ∆χ 2 global ∼ χ 2 1 is wrong with a significance that is grossly overestimated. The fact that the classical χ 2 method uses ∆χ 2 instead of ∆χ 2 global compensates the overestimation. For the rejection of this H0 model, there is still a small overestimation of the significance if we identify ∆χ 2 with a χ 2 1 distribution.

TABLE 2 χ 2 1 ∆χ 2 global ∆χ 2 P (X ≥ 1.0) 0.32 (1σ) 0.66 (0.45σ) 8.0×10 -2 (1.75σ) P (X ≥ 2.25) 0.13 (1.5σ) 0.37 (0.9σ) 4.2×10 -2 (2.05σ) P (X ≥ 4.0) 4.5×10 -2 (2σ) 0.17 (1.4σ) 1.6×10 -2 (2.4σ) P (X ≥ 6.25) 1.2×10 -2 (2.5σ) 5.5×10 -2 (1.9σ) 5.7×10 -3 (2.75σ) P (X ≥ 9.0) 2.7×10 -3 (3σ) 1.3×10 -2 (2.5σ) 1.8×10 -3 (3.1σ)
NOTES.-Same as table 1 for the rejection of the particular H0 model with Ωmh 2 = 0.0823, α = 0.9. Again the assumption ∆χ 2 global ∼ χ 2 1 is wrong with a significance that is grossly overestimated. The fact that the classical χ 2 method uses ∆χ 2 instead of ∆χ 2 global compensates the overestimation. For the rejection of this H0 model, the significance becomes underestimated if we identify ∆χ 2 with a χ 2 1 distribution.

When using zero-baryon models for H 0 we find the same qualitative results. We find that ∆χ 2 global is very different from a χ 2 1 variable, and that identifying ∆χ 2 with a χ 2 1 distribution is also wrong, leading to either an overestimation or an underestimation of the significance for the rejection of particular H 0 models.

The assumption ∆χ 2 global ∼ χ 2 1 is broken because the spaces of model correlation functions (ξ noBAO,θ ) θ∈Θ and (ξ β,θ ) β∈R,θ∈Θ are not affine. It is easy to see for example that no-BAO correlations are more degenerate than BAO correlations with respect to the three parameters. So for a given range of parameters, the space of BAO correlations is larger than the space of no-BAO correlations, which tends to increase the values of ∆χ 2 global and ∆χ 2 . In this case, spaces of correlations fail to be affine because of their limited extent, which itself is due to the limited range of parameters.

As we saw in section 5, one needs to consider the "worst-case" H 0 model to obtain the correct significance. For a realization value ∆χ 2 = x the p-value is given by

p(x) = max θ∈Θ P (∆χ 2 ≥ x | H 0 , θ) (49)
We saw in tables 1 and 2 that the significance can be either overestimated or underestimated when rejecting particular H 0 models, if we identify ∆χ 2 with a χ 2 1 variable. If it is overestimated for a H 0 model with parameter θ and for ∆χ 2 = x ≥ 0, it means

P (∆χ 2 ≥ x | H 0 , θ) > p χ 2 1 (x) = 2Φ(- √ x) (50)
The consequence is that the classical χ 2 method overestimates the significance of the BAO detection. Indeed when using equation ( 49) for determining the significance of the full H 0 rejection, we get

p(x) > p χ 2 1 (x) = 2Φ(- √ x) (51)
When considering varying covariance matrices, the estimate of the significance by the classical χ 2 method could even be more wrong (see section 4.3 and section 7.5).

∆l method

In this section we test the modified version of the χ 2 method that we proposed in section 5, which we called the ∆l method. We still consider a constant covariance matrix, i.e. the hypotheses

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C)
There are two modifications that we proposed compared to the classical χ 2 method. One of the modifications consists in replacing ∆χ 2 by ∆l ∆l = -2 max θ ln (L noBAO,θ ) -max θ ln (L BAO,θ ) (52)

We are still in the case of a constant covariance, so the two statistics ∆χ 2 and ∆l are equal. We will only see the effect of this change in section 7.5 when we consider varying covariance matrices.

We also modify the procedure for computing the significance, so that we obtain a correct value as in equation ( 22) or (49). For a realization value ∆l = x the p-value is given by

p(x) = max θ∈Θ P (∆l ≥ x | H 0 , θ) (53)
Let us define a cumulative distribution function corresponding to this p-value

F 0 (x) = p(x) = max θ∈Θ P (∆l ≥ x | H 0 , θ) (54)
So the method requires to precompute F 0 in order to obtain the p-value p(x) for a given measurement ∆l = x. Here we consider the range of parameters Ω m h 2 ∈ [0.1023, 0.1623] with grid step 0.015, α ∈ [0.8, 1.2] with grid step 0.05, and b 2 ∈ [4, 9] with grid step 0.25. The grid is not very fine because the computation time is proportional to the square of the grid size, so it increases very rapidly. We tested to refine the grid for each parameter and found agreements at a few percents, which is enough for our purpose here.

For each model on the grid we generate 10.000 realizations with equation ( 48) as before, and compute the statistic ∆l. The fact that there are only 10.000 realizations for each H 0 model does not enable to quantify p-values smaller than approximately 10 -4 , i.e. significances higher than 3.85σ. Note that unlike in the χ 2 method, the imprecision here is only computational and limited to high significance. In particular when working with a data measurement ξ, the imprecision only occurs when the significance is high and the BAO detection is already clear.

We perform this procedure both for the H 0 hypothesis and for the H 1 hypothesis. With H 0 realizations we compute the function F 0 , and we get the significance obtained with H 1 realizations using equation ( 54). Since H 1 is composite, the distribution of ∆l is not well-defined under H 1 . For example we cannot speak about the expected significance obtained under H 1 . Here we simply consider the expected significance for every H 1 model, that we average over all models. This is actually equivalent to the expected significance obtained under H 1 when adding a constant prior p(θ) in the hypothesis.

We obtain an average significance of 2.11σ with this procedure. On the other hand, when estimating the significance by ∆χ 2 .σ as in the classical χ 2 method, we obtain an average of 2.33σ (we take the convention that ∆χ 2 ≤ 0 corresponds to 0σ). Let us see the effect of the imprecision at high significance by only considering realizations under the limit of 3.85σ. In this case we obtain an average significance of 1.92σ with the modified procedure, and an average of ∆χ 2 equal to 2.0. Note that this does not mean that the χ 2 method is better. Indeed it uses the same statistic, so this only means that the significance is overestimated.

Finally we repeat the same computations using zerobaryon models for H 0 . In this case we expect a larger significance because zero-baryon models not only lack the BAO feature, but also have a different global shape than baryonic models. The expected significance is higher, so we use a higher number of realizations equal to 50.000 for every H 0 model. When using the rigorous procedure for estimating the significance, we obtain an average of 2.34σ. On the other hand, when using ∆χ 2 we obtain an average of 2.9. Here the large difference is mostly due to the imprecision at high significance. Indeed if we restrict to significances under the limit of 4.25σ corresponding to this number of realizations, we obtain an average significance of 2σ with our modified procedure, and an average of ∆χ 2 equal to 2.21.

BAO detection with varying covariance matrix

We finally consider the general case where the covariance matrix depends on cosmological parameters

H 0 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ noBAO,θ , C noBAO,θ ) H 1 : ∃ θ ∈ Θ s.t. ξ ∼ N (ξ BAO,θ , C BAO,θ )
Our goal is only to illustrate the effect of a varying covariance instead of a constant covariance. So we consider a simple example where the covariance matrix only depends on the amplitude parameter b

C noBAO,θ = C BAO,θ = b b 0 4 C (55)
with b 0 = 2.5 the simulation value.

We apply the ∆l method with the same procedure as in section 7.4.2, except we take into account variations of the covariance matrix in the likelihoods and when generating H 0 , H 1 realizations using equation ( 48). In this case, the ∆l statistic is different from the ∆χ 2 statistic, and is computed using equations ( 27), ( 28), ( 29), (30), and (31).

This time we obtain an average significance of 1.96σ under H 1 , which is a bit lower than for a constant covariance. When using the ∆χ 2 statistic we obtain an average of 1.59σ under H 1 with our modified procedure. This justifies our choice of replacing ∆χ 2 by ∆l, which can be thought as a generalized likelihood ratio (see section 5). Using the classical χ 2 method, we obtain an average of ∆χ 2 equal to 2.32. In this case the estimate given by the classical χ 2 method is very far from the correct significance of 1.59σ. As we already mentioned in section 4.3, the classical χ 2 method cannot be used in the case of a varying covariance matrix.

We verify that these conclusions are not due to the imprecision at high significance of our procedure. When considering only realizations under the limit of 3.85σ, our modified procedure gives average significances of 1.89σ for ∆l and 1.52σ for ∆χ 2 , and the average of ∆χ 2 is 2.20.

These results show how the BAO detection is dependent on the tested hypotheses H 0 , H 1 , and the choice of the statistic.

7.6. Effect of varying covariance matrix on cosmological parameters constraints Finally let us see the effect of a varying covariance matrix on parameter constraints. As we saw in section 6, we must have a prior p(θ) for the posterior p(θ | ξ) to be well-defined. We consider a constant p(θ) so that the constraints only come from the measurement ξ. Then the posterior p(θ | ξ) is given by the likelihood

p(θ | ξ) ∝ L BAO,θ (56)
Changing the covariance matrix is equivalent to changing the likelihood function. Let us see the effect for a given measurement ξ. We use for the illustration the expected correlation function of the lognormal simulations ξ = ξ BAO,θ with θ = (Ω m h 2 , α, b) = (0.1323, 1.0, 2.5).

We compute the posterior p(Ω m h 2 , α | ξ) after marginalizing over the amplitude b 2 with b 2 ∈ [4, 9]. We plot the results in figure 6 and 7 respectively for a constant covariance matrix and for a varying covariance matrix. We also plot two lines of constant apparent horizon at matter-radiation equality α Ω m h 2 , and constant apparent sound horizon α (Ω m h 2 ) 0.25 . These would be degeneracy lines if we focused respectively on small scales and on the BAO scale. As expected the degeneracy direction for the constraints lies in between the two lines. 1323, 1.0, 2.5) for the illustration. We plot the 1σ to 5σ confidence regions with the approximation that p is a 2dimensional Gaussian. They correspond respectively to -2 ln(p) = -2 ln(pmax) + 2. 29, 6.16, 11.81, 19.32, 28.74 (see section "Confidence Limits on Estimated Model Parameters" in [START_REF] Press | Numerical Recipes 3rd Edition: The Art of Scientific Computing[END_REF]). We see deviations to a Gaussian posterior because the contours are not totally elliptical and symmetrical. We also plot the lines of constant apparent horizon at matter-radiation equality α Ωmh 2 , and constant apparent sound horizon α (Ωmh 2 ) 0.25 . As expected the degeneracy direction of the constraints lies in between the two lines.

First we notice that the posterior p(Ω m h 2 , α | ξ) is much farther to a 2-dimensional Gaussian for a varying covariance matrix than for a fixed covariance matrix.

Constraints on each parameter Ω m h 2 and α are obtained after marginalizing on the other parameter. We obtain different constraints in the two cases, with a small shift in the maxima of the posteriors. For a constant covariance we obtain Ω m h 2 = 0.134 ± 0.015 and α = 0.995 ± 0.070. For a varying covariance we obtain Ω m h 2 = 0.126 ± 0.014 and α = 0.976 ± 0.070.

Constraints on individual parameters are dominated by the degeneracy direction of the 2-dimensional posterior. In both cases this direction is approximately as poorly constrained. However the orthogonal direction is better constrained in the case of varying covariance. Overall the 2-dimensional constraints on (Ω m h 2 , α) are better for this particular example of correlation function.

CONCLUSIONS

We have presented different methods for BAO detection, and for each of them we detailed the tested H 0 , H 1 hypotheses and the underlying assumptions. We show in table 3 a summary with pros and cons for each method. ) 4 C. In this case the posterior p(Ωmh 2 , α | ξ) is much farther from a 2-dimensional Gaussian. The degeneracy direction is as poorly constrained as in figure 6. However the orthogonal direction is better constrained. Overall the 2-dimensional constraints on (Ωmh 2 , α) are better for this particular realization.

A first type of methods is based on wavelet filtering. Their main advantage is that they are mildly modeldependent, and mainly sensitive to the BAO feature in the correlation function. Thus they are only weakly affected by modeling errors. The price to pay is that they are outperformed by some model-dependent methods when H 0 and H 1 are well modeled.

Other methods are fully model-dependent. They assume that the measurement ξ is Gaussian, which is well verified on our simulations. They also forbid a too small binning since they can be unstable and use the inverse of the covariance matrix C of ξ. However this does not cause too much loss of information.

Among these methods, the most often used is the classical χ 2 , based on the χ 2 statistic. We found that it only gives a rough estimate for the significance of the BAO detection, and more precisely an overestimation. This comes from the method assumption that spaces of model correlation functions are affine, which is not verified in practice. As a consequence, the significance of the rejection of some particular H 0 models is overestimated. Since the rejection of the full H 0 hypothesis (the BAO detection) is based on the "worst-case" H 0 model, its significance is also overestimated. Moreover the estimate of the significance can become more wrong with hypotheses where the covariance matrix is model-dependent.

We proposed to use the ∆l method, which is a modified version of the classical χ 2 method. We first modify the procedure for obtaining the correct significance. Using simulations, we found that correct significances are indeed lower than estimates of the classical χ 2 method. The price to pay is that the method becomes much more expensive computationally. As a result we cannot use as many H 0 realizations as we want, which causes imprecision at high significance. Yet this limitation is only computational and restricted to the case of high significance where the BAO detection is already clear.

The second modification consists in replacing the ∆χ 2 statistic by the ∆l statistic, which coincide for a constant covariance matrix but are different for a varying covariance matrix. We found that the ∆l statistic gives better results in the case of varying covariance. As we have seen with a simple example, taking into these variations can affect both the BAO detection and cosmological parameter constraints.

In the course of our study we also found that nowiggles models are rejected at a lower level than zerobaryon models. It agrees with the analysis in Cabré & Gaztañaga (2011), which uses no-wiggles models for H 0 and finds that the BAO peak is rarely detected above 3σ for current galaxy surveys. This comes from the fact that no-wiggles correlation functions only lack the BAO peak, whereas zero-baryon correlation functions have a global different shape. So there must be a clear distinction between testing the existence of the BAO peak, and testing the existence of baryons.

Let us summarize our main conclusions:

1. The choice of the hypotheses H 0 and H 1 is important since it affects both the BAO detection and cosmological parameter constraints. To be rigorous one should take into account variations of the covariance matrix. It should also be clear whether one tests the existence of baryons or only the existence of the BAO peak, because the expected re-sults are quite different.

2. We have presented a new method, the ∆l method, which has two main advantages over the classical χ 2 method. Unlike the latter it provides the correct significance, apart from imprecisions at high significance. It also provides better results in the case of varying covariance matrix. We have α Λ (η Λ ) = α and α S (η S ) = α. Since α Λ and α S are increasing functions, the conditions Λ( ξ) ≤ η Λ and S( ξ) ≤ η S are equivalent respectively to α Λ Λ( ξ) ≤ α and α S S( ξ) ≤ α. If we simplify the notations and write α Λ ( ξ) for α Λ Λ( ξ) and α S ( ξ) for α S S( ξ) , we obtain from equation (C1) that for any α

P α Λ ( ξ) ≤ α | H 1 ≥ P α S ( ξ) ≤ α | H 1 (C2)
Let us show that this implies

E α Λ ( ξ) | H 1 ≤ E α S ( ξ) | H 1 (C3)
In what follows, we always keep the condition H 1 in expectations and probabilities, so we omit it to simplify the notations. We write F Λ and F S the cumulative distribution functions given by F Λ (α) = P (α Λ ( ξ) ≤ α) and F S (α) = P (α S ( ξ) ≤ α). Equation (C2) implies for any α and p F Λ (α) ≥ F S (α) (C4)

F -1 Λ (p) ≤ F -1 S (p) (C5)
The expectation of α Λ ( ξ) is given by

E[α Λ ( ξ)] = α dF Λ (α) = 1 0 F -1 Λ (p) dp (C6)
where we made the change of variable p = F Λ (α). The same computation can be made for E[α S ( ξ)], and since F -1 Λ (p) ≤ F -1 S (p) we get E[α Λ ( ξ)] ≤ E[α S ( ξ)]. So the expected significance given as a p-value is minimized for the likelihood ratio.

If we measure the significance as a number of σ instead of a p-value, both inequalities (C4) and (C5) are reversed. So the inequality on the expected values is also reversed, and the expected number of σ is maximized for the likelihood ratio.

Fig. 1 .

 1 Fig. 1.-Different analyzing wavelets used for BAO detection. We show the mexican hat with parameters R = 113.6 h -1 Mpc, s = 20 h -1 Mpc (blue) and the BAOlet with parameters R = 116 h -1 Mpc, s = 36 h -1 Mpc (red). These parameters correspond to the maximum responses Z max w in the respective studies Tian et al. (2011) and Arnalte-Mur et al. (2011).

4. 1 .

 1 The χ 2 statistic For a measured correlation function ξ ∼ N (ξ m , C) the χ 2 statistic writes

Fig. 3 .

 3 Fig. 3.-Estimated pdf of χ 2 BAO (black) using the histogram on 2000 lognormal simulations and pdf of a χ 2 18 distribution (red). Error bars give the Poisson error in the estimate due to finite number of simulations.

Fig. 4 .

 4 Fig. 4.-Estimated pdf of Zw for the mexican hat filter with parameters R = 113.6 h -1 Mpc, s = 20 h -1 Mpc using the histogram on 2000 lognormal simulations (black), and pdf of a standard Gaussian centered on Zw (red). Error bars give the Poisson error in the estimate due to finite number of simulations.

7. 4 .

 4 BAO detection with constant covariance matrix 7.4.1. Classical χ 2 method

Fig. 5 .

 5 Fig. 5.-Estimated pdf of Zw for the BAOlet filter with parameters R = 116 h -1 Mpc, s = 36 h -1 Mpc using the histogram on 2000 lognormal simulations (black), and pdf of a standard Gaussian centered on Zw (red). Error bars give the Poisson error in the estimate due to finite number of simulations.

Fig. 6 .

 6 Fig.6.-Posterior p(Ωmh 2 , α | ξ) in the case of constant covariance matrix, with ξ = ξ BAO,θ and θ = (Ωmh 2 , α, b) = (0.1323, 1.0, 2.5) for the illustration. We plot the 1σ to 5σ confidence regions with the approximation that p is a 2dimensional Gaussian. They correspond respectively to -2 ln(p) = -2 ln(pmax) + 2.29, 6.16, 11.81, 19.32, 28.74 (see section "Confidence Limits on Estimated Model Parameters" in[START_REF] Press | Numerical Recipes 3rd Edition: The Art of Scientific Computing[END_REF]). We see deviations to a Gaussian posterior because the contours are not totally elliptical and symmetrical. We also plot the lines of constant apparent horizon at matter-radiation equality α Ωmh 2 , and constant apparent sound horizon α (Ωmh 2 ) 0.25 . As expected the degeneracy direction of the constraints lies in between the two lines.

Fig. 7 .

 7 Fig. 7.-Same as figure 6 but for varying covariance matrices C BAO,θ = ( b b 0

P

  Λ( ξ) ≤ η Λ | H 1 ≥ P S( ξ) ≤ η S | H 1 (C1)Now we consider the significances corresponding to realization values Λ( ξ) = x and S( ξ) = y that we write respectively α Λ and α S α Λ (x) = P (Λ( ξ) ≤ x | H 0 ) α S (y) = P (S( ξ) ≤ y | H 0 )

In statistics an hypothesis is said to be simple when the distribution of the random variable is completely specified. It is said to be composite when the distribution is not completely specified. For example a union of simple hypotheses gives a composite hypothesis.

We plan to apply the ∆l method for BAO detection in the LRG sample of SDSS DR7. For this we also plan to use more realistic hypotheses H 0 and H 1 , by modeling the variations of the covariance matrix C( ξ). A more realistic H 1 hypothesis would also give more realistic parameter constraints.

Part of this work was supported by the European

Research Council grant ERC-228261. We would like to thank Nicolas Clerc and Vicent Martínez for useful discussions and comments, as well as the anonymous referee for improving the quality of this paper.

-Variations of covariance matrix -Long computation time → -Better results than ∆χ 2 statistic imprecise for high significance for varying covariance matrix -Unstable for small binning -Gaussian hypothesis NOTES.-The most important points are in bold. We found that the Gaussian hypothesis is well verified in practice, and that using large bins is not a serious problem. A major difference is whether methods give correct estimate of the significance. Other important differences come from the tested hypotheses H 0 and H 1 : whether they are based on a full modeling of ξ, whether they allow variations of the covariance matrix.

APPENDIX

We consider a class of binned model correlation functions, ξ θ = (ξ θ (r i )) 1≤i≤n , with a k-dimensional parameter θ = (θ 1 , . . . , θ k ) ∈ Θ. We suppose that the estimator ξ of the correlation function is Gaussian with covariance matrix C and expectation inside the model space (i.e. ∃ θ 0 such that ξ ∼ N (ξ θ0 , C)). We look at the χ 2 θ statistic which has a dependence on θ

Now we make the important assumption that the space of model correlation function (ξ θ ) θ∈Θ is a k-dimensional affine subspace of R n . Then the best-fit χ 2 θ value follows a chi-square distribution with a number of degrees of freedom equal to n -k, i.e. the measurement dimension minus the parameter space dimension

Since C is a positive definite matrix, we can consider

), so that we can rewrite the χ 2 θ statistic as

This is the Karhunen-Loève transform which consists in whitening the measurement vector ξ. This means that the resulting vector X is a multivariate Gaussian variable with expected value 0 and covariance matrix equal to the identity. Indeed the covariance matrix of X is equal to

with I n the n x n identity matrix. Thus, in any orthonormal basis of R n , the n components of X are independent standard Gaussian variables. Let us write F Θ = (X θ ) θ∈Θ , which is a k-dimensional vectorial space, and F ⊥ Θ its orthogonal complement of dimension n -k. Let us write ( Ŷ1 , . . . , Ŷk+1 , . . . , Ŷn ) the components of X into an orthonormal basis, which has the first k vectors in F Θ and the last n -k vectors in F ⊥ Θ . Then the Ŷi 's are independent standard normal variables. Moreover χ 2 θ is minimized when X θ is the projection of X onto F Θ , and equals

This shows that the best-fit χ 2 θ follows a chi-square distribution with n -k degrees of freedom, i.e. min θ χ 2 θ ∼ χ 2 n-k .

DIFFERENCE OF BEST-FITS χ 2 IN NESTED MODELS

Here we consider two nested classes of model correlation functions, ξ θ with θ ∈ Θ 1 for the first class and θ ∈ Θ 2 for the second class. We suppose that Θ 1 is k-dimensional and that Θ 2 is (k + l)-dimensional with Θ 1 ⊂ Θ 2 .

We still suppose that the estimator ξ is Gaussian with covariance matrix C and expectation inside the restricted class (i.e. ∃ θ 0 ∈ Θ 1 such that ξ ∼ N (ξ θ0 , C)). We also keep the assumption that the spaces of model correlation functions (ξ θ ) θ∈Θ1 and (ξ θ ) θ∈Θ2 are affine subspaces of R n of respective dimensions k and k + l. Then the difference of best-fits between the two classes Θ 1 and Θ 2 follows a chi-square distribution with number of degrees of freedom equal to l, i.e. the difference in the number of parameters of the two classes.

This follows easily from appendix A. We consider again the Karhunen-Loève transforms X = C -1/2 ( ξ -ξ θ0 ) and

Let us write the model spaces F Θ1 = (X θ ) θ∈Θ1 and F Θ2 = (X θ ) θ∈Θ2 , and their orthogonal complements F ⊥ Θ1 and F ⊥ Θ2 . We can write ( Ŷ1 , . . . , Ŷk+1 , . . . , Ŷk+l+1 , . . . Ŷn ) the components of X into an orthonormal basis, which has the first k vectors in F Θ1 ∩F Θ2 , the next l components in F ⊥ Θ1 ∩F Θ2 and the last n-

The components Ŷi 's are independent standard normal variables. Moreover for each class of model, χ 2 θ is minimized when X θ is the projection of X onto the model space

So the best-fit difference is given by

This shows that the difference of best-fits χ 2 θ follows a chi-square distribution with l degrees of freedom, i.e. min θ∈Θ1 χ 2 θ -min θ∈Θ2 χ 2 θ ∼ χ 2 l .

OPTIMALITY OF THE LIKELIHOOD RATIO

Let us consider the likelihood ratio Λ( ξ) = L H0 ( ξ)/L H1 ( ξ) and another statistic S( ξ) for testing the hypotheses H 0 and H 1 . Let us suppose that H 1 is preferred over H 0 for low values of S( ξ) as for the likelihood ratio (if this is not the case we just consider -S).

We first consider statistical tests for a given significance of α. The test based on the likelihood ratio is

The test based on the statistic S is

• if S( ξ) ≤ η S then accept H 1

• if S( ξ) > η S then accept H 0 with α = P Λ( ξ) ≤ η Λ | H 0 and α = P S( ξ) ≤ η S | H 0 . The Neyman-Pearson lemma states that the power of the likelihood ratio test if larger than the power of any other test. This means that the probability of accepting H 1 if it is true is larger for the likelihood ratio test.


Kelly Ribeiro

Daniel Sadarnac
email: daniel.sadarnac@supelec.frcharif.karimi@supelec.frtanguy.phulpin@supelec.frlarbi.bendani@valeo.com

Charif Karimi

Tanguy Phulpin

Larbi Bendani

A Bidirectional Three-port Current-fed LC Parallel Resonant Converter

This paper presents a Three-port current-source bidirectional resonant converter. The converter achieves zerovoltage switching (ZVS) due to the presence of a parallel resonant tank circuit, a combination between parallel inductors and capacitors. Moreover, it has reduced output current ripple and reduced power losses in the transistors if compared with a voltage source topology. The conversion occurs with self-oscillating switching technique or with forced switching in a pre-determined frequency range without losing the ZVS conditions. Due to the high frequency transformer that allows a high voltage ratio between the ports, the converter applications includes different energy storage systems that require bidirectional power flow capability, such as hybrid vehicles. A design procedure to select the converter parameters and simulation results are presented to validate the converter operation principle. I.

INTRODUCTION

Nowadays, the demand for high performance isolated DC/DC converters in multiple-port applications has been increased due to the existence of different energy sources and loads present in hybrid systems.

In parallel, in order to increase converter's power density, resonant converters have become an interesting solution. Due to its soft switching strategies as zero-voltage switching (ZVS) and zero-current switching (ZCS), the frequency can be increased without generating higher power losses in the transistors. Also, due to the sinusoidal current or sinusoidal voltage characteristics, the transformer can benefit from lower power losses if compared to other converters.

Three-port resonant converters using a series LC or a series LLC circuit tanks have been presented in [START_REF] Kim | Multiport Resonant DC-DC Converter[END_REF], [START_REF] Krishnaswami | Constant Switching Frequency Series Resonant Three-port Bi-directional DC-DC Converter[END_REF]. A currentfed resonant converter without bidirectional capability is presented in [START_REF] Nareshkumar | Closed-loop control of current-fed fullcbridge parallel resonant soft-switched dc-dc boost converter for PV applications[END_REF]. Parallel resonance is usually applied to pushpull converter stages, as in [START_REF] Namadmalan | Bidirectional Current-fed resonant inverter for contactless energy transfer systems[END_REF], where a push-pull LC parallel converter is operated under soft-switching condition and does not require bulk capacitors.

In [START_REF] Krishnaswami | A Current-fed Three-port bidirectional DC-DC Converter[END_REF], a bidirectional three-port LC parallel converter is discussed for fuel-cell hybrid system applications. The converter structure is similar to the converter presented in this text; however ZVS was not achieved because of the control methodology implemented. In [START_REF] Moorthy | Zero Current Switching Currentfed parallel resonant push-pull (CFPRPP) Converter[END_REF], a push-pull converter uses a series-parallel resonance to operate under zero-current switching condition.

The converter presented in this paper has reduced current ripple, reduced conduction power loss and does not need bulk capacitors. In addition, to reduce the total number of magnetic components and the current ripple on the DC sides, different coupling factors between the three input inductors are investigated in the converter design.

The Three-port bidirectional current-source parallel resonant converter is shown in Fig. 1. The converter is a current-fed topology and in this application, instead of large capacitances, as required in common voltage-source converters, the inductor in series with the voltage source is large enough to smooth the current. The converter study is divided into 4 sections:

The section II presents the converter analysis and operation modes. In sections III and IV, the converter design is presented by taking in consideration different coupling factors between the three input series connected inductors. The converter behavior and characteristics are verified with simulations and the results are presented in section V and VI.

II. THE CONVERTER OPERATION ANALYSIS

In the circuit presented in Fig. 1, three-voltage sources are connected in series with individual inductors to reduce the current ripple. Due to the presence of these inductors, the input current ripple is low and the current can be considered like almost constant with a small variation around its average value. The LC circuits connected in parallel to the transformer terminals are used in the resonance and a RC load is connected at the second DC side. Full-bridges are required to keep the current continuity during the switching transitions and to allow the reverse power flow operation required in the bidirectional functionality.

The transistors are switched by pairs alternatively in the diagonal direction and an overlap is used to keep uninterrupted the current flow. When the three full-bridges are simultaneously switched at a same frequency, a square current waveform is injected by each port to the three windings of the high-frequency transformer connected in parallel with the LC tank circuits.

A simplified circuit is represented in Fig. 2. In this figure the full-bridges are replaced by current sources and the transformer is considered with ideal characteristics, such as, infinite magnetizing inductance and very small leakage inductances. The passive elements are all transferred to the primary side of the transformer and an equivalent parallel inductance and capacitance are represented.

A second representation is interested in order to consider the input inductance values in the analysis. The schematic from Fig. 3 considers the voltage sources connected in series with the input inductances and the DC voltage redressed between the points "a" and "b", Vab, in Fig. 1. An equivalent AC-load from the second port is also represented and it could be connected to any of the ports. In this schematic, the MOSFET's on-resistance are considered small without generating significant voltage drop. The total current flowing to the point A, also represented in Fig. 1 is:

i ab = i 1 + i 3 + i 2 = i C + i L + i Lm (1) i Lm = i T1 + i T3 + i T2 ≈ 0, L m → ∞ (2)
The sum of the currents flowing into the transformer windings is considered as zero, as can be seen in [START_REF] Kim | Multiport Resonant DC-DC Converter[END_REF]. The voltage drop on the series inductances is calculated as a function of the three current sources on the DC sides.

V 1 = L s1 di 1 dt + u; V 3 = L s3 di 3 dt + u; V 2 = -L s2 di 2 dt + u (3)
Combining ( 1) and ( 3), the resonant voltage is obtained from:

C du 2 dt 2 + u L + u ( 1 Ls 1 + 1 Ls 2 + 1 Ls 3 ) -( V 1 Ls 1 + V 2 Ls 2 + V 3 Ls 3 ) = 0 (4)
The solution of the differential equation depends on the initial conditions: The transistors are switched during the zerovoltage crossing, while the current still keeps its same signal just before the transition, so the current in the parallel inductor achieves the maximum value when the voltage u (t) passes by zero and the capacitor is completely discharged. The equations for the voltage and current in the parallel resonant cell are presented in ( 5) and [START_REF] Moorthy | Zero Current Switching Currentfed parallel resonant push-pull (CFPRPP) Converter[END_REF].

u(t) = I Lpeak L W r sin(Wt ) + K[1-cos(Wt)] (5) i L (t) = - K L W r sin(W r t)-I Lpeak cos(W r t) + Kt L (6) W r = √ L+L s L s L C ; K = ( V 1 L s1 + V 2 L s2 + V 3 L s3 ) L L s (L+L s ) (7) 1 L s = ( 1 L s1 + 1 L s2 + 1 L s3 ) (8)
The resultant resonant frequency of the three-port converter depends on the combination between the parallel elements, but also on the series connected inductors. In the case that the value of the combined DC side series inductors in (8) is large if compared with the inductance resultant from the parallel arrangement, the resonant frequency is mainly influenced by the parallel LC tank circuit, as can be seen in ( 9):

F r = 1 2π √ L+L s L s L C ≅ 1 2π √ 1 L C (9)
At steady state, the gain between the voltage sources and the load is unitary. In reality will be less than one due to the voltage drop on the MOSFET on-resistances and parasitic elements.

From the equivalent circuit in Fig. 3, the voltage-time balance in the DC side inductor can be used to determine the resonant peak voltage, in the transformer primary side. The same calculation can be made for Ls2 or Ls3. The voltage peak can be determined using the equivalence between the voltagetime balance principle over Ls1 in half of the period and the DC voltage source V1:

V 1 - 1 π ∫ u(t) dωt π 0 = 0 => V 1 = 1 π ∫ u(t) dωt π 0 (10) 1 π ∫ u(t) dωt π 0 = 1 π [2V peak + Kπ] = V 1 V peak = π(V 1 +K) 2 (11)
If the DC side input inductances, Ls1, Ls2 and Ls3 are very large when compared to the parallel inductance L, the voltage peak is minimized and is determined in (12):

V peak = π 2 [V 1 + ( V 1 L s1 + V 2 L s2 + V 3 L s3 ) L] ≅ V 1 π 2 (12)
1. The ZVS mode The Three-port parallel resonant converter can be operated under ZVS conditions when the switching frequency is equal or smaller than the resonant frequency. At the resonant frequency, the parallel equivalent impedance is increased and the time interval during the switching transition is minimized. At this condition the circulating energy is reduced because the minimum power loss is produced in the parallel resonant elements to obtain the soft switching mode.

At frequencies smaller than the resonant frequency, the ZVS is still obtained, but the power loss is increased due to the circulating current during the switching transition times. In addition, as can be verified in (10), smaller the switching frequency, higher will be the voltage peak applied in the transistors.

The converter operation of a single full-bridge is illustrated in Fig. 4. The parasitic capacitors of the MOSFET are also indicated because they also have an impact in the resonance, being in parallel with the capacitance C1 in case that parasitic inductances are minimized.

During the first half-cycle, the switches S1 and S4 conduct. The current and the voltage are positive; when the voltage cross zero, the opposite switches (S2 and S3) can be turned-on at zero-voltage switching. Only after an overlap time, the transistors S1 and S4 are turned-off to prevent an open circuit. The reactive energy oscillates in the resonating cell and then is returned to the voltage sources, all during a half-cycle allowing the switching when the capacitors are fully discharged.

The current and the voltage waveforms are illustrated in Fig. 5 This behavior is verified in the other full-bridge sides since the same voltage u (t) is imposed by the parallel circuit tank in theirs terminals due to the transformer characteristics, which must have minimized leakage inductances and a high magnetizing inductance.

III. THE DC SIDE INDUCTORS

One of the characteristics in the three-port current-fed resonant converter is the possibility of implementing different arrangements with the DC side inductors. Their values are mainly determined by the average current in each of the ports and we can see from (11) that it plays an important role in the resonant voltage peak.

In order to obtain a good converter's performance, MOSFETs with low on-resistance should be used; and since their resistances are directly associated with their nominal voltages, an effort should be done to keep a minimized voltage peak during the converter operation.

Increase the DC side inductances it is one of the solutions to achieve this objective, but will have an impact on the converter volume. A second option will consist in the use of a single magnetic core. Using a coupling factor between the three DC side inductors will bring some advantages to the configuration, especially in the case that the DC side inductances are smaller than the parallel ones. The first advantage is the reduction of the total number of magnetic elements; and in addition, the reduction of the peak voltage, allowing the use of a MOSFET with smaller conduction on-resistance without needing to actually increasing the DC inductor values.

The common coupling factor kc between the three series connected inductances allow to write the following equations:

V 1 -u = L s1 di 1 dt -k c L s2 di 2 dt -k c L s3 di 3 dt (12) V 2 -u = L s2 di 2 dt -k c L s1 di 1 dt -k c L s3 di 3 dt (13) V 3 -u = L s3 di 3 dt -k c L s1 di 1 dt -k c L s2 di 2 dt (14)
In order to simplify, the DC sides inductances remained to the primary side in Fig. 3 are considering with the same value, Ls1, and the solution of the resonant voltage is then obtained from (15). The equivalent voltage and the resonant frequency are presented in ( 16) and (17).

C du 2 dt 2 + u L + 3 u L s1 (1-2k c ) - V 1 +V 2 +V 3 L s1 (1-2k c ) = 0 (15) u(t) = I Lpeak L W r 2 sin(Wt ) + K 2 [1-cos(Wt)] (16) W r 2 = √ 3L+L s1 (1-2k c ) (1-2k c )L s1 L C ; K 2 = ( V 1 L s1 + V 2 L s2 + V 3 L s3 ) L L s1 (1-2k c ) [L+L s1 (1-2k c )] (17)
The voltage-time principle balance when applied over Ls1 during half of the period corresponds to the DC side voltage V1. The peak voltage is now calculated as:

V peak = π V 1 2 + ( V 1 L s1 + V 2 L s2 + V 3 L s3 ) L L s1 (1-2k c ) 2 [L+L s1 (1-2k c )] (18)
The coupling factor kc represents a symmetrical contribution from two inductances to the third one, and is not higher than one. When this factor is small, the voltage peak in (18) becomes the result presented in (11). Fig. 6 presents an example to show the impact in the MOSFET's peak voltage for different ratios between the parallel and the series DC side inductances. We see that, even when the series inductance is smaller than the parallel inductance, is still possible to find a compromise and keep a security margin for the equivalent voltage applied in the MOSFET during its turned-off stage.

IV. THE CONVERTER DESIGN

The converter design should start at first with the selection of the desirable switching frequency and the nominal operating voltage. Second, after selecting the resonant frequency, for different quality factors between the parallel inductance and the load, the parallel inductance can be pre-determined. For the parallel resonant converter is interesting to obtain the parallel impedance higher than the load in order to obtain a good quality factor to the circuit. Then, using different ratios between the parallel and the DC series side inductance the last one is evaluated. A compromise should be done in order to obtain a security margin for the mosfet nominal voltage and a reduced volume. A minimum value must be ensured in order to keep low the current ripple. Because of the sinusoidal voltage applied on the transistors, the power transistor selection will have an important influence on the converter performance, since the MOSFETs for high voltage applications have an increased conduction resistance. This resistance should be as low as possible, not only to reduce the conduction power losses, but also to maintain the voltage gain close to one.

The following table shows the converter design for a typical application example.

The converter control

The converter control can be implemented in two different ways. The first one consists in the self-oscillating switching technique presented in [START_REF] Namadmalan | Bidirectional Current-fed resonant inverter for contactless energy transfer systems[END_REF] for a current-fed push-pull converter. The technique consists in to measure the voltage on the transformer terminals and detect the zero-crossing to ensure that the converter will be switched at ZVS by tracking its natural resonant frequency. This technique has the advantage that the power losses are reduced because the converter is switched with the minimum required overlap transition time. In addition, because of the self-tracking, the variation on the capacitor and inductor values are naturally considered. However, the control is not robust and load variations can provoke some instability. The second technique is to force the converter operation in a frequency smaller than the resonant frequency; in this frequency range the current will be delayed from the voltage because the equivalent impedance from the parallel inductance is smaller than the parallel capacitance, allowing the complete discharge of the capacitors before the switching transition. However, smaller the switching frequency, higher will be the voltage peak and the overlap transition time, increasing also the circulating energy and the power losses.

In order to validate the converter operation, the control proposed here is a combination of both techniques and is similar to the control technique presented in [START_REF] Namadmalan | Bidirectional Current-fed resonant inverter for contactless energy transfer systems[END_REF] for the pushpull converter. It consists in to measure the voltage at the transformer terminals to detect its zero-crossing, while keeping a forced switching signal at a pre-determined frequency, as can be verified in Fig. 7.

V. SIMULATION RESULTS

Using the design requirements from Table I and the design steps from section IV, the converter presented in this section considers the characteristics from Table II. Parasitic capacitances in the MOSFETs and also a small leakage inductance in the transformer primary side were introduced. The input DC side inductances are reduced if compared to the parallel elements and the designed resonant frequency is 200 kHz, however, due to the parasitic elements, this frequency was modified, but still kept closed of the desirable value.

As represented in Fig. 1 and according to Table 1, 480V is connected in the voltage source V1 and 20V in the voltage source V3. A MOSFET with a nominal voltage of 1200V is used in this simulation and a margin security of 350V is obtained with the design presented in Table 2. A coupling factor it is considered between the DC side inductances. In Fig. 8, the AC current in each winding of the converter is presented. The voltage source V1 connected to the first port provides 7 kW to the loads connected at the second and the third ports. In Fig. 9, the currents in the transformer terminals after the resonant cells are shown. Note that the transistors are switched when the sinusoidal voltage in the transformer winding pass by zero. The current waveforms are very similar to the expected square current waveforms. The current ripple obtained at steady-state in the DC side can be seen in Fig. 10.

VI. CONCLUSION

In this paper a modified Three-port Bidirectional Currentfed resonant converter is presented. The converter has reduced current ripple and operates with zero-voltage switching due to the presence of a parallel tank circuit connected with the transformer windings. The coupling factor between the three DC series inductances is also studied in order to reduce the number of magnetic components and the converter conduction power losses by making possible the use of transistors with a smaller on-resistance.

Fig. 1 .

 1 Fig. 1. The Three-port parallel resonant current-fed converter.

Fig. 2 .

 2 Fig. 2. The equivalent circuit with large input inductances. 1/L = (1/L1 + 1/L2 + 1/L3); C = C1 + C2 + C3.

  for a single full-bridge. During the time interval t1, the transistors S1 and S4 conduct while S2 and S3 are OFF. The input current i1 and the voltage Vab are positive, while the currents in the capacitors Cs2 and Cs3 are decreasing due to the resonant circulating energy. During t2, while S1 and S4 are still ON, the transistors S2 and S3 can be turned-on under ZVS condition, because Cs2 and Cs3 are fully discharged and their currents become zero again. At interval t3, S1 and S4 are turned-off only after S2 and S3 had been turned-on to prevent an open-circuit. The voltage Vab and the current i1 are negatives. The capacitors Cs1 and Cs4 are being charged and discharged during this interval. When they are fully discharged, during t4, the switches S2 and S3 are turned-off after the switches S1 and S4 had been turned-on. The capacitances are discharged and the parallel inductor conducts the maximum current while the voltage is zero.

Fig. 3 .

 3 Fig. 3. The equivalent circuit for analysis considering the input inductances values. 1/L = (1/L1 + 1/L2 + 1/L3); C = C1 + C2 + C3.

Fig. 4 .

 4 Fig. 4. The converter operation modes for a single full-bridge. a) S1 and S4 are ON and Vab =0. b) The capacitors Cs2, Cs3 and C1 are charged and discharged. The remaining energy goes back to V1. c) While Vab =0, S2 and S3 are turned-on under ZVS. The remaining energy is stored in L1, while the current is distributed between both arms. d) S1 and S4 are turned-OFF under ZVS. Cs1 and Cs4 are progressively charged while Vab increases. The circulating energy is now absorbed from V1.

Fig. 5 .

 5 Fig. 5. The converter operation modes for a single full-bridge.

Fig. 6 .

 6 Fig. 6. Peak voltage under the power transistors for different DC sides inductances values and coupling factor.

Fig. 7 .

 7 Fig.7. Command generation for the three-port resonant converter.

TABLE II CONVERTER

 II DESIGN SPECIFICATION:

		Appearance
		Port 1	Port 2	Port 3
	Nom.Volt (V)	480	480	20
	N	24	24	1
	Lsx	42 uH	42 uH	73 nH
	Cx	4nF	4nF	2.3 uF
	Lx	165 uH	165 uH	290 nH