C Arnaboldi

J W Beeman

O Cremonesi

L Gironi

M Pavan

G Pessina

S Pirro

E Previtali

CdWO 4 scintillating bolometer for Double Beta Decay: Light and Heat anticorrelation, light yield and quenching factors

Keywords: Double Beta Decay, Bolometers, CdWO 4, Quenching Factor

We report the performances of a 0.51 kg CdWO 4 scintillating bolometer to be used for future Double Beta Decay Experiments. The simultaneous read-out of the heat and the scintillation light allows to discriminate between different interacting particles aiming at the disentanglement and the reduction of background contribution, key issue for next generation experiments. We will describe the observed anticorrelation between the heat and the light signal and we will show how this feature can be used in order to increase the energy resolution of the bolometer over the entire energy spectrum, improving up to a factor 2.6 on the 2615 keV line of 208 Tl. The detector was tested in a 433 h background measurement that permitted to estimate extremely low internal trace contaminations of 232 Th and 238 U. The light yield of γ/β, α's and neutrons is presented. Furthermore we developed a method in order to correctly evaluate the absolute thermal quenching factor of α-particles in scintillating bolometers.

Introduction

Double Beta Decay (DBD) searches became of critical importance after the discovery of the neutrino oscillations and plenty of experiments are now in the construction phase and many others are in R&D phase [1,2,3,4]. The main challenges for all the different experimental techniques are the same [5]: i) increase the active mass, ii) decrease the background, and iii) increase the energy resolution.

Thermal bolometers are ideal detectors for this survey: crystals can be grown with different interesting DBD-emitters and, fundamental for next generation experiments, they show an excellent energy resolution.

The CUORICINO Experiment [6], constituted by an array of 62 TeO 2 crystal bolometers, demonstrated not only the power of this technique but also that the main source of background for these detectors arises from surface contaminations of radioactive α-emitters. Moreover simulations show that this contribution will largely dominate the expected background of the CUORE Experiment [7,8] in the region of interest, since there is no possibility to separate this background from the two DBD electrons. The natural way to discriminate this background is to use a scintillating bolometer [9]. In such a device the simultaneous and independent read out of the heat and the scintillation light permits to discriminate events due to γ/β, α and neutrons thanks to their different scintillation yield. Moreover if the crystal is based on a DBD emitter whose transition energy exceeds the 2615 keV γ-line of 208 Tl then the environmental background due to natural γ's will decrease abruptly. CdWO 4 is an ideal candidate for such kind of detector:

• it is a well established scintillator

• 116 Cd (7.5 % i.a) has a DBD transition at 2805 keV

• the light yield (LY) is rather large

• the radiopurity of this compound is "naturally" high Due to these favourable features this crystal compound was already used to perform a DBD experiment [10] using "standard" Photomultipliers. A large mass experiment, based on enriched CdWO 4 crystals, readout by Photomultipliers, was also proposed [11]. But this technique, limited by the modest achievable energy resolution, key point of future experiments, is no more pursued. 116 Cd is also used in different DBD Experiments. The Nemo Experiment [12] uses 405 g of enriched 116 Cd thin foils in a calorimetric-tracker detector (mainly for the study of the 2 neutrino decay spectrum). The COBRA collaboration is developing small CdZnTe semiconductor detectors [13] for a future DBD Experiment.

Experimental details

The dimensions of the CdWO 4 crystal tested are 4 cm diameter, 5 cm height. All the surfaces of the crystal are polished at optical grade. The crystal was tested as a standard scintillator at room temperature. It was wrapped with a reflecting sheet (3M Radiant Mirror film VM2000) and coupled with optical grease to a Photomultiplier (Hamamatzu R6233). The energy resolution evaluated on the 137 Cs line is 12.3 % FWHM.

Our detector setup is schematized in Fig. 1. It is held by means of four L-shaped Teflon pieces fixed to the two cylindrical Cu frames. The frames are held together through two Cu columns. The crystal is surrounded (without being in thermal contact) by a 43 mm diameter reflecting sheet (3M VM2002). The Light Detector (LD) [14] is constituted by a 36 mm diameter, 1 mm thick pure Ge crystal absorber. The surface of the Ge facing the crystal is further "darkened" through the deposition of a 600 Å layer of SiO 2 in order to increase the light abortion. Furthermore since the diameter of the LD is smaller with respect to the one of the reflecting cavity, a "reflecting ring" has been mounted in order to decrease light losses. On the opposite face of the crystal a reflecting sheet is mounted. The temperature sensor of the CdWO 4 crystal is a 3x3x1 mm 3 neutron transmutation doped Germanium, the same used in the CUORICINO experiment. The temperature sensor of the LD has smaller volume (3x1.5x0.4 mm 3 ) in order to decrease its heat capacity, increasing therefore its thermal signal. A resistor of ∼300 kΩ, realized with a heavily doped meander on a 3.5 mm 3 silicon chip, is attached to each absorber and acts as a heater to stabilize the gain of the bolometer [15,16]. The detectors were operated deep underground in the Gran Sasso National Laboratories in the CUORE R&D test cryostat. The details of the electronics and the cryogenic facility can be found elsewhere [17,18,19].

The heat and light pulses, produced by a particle interacting in the CdWO 4 crystal and transduced in a voltage pulse by the NTD thermistors, are amplified and fed into a 16 bit NI 6225 USB ADC unit. The entire waveform of each triggered voltage pulse is sampled and acquired. The time window has a width of 256 ms sampled with 512 points. The trigger of the CdWO 4 is software generated while the LD is automatically acquired in coincidence with the former. The amplitude and the shape of the voltage pulse is then determined by the off line analysis that makes use of the Optimal Filter technique. The energy calibration of the CdWO 4 crystal is performed using γ sources placed outside the cryostat. The Heat axis is calibrated attributing to each identified γ peak the nominal energy of the line, as if all the energy is converted into heat. Consequently this calibration does not provide an absolute evaluation of the heat deposited in the crystal. The dependency of amplitude from energy is parameterized with a second order polynomial in log(V) where V is the heat pulse amplitude. The three coefficients of the polynomial are fitted on calibration data (the heat pulse amplitude corresponding to the known γ lines visible in the spectrum).

The choice of such a function was established by means of simulation studies based on a thermal model of the detectors. The energy calibration of the LD is obtained thanks to a weak 55 Fe source placed close to the Ge that illuminates homogeneously the face opposed to the CdWO 4 crystal. During the LD calibration its trigger is set independent from the one of the CdWO 4 . The LD is calibrated using a simple linear function. The FWHM energy resolution of the LD, evaluated on the X doublet at 5.90 and 6.49 keV, is 480 eV FWHM. In the range of 0÷50 keV, i.e. the energy interval of the light signals, the response of the LD can be definitely assumed to be linear.

Three sets of data have been collected with this device: two calibrations using 232 Th and 40 K sources, a long background measurement (433 h) and a neutron measurement (8 h) done exposing the detector to an Am-Be source.

Light-Heat scatter plot

In this field the usual way to present the results is to draw the Light vs. Heat scatter plot. Here each event is identified by a point with abscissa equal to the heat signal (recorded by the CdWO 4 bolometer), and ordinate equal to the light signal (contemporary recorded by the LD). In the scatter plot, γ/β, α and neutrons give rise to separate bands, in virtue of their characteristic Light to Heat ratio (see Fig. 2). This feature is the result of the different LY's characterizing these particles: γ/β events belong to a distribution (see Fig. 3) characterized by the θ βγ angle, while the α events are characterized by θ α . The scintillation Quenching Factor (QF) is defined as the ratio of the scintillating yield of an interacting particle (α, neutron, nucleus) with respect to the LY of a γ/β event at the same energy. From Fig. 3 we have QF α =tan(θ α )/tan(θ βγ ). But this holds only in first approximation, as will be exposed in Sec. 5 and Sec. 6.

Within each band, monochromatic events (i.e. those corresponding to the complete absorption of a monochromatic particle in the CdWO 4 crystal) appear as sections of lines with negative slope, showing a strong anticorrelation between Heat and Light. This is evident in the Light vs. Heat scatter plot obtained with a source of 232 Th (Fig. 2) where different γ-lines are clearly visible. Despite the much lower statistics, the same feature is perceptible also for the α lines that appear in the background spectrum (Fig. 7).

A simple model accounts for the observed pattern. The energy E of a monochromatic particle, fully absorbed in the CdWO 4 crystal, is divided into two channels: a fraction (L) is spent in the production of Light (photons) and a fraction (H) is spent in the production of Heat (phonons). In absence of "blind" channels (i.e. channels in which the energy deposited into the crystal is stored in some system that do not takes part in signal formation) the energy conservation requires:

E = E Heat + E Light = (1 -k)E + kE = H + L (1)
being k the fraction of the total energy that escapes the crystal in form of light. The last equation states not only the energy conservation, but also a subtle distinction between the Heat (what is measured) and the energy definition. Monochromatic events should produce a spot in the scatter plot, with a size determined by the intrinsic energy resolution of the heat and light detectors. The observed spread along a negative slope line, however, can be accounted for assuming the existence of fluctuations in the H and L signal amplitudes. Since the energy has to be conserved the amplitude of the fluctuations (δL and δH) must compensate each other: δL=-δH. In other words the fluctuations are anti-correlated between each other and produce the observed pattern. They are visible because their amplitude is much larger than the intrinsic resolution of the detectors.

We can devise two mechanisms that introduce a fluctuation in the L/H ratio: i) the statistical Poissonian fluctuation in the emitted light: the more the energy spent in light production, the less the energy spent in heat; ii) any variation of the LY as those due to position dependent effects (inhomogeneities and defects of the crystals that modify the light emission or self absorption of scintillation photons). Regardless the mechanism, we will demonstrate, in the next section, how this anticorrelation can be used in order to increase (correct) the energy resolution of the detector over the entire energy spectrum. Fig. 3 tries to summarize the model here discussed. The dark spots represent monochromatic energy deposition in the CdWO 4 crystal, as those observed during a γ calibration. For each of them a negative slope (iso-energy) line is tracked. The intercept on the Heat axis (H Light=0 ) is the Heat that would correspond to a full heat conversion of the deposited energy (in absence of light emission or in the case in which all the light is absorbed by the crystal itself). In this simple model γ/β's and α's releasing the same energy within the CdWO 4 crystal lie on the same iso-energy line. But due to the larger LY, γ/β events will convert less heat within the crystal with respect to an α particle releasing the same energy: part of the energy escapes the crystal in form of photons. This feature can be easily deduced from Fig. 3: the α 1 particle release the same total energy of the γ 2 particle (they belong to the same iso-energy line) but they show a different position on the Heat axis. In fact the α 1 shows the same heat of the γ 1 particle that releases a larger energy within (but not "into") the crystal. The "usual" picture (in which γ/β and α of the same energy show the same heat signal) corresponds to the limit in which the light emitted is negligible (or completely re-absorbed by the crystal): in this case the distinction between H and E becomes meaningless.

H e a t θ θ βγ θ α γ 1 γ 2 α 1 γ α 1 2 α 1 γ 1 γ 1 α 1 α 1 L=0 α1 θ fit θ Light Heat Iso-energy Line γ1 L=0 =E( ) E( ) H( )=H( ) = E( ) E( ) =E( ) H H
It is important to remark that up to now it was assumed to be able to measure the absolute values of energy converted into heat and light.

Actually, the measured experimental heat and light signals can be written as

H = α(1 -k)E L = β kE (2)
being α and β the absolute calibration factors for the Heat and Light axis. The factor , instead, takes into account the overall light collection efficiency. Under the condition that α and β are not depending on the energy, and considering that the conservation of energy requires that δ

[kE]=-δ[(1 -k)E], Eq. 2 implies that δL δH = - β α = -tan(θ f it ) (3)
The last equation states that the effect of miss-calibrations of the Heat and Light axis simply implies a variation of the slope of monochromatic lines, which is exactly what we observe experimentally. Also the light collection efficiency plays a fundamental role in the evaluation of the slope, as discussed in the next section.

It is important to note that, unlike the Heat channel, the LD energy calibration is more delicate. The energy calibration on this detector is performed using ionizing X-Rays from the 55 Fe source. A priori the thermal signal that arises by the absorption of 1 keV of photons at the boundaries of the Ge crystal (moreover covered by the SiO 2 layer) could give a different thermal signal with respect to a 1 keV ionizing energy. In any case this effect still preserves proportionality, so that it can be considered within the factor β.

Light-Heat anticorrelation: energy resolution

The energy anticorrelation between the two signals of a double readout system was already observed in ionization/scintillation detectors [20] as well as in heat/scintillation bolometers [21] with a device very similar to the one presented here. In the former case this anticorrelation was demonstrated to improve the energy resolution by a factor close to ∼ 25 % (evaluated @570 keV). In paper [21] two different γ-lines were studied at rather low energy. An evident anticorrelation was found but, however, the correlation factor of the two lines was found to strongly depend on the energy.

The energy correction procedure is graphically explained in Fig. 3. Here blue (dark) and red (light) spots represent γ/β and α produced by monochromatic events. In order to obtain a spectrum with improved energy resolution we have to combine the heat and light values of each single event in an appropriate way. This can be obtained in a natural way performing a SU(1) rotation in the scatter plane an projecting the points on the Heat θ axis:

Heat θ = Heat cosθ + Light sinθ Ligh θ = -Heat sinθ + Light cosθ (4)
The value of θ can be evaluated in two different ways: i) it can be "optimized" in order to obtain the best energy resolution on the Heat θ axis (energy minimization); ii) by fitting the monochromatic spots in the scatter plot with a negative slope line, evaluating the θ f it angle for each single distribution (for construction we have

θ = π/2 -θ f it )
It is clear from this scheme that some assumptions should be satisfied in order for this procedure to work properly: For each angle a spectrum is produced and the corresponding peaks are fitted. The curves through the data points are only to guide eyes. For two peaks the corresponding error on the FWHM is not plotted for better clarity.

• the slope of each monochromatic spot (i.e. the θ f it angle) or the θ angle minimizing the FWHM has to be the same within all the lines;

• the LY of the different class of interacting particles (α's and γ/β) has to be independent from energy.

The first condition ensures that the rotation minimizes the FWHM on the whole energy spectrum. It has to be noted that if the first condition is not verified, this technique still provides a useful tool for resolution optimization. In case of different θ f it values the rotation angle will be chosen in order to reach the best performances in the region of more interest for the physics (e.g. the region where the DBD peak should appear).

The second condition ensures that the projected spectrum keeps linearity and energy calibration. This point is less trivial and will be treated in more detail in Sec. 5.

We applied both the methods in order to evaluate the rotation angle. Using the energy minimization method, the value of the rotation θ angle is allowed to vary over a wide range (0-80 degrees) and for each angle a spectrum is produced projecting the data on the Heat θ axis. Finally the FWHM energy resolution of the more intense peaks is evaluated by a fit (we use an asymmetric Gaussian reporting then the average FWHM).

We applied this method not only to the calibration spectrum but also to the background measurement. In this way we were able to study also the the 88 keV γ-line of 109 Cd (accidentally contamination of our CdWO 4 sample) and the α peaks due to internal trace contaminations: in this latter case we observe a line at the full energy of the decay (Q=α+nuclear recoil). The internal α-lines observed are 234 U-α(Q=4858 keV), 180 W-α (Q=2516 keV) and 238 U-α (Q=4270 keV).

In Fig. 4 we report the FWHM energy resolutions evaluated on the most intense lines obtained varying the rotation angle θ and projecting on Heat θ . The energy resolutions of the peaks are evaluated for each single rotated spectrum. The minimum of the energy resolution occurs within θ=(67.5 ± 0.5) deg for all the lines, except for the 88 keV Electron Capture internal line of 109 Cd.

In fact the minimum of the energy resolution for this line, 1.03 ± 0.06 keV FWHM is obtained at θ= 60 deg while at 67.5 deg it becomes 1.11 ± 0.06 keV. These two values seem rather similar, but from Fig. 4 it turns out that the resolution curve for this peak shows a different behaviour with respect to all the other lines. We do not have an explanation for this enhancement in the anticorrelation at low energies, but the same effect was already observed in [21].

Each monochromatic events distribution in the scatter plot was also fitted with a linear function, as previously discussed.

The θ f it angles as well as the most relevant parameters of this analysis are summarized in Tab. 1. The values of the θ f it angles for the α's and γ/β's lines reported here are compatible within the experimental errors within the two groups (with the exception of the 88 keV line). However, with respect to the γ/β's, the angular coefficient calculated by fitting linearly the single lines gives a result that is systematically slightly larger with respect to the value (67.5±0.5) obtained through the energy minimization. This systematic slightly difference could arise from the fact that while the bidimensional fit is model independent, the energy minimization depends on the fitting procedure of the peaks. Moreover it should be pointed out that the quoted error on θ f it don't include systematic effects (introduced, for example, by the continuum background) whose values could dominate the total error.

The values of the θ f it angles reported in Tab. 1 seems to show a slightly difference between α and γ/β. However it has to be pointed out that the quoted error on the θ f it for the α particles has a (further) systematic error that cannot be easily evaluated. The point is that the α events were collected during the long background measurement. During this period the CdWO 4 crystal was slowly cooling down; this happen very often with cryogenic detectors and is due to residual heat leaks of the materials constituting the overall setup. The base (stable) temperature is normally reached after 2-3 weeks. In our case the background measurement started immediately after the down cooling so that the initial measurements suffered by this (large) drift (the LD, on the contrary, was rather stable during the same period). This drift, as explained in Sec. 2, is usually corrected through the Heater pulse. If the temperature drift is rather small ( 50 µK) then the correction is independent from the energy. On the contrary, if the drift is larger (in this case the drift was ≈ 0.3 mK ) then the shape of the signal changes appreciably. In such cases the correction starts to depend on the energy. As a consequence a small drift in function of the time (or the temperature) becomes appreciable. Moreover this drift depends on the energy. So, in order to evaluate the second, third and last column in Tab. 1 each α peak was "self-stabilized". This can introduce an uncontrolled systematics since we are somehow forcing a distribution to obtain a minimum. On the other hand, a check was made on the weak 2615 keV γ-line present in the background spectrum. As in the case of the α's, we had to self-correct the drift. The obtained θ f it value is 71.5±0.4, rather consistent with the α background. For the following we will assume θ γ/β f it = θ α f it . The 232 Th and 40 K calibrations, as well as the 88 keV analysis, on the other hand, were made at the end of the measurement in a rather stable working temperature.

As a final remark it has to be noticed that the projection of the events on the Heat θ axis introduces an Heat (energy) normalization factor K that will be different for different classes of particles, namely K γβ,α =cosθ(1+tanθ tanθ γβ,α ). This feature is rather evident in Fig. 3: in the Heat axis we have H(α1)=H(γ1) while in the rotated Heat θ axis we have H(α1)=H(γ2). This has two consequences: the first, obvious, is that the intercalibration between α and γ/β is different within the Heat spectrum and the Heat θ spectrum. The second is that the same consideration holds within the same class of events: if the LY of α particles (for example) depends slightly on energy, then the linearity of the projected spectrum will change.

Light Yield and Scintillation Quenching Factor

In Tab. 1 we also present the values of the LY, determined on the different monochromatic lines. We define the LY as the energy released in the LD (in keV) for a nominal energy deposition of 1 MeV in the scintillating crystal. The table shows that the LY of the γ/β events is, within the error, the same for all the observed lines (including the 88 keV line that shows a different resolution curve in Fig. 4). On the other hand a small, but rather clear, energy dependence is evident for α particles. This effect can also be observed in Fig. 6 in which it is evident that the α particles belong to a curve rather than a straight line. Moreover, as already stated, the observed α-lines are internal so that the effect cannot be ascribed to surface effects. A theoretical explanation for this behaviour and a detailed discussion of the energy dependence of LY and QF in different scintillating crystals can be found in [22]. Essentially, the observed behaviour reflects the fact that an α particle scintillates less with respect to an electron because of its larger dE/dx which can induce saturation effects in the scintillator (the Birks law [23]). The energy dependence of the LY and QF is a consequence of the energy dependence of the stopping power. Electrons, due to their low stopping power, do not suffer of saturation effects and their LY is, consequently, energy independent. For this reason the definition of QF α =tan(θ α )/tan(θ βγ ), given in Sec. 3, represents only an approximation.

Using the values in the table we can evaluate the scintillation Quenching Factors of α particles (QF=LY(α)/LY(β)): QF180 W =0.160±.006, QF238 U =0.192±.006, QF234 U =0.201±.006.

Within this framework, furthermore, it is possible to evaluate the total amount of light k that escapes the crystal.

Using Eq. 3 and considering that, for definition, LY≡L/E=β k, we obtain:

LY = tan(θ f it )αk (5)
Now we have to note that the energy calibration we adopted for the Heat axis (as discussed in Sec. 2) implies H≡E. This means (see Eq. 2) that α=1/(1-k γ ), being k γ the absolute LY for γ/β events. Using this last relation in the last equation we finally get for γ/β

k γ = LY γ tan(θ f it ) + LY γ (6)
Using the values of Tab. 1 we get k γ = 4.6 % As a final remark, we point out that the large spread in the light channel (6.8 % at the 2615 keV scintillation signal), as observed in the scatter plots, cannot be dominated by the fluctuation of the Poissonian statistics of the absorbed photons. In fact we have that the mean energy absorbed by the LD for a 2615 keV interaction in CdWO 4 is of the order of 2.615•17.6= 46 keV. Assuming ≈ 3 eV/photon we get (assuming a Fano fac-tor=1) the variation in the light channel to be of the order of 1/ √ 46000/3=0.8 %. Moreover in other tests performed with different CdWO 4 crystal samples (showing "only" a slightly larger LY) we observed the same anticorrelation (in terms of θ f it ) but with a much smaller spread. Thus the magnitude of this spread is probably dominated by fluctuation in the overall light collection efficiency, mainly due to variation of the LY due to position dependent effects. In particular we believe that the grade of the surfaces of the crystal (in terms of diffusion of the scintillation light) plays an important role in this mechanism.

Heat absolute scale and Heat Quenching Factor

From the previous sections is evident that for scintillating bolometers the Heat/Energy scale is different for different kinds of interacting particles (or, better, for particles characterized by different LY's). In other words, with our convention for the Heat axis calibration, H measures the total energy of a γ/β particle while it doesn't for an α: the energy calibration of αparticles has to be dealt in a separate way.

This "displacement" for the α lines is evident in our experimental data. In particular the 180 W-α line (as can be seen in Fig. 2), whose nominal Q-value is 2516 keV [24], appears (in the γ-calibrated spectra) at 2627±2 keV. This line is rather peculiar since it definitely proves that the mis-calibration cannot be ascribed to a wrong extrapolation of the detector amplitude vs. energy calibration (the highest γ-line used for the calibration of the Heat signal is indeed the 208 Tl line at 2615 keV).

The observed shift in keV, ∆E = E α -E γ/β , between the experimentally reconstructed energy and its nominal value is given -for the three mentioned α's -by:

∆E180 W = 110±3 ∆E238 U = 172±3.5 ∆E234 U = 180±2 (7)
This shift can be easily recognized by looking at Fig. 3. Let us consider the case in which the crystal does not scintillate: in this case the two particles, α1 and γ1 will show, as example, E(α1) < E(γ1). If the crystal scintillates and the light escapes the crystal, then the two events will move left on the two different iso-energy lines up to reach the case in which H(α1)=H(γ1). This means that the energy shift can be written as ∆E = H L=0 γ1 -H L=0 α1 . By simple geometric consideration we get

E Light=0 γ/β = E + Light γ/β tan(θ f it ) = E 1 + LY γ/β tan(θ f it ) (8)
So, finally, we obtain

∆E = E [MeV] • LY γ/β -LY α tan(θ f it ) [keV] (9)
This formula represents the energy displacement that characterizes two isoenergetic particles impinging a generic scintillating crystal, induced by their different scintillation yield (in our case LY γ/β and LY α ).

Inserting the experimental values we get:

∆E180 W = 107±3 ∆E238 U = 175±6 ∆E234 U = 195±7 (10)
These last values are in very good agreement with the experimental ones in Eq. 7, confirming our model.

To conclude this section, we need to face an important question about the definition of the Heat Quenching Factor (HQF). The conventional way to treat differences in the response of a detector with respect to different type of interacting particles is to introduce relative quenching factors. This is exactly what is done for the scintillation light. In the same way, following what is often done in literature [25,26,27,28], it is possible to introduce the HQF as the ratio between the Heat signal of an α particle and that of a γ/β of the same energy, fully absorbed in the detector. With this definition, for our detector we can quote the HQF measured for the 180 W α decay as: HQF180 W =2627/2516=1.044. It is straightforward to show that in our framework we can easily evaluate HQF=1+

LY γ/β -LY α tan(θ f it ) = 1 + LY γ/β tan(θ f it ) (1 -QF α ).
Even if the definition of the (scintillation) QF and the HQF is identical, due to the different physical mechanism that generates them, they show a substantial difference. While the former depends, substantially, only on the "nature" of the scintillator, the HQF, instead, depends on the amount of light (energy) that escapes the crystal: it depends on the size and/or on the quality of the crystal. It is also obvious that since some (minor) amount of light that escapes the crystal can be reflected by the reflecting foil and can be absorbed by the crystal itself, then the HQF will depend also by the setup itself.

With the 3x3x2 cm 3 CdWO 4 crystal measured in [9], we found the position of the 180 W line at 2670 keV, giving therefore HQF=1.061. Now one relevant question that could be addressed is: do we have blind channels in our detector? We don't have a general method to investigate this point, but the results presented so far are still valid even if we assume that a fixed fraction of the total energy is lost in blind channels (we used only Heat/Light ratios in our framework). Then we can investigate if this fraction is different between a γ/β and an α. In other words we could ask ourselves if the sum of the heat and light signals is the same for particle converting the same total energy in the scintillator. We provide therefore an alternative definition of HQF as the ratio of E Light=0 . Provided that our model that compensates the energy losses due to scintillation is correct, using the values of Eq. 10 we can evaluate the (corrected) HQF: HQF( 180 W)=1.001± .002 HQF( 238 U)=1.000 ± .002 and HQF( 234 U)=0.997 ± .002. These HQF's are compatible with 1 meaning that there are no intrinsic differences in the heat signal generated by an α with respect to a signal generated by a γ/β1 . Alternatively we can state that any eventually present blind channel have the same behaviour for the two. The liquid standard also contains a small amount of 234 U. The alpha curve is completely separated from the γ/β line. Moreover it can be noted that the α events are not belonging to a straight line, showing a decrease of the QF with the energy, as discussed in Sec. 5.

Background Rejection: α's and neutrons

The detector was also exposed to a neutron source. It consists of a 185 kBq Am-Be source with a neutron production rate of ≈10 n/s. The neutron spectrum has its maximum at ∼5 MeV with an high energy tail reaching 10 MeV. The source was placed inside the shielding of the cryostat, in the same place where the γ sources are inserted. The scatter plot is shown in Fig. 5: neutron direct interactions are clearly visible, especially in the inset. The γ/β events extends well above 2615 keV due to (n,γ) reaction on the surrounding materials but, mainly, by the source itself: for each neutron produced there is 60 % probability to produce an excited state of 12 C that emits a γ of 4.44 MeV.

Moreover it has to be remarked (see also Sec. 8) that 113 Cd has a huge neutron capture cross section, with a Q-value larger than 9 MeV, so that "mixed events" are possible (a neutron scatters on the Oxygen of the crystal and then is absorbed by the 113 Cd with subsequent de-excitation of the nucleus). From Fig. 5 we evaluate the neutron scintillation QF with respect to γ/β to be (0.14±0.03). The error is dominated by the systematics induced by the extreme weakness of the scintillation light in the region 0÷300 keV, where most of the events are recorded.

In a subsequent test, an α source was mounted close to the detector in order to evaluate its rejection capability. The custom source was build up using 2 µL of a 0.1 % 238 U liquid Standard. The U "droplets" were dried on a 1 cm 2 Al tape and then covered with a 6 µm Aluminized Mylar foil in order to "smear" the alpha energy down to the region of interest. The source was then faced to the crystal on the face opposed to the LD. In Fig. 6 we present the scatter plot obtained with the above mentioned source during a 93 hours 232 Th calibration measurement. Fixing at 4σ the acceptance on the light signal of the 2615 keV γ-line, we evaluate a rejection factor for α particles >35σ at the same energy. As far as the fast neutron interaction is concerned , the rejection factor will be obviously larger, but the most dangerous contribution will arise from thermal neutron absorption through 113 Cd. But, at least in principle, thermal neutron can be effectively shielded while fast neutron, via µ-spallation within the shielding or close to it, cannot.

ACCEPTED MANUSCRIPT

Background measurement: internal contaminations

In Fig. 7 we present the Light vs. Heat scatter plot obtained in the 433 h live time background measurement. It can be easily recognized that we observe only two transitions belonging to the Uranium chain, namely 238 U and 234 U, corresponding to an internal contamination of 3.1•10 -12 g/g and 5.7•10 -17 g/g, respectively. The Uranium chain is broken at 234 U, while a contamination of 210 Po is clearly seen. As it often happens, it is rather difficult to determine if this Po contamination arise from its parent, 210 Pb. Moreover the absence of a clear peak indicates that the contamination is partially penetrating (few µm) the surface of the crystal (or its surroundings). No events can be ascribed to the Thorium chain, giving a limit of 9•10 -13 g/g (95% CL) for 232 Th, one of the most "dangerous" contaminant for DBD searches. But the main (unexpected) feature resulting from Fig. 7 is the γ/β background above 2615 keV. These 11 events cannot arise from external γ's due to 208 Tl or 214 Bi since the observed background lines (583 keV, 2615 keV and 1764 keV) are too weak to allow such contribution at high energy. The same holds for internal contaminations belonging to Th and U chain. Unrecognized pile-up of the 113 Cd-β decay (Q= 320 keV) with a 2615 keV γ cannot contribute at this level (we expect ∼ 2 ev/year). Internal contaminations due to rare high-energy β emitters like 106 Rh (Q=3541 keV) are extremely difficult to evaluate. But, in any case, we observe γ/β events up to 4.5 MeV, extremely large to be associated with standard "known" β emitters. On the other hand we previously tested different CdWO 4 scintillating crystals [30], obtaining, with larger statistics, no events above 2615 keV. The main difference with respect to [30] is that in the present work the environmental neutron shielding (consisting of 7 cm of polyethylene and 1 cm of CB 4 ) was removed from the top of the cryostat, leaving ∼0.6 m 2 opening, corresponding to ≈ 8 % of the total neutron shielding. We think that the thermal neutron absorption by 113 Cd is the most probable mechanism in order to explain the anomalous background above 2615 keV.

Conclusions

For the first time a "large" scintillating bolometer to be used for future DBD searches was fully characterized in details in terms of energy resolution, internal contaminations and particle identification capabilities. It was shown how the use of the anticorrelation between light and heat improves the energy resolution by a factor 1.26 @ 88 keV and up to a factor 2.6 @ 2615 keV. We developed, for the first time, a procedure that evaluates the Heat Quenching Factor corrected for the loss of the scintillation light.

This CdWO 4 , grown without any precaution in terms of radiopurity, shows extremely low trace contaminations in U and Th. Nonetheless a future experiment based on this compound should use isotopic enrichment, for two reasons: i) increase the mass of 116 Cd; ii) decrease the fraction of 113 Cd in order to avoid β decay and, more important, decrease (n,γ) reaction due to 113 Cd.

As a final remark we point out that such a crystal, depleted in 113 Cd, could be an extremely interesting compound not only for DBD experiments but also for Dark Matter searches, due to its "uncommon" radiopurity. Our LD, obviously, does not permit to discriminate charged particles with respect to neutral particles below 40 keV. On the other hand, the LY we measured is comparable with that one of the scintillators (CaWO 4 , ZnWO 4 ) used in the CRESST [31] Experiment , opening, at least in principle, the possibility of having a cryogenic experiment both for DBD and DM Searches.

Acknowledgements

The results reported here have been obtained in the framework of the Bolux R&D Experiment funded by INFN, aiming at the optimization of a cryogenic DBD Experiment for a next generation experiment. Thanks are due to E. Tatananni, A. Rotilio, A. Corsi and B. Romualdi for continuous and constructive help in the overall setup construction. Finally, we are especially grateful to Maurizio Perego for his invaluable help in the development and improvement of the Data Acquisition software.

Figure 1 :

 1 Figure 1: Setup of the detector. The two Cu columns, holding the two frames, are not visible from the chosen perspective.

Figure

  Figure Scatter plot Light vs. Heat obtained in a 96 h calibration using an external 232 Th source. In the inset the highlight of the 208 Tl line. In the circle two events due to the internal α-decay of 180 W. The energy spectrum below 400 keV is completely dominated by the 113 Cd β-decay, with a rate close to 0.4 Hz.

Figure 3 :

 3 Figure 3: Light vs. Heat scheme for the interpretation of the energy correction method. The blue (dark) spots represent the γ/β monochromatic events, while the red (light) represents an α.

Figure 4 :

 4 Figure4: FWHM energy resolutions evaluated for different θ angles. For each angle a spectrum is produced and the corresponding peaks are fitted. The curves through the data points are only to guide eyes. For two peaks the corresponding error on the FWHM is not plotted for better clarity.

Figure

  Figure Scatter plot Heat vs. Light obtained in a 8 hour measurement with an Am-Be neutron source. In the inset the low energy region. The energy spectrum extends up to 9 MeV. The 4.44 MeV "line" due to the 12 C * is visible. The intrinsic energy spread of this peak is rather large (≈ 1 %) due to Doppler shift induced by the "in flight" decay of 12 C * in the source.

Figure 6 :

 6 Figure 6: Scatter plot Light vs. Heat obtained in a 93 hours232 Th calibration measurement with a smeared α source facing the CdWO 4 . The liquid standard also contains a small amount of 234 U. The alpha curve is completely separated from the γ/β line. Moreover it can be noted that the α events are not belonging to a straight line, showing a decrease of the QF with the energy, as discussed in Sec. 5.

Figure 7 :

 7 Figure 7: Scatter plot Light vs. Heat obtained in a 433 hours live time background measurement. Part of the 210 Po "bump" is characterized by smaller light emission. This could arise from a worse light collection efficiency, especially from the lateral part of the crystal. The energy scale of the α is corrected according with what exposed in Sec. 6.

Table 1 :

 1 Tablewiththe main parameters of the detector. The theoretical FWHM energy resolution of the 4 crystal, evaluated through Optimal Filtering, is 0.8 keV. The LY of the α-particles (last three lines) is evaluated for the overall Q-value of the decay. The correct value (assuming negligible scintillation for the nucleus recoil) is ≈ 2% larger. †See text.

	Line	FWHM θ=0 FWHM θ=67.5	LY	π/2 -θ f it
	[keV]	[keV]	[keV]	[keV/MeV]	[deg]
	88	1.4±0.08	1.11±0.06	17.3±0.2	61.6±0.9
	583	5.2±0.4	2.1±0.2	17.0±0.35	70.8±0.4
	911	8.4±0.5	3.05±0.2	17.5±0.2	70.6±0.3
	1461	11.7±0.5	3.27±0.2	17.45±0.15	70.4±0.3
	2615	16.5±0.5	6.25±0.22	17.6±0.1	69.8±0.2
	2516	5.2±1.6 †	4.3±1.3 †	2.74±0.08 73.0±0.2 †
	4270	7.0±1.5 †	4.4±1.0 †	3.29±0.07 72.3±0.2 †
	4858	7.8±1.2 †	4.8±0.7 †	3.44±0.08 72.7±0.2 †

Preprint submitted to Atroparticle PhysicsJuly 1, 2010

We note here that although we always discuss of α's it is true that we are always considering the system (α + nuclear recoil), where however the nuclear recoil transports only a small fraction of the total energy.


Philippe Angot
email: angot@cmi.univ-mrs.fr

Pierre Fabrie
email: pierre.fabrie@math.u-bordeaux1.fr

Roger Temam

CONVERGENCE RESULTS FOR THE VECTOR PENALTY-PROJECTION AND TWO-STEP ARTIFICIAL COMPRESSIBILITY METHODS

Keywords: 2000 Mathematics Subject Classification. Primary: 35Q30, 76D05, 76N10, 35A35; Secondary: 65M12, 65N12 Artificial compressibility, Navier-Stokes equations, vector penaltyprojection, pseudo-compressibility, penalty method

Convergence results for the vector penalty-projection and two-step artificial compressibility methods

1. Introduction and setting of the problem. The artificial compressibility method was introduced by Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] and Temam [START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF] for the solution of the unsteady incompressible Stokes or Navier-Stokes equations; see also [START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF] for the theoretical analysis. Then, some other numerical schemes to efficiently compute the solutions of Navier-Stokes problems can be viewed as discretizations of perturbed systems of the type of penalization [START_REF] Shen | On error estimates of the penalty method for unsteady Navier-Stokes equations[END_REF] or pseudo-compressibility. This is the case with the famous projection methods from Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and Temam [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires I[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF] and their variants [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], see e.g. [START_REF] Shen | On a new pseudocompressibility method for the incompressible Navier-Stokes equations[END_REF].

Here, we present a new approximation method for the Navier-Stokes equations modeling incompressible viscous flows in a bounded regular open set Ω endowed with Dirichlet boundary conditions on Γ = ∂Ω (Lipschitz-continuous). With a given source term f , the Navier-Stokes system reads:

2 PHILIPPE ANGOT AND PIERRE FABRIE        ∂v ∂t + (v • ∇)v - 1 R e ∆v + ∇p = f, div v = 0, v(0) = v 0 , v |Γ = 0,
where R e denotes the Reynolds number.

According to the identity, -∆ϕ = curl curl ϕ-∇ div ϕ, we consider the following approximate method to obtain a solution of the above Navier-Stokes system, with the parameters r ≥ 0, γ > 0 and, ε > 0

                                 ∂ v ε ∂t + (v ε • ∇) v ε + 1 2 (div v ε ) v ε + 1 R e curl curl v ε - 1 R e ∇ div v ε -r∇div v ε + ∇p ε = f ∂ v ε ∂t + (v ε • ∇) v ε + 1 2 (div v ε ) v ε + 1 R e curl curl v ε - 1 R e ∇ div v ε -r∇ div v ε - 1 εR e ∇ (div v ε + div v ε ) = 0 v ε = v ε + v ε γ ∂p ε ∂t + γp ε + 1 ε div v ε + r div v ε = 0.
We associate to the previous system the following boundary conditions and initial data,

v ε (0) = v 0 , v ε (0) = 0, p ε (0) = p 0 , v ε • ν Γ = 0, v ε • ν Γ = 0, (curl v ε ) ∧ ν Γ = 0, (curl v ε ) ∧ ν Γ = 0,
where ν denotes the outward unit normal vector on Γ. To vanish, at the limit process, the two tangential component of the velocity fields, v ε ∧ ν and v ε ∧ ν, we use a penalization method which will be detailed below. This method is close to the artificial compressibility method of Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] and Temam [START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF], but presents one important difference. It is a two-step splitting method. The first equation of the previous system gives a predicted velocity v ε and the second one is the approximate projection of v ε on the free-divergence vector fields. This equation may be seen as an approximate method to solve the well-posed problem (see appendix A) :

div v ε = -div v ε , curl v ε = 0, v ε .ν |Γ = 0.
Remark 1. This approximate method is issued from the Vector Penalty-Projection (VPP r,ε ) methods for the numerical solution of unsteady incompressible viscous flows introduced in [START_REF] Ph | Vector penalty-projection method for the solution of unsteady incompressible flows[END_REF] and [START_REF] Ph | A new fast method to compute saddle-points in constrained optimization and applications[END_REF]. A fast version of these methods, the so-called (VPP ε ) method, is recently proposed also for the numerical solution of the nonhomogeneous Navier-Stokes equations in [START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-stokes problems[END_REF]. It is shown to be very efficient to compute multiphase flows, i.e. fast, cheap, and robust whatever the density, viscosity or permeability jumps.

Even for r = 0, the resulting method which corresponds to a two-step pseudocompressibility method, is different from the original artificial compressibility method of Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] and Temam [START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF][START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF].

The new important point is the penalty term 1 ε (∇div v ε + ∇div v ε ) that appears in the velocity correction step which allows us a direct estimate on the divergence of the velocity. Moreover, this system is quite easy to solve and presents good stability properties, see [START_REF] Ph | Vector penalty-projection method for the solution of unsteady incompressible flows[END_REF][START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A new fast method to compute saddle-points in constrained optimization and applications[END_REF]. The velocity v ε and the pressure p ε satisfy the equations:

                           ∂v ε ∂t + (v ε • ∇)v ε + 1 2 (div v ε ) v ε + 1 R e curl curl v ε - 1 R e ∇ div v ε -r ∇div v ε - 1 εR e ∇div v ε + ∇p ε = f γ ∂p ε ∂t + γ p ε + 1 ε div v ε + r div v ε = 0, v ε (0) = v 0 , p ε (0) = p 0 , v ε • ν Γ = 0, (curl v ε ) ∧ ν Γ = 0.
The vanishing of the tangential component at the limit process, is fullfilled by a penalization method, which implies that this boundary condition is satisfied at the order ε for the approximate solution.

L p (Ω) = (L p (Ω)) d H div (Ω) = {v ∈ (L 2 (Ω)) d , div v ∈ L 2 (Ω)} H = {v ∈ (L 2 (Ω)) d , div v = 0, v • ν Γ = 0} H 1 ν (Ω) = {v ∈ (H 1 (Ω)) d , v • ν Γ = 0} G = {v ∈ (L 2 (Ω)) d , ∃q ∈ H 1 (Ω), v = ∇q} 1.2. Mathematical recalls.
Proposition 1. Under the previous hypothesis, one has the following properties:

L 2 (Ω) = H ⊕ G, Ker (curl ) = G.
Moreover, there exists one constant C > 0 depending only on Ω such that:

u 2 H 1 = u 2 L 2 + ∇u 2 L 2 ≤ C u 2 L 2 + div u 2 L 2 + curl u 2 L 2 , ∀u ∈ H 1 ν (Ω). (1)
Besides, if we suppose that the open set Ω is simply-connected, there exist two constants λ 0 and λ 1 depending only on Ω such that:

u 2 L 2 ≤ λ 0 div u 2 L 2 + curl u 2 L 2 , ∀u ∈ H 1 ν (Ω), u 2 L 2 + ∇u 2 L 2 ≤ λ 1 div u 2 L 2 + curl u 2 L 2 , ∀u ∈ H 1 ν (Ω)
and we have:

Ker (curl ) ∩ H = {0}.
Proof. All these results may be found in [START_REF] Girault | Finite Element Methods for the Navier-Stokes Equations[END_REF] and [START_REF] Foias | Remarques sur les équations de Navier-Stokes stationnaires et les phénomènes successifs de bifurcation[END_REF].

For a Banach space E we introduce the Nikolskii space defined for 1 ≤ q < +∞, 0 < σ < 1:

N σ q (]0, T [; E) = f ∈ L q (]0, T [; E), sup 0<h<T f (• + h) -f (•) L q (]0,T -h[;E) h σ < +∞ ,
endowed with the following norm:

f N σ q = f q L q (]0,T [;E) + sup 0<h<T 1 h σ f (• + h) -f (•) L q (]0,T -h[;E) q 1 q .
Let us recall the following property see for example [START_REF] Boyer | Éléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF] page 105, Proposition 2. Let H be an Hilbert space and f a function given in

L 2 (]0, T [; H) such that, for some 0 < σ < 1, R |τ | 2σ F( f )(τ ) 2 H dτ ≤ C 2 ,
where f denotes the extension by 0 of the function f outside [0, T ]. Then f ∈ N σ 2 (]0, T [; H) and we have

f N σ 2 ≤ M σ (1 + C),
where M σ is a constant depending only on σ.

We now recall the important compactness theorem, see for example [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF] Theorem 1. Then, for all 1 ≤ q ≤ +∞ and 0 < σ < 1, the imbedding

L q (]0, T [; B 0 ) ∩ N σ q (]0, T [; B 2 ) ֒→ L q (]0, T [; B 1 )
is compact.

2. Main result. We associate to the previous approximate system, the variational problem where the tangential components of the velocities v ε and v ε are penalized. This problem is studied in the next section.

Find ( v ε , v ε , p ε ) in L ∞ (]0, T [; L 2 (Ω))∩L 2 (]0, T [; H 1 ν (Ω)) 2 × L ∞ (]0, T [; L 2 0 (Ω)) satisfying in D ′ (]0, T [),                                                                                                          Ω ∂ v ε ∂t • ϕ dω + Ω (v ε • ∇) v ε + 1 2 (div v ε ) v ε • ϕ dω + 1 R e Ω curl v ε • curl ϕ dω + 1 R e Ω div v ε div ϕ dω + r Ω div v ε div ϕ dω - Ω p ε div ϕ dω + 1 ε Γ ( v ε ∧ ν) • (ϕ ∧ ν) dσ = Ω f • ϕ dx, Ω ∂ v ε ∂t • ψ dω + Ω (v ε • ∇) v ε + 1 2 (div v ε ) v ε • ψ dω + 1 R e Ω curl v ε • curl ψ dω + 1 R e Ω div v ε div ψ dω + r Ω div v ε div ψ dω + 1 εR e Ω div v ε + div v ε div ψ dω + 1 ε Γ ( v ε ∧ ν) • (ψ ∧ ν) dσ = 0, v ε = v ε + v ε , γ Ω ∂p ε ∂t π dω + γ Ω p ε π dω + 1 ε Ω π div v ε dω + r Ω π div v ε dω = 0, ∀(ϕ, ψ, π) ∈ (H 1 ν (Ω)) 2 × L 2 0 (Ω), v ε (0) = v 0 , v ε (0) = 0, p ε (0) = p 0 . (2)
Then the velocity v ε and the pressure p ε satisfy in

D ′ (]0, T [),                                        Ω ∂v ε ∂t • ϕ dω + Ω (v ε • ∇)v ε + 1 2 (div v ε ) v ε • ϕ dω + 1 R e Ω curl v ε • curl ϕ dω + 1 R e Ω div v ε div ϕ dω +r Ω div v ε div ϕ dω + 1 εR e Ω div v ε div ϕ dω - Ω p ε div ϕ dω + 1 ε Γ (v ε ∧ ν) • (ϕ ∧ ν) dσ = Ω f • ϕ dω, ∀ϕ ∈ H 1 ν (Ω), v ε (0) = v 0 (3)
Remark 2. In order to establish the strong convergence of the sequence (v ε ) ε>0 when ε → 0, we use in Section 4 the Leray's orthogonal decomposition in the bounded domain. The curl-free component vanishes with the penalty term introduced by our method, whereas the divergence-free component strongly converges thanks to an estimate of a fractional derivative in time, see [START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF]. However, this requires to consider velocity fields having only their normal component which is zero on the boundary. Since at the limit process, we aim at solving the Navier-Stokes problem with homogeneous Dirichlet boundary condition, we also penalize the tangential part of the velocity fields.

We prove in section 3 the following results.

Lemma 2.1. Let us suppose that f belongs to L 2 (]0, T [; L 2 (Ω)). Then, there exists at least a solution to the system (2). This solution is unique in two space dimension.

For the dimension d ≤ 3, this solution satisfies the following energy inequality:

1 2 d dt r ε v ε 2 L 2 + ε v ε 2 L 2 + v ε 2 L 2 + γε p ε 2 L 2 + γε p ε 2 L 2 + r ε 2R e curl v ε 2 L 2 + ε 2R e curl v ε 2 L 2 + 1 2R e curl v ε 2 L 2 + 1 2R e div v ε 2 L 2 + r ε 2R e div v ε 2 L 2 + ε 2R e div v ε 2 L 2 + εr 2 div v ε 2 L 2 + εr div v ε 2 L 2 + 1 2εR e div v ε 2 L 2 + r |( v ε ∧ ν)| 2 L 2 (Γ) + |( v ε ∧ ν)| 2 L 2 (Γ) + 1 ε |(v ε ∧ ν)| 2 L 2 (Γ) ≤ λ R e 2 (1 + rε) f 2 L 2
For the dimension d = 2, one has the following energy equality:

1 2 d dt r ε v ε 2 L 2 + v ε 2 L 2 + γε p ε 2 L 2 + r ε R e curl v ε 2 L 2 + 1 R e curl v ε 2 L 2 + γ ε p ε 2 L 2 + r ε R e div v ε 2 L 2 + 1 R e div v ε 2 L 2 + εr 2 div v ε 2 L 2 + 1 εR e div v ε 2 L 2 + |( v ε ∧ ν)| 2 L 2 Γ + 1 ε |(v ε ∧ ν)| 2 L 2 (Γ) = r ε Ω f v ε dω + Ω f v ε dω.
This result is quite classical and we only give the sketch of proof in the section 3.

In fact, we can precise the previous energy inequality if we suppose that the data f belongs to L ∞ (R + ; L 2 (Ω)). This shows the absolute stability of the approximate method.

Theorem 2.2. Suppose that the data f satisfies

f ∈ L ∞ (R + , L 2 (Ω)),
then, there exists a constant α independent of the data, such that

r ε v ε (t) 2 L 2 + ε v ε (t) 2 L 2 + v ε (t) 2 L 2 + γε p ε (t) 2 L 2 ≤e -αt (1 + r ε) v 0 2 L 2 + γ ε p 0 2 L 2 + λ R e α (1 + r ε) f 2 L ∞ , L 2 , ∀t ∈ R + .
The main goal of this paper is to prove the following convergence theorem: 3) that converges to a solution (v, p) to the system of Navier-Stokes equations with homogeneous Dirichlet boundary conditions.

Theorem 2.3. For d ≤ 3, there exists a subsequence (v ε k , p ε k ) k solution of (
For d = 2, the solution (v, p) is unique, and for all sequences

ε k , (v ε k , p ε k ) ε k converges to (v, p). Moreover, for all sequences ε k , (v ε k ) k converges strongly to v in L 2 (0, T ; H 1 ν (Ω))
. We now give an interpretation of the pressure and precise its convergence. Let us define

q ε = p ε - 1 + ε εR e + r div v ε .
The scalar function q ε appears to be the effective approximate pressure, and we have Theorem 2.4. The function ∇q ε k satisfies

• if d = 3, ∇q ε k converges weakly to ∇p in H -1 ((0, T ) × Ω) 3 • if d = 2, ∇q ε k converges strongly to ∇p in H -1 ((0, T ) × Ω) 2
These convergence results for both velocity and pressure are proved in Section 4.

3. Energy estimates. We first establish the following existence result.

Proposition 3. For v ε 0 , p 0 given in L 2 (Ω) × L 2 0 (Ω)
, there exists at least a solution of the system (2) satisfying for d = 3:

v ε ∈ L ∞ (]0, T [; L 2 (Ω) ∩ L 2 (]0, T [; H 1 ν (Ω)), ∂ v ε ∂t ∈ L 4 3 (]0, T [; (H 1 ν (Ω)) ′ ) v ε ∈ L ∞ (]0, T [; L 2 (Ω) ∩ L 2 (]0, T [; H 1 ν (Ω)), ∂ v ε ∂t ∈ L 4 3 (]0, T [; (H 1 ν (Ω)) ′ ) p ε ∈ L ∞ (]0, T [; L 2 0 (Ω)), ∂p ε ∂t ∈ L 2 (]0, T [; L 2 0 (Ω)) v ε (0) = v ε o in (H 1 ν (Ω)) ′ , v ε (0) = 0 in (H 1 ν (Ω)) ′ , p ε (0) = p 0 in L 2 0 (Ω).
If d = 2, the unique solution of (2) satisfies the following regularity results:

v ε ∈ L ∞ (]0, T [; L 2 (Ω) ∩ L 2 (]0, T [; H 1 ν (Ω)), ∂ v ε ∂t ∈ L 2 (]0, T [; (H 1 ν (Ω)) ′ )) + L 4 3 (]0, T [; L 4 3 (Ω)) v ε ∈ L ∞ (]0, T [; L 2 (Ω) ∩ L 2 (]0, T [; H 1 ν (Ω)), ∂ v ε ∂t ∈ L 2 (]0, T [; (H 1 ν (Ω)) ′ )) + L 4 3 (]0, T [; L 4 3 (Ω)) p ε ∈ L ∞ (]0, T [; L 2 0 (Ω)), ∂p ε ∂t ∈ L 2 (]0, T [; L 2 0 (Ω)) v ε (0) = v ε o in L 2 (Ω), v ε (0) = 0 in L 2 (Ω), p ε (0) = p 0 in L 2 0 (Ω).
Remark 3. In the three-dimensional case, the equalities

v ε (0) = v ε o in (H 1 ν (Ω)) ′ , v ε (0) = 0 in (H 1 ν (Ω)) ′ , are valid in the trace sense.
Proof. For fixed parameters ε > 0, r ≥ 0 and γ > 0, we build approximate solutions by a classical Galerkin process.

Let us introduce the self-adjoint operator A = curl curl -∇ div defined on the domain

H 1 ν (Ω) ∩ (H 2 (Ω)) d .
Then, for the approximation of the two fields of velocity, we use as special basis the eigenfunctions of this operator associated with the following boundary conditions:

u • ν Γ = 0, (curl u) ∧ ν Γ = 0.
For the pressure, one can use as special basis the eigenfunctions of the self-adjoint operator A = -∆ with domain H 2 (Ω) associated to the Neumann boundary conditions. This approximate finite dimensional system is then a classical ordinary differential equation which has a unique solution. Next, to perform the limit we use the same strategy as for the classical Navier-Stokes equations i.e. a priori estimates and compactness results using an estimate on the temporal derivative, see for example [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], [START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF], [START_REF] Boyer | Éléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF]. Now, we will focus our attention on the estimates on the time derivative according to the dimension d. Let us begin with the three-dimensional case. We have to estimate the two nonlinear terms v ε • ∇w, w div v ε , with either w = v ε or w = v ε and the pressure p ε .

Suppose first that d = 3. By Sobolev embedding, the two nonlinear terms of the form v ε • ∇w and w div v ε belong to L 4 3 (0, T ; H 1 ν (Ω) ′ ) since we have, for example for all ϕ ∈ H 1 ν (Ω) :

Ω (v ε • ∇)w • ϕ dω + 1 2 Ω (div v ε ) w • ϕ dω ≤ C v L 3 w H 1 ν ϕ L 6 + w L 3 v ε H 1 ν ϕ L 6 ≤ C v ε 1 2 L 2 v ε 1 2 H 1 ν w H 1 ν + w 1 2 L 2 w 1 2 H 1 ν v ε H 1 ν ϕ H 1 ν .
The bounds of the linear terms are straightforward and we have

Ω curl v ε • curl ϕ dω + Ω div v ε div ϕ dω + r Ω div v ε div ϕ dω ≤ C v ε H 1 ν ϕ H 1 ν Ω p ε div ϕ dω ≤ C p ε L 2 ϕ H 1 ν Ω div v ε + div v ε div ϕ dω ≤ C v ε H 1 ν + v ε H 1 ν ϕ H 1 ν . and, with standard trace theorems Γ (v ε ∧ ν) • (ϕ ∧ ν) dω ≤ C v ε H 1 ν ϕ H 1 ν . Thus it follows, from equation (2) that ∂ v ε ∂t ∈ L 4 3 (]0, T [; (H 1 ν (Ω)) ′ ), ∂ v ε ∂t ∈ L 4 3 (]0, T [; (H 1 ν (Ω)) ′ ) ∂p ε ∂t ∈ L 2 (]0, T [; L 2 0 (Ω))
These estimates show that the velocities ( v ε , v ε ) are equal almost everywhere to continuous functions with values in (H 1 ν (Ω)) ′ . Besides, the pressure p ε is equal almost everywhere to a continuous function with value in L 2 0 (Ω). For the two-dimensional case, the situation is quite different. We observe first that the velocity fields v ε and v ε belong to L 4 (]0, T [; L 4 (Ω)), so that:

v ε ∇) v ε + 1 2 div v ε v ε ∈ L 4 3 (]0, T [; L 4 3 (Ω)), v ε ∇) v ε + 1 2 div v ε v ε ∈ L 4 3 (]0, T [; L 4 3 (Ω)).
So it follows from the equation ( 2) that

∂ v ε ∂t ∈ L 2 (]0, T [; (H 1 ν (Ω)) ′ ) + L 4 3 (]0, T [; L 4 3 (Ω)) ∂ v ε ∂t ∈ L 2 (]0, T [; (H 1 ν (Ω)) ′ ) + L 4 3 (]0, T [; L 4 3 (Ω)) ∂p ε ∂t ∈ L 2 (]0, T [; L 2 0 (Ω))
We now observe that the two velocity fields v ε and v ε belong to

L 4 (]0, T [; L 4 (Ω)) ∩ L 2 (]0, T [; (H 1 ν (Ω)) ′ ) which is the dual space of L 2 (]0, T [; (H 1 ν (Ω)) ′ ) + L 4 3 (]0, T [; L 4 3 (Ω)).
Thus the functions ( v ε , v ε ) are equal almost everywhere to continuous functions with values in L 2 (Ω).

This ends the proof of proposition 3.

3.1. Stability. In the case of three-dimensional vector spaces, we do not have an equality for the conservation of the energy, we have only an inequality. Nevertheless for two-dimensional vector spaces, the weak solutions satisfy the energy equality.

Proof. Through classical computations one obtains, with equations ( 2) and (3):

1 2 d dt v ε 2 L 2 + 1 R e curl v ε 2 L 2 + 1 R e div v ε 2
L 2 + r div v ε 2 L 2 + 1 ε |( v ε ∧ ν)| 2 L 2 (Γ) - Ω p ε div v ε dω = Ω f • v ε dω, (4) 1 2
d dt v ε 2 L 2 + 1 R e curl v ε 2 L 2 + 1 R e div v ε 2
L 2 + r div v ε 2 L 2 + 1 ε |( v ε ∧ ν)| 2 L 2 (Γ) + 1 εR e 0 div v ε div v ε dω = 0, (5) 1 2
d dt v ε 2 L 2 + 1 R e curl v ε 2 L 2 + 1 R e div v ε 2 L 2 + 1 εR e div v ε 2
L 2 + r div v ε 2 L 2 + 1 ε |(v ε ∧ ν)| 2 L 2 (Γ) - Ω p ε div v ε dω = Ω f • v ε dω, (6)
γ ε 2

d dt p ε 2 L 2 + γ ε p ε 2 L 2 + Ω p ε div v ε dω + r ε Ω p ε div v ε dω = 0. (7)
Multiplying (4) by r ε and ( 5) by ε and summing with ( 6) and ( 7), one obtains:

1 2 d dt r ε v ε 2 L 2 + ε v ε 2 L 2 + v ε 2 L 2 + γε p ε 2 L 2 + r ε R e curl v ε 2 L 2 + ε R e curl v ε 2 L 2 + 1 R e curl v ε 2 L 2 + γ ε p ε 2 L 2 + r ε R e div v ε 2 L 2 + ε R e div v ε 2 L 2 + 1 R e div v ε 2 L 2 + r|( v ε ∧ ν)| 2 L 2 (Γ) + |( v ε ∧ ν)| 2 L 2 (Γ) + 1 ε |(v ε ∧ ν)| 2 L 2 (Γ) + εr 2 div v ε 2 L 2 + rε div v ε 2 L 2 + 1 εR e div v ε 2 L 2 = - 1 R e Ω div v ε div v ε dω + r ε Ω f • v ε dω + Ω f • v ε dω,
Let us now give some bounds of the right-hand side terms.

The term 1 R e (div v ε , div v ε ) is bounded by 1 2 ε R e div v ε 2
L 2 + ε 2 R e div v ε 2 L 2 .
According to the estimate of the L 2 norm in H 1 ν (Ω) given by equation ( 1), we bound the source terms in the following way:

Ω f • v d ω ≤ f L 2 v L 2 ≤ √ λ f L 2 div v 2 L 2 + curl v 2 L 2 1 2 ≤ 1 2 R e div v 2 L 2 + curl v 2 L 2 + R e λ 2 f 2 L 2 .
Using these bounds, we get from the previous equation the following fundamental estimate:

1 2 d dt r ε v ε 2 L 2 + ε v ε 2 L 2 + v ε 2 L 2 + γε p ε 2 L 2 + γε p ε 2 L 2 + r ε 2R e curl v ε 2
L 2 + ε R e curl v ε 2 L 2 + 1 2R e curl v ε 2 L 2 + 1 2R e div v ε 2 L 2 + r ε 2R e div v ε 2 L 2 + ε 2R e div v ε 2
L 2 + εr 2 div v ε 2 L 2 + εr div v ε 2 L 2 + 1 2εR e div v ε 2 L 2 + r |( v ε ∧ ν)| 2 L 2 (Γ) + |( v ε ∧ ν)| 2 L 2 (Γ) + 1 ε |(v ε ∧ ν)| 2 L 2 (Γ) ≤ λ R e 2 (1 + ε r) f 2 L 2 . (8)
After integration in time, we deduce from the previous estimate that there exists a continuous function g defined on [0, T ] such that: for all t > 0,

1 2 r ε v ε (t) 2 + ε v ε (t) 2 + v ε (t) 2 L 2 + γε p ε (t) 2 L 2 + γ ε t 0 p ε (s) 2 L 2 ds + rε 2R e t 0 curl v ε (s) 2 L 2 ds + ε R e t 0 curl v ε (s) 2 L 2 ds + 1 2R e t 0 curl v ε (s) 2 L 2 ds + rε 2R e t 0 div v ε (s) 2 L 2 ds + ε 2R e t 0 div v ε (s) 2 L 2 ds + 1 2R e t 0 div v ε (s) 2 L 2 ds + r t 0 |( v ε (s) ∧ ν)| 2 L 2 (Γ) ds + t 0 |( v ε (s) ∧ ν)| 2 L 2 (Γ) ds + 1 ε t 0 |(v ε (s) ∧ ν)| 2 L 2 (Γ) ds + εr 2 t 0 div v ε (s) 2 L 2 ds + εr t 0 div v ε (s) 2 L 2 ds + 1 2εR e t 0 div v ε (s) 2 L 2 ds ≤ g(t), (9)
with

g(t) = r ε v ε (0) 2 +ε v ε (0) 2 + v ε (0) 2 L 2 +γε p ε (0) 2 L 2 + λR e 2 (1+rε) t 0 f (s) 2 L 2 ds.
This inequality is the key point to establish the convergence result.

To improve the convergence result in the two-dimensional case, we use the following energy equality derived as above without using (5).

1 2 d dt r ε v ε 2 L 2 + v ε 2 L 2 + γε p ε 2 L 2 + r ε R e curl v ε 2 L 2 + 1 R e curl v ε 2 L 2 + γ ε p ε 2 L 2 + r ε R e div v ε 2 L 2 + 1 R e div v ε 2 L 2 + εr 2 div v ε 2 L 2 + 1 εR e div v ε 2 L 2 + r |( v ε ∧ ν)| 2 L 2 (Γ) + 1 ε |(v ε ∧ ν)| 2 L 2 (Γ) = r ε Ω f • v ε dω + Ω f • v ε dω. (10)
This concludes the proof of lemma 2.1.

3.2.

Uniform stability for the approximate solution. In this section, we deal with the stability of the proposed approximation method. We notice that this property is valid for all solutions satisfying the energy inequality (8), as it is the case when they are built by a finite dimensional approximation method such as the Galerkin method for example.

Proof. Let us write

χ ε (t) = ε r v ε (t) 2 L 2 + ε v ε (t) 2 L 2 + v ε (t) 2 L 2 + γε p ε (t) 2 L 2 .
We note λ > 0 the smallest eigenvalue of the self-adjoint operator A = curl curl -∇ div with the domain

D = H 1 ν (Ω) ∩ (H 2 (Ω)) d ,
and we introduce α = min λ Re , 1 . Classically, the inequality (8) leads to the differential inequality

d dt χ ε (t) + αχ ε (t) ≤ λ R e (1 + r ε) f (t) 2 L 2 , which implies the following uniform bound χ ε (t) ≤ e -αt χ ε (0) + λ R e α (1 + r ε) f L ∞ ,L 2 .
This concludes the proof of theorem 2.2.

4. Convergence analysis and compactness results.

4.1.

Compactness results for the velocity. Let us introduce the Leray projection w ε of a velocity field v ε (t) ∈ H 1 ν (Ω) defined as follows

v ε = w ε + ∇q ε , div w ε = 0, w ε • ν |Γ = 0, ∇q ε • ν |Γ = 0, Ω q ε dω = 0.
By the estimate (9), we see that the irrotational part of v ε goes to zero with ε. Thus it remains to bring to the fore the behavior of the free divergence part w ε and to obtain an estimate on a fractional time derivative of this term. We detail the different steps of this strategy.

From the regularity of the Leray projector (see R. Temam [START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF] page 18), one has:

w ε L ∞ ,L 2 ≤ c v ε L ∞ ,L 2 , w ε L 2 ,H 1 ≤ c v ε L 2 ,H 1 . ( 11
)
Moreover, we can easily prove the following lemma.

Lemma 4.1. There exists two constants depending only on T and Ω such that:

∇q ε L 2 ,H 1 ≤ c √ ε, ∇q ε L ∞ ,L 2 ≤ c (12)
Proof: The function q ε belongs to H 2 (Ω) and satisfies

-∆q ε (t) = -div v ε (t), ∇q ε (t) • ν |Γ = 0.
This implies using the estimate ( 9)

∆q ε (t) L 2 ,L 2 = div v ε L 2 ,L 2 ≤ C √ ε.
Besides, we have

Ω ∇q ε • ∇q ε dω = - Ω ∆q ε q ε dω,
so that, with Poincaré-Neumann inequality, we get

∇q ε (t) 2 L 2 ,L 2 ≤ C ∆q ε (t) L 2 ,L 2 q ε L 2 ,L 2 , ≤ C √ ε ∇q ε (t) L 2 ,L 2 .
The regularity properties of the Neumann problem give

∇q ε L 2 ,H 1 ≤ C q ε L 2 ,H 2 ≤ C ∆q ε L 2 ,L 2 ≤ C √ ε. (13)
Moreover, by orthogonality of the Leray projector in L 2 , one has

∇q ε L ∞ ,L 2 ≤ v ε L ∞ ,L 2 ≤ C. ( 14
)
This concludes the proof of the lemma 4.1.

So by interpolation and using estimates ( 13)-( 14), we have proved the result below.

Corollary 1. The function q ε satisfies:

∇q ε strongly converges to 0 in L p (]0, T [; L 2 (Ω)) , ∀p, 1 ≤ p < +∞.
Now we have to write the equation satisfied by w ε . As the Leray projection is orthogonal in L 2 (Ω), this equation reads

                     ∀ϕ ∈ H 1 ν (Ω), div ϕ = 0, Ω ∂w ε ∂t • ϕ dω + Ω (v ε • ∇)v ε + 1 2 (div v ε )v ε • ϕ dω + 1 R e Ω curl v ε • curl ϕ dω + 1 ε Γ (v ε ∧ ν) • (ϕ ∧ ν) dγ = Ω f • ϕ dω in L 1 (0, T ). ( 15
)
Now we introduce the extension by 0 of w ε (resp. v ε ) outside [0, T ] denoted, only in this part, by w ε (resp. v ε ) and we take the Fourier transform in time of the equation ( 15) to obtain

                     ∀ϕ ∈ H 1 ν (Ω), div ϕ = 0, iτ Ω F( w ε )(τ ) • ϕ dω + Ω F ( v ε • ∇) v ε + 1 2 (div v ε ) v ε (τ ) • ϕ dω + 1 R e Ω curl F( v ε )(τ ) • curl ϕ dω + 1 ε Γ F( v ε ∧ ν)(τ ) • (ϕ ∧ ν) dγ = Ω F( f )(τ ) • ϕ dω + 1 √ 2π Ω v ε (0) • ϕ dω - e -iτ T √ 2π Ω v ε (T ) • ϕ dω.
Following Boyer-Fabrie [5, page 253], we take ϕ = F( w ε )(τ ) as test function in the previous equation to obtain for all τ ∈ R:

iτ Ω |F( w ε )(τ )| 2 dω = - Ω F( v ε • ∇) v ε (τ ) • F( w ε )(τ ) dω - 1 2 Ω F(div v ε ) v ε (τ ) • F( w ε )(τ ) dω - 1 R e Ω curl F( v ε )(τ ) • curl F( w ε )(τ ) dω - 1 ε Γ F( v ε ∧ ν)(τ ) • (F( w ε )(τ ) ∧ ν) dγ + Ω F( f )(τ ) • F( w ε )(τ ) dω + 1 √ 2π Ω v ε (0) • F( w ε )(τ ) dω - e -iτ T √ 2π Ω v ε (T ) • F( w ε )(τ ) dω.
As we look for an estimate independent of ε, we have to pay a special attention to the imaginary part of the penalty term:

A ε = - 1 ε Γ F( v ε ∧ ν)(τ ) • (F( w ε )(τ ) ∧ ν) dγ. ( 16
)
By writing w ε = v ε -∇q ε , we have:

1 ε Γ F( v ε ∧ ν)(τ )•(F( w ε )(τ ) ∧ ν) dγ = 1 ε Γ |F( v ε )(τ ) ∧ ν| 2 dγ - 1 ε Γ F(∇ q ε )(τ ) ∧ ν • F( v ε )(τ ) ∧ ν) dγ.
So, the imaginary part of A ε is bounded as follows:

1 ε Γ F(∇ q ε )(τ ) ∧ ν • F( v ε )(τ ) ∧ ν) dγ ≤ 1 ε |F(∇ q ε )(τ )| L 2 (Γ) |F( v ε )(τ ) ∧ ν| L 2 (Γ) , ≤ C 1 ε F(∇ q ε )(τ ) H 1 |F( v ε )(τ ) ∧ ν| L 2 (Γ) . (17)
From estimates ( 12) and ( 9), we have

∇q ε L2,H 1 ≤ C √ ε, |v ε ∧ ν| L 2 ,L 2 (Γ) ≤ C √ ε.
So there exists a function f 4 ε (τ ) ∈ L 1 (R) bounded independently on ε such that:

1 ε Γ F(∇ q ε )(τ ) ∧ ν • F( v ε )(τ ) ∧ ν) dγ ≤ f 4 ε (τ ) (18)
Now we can derive the estimate of |τ | Ω |F( w ε )(τ )| 2 dω, and we have

|τ | Ω |F( w ε )(τ )| 2 dω ≤ Ω F( v ε • ∇) v ε (τ ) • F( w ε )(τ )dω + 1 2 Ω F(div v ε ) v ε (τ ) • F( w ε )(τ )dω + 1 R e Ω curl F( v ε )(τ ) • curl F( w ε )(τ )dω + 1 ε Γ F(∇ q ε )(τ ) ∧ ν • F( v ε )(τ ) ∧ ν) dγ + Ω F(f )(τ ) • F( w ε )(τ ) dω + 1 √ 2π Ω v ε (0) • F( w ε )(τ )dω + 1 √ 2π Ω v ε (T ) • F( w ε )(τ )dω ≤ f 1 ε (τ ) + f 2 ε (τ ) + f 3 ε (τ ) + f 4 ε (τ ) + f 5 ε (τ ) + f 6 ε (τ ) + f 7 ε (τ ). (19)
We now estimate each term of the right-hand side of the previous inequality for d ≤ 3.

Term f 1 ε = Ω F( v ε • ∇) v ε (τ ) • F( w ε )(τ )dω
According to the energy estimate (9), the function v ε is bounded in L 2 (0, T ; H 1 ν (Ω)) and hence, by Sobolev injection, it is bounded in L 6 5 (0, T ; L 6 (Ω)). So by Hausdorff-Young theorem, (F(

w ε )) ε is bounded in L 6 (R; L 6 (Ω)). ( 20
)
We also have the inequality:

v ε • ∇ v ε L 6 5 ≤ ∇v ε L 2 v ε L 3 ≤ C ∇v ε 3 2 L 2 v ε 1 2 L 2 ,
which implies, according to [START_REF] Girault | Finite Element Methods for the Navier-Stokes Equations[END_REF], that v ε • ∇v ε is bounded in L

F( v ε • ∇) v ε is bounded in L 6 (R; L 6 5 (Ω)). Then with Hölder inequality, (f 1 ε ) ε is bounded in L 3 (R). ( 21
)
Term f 2 ε = 1 2 Ω F(div v ε ) v ε (τ ) • F( w ε )(τ )dω
The same arguments show that

(f 2 ε ) ε is bounded in L 3 (R). ( 22
)
Term f 3 ε = 1 R e Ω curl F( v ε )(τ ) • curl F( w ε )(τ )dω
According to the regularity of the Leray projection recalled above and estimate (9), one has:

(f 3 ε ) ε is bounded in L 1 (R). ( 23
)
Terms f 4 ε = 1 ε Γ F(∇ q ε )(τ ) ∧ ν • F( v ε )(τ ) ∧ ν) dγ
As we have seen by the estimate [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires I[END_REF],

(f 4 ε ) ε is bounded in L 1 (R). ( 24
)
Term f 5 ε = Ω F( f )(τ ) • F( w ε )(τ ) dω
By hypothesis, f is a given function in L 2 (R; L 2 (Ω)) and from estimate (9), we get that v ε is bounded in L 2 (R; L 2 (Ω)), so as w ε is the Leray projection of v ε , the function w ε is also bounded in L 2 (R; L 2 (Ω)). Finally, we obtain

(f 5 ε ) ε is bounded in L 1 (R). ( 25
)
Terms f 6 ε = 1 √ 2π Ω v ε (0) • F( w ε )(τ )dω and f 7 ε = 1 √ 2π Ω v ε (T ) • F( w ε )(τ )dω
These two terms come from the Dirac measure when we derive discontinuous functions. Let us consider f 7 ε .

f 7 ε (τ ) ≤ 1 √ 2π v ε (T ) L 2 F( w ε )(τ ) L 2 .
According to [START_REF] Girault | Finite Element Methods for the Navier-Stokes Equations[END_REF], this term is bounded in L ∞ (R). Moreover, the set of functions

F( w ε )(τ ) ε is bounded in L 2 (R; L 2 (Ω)) so, f 7 ε ε is bounded in L 2 (R).
We treat in the same way the term f 6 ε , and thus:

(f 6 ε ) ε and (f 7 ε ) ε are bounded in L 2 (R) (26)
We are now able to show that the set of functions ( w ε ) ε is bounded in an appropriate Nikolskii space.

For all γ < 1, there exists a constant d such that:

|τ | 1-γ ≤ d 1 + |τ | 1 + |τ | γ , so that, |τ | 1-γ F( w ε )(τ ) 2 L 2 ≤ d F( w ε )(τ ) 2 L 2 + |τ | 1 + |τ | γ F( w ε )(τ ) 2 L 2 Let us denote f 0 ε (τ ) = F( w ε )(τ ) 2 L 2
, which belongs to L 1 (R), the previous inequality reads with [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF]:

|τ | 1-γ F( w ε )(τ ) 2 L 2 ≤ d f 0 ε (τ ) + d 1 + |τ | γ f 1 ε (τ ) + f 2 ε (τ ) + f 3 ε (τ ) + f 4 ε (τ ) + f 5 ε (τ ) + f 6 ε (τ ) + f 7 ε (τ ) ≤ h(τ )
If we suppose that the function τ → 1 1 + |τ | γ belongs to L ∞ (R) ∩ L 2 (R), then the function τ → h(τ ) belongs to L 1 (R). This condition is satisfied for γ ∈] 2 3 , 1[. So we have proved : Lemma 4.2. Let us suppose that σ ∈]0, 1 6 [, then there exists a constant C such that

R |τ | 2σ F(w ε )(τ ) 2 L 2 dτ ≤ C. (27)
Then, from lemma 4.1 and 4.2 we deduce the following key result:

Theorem 4.3. There exists a sequence (ε k ) k which converges to zero and a function v ∈ L 2 (]0, T [; L 2 (Ω)) satisfying div v = 0 such that:

(v ε k ) k → v in L 2 (]0, T [; L 2 (Ω)) strongly.
Proof. The function v ε is the sum of two terms ∇q ε and w ε . From corollary 1 the first term converges strongly to 0 in L 2 (]0, T [; L 2 (Ω)). Now, from Aubin-Lions-Simon Theorem, it follows from lemma 4.2, that there exists a sequence (ε k ) k such that:

(w ε k ) k -→ v in L 2 (]0, T [; L 2 (Ω)) strongly .
Moreover, since div w ε k = 0, we have div v = 0.

4.2.

Convergence of the method. We first give a general convergence theorem for a subsequence solution of the approximate scheme (3), to a weak solution of the initial Navier-Stokes problem, in the case d ≤ 3. For the two-dimensional case, since the weak solution of the Navier-Stokes equation is unique, the whole sequence of approximate solution v ε converges to v. Moreover, in this case, we prove that the convergence is strong.

- T 0 Ω v ε • ϕ dω θ ′ (τ ) dτ + T 0 Ω (v ε • ∇)v ε + 1 2 (div v ε ) v ε • ϕ dω θ(τ ) dτ + 1 R e T 0 Ω curl v ε • curl ϕ dω θ(τ ) dτ = T 0 Ω f ϕ dω θ(τ ) dτ + Ω v 0 ϕ dω θ(0). (28)
According to estimate [START_REF] Girault | Finite Element Methods for the Navier-Stokes Equations[END_REF], there exists a sequence ε k such that

v ε k → v in L 2 (0, T ; L 2 (Ω)) strongly, div v ε k → div v = 0 in L 2 (0, T ; L 2 (Ω)) strongly, curl v ε k ⇀ curl v in L 2 (0, T ; L 2 (Ω)) weakly. Since u 2 L 2 + div u 2 L 2 + curl u 2 L 2 1
2 is a norm equivalent to the H 1 -norm on H 1 ν (Ω), we have

∇v ε k ⇀ ∇v in L 2 (0, T ; L 2 (Ω))
weakly. and so, we can take the limit on the term (v ε • ∇)v ε as

v ε k • ∇v ε k → v • ∇v in L 1 (0, T ; L 1 (Ω)) div v ε k v ε k → 0 in L 1 (0, T ; L 1 (Ω)).
From estimate (9), we have for the tangential traces

t 0 |(v ε ∧ ν)(τ )| 2 L 2 (Γ) dτ ≤ ε g(t),
and since for any function v in

H 1 (Ω) d , |v| L 2 (Γ) ≤ C v 1 2 L 2 v 1 2 H 1 , we obtain that v ε k → v in L 2 ]0, T [; (L 2 (Γ)) d strongly. This implies (v ∧ ν) |Γ = 0, and so, since by construction (v • ν) |Γ = 0, v belongs to H 1 0 (Ω) d Finally, at the limit process we obtain                    - T 0 Ω v • ϕ dω θ ′ (τ ) dτ + T 0 Ω (v • ∇)v • ϕ dω θ(τ ) dτ + 1 R e T 0 Ω curl v • curl ϕ dω θ(τ ) dτ = T 0 Ω f • ϕ dω θ(τ ) dτ + Ω v 0 ϕ dω θ(0), div v = 0, v |Γ = 0.
From the identity

Ω ∇v : ∇ϕ dω = Ω curl v • curl ϕ dω + Ω div v • div ϕ dω,
which is valid for all functions (v, ϕ) ∈ (H 1 0 (Ω)) d × (H 1 0 (Ω)) d , the previous equality shows that the limit function v satisfies the classical Navier-Stokes equations in a weak sense. 4.2.2. The special case d = 2. The key point to establish the strong convergence in the two-dimensional case, lies on an idea of R. Temam [START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF]. It is based on the fact that in this case, the solution of the approximate problem and the solution of the Navier-Stokes equation verify the equality of energy. The idea is to bring to the fore an energy equation satisfied by the difference between the approximate solution and the exact solution.

We first observe that, according to the equality -∆ = curl curl -∇ div the classical weak solution v of the Navier-Stokes equation satisfies for any test function in (H 1 0 (Ω)) d , with free-divergence:

Ω ∂v ∂t • ϕ dω + Ω (v • ∇)v • ϕ dω + 1 R e Ω curl v • curl ϕ dω = Ω f • ϕ dω.
The equation satisfied by the error v

ε -v = u ε with a free-divergence test function ϕ in H 1 0 (Ω) d ⊂ H 1 ν (Ω).        Ω ∂u ε ∂t • ϕ dω + Ω (u ε • ∇)u ε + 1 2 (div u ε ) u ε • ϕ dω + Ω (v • ∇)u ε • ϕ dω + Ω (u ε • ∇)v + 1 2 (div u ε ) v • ϕ dω + 1 R e Ω curl u ε • curl ϕ dω = 0.
After integration in time, this equation gives

                 Ω u ε (t) • ϕ dω + t 0 Ω (u ε • ∇)u ε + 1 2 (div u ε ) u ε • ϕ dω dτ + t 0 Ω (v • ∇)u ε • ϕ dω dτ + t 0 Ω (u ε • ∇)v + 1 2 (div u ε ) v • ϕ dω dτ + 1 R e t 0 Ω curl u ε • curl ϕ dω dτ = Ω u ε (0) • ϕ dω.
Taking the limit when ε goes to 0, one obtains with the convergences properties stated in the previous section:

Lemma 4.4. ∀ϕ ∈ (H 1 0 (Ω)) d , div ϕ = 0, lim ε→0 Ω u ε (t) • ϕ dω = 0. ( 29
)
Following R. Temam [START_REF] Temam | Navier-Stokes Equations; Theory and Numerical Analysis[END_REF], we introduce

χ ε (t) = 1 2 ε r v ε (t) 2 L 2 + γε p ε (t) 2 L 2 + (v ε -v)(t) 2 L 2 + r ε R e t 0 curl v ε (τ ) 2 L 2 dτ + γε t 0 p ε (τ ) 2 L 2 dτ + 1 R e t 0 curl (v ε -v)(τ ) 2 L 2 dτ + 1 R e t 0 div (v ε -v)(τ ) 2 L 2 dτ +r t 0 |( v ε ∧ ν)(τ )| 2 L 2 (Γ) dτ + 1 ε t 0 |(v ε ∧ ν)(τ )| 2 L 2 (Γ) dτ + εr R e t 0 div v ε (τ ) 2 L 2 dτ + 1 εR e t 0 div v ε (τ ) 2 L 2 dτ.
Here v is the unique solution of the two-dimensional Navier-Stokes equations. Due to the energy equality [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], χ ε (t) satisfies

χ ε (t) = 1 2 ε r v ε (0) 2 L 2 + v ε (0) 2 L 2 + γ ε p ε (0) 2 L 2 + rε t 0 Ω f (τ ) • v ε (τ )dω dτ + t 0 Ω f (τ ) • v ε (τ )dω dτ - Ω v ε (t) • v(t)dω - 2 R e t 0 Ω curl v ε (τ ) • curl v(τ )dω dτ - 2 R e t 0 Ω div v ε (τ )div v(τ )dω dτ + 1 2 v(t) 2 L 2 + 1 R e t 0 curl v(τ ) 2 L 2 dτ + 1 R e t 0 div v(τ ) 2 L 2 dτ.
By weak convergence in L 2 (]0, t[; H 1 0 (Ω)) of the sequence (v ε ) ε and from the lemma 4.4, we observe that :

lim ε→0 1 2 v ε (0) 2 L 2 - Ω v ε (t) • v(t)dω + t 0 Ω f (τ ) • v ε (τ ) dω dτ - 2 R e t 0 Ω curl v ε (τ ) • curl v(τ )dω dτ - 2 R e t 0 Ω div v ε (τ ) • div v(τ )dω dτ = 1 2 v(0) 2 L 2 + t 0 Ω f (τ ) • v(τ ) dω dτ -v(t) 2 L 2 - 2 R e t 0 curl v(τ ) 2 L 2 dτ - 2 R e t 0 div v(τ ) 2 L 2 dτ.
In the two-dimensional case, the unique solution of the Navier-Stokes equation satisfies the following energy equality We are now able to precise the convergence for the effective pressure and establish the theorem 2.4

Let us write the equation satisfied by the velocity v ε and the pressure p ε in the distribution sense. We have

∂v ε ∂t + (v ε • ∇)v ε + 1 2 (div v ε )v ε + 1 R e curl curl v ε + ∇ p ε - 1 + ε εR e + r div v ε = f.
Introducing the effective pressure

q ε = p ε - 1 + ε εR e + r div v ε ,
this equation reads

∇q ε = f - ∂v ε ∂t + (v ε • ∇)v ε + 1 2 (div v ε )v ε + 1 R e curl curl v ε
and the proof follows, from the previous steps.

5.

Appendix. Let us consider the following problem:

Proposition 4. For a fixed ε > 0 and a couple of functions (f, g) given in L 2 0 (Ω) × L 2 (Ω) , there exists a unique solution v ε ∈ {w ∈ H div (Ω), w • ν |Γ = 0} solution of:

εv ε -∇div v ε = ∇f + εg. ( 30
)
Moreover, if (u, u 1 ) ∈ H div (Ω) × H is solution of

div u = -f, curl u = 0, u • ν |Γ = 0, (31)
div u 1 = 0, curl u 1 = curl g, u 1 • ν |Γ = 0, and we have the following estimate

v ε -u -u 1 H 1 ≤ ε u + u 1 -g L 2 .
Proof.

Step 1: Existence of v ε Let us note H div,0 (Ω) = {v ∈ L 2 (Ω), div v ∈ L 2 , v • ν |Γ = 0}. The existence of a unique solution to the equation ( 30) is obtained by a straightforward application of the Lax-Milgram theorem with the bilinear form defined on H div,0 (Ω) by ε(u, v) + (div u, div v), and the right-hand side: -(f, div v) + ε(g, v).

Step 2: existence of u and u 1

The existence of u satisfying (31) comes from the resolution of the following Neumann problem

-∆q = f, ∇q • ν |Γ = 0,
and we set u = ∇q, with q ∈ H 1 (Ω)/R ∩ H 2 (Ω). The existence of u 1 ∈ H is the consequence of the Leray projection applied to g by writing g = u 1 + ∇p, div u 1 = 0, u 1 • ν |Γ = 0. Now, writting v ε = u + u 1 + u ε , we get that u ε ∈ H div,0 satisfies εu ε -∇div u ε = -εu -εu 1 + εg, curl u ε = 0,

6 5

 6 (Ω)). So, by Hausdorff-Young theorem, the family of functions

4. 2 . 1 .

 21 The general case d ≤ 3. Let θ an element of C ∞ (0, T ), satisfying θ(T ) = 0, and ϕ a free-divergence vector field in H 1 0 (Ω) d ∩ H 2 (Ω)) d . An integration by parts gives from the equation[START_REF] Ph | A new fast method to compute saddle-points in constrained optimization and applications[END_REF]

ff

  (τ ) • v(τ ) dω dτ. (τ ) • v ε (τ ) dω dτ = 0. So, we have proved that lim ε→0 χ ε (t) = 0.In other words, we have established the following resultv ε → v in C 0 ([0, T ]; L 2 (Ω)) ∩ L 2 (]0, T [; H 1 ν (Ω)) v ∈ C 0 ([0, T ]; H) ∩ L 2 (]0, T [; H 1 0 (Ω))This concludes the proof of theorem 2.3.

  1. Aubin-Lions-Simon Let B 0 , B 1 , B 2 three Banach spaces with B 0 ⊂ B 1 ⊂ B 2 with continuous imbedding. Suppose moreover that the injection of B 0 in B 1 is compact.