G Angloher

M Bauer

I Bavykina

A Bento

A Brown

C Bucci

C Ciemniak

C Coppi

G Deuter

F Von Feilitzsch

D Hauff

S Henry

P Huff

J Imber

S Ingleby

C Isaila

J Jochum

M Kiefer

M Kimmerle

H Kraus

J.-C Lanfranchi

R F Lang

B Majorovits

M Malek

R Mcgowan

V B Mikhailik

E Pantic

F Petricca

S Pfister

W Potzel

F Pröbst

W Rau

S Roth

K Rottler

C Sailer

K Schäffner

J Schmaler

S Scholl

W Seidel

L Stodolsky

A J B Tolhurst

I Usherov

W Westphal

Commissioning Run of the CRESST-II Dark Matter Search

Keywords: Dark Matter, WIMP, Low-temperature detectors, CaWO 4

The CRESST cryogenic direct dark matter search at Gran Sasso, searching for WIMPs via nuclear recoil, has been upgraded to CRESST-II by several changes and improvements. The upgrade includes a new detector support structure capable of accommodating 33 modules, the associated multichannel readout with 66 SQUID channels, a neutron shield, a calibration source lift, and the installation of a muon veto. We present the results of a commissioning run carried out in 2007.

The basic element of CRESST-II is a detector module consisting of a large (∼300 g)CaWO 4 crystal and a very sensitive smaller (∼ 2 g) light detector to detect the

scintillation light from the CaWO 4 .The large crystal gives an accurate total energy measurement. The light detector permits a determination of the light yield for an event, allowing an effective separation of nuclear recoils from electron-photon backgrounds. Furthermore, information from light-quenching factor studies allows the definition of a region of the energy-light yield plane which corresponds to tungsten recoils. A neutron test is reported which supports the principle of using the light yield to identify the recoiling nucleus.

Data obtained with two detector modules for a total exposure of 48 kg-days are presented. Judging by the rate of events in the "all nuclear recoils" acceptance region the apparatus shows a factor ∼ten improvement with respect to previous results, which we attribute principally to the presence of the neutron shield. In the "tungsten recoils" acceptance region three events are found, corresponding to a rate of 0.063 per kg-day. Standard assumptions on the dark matter flux, coherent or spin independent interactions, then yield a limit for WIMP-nucleon scattering of 4.8 × 10 -7 pb, at M WIMP ∼ 50 GeV.

Introduction

Evidence continues to grow that the majority of the mass of our galaxy and in the universe is made up of non-baryonic dark matter [1]. Precision measurements of the cosmic microwave background have given accurate figures for the density of matter and energy in the universe as a whole, suggesting that about a fourth of the energy density of the universe is in the form of dark matter [2]. Measurements of the rotation curves of other spiral galaxies, indicate large amounts of dark matter, which would suggest it also dominates our own galaxy. Gravitational lensing of light by galactic clusters indicates large amounts of dark matter, and dark matter is necessary for a reasonable description of the formation of galaxies.

Although some of these observations might also be explained by alternative theories such as modified Newtonian dynamics (MOND) [3], recent studies of the Bullet Cluster appear to favour the dark matter hypothesis [4]. However the need for direct detection of dark matter evidently remains strong.

A plausible origin for dark matter comes from particle physics in the form of WIMPs (weakly interacting massive particles). These would be stable massive particles with an interaction cross section typical for weak interaction processes. A relic density sufficient to make a significant contribution to the energy density of the universe arises naturally in this way. Supersymmetry provides a well motivated candidate in the form of the neutralino and this has been extensively studied theoretically.

According to the WIMP hypothesis a galaxy has a large gravitationally bound halo of such particles, making up most of the mass. The local mass density may be estimated from halo models and is found to be on the order of 0.3 GeVcm -3 . The range of the coherent or spin-independent WIMP-nucleon scattering crosssections predicted by minimal supersymmetric models extends from below 10 -10 to above 10 -7 pb [5,6]; this makes the direct detection of such particles a difficult but in principle achievable task.

The CRESST project is a dark matter search aiming to detect the nuclear scattering of WIMPs by the use of cryogenic detectors. Since only small recoil energies (∼ keV's) are anticipated in WIMP-nucleus scattering, cryogenic detectors with their high sensitivity are well suited to the problem.

CRESST-II is an upgrade of the original CRESST appartus where the detectors are arranged in modules, each consisting of two detectors, a large detector comprising the target mass and a similar but much smaller light detector measuring the light yield. While the large detector gives an accurate total energy measurement, the light yield measurement permits an effective rejection of events which are not nuclear recoils. This arises from the property of electronphoton events, constituting the dominant background, to give a much higher light output than nuclear recoils.

The large detector consists of a CaWO 4 crystal of ∼ 300 g with a tungsten superconducting-to-normal thermometer deposited on the surface. The energy of an event is measured in this detector. The energy deposited in the crystal is quickly converted into a gas of phonons which are then absorbed in the superconducting thermometer. The separate light detector, based on the same principle, measures the simultaneously emitted scintillation light. A comparison of these two signals then allows an effective background discrimination.

Results from an earlier prototype run with two modules were presented in Ref. [7], where taking the better of the two modules and assuming coherent WIMP scattering on the tungsten led to the limit 1.6 × 10 -6 pb. The improvements discussed here are aimed at ultimately increasing this sensitivity by up to two orders of magnitude. Such an increase would result, for example, with background free data and a threshold of 12 keV from 10 kg for 100 days giving 1000 kg-days. Or if the threshold could be reduced to 5 keV, with 300 kg days.

During 2007 an extended commissioning run of the new setup was carried out. Although this was primarily for optimization purposes, the results also represent an interesting further improvement. We report on them and the modifications for CRESST-II in the present paper. 3 He-4 He dilution refrigerator, low background detector environment, and the shielding as upgraded for CRESST-II. The detector carousel (CA) is connected to the mixing chamber of the cryostat (CR) by a long copper cold finger (CF) in order to reduce background originating from the dilution refrigerator. The gas-tight radon box (RB) encloses the low background copper (CU) and low background lead shielding (PB). It is covered by a plastic scintillator muon-veto (MV) and a 45 cm thick polyethylene neutron moderator (PE). Additional granular PE is placed between the baffles in the upper part of the cryostat to close the line of sight for neutrons coming from above.

Experimental Setup

The detector volume of the CRESST installation at Gran Sasso is surrounded by passive shielding consisting of 14 cm thick low background copper and 20 cm of lead. This inner shielding is entirely enclosed within a gas-tight radon box which is continuously flushed with N 2 gas. This part of the setup is as used in previous phases of the experiment and is described in Ref. [8]. The changes made in the upgrade to CRESST-II include the installation of a neutron shield and a muon veto. Figure 1 shows the setup after the upgrade.

In the detector volume itself, a new detector support structure ("carousel") capable of holding the 33 modules has been mounted and a new 66-channel SQUID readout with associated wiring and data acquisition electronics has been installed. The various changes will be discussed in this section.

Detector Support Structure ("Carousel")

The new support structure for holding the detector modules is made from copper suitable for milliKelvin (mK) operation and electropolished to reduce surface contamination. The structure, which is depicted in Fig. 2, can accommodate a total of 33 detector modules, thus allowing an increase of the target mass by an order of magnitude (to ∼ 10 kg) relative to CRESST-I. Individual handling and dismounting of the modules is possible. The structure rests on custom-made CuSn6 springs [START_REF] Majorovits | submitted to Applied Radiation and Isotopes[END_REF] to reduce interference from vibrations. CuSn6 was chosen since it exhibits both intrinsic radiopurity and elasticity at mK temperatures. The thermal coupling of the carousel to the cold finger is intentionally set to be weak, with a relaxation time of about half an hour. This filters out high-frequency temperature variations of the cryostat, giving a high degree of temperature stability at the detectors.

Sixty-six channel SQUID readout

With the full complement of 33 modules, each with two superconducting thermometers, there will be 66 readout channels in CRESST-II. In the basic CRESST design these channels are read out with SQUID's. Figure 3 shows one such readout circuit, simplified for clarity. The actual implementation is discussed in detail in Ref. [START_REF] Henry | The 66-channel SQUID readout for CRESST-II[END_REF]. A constant bias current is shared between two branches, one containing the tungsten thermometer, the other the input coil for a SQUID in series with two reference resistors. The reference resistors assure proper branching of the current; two are used here to make the circuit more balanced and to minimize the level of crosstalk. Any change in the resistance of the thermometer changes the current through the SQUID input coil. The SQUID system then transforms this signal to a voltage pulse, which is processed by conventional follow-on electronics. The SQUIDs are DC SQUIDs using AC modulated flux-locked loops. Full details of the system are given in [START_REF] Henry | The 66-channel SQUID readout for CRESST-II[END_REF].

To handle the large number of detectors, a multichannel system was installed, consisting of 66 SQUID units. These are located at the bottom of the main helium bath, at 4 K. Woven twisted pair cables are used to connect the SQUIDs to the detectors in the mK environment, while another set of cables leads to connector boxes at room temperature, which then link to the control electronics and the data aquisition system.

In addition to the readout-line pair for the SQUID, each thermometer has a heater-line pair for temperature control and calibration pulses. Thus there are two pairs of lines for each thermometer and so a total of eight lines per module. With the many lines into the mK environment close attention must be paid in the design and fabrication of flanges and feed-throughs [10][11].

Neutron Shield

A significant new element in the setup is a neutron shield, consisting of a 45 cm thickness of polyethylene (PE) plates. The shield is mounted on rails to allow access to the inner parts of the experiment. The polyethylene moderates ambient neutrons to thermal energies where they are then captured on the hydrogen. The neutron capture range in CH 2 materials is 4.3 cm for MeV energies [12], so the shield should eliminate essentially all incoming neutrons.

Calibration Source Lift

During a run, CRESST detectors are regularly calibrated with a radioactive source supplying gamma quanta of known energy. In view of the many and closely spaced detectors, in CRESST-II it is necessary to be able to move and position the radioactive source during calibration runs.

Therefore the apparatus has been equipped with a calibration source lift. The lift consists of a plastic tube through which the source is moved by compressed air, the position being determined by strings attached from both ends of the tube. It spirals around the cold box inside the shielding and allows calibration of individual detector modules without opening the shielding. It is designed so that a rate of ∼ 0.5 Hz can be generated at any detector module with a 0.7 MBq 57 Co source. The system was used successfully several times during the commissioning run.

Muon Veto

With the aim of reducing the effects of muons entering the apparatus a muon veto has been installed. It consists of 20 plastic scintillator panels of area about 1-1.2 m 2 surrounding the radon box and the Cu/Pb shielding. By using three different types of panels, a total solid angle coverage of 98.7 % is achieved. A circular opening of 0.27 m 2 on top was unavoidable due to the entry for the cryostat (see Fig. 1). In a dedicated setup [13], the characteristics of each panel type were measured and calibrated. The photomultiplier signals are sent out of the Faraday cage optically to avoid interference with the cryogenic detectors. The veto triggered at ∼ 5 Hz during the run.

Detectors

Module housing

An important feature of the CRESST detector modules is their ability to reject events which are not nuclear recoils via a simultaneous measurement of a heat/phonon signal in the large detector and a scintillation signal [14] in the light detector. The crystal and the light detector are enclosed together in a reflective and scintillating housing (Fig. 4) of 3M foil, as described in Ref. [7]. This serves two purposes. Firstly, it increases the amount of light collected by the light detector. Secondly, its scintillation permits rejection of a potentially dangerous background from surface radioactivity on the crystal [15]. An important such possibility is polonium, a daughter of the radon decay chain with an alpha decay. If the alpha is lost and the recoiling lead nucleus stays in the crystal, this can mimic a tungsten recoil, i.e. a WIMP interaction. However, it was found that the scintillation of the alpha particle on the reflecting foil gives a light signal which discriminates events of this type and allows for their separation.

Thermometers (TES's)

The basic sensor in CRESST is a superconducting thin film thermometer, with the temperature stabilized at a suitable operating point in the middle of the superconducting-to-normal transition of the film . A small variation of the film resistance in such an arrangment gives a very sensitive temperature measurement; this system is often called a Transition Edge Sensor (TES) .

In the work described here tungsten thin films are used. Tungsten in its crystalline α phase becomes superconducting at T c ≈ 15 mK. The films are equipped with heaters with the dual function of stabilizing the detectors' operating points and the injection of heat pulses for calibration. The thermometer layouts for the CaWO 4 target and the light detector are similar. However, the light detector has additional Al/W phonon collectors for better energy collection. This is as described in ref. [7], but here the light detectors use a small portion of the thermal-link gold structure as the heater. This light detector layout is shown in Fig. 5.

For temperature regulation heater pulses of a fixed amplitude, so-called control pulses are injected through the heater lines every three seconds. Online pulse height evaluation of the resulting signal is used for regulation of the heater current which keeps the thermometer's temperature at the operating point. This is somewhat different from older work where temperature stabilization was on the baseline, and is found to provide excellent long term stability of operation. In addition to these control pulses for temperature stabilization, calibration pulses delivering a known energy are injected every 30 seconds. These are used in determining the energy to be identified with a pulse, as will be explained in section 4.2. The calibration and stabilization procedures are the same for both target and light detector thermometers.

Target Detector

The CRESST-II detectors use large CaWO 4 single crystals as the target or "absorber" mass. CaWO 4 has been selected for its high light output and the presence of the heavy nucleus tungsten, which gives a large factor ∼ A 2 in the presumed coherent (spin-independent) scattering of the WIMP. The crystals are cylindrical, of 4 cm height and 4 cm diameter, and about 300 g in weight.

Ten such crystals were installed for the commissioning run; two are used here for the dark matter analysis.

Light Detector

The scintillation light is measured by a separate cryogenic detector. This is made from a sapphire wafer of 40 mm diameter and 0.4 mm thickness, with an epitaxially grown silicon layer on one side for photon absorption. The weight is 2.3 g. This silicon-on-sapphire (SOS) material exploits the excellent phonon transport properties of sapphire. The resulting detector with the sensor as in Fig. 5 is sensitive down to about 20 eV, or about seven optical photons. Further progress in using the light yield to identify which nucleus is recoiling, as discussed below, depends largely on the improvement of the resolution of the light detector.

Data Taking and Analysis

Only very basic quality cuts are performed on the data. Pulses are rejected if they occur in one of the rare periods when the temperature of the cryostat is not stable, as inferred from the pulse height of the control pulses. Pileup is rejected by simple cuts on the baseline of the pulse. Due to the low trigger rate this does not introduce a significant dead time into the analysis.

Signal Processing and Pulse Fitting

For recording thermometer pulses, the output voltage of the SQUID electronics is split into two branches. In one branch the pulse is shaped and AC-coupled to a trigger unit, while in the other branch the signal is passed through an 8-pole anti-aliasing low-pass filter and then DC-coupled to a 16-bit transient digitizer.

The time bin of the transient digitizer was chosen to be 40 µs, providing about 30 samples for the rising part of the pulse in the phonon channel. The record length of 4096 time bins includes a pre-trigger region of 1024 samples to record the baseline before the event, while the remaining 3072 post-trigger samples contain the pulse itself. The phonon and light channels of each module are read out together, whenever one or both trigger.

The pulses are fitted offline using a template fit, with a template constructed from pulses collected during a gamma calibration run. There is an individual template for each thermometer. Both pulses from a module are fitted simultaneously. Free parameters in the fit are the two baselines, the two pulse height amplitudes, and a common time shift relative to the trigger.

Energy Determination

Finally, the assignment of an energy to a pulse is made on the basis of the resulting amplitude or pulse height parameter. Since the pulses from the 57 Co calilbration with ∼ 122 keV and ∼ 136 keV photons establish a relation between pulse height and energy for each detector, a first simple determination of the pulse height-energy relation for a detector would be provided by a linear extrapolation from 122 keV to lower energies, and in early work this gave good results. Direct calibration with lower energy gammas is not possible since the calibration source is outside the coldbox and the photon must penetrate about 12 mm of copper.

A finer and more detailed determination of the pulse height-energy relation for a given detector is possible, however, by using the heater calibration pulses mentioned in section 3.2. For each thermometer, heater pulses delivering a known energy via ohmic heating are injected every 30 seconds. These are applied over the whole energy range of interest, so the thermometer response can be mapped down to low energies and the assumption of a simple linear extrapolation can be avoided. Both target and light detectors are calibrated by this procedure. The accuracy and resolution of this method is illustrated in Fig. 6, which shows the energy spectrum of the target detector "Verena" in the commissioning run. The energy resolution may be gauged from the prominent feature starting around 46 keV and the line near 8 keV, which can be identified as a copper fluorescence line. The large feature is due to the gamma transition of a 210 Pb impurity in the crystal, with a broadened right shoulder from an accompanying 17 keV β spectrum. The steep left edge appears precisely at 46.5 keV as expected. The peak for the copper fluorescence line is at 8.05 keV, also as expected, and the energy resolution is ∼ 300 eV. These two features confirm the excellent accuracy of the calibration method with heater pulses. Further details of the energy spectrum can also be understood and will be discussed in a forthcoming publication.

Energy-Light Plots and Quenching Factor

Events from a given module are plotted as points in the energy-light yield plane, as in the following Figs. 7 and8. An energy is assigned to an event as described in the previous section. The light yield is given by the ratio of the energy or pulse height in the light detector to that in the target detector, normalized such that it is 1 for a 122 keV calibration photon absorbed in the large detector. While the light yield for a 122 keV photon is therefore one by definition, the fact that the bands are roughly flat on the plots, and the average light yield is always near one, is not a matter of definition and reflects the accuracy of the thermometer calibrations and the approximate linearity1 of the light production process [17].

It will be seen that when a source inducing nuclear recoils is present, as with neutrons in Fig. 7, there are two distinct bands, one of lower light yield from nuclear recoils and one of higher light yield from electron-photon events. The factor of reduction for the light yield of a nuclear recoil relative to that of the electron-photon event of the same energy is called (in our usage) the quenching factor. With the electron-photon band centered on 1, the center of the nuclear recoil band on Fig. 7 is simply at the inverse of the quenching factor. Fig. 7. Low energy event distribution in the (Verena/SOS21) module in the presence of a neutron source during the commissioning run. Below the red curve 90 % of the oxygen recoils are expected, as calculated from the known quenching factor for oxygen and the energy resolution of the light detector. The good agreement with the prediction supports the use of this method to identify tungsten recoils in the dark matter analysis. The events above 47 keV between the electron and nuclear recoil bands can be understood as inelastic neutron scatterings where the tungsten nucleus is excited to a gamma-emitting level.

Since the quenching factor varies widely, up to ∼ 40 for heavy elements, the level of light output can be used not only to distinguish nuclear recoils from electron-photon events, but also to distinguish tungsten recoils from those of the lighter nuclei. It is hoped that it may be ultimately possible to determine which particular nucleus, in a compound material like CaWO 4 , is recoiling. This would be very useful in verifying a positive dark matter signal, since a simple behavior for the signal with respect to changing the target nucleus is expected.

Since an understanding of the quenching factor is thus of great interest, a number of investigations have been carried out within the CRESST collaboration [15], [18], [19], [20], [21]. The results are summarized and discussed in [17].

Neutron Test

To directly test the response of the system to nuclear recoils a neutron test was carried out during the commissioning run. The test can be used to check the principle of identifying the recoiling nucleus via the light yield and to also determine the dispersion in the light output, which is important in applying the method.

Above ∼ 10 keV recoil energy, calculations show that the spectrum for neutron scattering on CaWO 4 is totally dominated by oxygen recoils [START_REF] Wulandari | The dominance of oxygen recoils above 10 keV is principally due to kinematics; the test neutrons have ∼MeV energy and the recoil energy of a nucleus is ∼ 2E neutron[END_REF]. Therefore one anticipates that the events of the nuclear recoil band in this test should lie roughly on or below the line corresponding to the oxygen quenching factor ∼ 9 known from the previous quenching factor studies.

A neutron source was placed inside the polyethylene shield and Fig. 7 shows the results for one of the two modules (Verena/SOS21). Using the quenching factor of 9 and the measured energy-dependent resolution of the light detector, 90 % of the neutron-induced recoils should be below the red line.

The neutron test data appears to agree well with this prediction, and so supports our arguments in the dark matter analysis that an acceptance region for nuclear or tungsten recoils based on the same principle can be defined.

Fig. 8. Low-energy event distribution measured with two 300 g CaWO 4 detector modules during the commissioning run. The vertical axis represents the light yield (see text), and the horizontal axis the total energy, as measured by the phonon channel. Below the dashed black curve 90 % of all nuclear recoils, and below the solid red curve 90 % of the tungsten recoils are expected. The heavy black dots show the events in the "tungsten recoils" acceptance region. The intense regions of the electron recoil bands correspond to the lines of Fig. 6.

Dark Matter Limits

For the dark matter analysis we use data taken between March 27th and July 23rd 2007 with the two detector modules (Zora/SOS23) and (Verena/SOS21). The data are shown in Fig. 8. The cumulative exposure was 47.9 kg -days. For the tungsten this corresponds to an exposure of 30.6 kg -days.

We perform an analysis on the assumption of coherent or spin-independent scattering for the WIMP. This process should strongly favor tungsten recoils due to the ∼ A 2 factor in the WIMP-nucleus cross section. We define an acceptance region on the plots based on a) the quenching factor for the light yield and b) the maximum energy expected for tungsten recoil. A similar region for "all nuclear recoils" can be defined using the quenching factors for Ca and O. The quenching factor boundaries are shown on the plots of Fig. 8 by the curves. Below the upper curve 90 % of all nuclear recoils are expected, and below the lower curve 90 % of the tungsten recoils are expected. The energy boundaries are shown by the vertical lines. The upper limit at 40 keV is set by form-factor [START_REF] Helm | [END_REF] effects, which effectively limit the energy transfer to the tungsten nucleus. The lower limit is set at 10 keV, where "leakage" from the electron/photon events becomes evident and so recoil discrimination becomes inefficient.

Three candidate events (heavy black dots) are observed in Fig. 8 for the "tungsten recoils" acceptance region. The individual events are at 16.89 keV (Verena/SOS21) and at 18.03 keV and 33.09 keV (Zora/SOS23). The corresponding rate is 0.063 per kg-day. From this rate we can derive a limit on the coherent spin-independent WIMP-nucleon scattering cross section. Using standard assumptions on the dark matter halo [24,25] (WIMP mass density of 0.3 GeVcm -3 ), an upper limit for the coherent or spin-independent WIMPnucleon scattering cross-section may be obtained using the maximum energy gap method [26]. This limit is plotted as the red curve in Fig. 9. The minimum of the curve, for a WIMP mass of ∼50 GeV, is at 4.8 × 10 -7 pb. Fig. 9. Coherent or spin-independent scattering cross section exclusion limit derived from the data of Fig. 8 using the maximum energy gap method. For comparison the limits from other experiments [27], [28], [29], [30], [31] and the range predicted by some supersymmetry models [32] are also shown.

Backgrounds

To judge the improvement with respect to our previous setup of Ref. [7] we take the "all nuclear recoils" region and reset the low energy boundary of the acceptance region to 12 keV as it was there. In this acceptance region the earlier work had 16 events, for a rate of 0.87 per kg-day. As noted in Ref. [7], this rate was compatible with the estimates of [33] for neutron background. The present commissioning run has 4 events in the same acceptance region, for a rate of 0.083 per kg-day, so we may say there is a factor ∼ 10 improvement. We attribute this to the presence of the neutron shield.

On the other hand, according to the arguments in section 5 the "tungsten recoils" should be insensitive to neutron background. Indeed, the three events here from 48 kg-days are quite compatible with the zero events from 10 kgdays in Ref. [7]; the present rate of 0.063 per kg-day predicts only 0.6 events for 10 kg-days. The fact that the neutron shield appears to have a great effect on "all nuclear recoils" but little effect on "tungsten recoils" supports our argument that the light yield successfully distinguishes tungsten and other nuclear recoils.

The question arises, however, as to the nature of the few observed tungsten or nuclear recoil candidates. One possibility is evidently remaining neutrons. However with many absorption lengths for the shield very few should penetrate and indeed simulations [34] would give a rate of only ∼ 10 -5 per kg-day, much less than the few events. A possible point here is that during the run a weak spot in the neutron shielding above the muon veto was identified and patched only after data taking was completed.

Another possible source for neutrons are some non-operational modules which were present during the run. These can act as non-vetoable neutron sources for the operational modules. However an estimate for the background events from an inactive detector module is well below 1.4 × 10 -5 per kg-day.

Apart from neutrons, another possibility arises from incomplete coverage of the inner surfaces of the detector module with scintillator, which could lead to unvetoed nuclear recoils from surface α-decays as discussed in section 3.1. In further work it is intended to paint some possibly uncovered areas with scintillating material.

Concerning muons, estimates [34] of muon-induced neutrons in the setup result in only 2.8 × 10 -3 per kg-day. No muon signals were found in coincidence with nuclear recoil events.

There thus appears to be no conclusive explanation for the few candidate events from conventional radioactive or particle sources. We hope to clarify some of these points in further work.

Conclusions

CRESST-II successfully completed its commissioning run in 2007. New elements of the apparatus were successfully installed and operated. A neutron test demonstrated the ability to detect nuclear recoils with a light yield consistent with that found from quenching factor studies, supporting the principle of identifying the recoil nucleus via the light yield.

Data were taken with two detector modules for a total of 48 kg-days. Three candidate events of uncertain origin are present in the acceptance region for tungsten recoils, yielding a rate of 0.063 per kg-day. A factor ∼ 10 improved performance is found with respect to previous work for the "all nuclear recoils" acceptance region.

A limit on coherent WIMP-nucleon scattering, is obtained, which at its lowest value, for M WIMP ≈ 50GeV, is 4.8 × 10 -7 pb.

Fig. 1 .

 1 Fig.1.3 He-4 He dilution refrigerator, low background detector environment, and the shielding as upgraded for CRESST-II. The detector carousel (CA) is connected to the mixing chamber of the cryostat (CR) by a long copper cold finger (CF) in order to reduce background originating from the dilution refrigerator. The gas-tight radon box (RB) encloses the low background copper (CU) and low background lead shielding (PB). It is covered by a plastic scintillator muon-veto (MV) and a 45 cm thick polyethylene neutron moderator (PE). Additional granular PE is placed between the baffles in the upper part of the cryostat to close the line of sight for neutrons coming from above.

Fig. 2 .

 2 Fig. 2. Detector carousel of ultrapure copper, electropolished to reduce surface contamination. The structure can accommodate 33 detector modules (i.e. 10 kg of target mass) which can be mounted or dismounted individually. The carousel is mounted at the lower end of the cold finger. The whole structure is cooled to ∼ 10 mK.

Fig. 3 .

 3 Fig. 3. CRESST-II readout circuit. The tungsten film thermometer is in parallel with two reference resistors and the input coil of a DC SQUID. The circuit is biased by a constant current.

Fig. 4 .

 4 Fig.4. Schematic arrangement of a CRESST detector module. The module consists of two independent detectors: one with the CaWO 4 target crystal, providing a total energy measurement, and one with a silicon-on-sapphire (SOS) wafer for measuring the scintillation light from the target crystal. Both detectors are enclosed in a reflective, scintillating housing.

Fig. 5 .

 5 Fig.5. Layout and connection scheme of a light detector used in the commissioning run. Aluminum/tungsten phonon collectors surround the tungsten thermometer or TES. A small portion of the gold thermal link is used as the heater. Connections shown are the aluminum bond wires for the SQUID and heater circuits, and a gold bond wire for thermal coupling, leading to an electrically insulated pad on the detector holder.

Fig. 6 .

 6 Fig.6. Energy spectrum of the 300g detector "Verena" in the commissioning run (48 days). The prominent feature at 46.5 keV is due to a gamma from 210 P b decay in the crystal with its associated 17 keV β spectrum. The Cu fluorescence line at 8.05 keV is also clearly visible, and has resolution ∼ 300 eV. The close agreement of these energies with known values confirms the excellent accuracy of the heater pulse calibration method.

Upon closer examination it appears that the centroid of the electron-photon band falls slightly at lower energies. Such a decrease for γ rays in CaWO 4 has been seen in Ref.[16].

Acknowledgements

This article is dedicated to the memory of our friend and colleague Wolfgang Westphal, who significantly contributed to the success of CRESST.

This work was partially supported by funds of the DFG (SFB 375, Transregio 27:"Neutrinos and Beyond"), the Munich Cluster of Excellence ("Origin and Structure of the Universe"), the EU networks for Cryogenic Detectors (ERB-FMRXCT980167) and for Applied Cryogenic Detectors (HPRN-CT2002-00322), and the Maier-Leibnitz-Laboratorium (Garching). We also thank Hans Ulrich Friebel for valuable technical support. Support was provided by the Science and Technology Facilities Council.


Hugo Leclerc
email: hugo.leclerc@lmt.ens-cachan.fr

Jean-Noël Périé
email: jean-noel.perie@iut-tlse3.fr

Stéphane Roux
email: stephane.roux@lmt.ens-cachan.fr

François Hild
email: francois.hild@lmt.ens-cachan.fr

J.-N Périé

Voxel-Scale Digital Volume Correlation

Keywords: Digital Volume Correlation, Equilibrium Gap Method, Regularization

Among various correlation techniques to find the displacement field of a volume imaged by X-Ray tomography at several deformation states, a new approach is proposed where the displacement is measured down to the voxel scale and determined from a mechanically regularized system using the equilibrium gap method, and an additional boundary regularization.

It is shown that even if the underlying material behavior is not very well known, this approach leads to extremely small correlation residuals. An excellent stability of the estimated displacement field for noisy (reconstructed) volumes is also observed.

1 Introduction X-ray computed microtomography (X-CMT) is a very powerful means for visualizing in a non-destructive way different (opaque) materials [START_REF] Baruchel | X-Ray Tomography in Material Sciences[END_REF][START_REF] Maire | On the application of X-ray microtomography in the field of materials science[END_REF][START_REF] Bernard | 1st Conference on 3D-Imaging of Materials and Systems 2008[END_REF][START_REF] Stock | MicroComputed Tomography: Methodology and Applications[END_REF][START_REF] Stock | Recent advances in X-Ray microtomography applied to materials[END_REF][START_REF] Banhart | Advanced Tomographic Methods in Materials Research and Engineering[END_REF]. For instance, yielding [START_REF] Bastawros | Compressive deformation and yielding mechanisms in cellular Al alloys determined using X-ray tomography and surface strain mapping[END_REF] and damage mechanisms [START_REF] Babout | Characterisation by X-Ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF][START_REF] Bontaz-Carion | X-Ray Microtomography Analysis of Dynamic Damage in Tantalum[END_REF][START_REF] Schilling | X-ray computed microtomography of internal damage in fiber reinforced polymer matrix composites[END_REF], cracks [START_REF] Ludwig | Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography[END_REF][START_REF] Sinclair | The effect of fibre fractures in the bridging zone of fatigue cracked Ti6Al4V/SiC fibre composites[END_REF][START_REF] Ferrié | Fatigue crack propagation: In situ visualization using X-ray microtomography and 3D simulation using the extended finite element method[END_REF] are observed in situ. Another feature is to have access to microstructural features [START_REF] Wolfsdorf | The morphology of high volume fraction solid-liquid mixtures: An application of microstructural tomography[END_REF][START_REF] Ludwig | In-situ Three-Dimensional Microstructural Investigation of Solidification of an Al-Cu Alloy by Ultrafast X-ray Microtomography[END_REF][START_REF] Maire | Characterization of the morphology of cellular ceramics by 3D image processing of X-ray tomography[END_REF][START_REF] Viot | Impact test deformations of polypropylene foam samples followed by microtomography[END_REF], some of which may directly influence the observed phenomenon (e.g., fatigue cracking [START_REF] Ludwig | Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography[END_REF]). From the reconstructed volumes, it is possible to mesh the microstructure [START_REF] Youssef | Finite element modelling of the actual structure of cellular materials determined by X-ray tomography[END_REF][START_REF] Maire | X-ray tomography applied to the characterization of cellular materials. Related finite element modeling problems[END_REF] in order to perform a finite-element simulation, tailored to the specimen microstructure, so as to describe its mechanical behavior.

When considering in situ experiments for which scans are obtained at different loading steps, it is possible to evaluate displacement fields with a coarse spatial resolution by resorting to marker tracking [START_REF] Nielsen | Measurements of plastic displacement gradient components in three dimensions using marker particles and synchrotron X-ray absorption microtomography[END_REF][START_REF] Toda | Assessment of the fatigue crack closure phenomenon in damage-tolerant aluminium alloy by in-situ high-resolution synchrotron X-ray microtomography[END_REF][START_REF] Withers | Crack opening displacements during fatigue crack growth in Ti-SiC fibre metal matrix composites by X-ray tomography[END_REF]. Another technique consists in local correlations [START_REF] Bay | Digital volume correlation: three-dimensional strain mapping using X-ray tomography[END_REF][START_REF] Bornert | Mesure tridimensionnelle de champs cinématiques par imagerie volumique pour l'analyse des matériaux et des structures[END_REF] in which small interrogation volumes in two scans are registered [START_REF] Mckinley | Trabecular bone strain changes associated with subchondral stiffening of the proximal tibia[END_REF].

Global correlation techniques provide an alternative route [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF][START_REF] Réthoré | Extended three-dimensional digital image correlation (X3D-DIC)[END_REF]. The latter was used to analyze the correlation between local density and strain levels [START_REF] Hild | Three dimensional analysis of a compression test on stone wool[END_REF], and stress intensity profiles in a cracked sample [START_REF] Limodin | Crack closure and stress intensity factor measurements in nodular graphite cast iron using 3D correlation of laboratory X ray microtomography images[END_REF][START_REF] Limodin | Influence of closure on the 3D propagation of fatigue cracks in a nodular cast iron investigated by X-ray tomography and 3D Volume Correlation[END_REF][START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]. When measurement uncertainties are compared, the correlation approaches usually outperform the marker tracking technique [START_REF] Limodin | Crack closure and stress intensity factor measurements in nodular graphite cast iron using 3D correlation of laboratory X ray microtomography images[END_REF].

The crucial question common to all correlation techniques is to select the appropriate discretization of the sought displacement field. Large interrogation windows (for local analyses) or coarse meshes (for global approaches) allow for a low uncertainty in the estimated local mean displacement at the expense of a poor spatial resolution [START_REF] Bergonnier | Digital image correlation used for mechanical tests on crimped glass wool samples[END_REF][START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF][START_REF] Bornert | Assessment of Digital Image Correlation measurement errors: Methodology and results[END_REF]. Conversely, fine meshes (or small interrogation windows) induce a better spatial resolution, but because of the fewer pixels or voxels considered, displacement uncertainties are larger. Hence, a compromise has to be found. The proper balance between spatial resolution and displacement uncertainty is a question that has to be addressed for each analysis since both elements are influenced by the image texture and complexity of the displacement field [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF]. Three dimensional X-CMT images display specific features in this connection. Reconstruction artifacts are inevitably present in the image [START_REF] Davis | Artefacts in X-ray microtomography of materials[END_REF][START_REF] Ketcham | New algorithms for ring artefact removal[END_REF] and hence small element sizes can hardly be used.

The present paper aims at breaking the above limitation, and a voxel-wise spatial resolution is aimed for, which may appear as an unreachable goal. The key to allow for a kinematic flexibility at the voxel scale (i.e., two neighboring voxels may have very different displacements) is to introduce a mechanical regularization based on the hypothesis that the local behavior is, say, elastic for the sake of simplicity. It is proposed to couple digital volume correlation [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF] and a regularization based on the equilibrium gap method [START_REF] Roux | Digital Image Mechanical Identification (DIMI)[END_REF]. The interest of the proposed approach is to be able to provide a much more resolved displacement field with an excellent quality of the registration.

In Section 2, the principle of regularized DVC is introduced. Section 3 deals with artificial cases generated to test the methodology and its performance. Section 4 is devoted to the analysis of an in situ tensile test on a spheroidal graphite cast iron specimen whose tomography is performed on a synchrotron beam-line.

Voxel-Based Digital Volume Correlation

In the following, a voxel-based digital volume correlation (V-DVC) technique is presented. This is made possible by regularizing the (ill-posed) measurement procedure by a mechanical analysis that is performed in a fully integrated way. Such type of procedures has been implemented to deal with 2D pictures. The first one was based upon updating finite element models [START_REF] Leclerc | Integrated Digital Image Correlation for the Identification of Mechanical Properties[END_REF]. The second one [START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF] resorted to the equilibrium gap method [START_REF] Claire | A finite element formulation to identify damage fields: The equilibrium gap method[END_REF]. The latter will be used in the following. However, none of them was used to measure displacements at the voxel scale as discussed herein.

Global Approach to DVC

Local approaches to DVC [START_REF] Bay | Digital volume correlation: three-dimensional strain mapping using X-ray tomography[END_REF][START_REF] Bornert | Mesure tridimensionnelle de champs cinématiques par imagerie volumique pour l'analyse des matériaux et des structures[END_REF] consider sub-volumes that are registered [START_REF] Mckinley | Trabecular bone strain changes associated with subchondral stiffening of the proximal tibia[END_REF] in two scans. In the present case, a global approach [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF][START_REF] Réthoré | Extended three-dimensional digital image correlation (X3D-DIC)[END_REF] is chosen. It allows for a direct coupling with mechanical simulations that will be used to regularize the correlation procedure.

Correlation Procedure

The registration of two gray level volumes f and g (f is the reference, and g the deformed one) is based upon the conservation of the image texture, which for X-CMT is the X-Ray absorption (whose gray-level encoding gives the reconstructed 3D images)

f (x) = g(x + u(x)) (1)
where u denotes the unknown displacement field. It consists in finding the best displacement field by minimizing the square of the correlation residual φc

φc(x) = |g(x + u(x)) -f (x)| (2)
The minimization of φ 2 c is a non-linear and ill-posed problem. In particular, if no additional information is available, it is impossible to determine the displacement for each voxel independently since there are three unknowns for a given (scalar) gray level.

For these reasons, a weak formulation is preferred. After integration over the whole Region Of Interest (ROI) Ω, the global correlation residual Φ 2 c reads

Φ 2 c = 1 2 Ω (g(x + u(x)) -f (x)) 2 dx (3)
in which the displacement field is expressed in a (chosen) basis

u(x) = n unψn(x) (4)
where ψn are (chosen) vector functions, and un the associated degrees of freedom. The measurement problem then consists in minimizing Φ 2 c with respect to the unknowns un. A Newton iterative procedure is followed to circumvent the non-linear aspect of the minimization problem. Let u i denote the displacement at iteration i, and {u} i the vector containing all the unknown degrees of freedom. By assuming small increments du = u i+1u i of the solution, a Taylor expansion is used to linearize g(x + u(x)) ≈ g(x) + u(x) • ∇f (x) and then, ∂Φ 2 c /∂{u} i is recast in a matrix-vector product as

∂Φ 2 c ∂{u} i = [M] i {du} -{b} i = {0} (5)
with

M i mn = Ω ∇f (x) • ψm(x)∇f (x) • ψn(x)dx (6)
and

b i m = Ω ∇f (x) • ψm(x)(f (x) -g(x + u i (x)))dx (7)
At this level of generality, many choices can be made to measure displacement fields.

When dealing with volumes, 8-noded elements with classical [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF] or enriched [START_REF] Réthoré | Extended three-dimensional digital image correlation (X3D-DIC)[END_REF] kinematic bases are used. The former is referred to as C8-DVC, and the latter XC8-DVC.

In the sequel, C8-DVC will be used since no displacement discontinuity is sought.

Consequently, the region of interest is discretized with finite elements for which the displacement field is interpolated with trilinear shape functions. The equivalent of a sub-volume in local approaches to DVC is the element for global DVC. The size of each element is defined by the length of any of its edge.

Resolution Analysis

To illustrate the ill-posedness of the measurement problem, let us perform the following resolution analysis [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF][START_REF] Roux | Stress intensity factor measurements from digital image correlation: postprocessing and integrated approaches[END_REF]. A reference volume is considered and the sensitivity of the measured degrees of freedom un to noises associated with the image acquisition and the reconstruction process is investigated. Both reference and deformed volumes are affected by the same noise. However, the noise-free reference being unknown, this is equivalent to considering a noise in the difference (g -f ) of variance 2σ 2 . It is therefore assumed that the deformed volume g is polluted by a random white noise η, of zero mean, variance 2σ 2 , and normalized covariance ρ(v, w), where v and w denote the indices of the voxels. The [M] matrix is unaffected by this noise, but only the second member vector {b} is modified by an increment

δbm = v (∇f • ψm)(v)η(v) (8)
whose average is thus equal to 0, and its covariance reads

C(δbm, δbn) = 2σ 2 v w (∇f.ψm)(v)ρ(v, w)(∇f • ψn)(w) (9)
The impact of this quantity on the uncertainty of un is sought. By linearity [see Equation (5)], the mean value of δu i vanishes too, and its covariance becomes

C(δu i , δu j ) = 2σ 2 M -1 im v w (∇f • ψm)(v)ρ(v, w)(∇f • ψn)(w) M -1 nj (10)
Equation ( 10) is a general result for any type of noise, kinematic decomposition, and texture. It shows that the key quantity to evaluate the covariance matrix is the local sensitivity (∇f.ψm), which is localized in all elements whose connectivity involves the considered degree of freedom um.

If the noise can be considered as white, i.e., the normalized covariance in its discrete definition is equal to δvw, where δ is Kronecker's delta, the previous result becomes

C(δu i , δu j ) = 2σ 2 M -1 ij (11)
Further, if the correlation length (i.e., the characteristic distance over which the normalized pair correlation function decays to 0) of the texture is significantly larger than the voxel size, which is desirable if the noise sensitivity is to be limited, but not too large, so that the matrix [M] is invertible, the standard displacement uncertainty

σu reads σu = √ 6Apσ (∇f ) 2 3/2 (12)
where • • • denotes the volume average, 3 the number of voxels in the considered finite element, p the physical size of one voxel, A a dimensionless constant dependent on the interpolation function [START_REF] Roux | Stress intensity factor measurements from digital image correlation: postprocessing and integrated approaches[END_REF][START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF].

To illustrate Equation ( 12), let us consider the volume f to be studied in Section 4.

A Gaussian white noise ( √ 2σ = 1 gray level) is added to the latter and a corrupted volume g is thus obtained without displacement. A C8-DVC correlation is performed to register the volumes for different element sizes. The measured displacement is expected to be identically zero. Figure 1 shows the standard deviation of the measured displacements. In all the results reported herein, it was checked that the bias (or average displacement error) was at least one order of magnitude lower than the standard deviation. Consequently, only the latter is shown. This result is in accordance with the resolution analysis performed above. As the element size increases, the displacement uncertainty decreases. When extrapolated to small element sizes, the standard displacement uncertainty reaches levels that make any measurement meaningless and eventually impossible (i.e., matrix [M] is no longer invertible).

A power law interpolation of the standard displacement uncertainty σu = B/ α describes very well the results for which a significant amount of degrees of freedom does not belong to any edge of the ROI. For the latter, the sensitivity to noise is higher since part of the information is missing (e.g., one half for faces, three fourth for edges, and seven eighth for corners). For element sizes less than 32 voxels, a value α = 1.42 is found, which is in good agreement with the predicted value (1.5) given in Equation ( 12).

Regularized DVC

The previous results show that the smaller the element size , the larger the standard displacement uncertainty. Consequently, if a voxel-scale determination of the displacement field is sought, additional information is needed to regularize the measurement problem. This will be achieved by using the equilibrium gap method.

Equilibrium Gap

To enforce mechanical admissibility in an FE sense, the equilibrium gap is first introduced. If linear elasticity applies, the equilibrium equations read

[K]{u} = {f } (13)
where [K] is the stiffness matrix, and {f } the vector of nodal forces. When the displacement {u} is prescribed and if the (unknown) stiffness matrix is not the true one, load residuals {fr} will arise

{fr} = [K]{u} -{f } (14)
In the absence of body forces, interior nodes are free from any external load. Consequently, the minimization of the equilibrium gap consists in minimizing the following quantity

Φ 2 m = 1 2 {u} t [K] t [K]{u} ( 15
)
where t is the transposition operator, and 2Φ 2 m corresponds to the sum of the squared norm of all equilibrium gaps at each interior node. It is worth noting that the unknowns are related to the components of the stiffness matrix. As such, the determination of stiffness fields requires a discretization that is coarser than that of the kinematic mesh [START_REF] Roux | From image analysis to damage constitutive law identification[END_REF].

However, if the local stiffness can be related to the gray level, then the problem can be regularized, as was shown when identifying a damage law [START_REF] Roux | Digital Image Mechanical Identification (DIMI)[END_REF]. The same type of procedure will be followed in the sequel, pushing the procedure down to its ultimate limit, namely, the voxel size.

Regularized Correlation Procedure

To solve the coupled minimization of correlation residuals (Φ 2 c ) and equilibrium gap (Φ 2 m ), a weighted sum of both residuals is defined as

Φ 2 t = (1 -β) Φ 2 m + β Φ 2 c ( 16
)
where β is a dimensionless parameter in the interval [0, 1], Φ 2 m the normalized equilibrium gap, and Φ 2 c the normalized correlation residual.

To normalize both residuals, let us consider a unitary displacement field

u 1 = u u (17)
With the previous definition, both residuals read

ϕm 2 ≡ Φ 2 m (u 1 ) = 1 2 {u 1 } t [K] t [K]{u 1 } ϕc 2 ≡ Φ 2 c (u 1 ) = 1 2 {u 1 } t [M]{u 1 } (18)
The latter values are thus used to make their corresponding residuals dimensionless and comparable. It is worth noting that any displacement field u 1 may be used in the present normalization technique. In the following, a simple "plane wave" displacement is chosen in the x direction with a 16-pixel wavelength (short enough to fit in small volumes, large enough to avoid discretization artifacts), meaning that with β = 1 2 , the "cut-off wavelength" is equal to 16 pixels. The normalized residuals then read

Φm = Φm ϕm and Φc = Φc ϕc (19)
When the material parameters are known, the minimization of Φ 2 t with respect to the unknown degrees of freedom can be performed [START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF]. It corresponds to an integrated way of measuring displacement fields since the displacement field will satisfy equilibrium equations in addition to the texture conservation.

Regularization of the boundary conditions

The above introduced equilibrium gap functional exploits the residual body forces at interior nodes, which are to be as small as possible. However, along the boundary of the ROI neither static nor kinematic information is available. Any displacement field prescribed on the boundary gives rise to one displacement field for which Φm = 0.

Hence, without boundary conditions, the mechanical regularization vanishes on the former. Similarly, close to the edges, the DVC algorithm gives rise to a higher level of uncertainty than in the bulk, because of the fewer number of neighboring elements.

Hence, although regularization is achieved in the bulk, boundaries remain a weak point in the analysis, and as smaller and smaller element sizes down to the elementary voxel are sought, ill-posedness of the problem remains. This observation calls for a third regularization for boundary nodes.

Since no mechanical information is available, it is proposed to introduce a penalization of short wavelength displacement fluctuations along the faces of the region of interest. The third objective function to be considered should however vanish for any rigid body motion. It is therefore proposed to introduce

Φ 2 b = 1 2 {u} t [L] t [L]{u} ( 20
)
where L is a discretized Laplacian operator acting on the faces of the ROI. For any translation and rotation, the above functional vanishes. Let us note that in two dimensions, where the same procedure can be applied, the boundary regularization takes the form of the quadratic norm of the second derivative of the displacement field along the edges. In the special case where the boundary of the ROI meets a free surface, a complete static information (zero tractions: σ σ σ.n = 0) is available and should be used instead of the above form.

As for the above cases, the boundary objective function is normalized based on a similar reference displacement field

Φ b = Φ b ϕ b ( 21
)
where

ϕ b 2 ≡ Φ 2 b (u 1 ) = (1/2){u 1 } t [L] t [L]{u 1 }.
Finally, the total objective function to be minimized reads

Φ 2 t = (1 -β -γ) Φ 2 m + β Φ 2 c + γ Φ 2 b ( 22
)
where β > 0, γ > 0 and β + γ < 1.

Justification of the regularization

To understand that the above combination of the correlation and mechanical residuals will allow for a regularization, it is necessary to resort to a spectral analysis of both functionals. Let us consider a displacement field in the form of a plane wave du(x) = du 0 exp(ik.x) with a vanishing amplitude |du 0 | as a perturbation of a reference displacement field, ueq ideally assumed to minimize both functionals. The correlation residual amounts to

Φ 2 c (ueq + du) = Φ 2 c (ueq) + 1 2 {du} t [M]{du} (23)
On average, the residual perturbation dΦ 2 c scales as

|du 0 | 2 [M] |Ω| The important
feature to note is that this term is independent of the wavevector k.

The mechanical residual perturbation is also a quadratic form in du 0 . However, it scales as the fourth power of the wavenumber

dΦ 2 m ∝ |k| 4 |du 0 | 2 |Ω| (24)
A similar spectral analysis performed on the boundary term shows that the scaling of Φ b for a plane wave displacement field is similar to that of the mechanical regularization

dΦ 2 b ∝ |k| 4 |du 0 | 2 |∂Ω| ( 25
)
where ∂Ω is the outer boundary of the considered ROI. Hence, Φm and Φ b are natural complements.

Balancing the expressions of Φc and Φm or Φ b shows that there are two characteristic length scales, ξ and ξ b , such that

ξ ∝ 4 1 -β -γ β and ξ b ∝ 4 γ β ( 26
)
where ξ b gives the depth of the boundary layer that is affected by the boundary regularization. Concerning the bulk, below the ξ length scale, Φm dominates over Φc, and hence the small scale displacement is mechanically admissible, whereas at large length scales the correlation residual is dominant. These two length scales play the same role in terms of regularization as the mesh size in the C8-DVC formulation. However, rather than choosing finite element shape functions, which are convenient but not mechanically significant, the introduced formulation implements, in addition to a conventional boundary regularization, a physically motivated short scale regularization where the information is available (i.e., in the bulk). Adjusting the value of the parameters β and γ is a way of tuning the regularization scales.

Artificial Cases

Two artificial cases are studied, namely, the first one whose texture has a low frequency content. The second one is directly related to the practical case that will be studied in Section 4.

First case

The deformation of an artificial isotropic material with slowly varying Young's modulus E is first analyzed. As shown in Figure 2, Dirichlet boundary conditions are prescribed to solve a regular 3D finite element problem, which allows then to generate a deformed image. In the present case uniform displacements corresponding to an average uniaxial strain of 1% are prescribed on two opposite faces, and the other external boundaries are traction-free. This synthetic case is considered in order to validate the method. The dimensionless correlation residual, η, is defined as

η ≡ 1 max f -min f 1 |Ω| x (f (x + u(x)) -g(x)) 2 (27)
Figure 3 shows the change of the dimensionless correlation residual for two different DVC techniques as a function of a regularization length scale. The first curve is a reference solution where a global C8-DVC is chosen. In that case, the regularization length scale ξ is taken as the element size . The increase of the residual is due to a systematic error caused by the fact that coarser and coarser elements cannot capture accurately the actual displacement field. It is to be stressed that for < 4 voxels the C8-DVC algorithm becomes unstable (see Figure 1). The second curve corresponds to the present mechanical regularization using homogeneous elastic properties. In that case, ξ b is set to 16 voxels. The correlation residual typically is cut down by a factor of 4. Moreover, ξ can be reduced to values as small as 2 voxels with no convergence problems.

The robustness of the proposed method is better than that achieved with a standard C8-DVC procedure. It is worth noting that the correlation residual is lower even if the material properties are not the actual ones.

Second case

The texture that is used in this example (Figure 4) was obtained from a X-Ray tomography (for a more detailed description, see Section 4). The goal of this example is to study a case that is more realistic in terms of texture.

The effect of noise is shown in Figure 5 with the proposed (V-DVC) approach and with C8-DVC. The same procedure as in Section 2 is followed. The standard displacement uncertainty is lowered by one order of magnitude when comparing C8-DVC and V-DVC. This result shows that the proposed regularization is very effective, even for low cut-off lengths. Since the V-DVC approach also leads to smaller correlation residuals, it means that V-DVC is not a tradeoff between stability and accuracy. In the present case, a small ROI (i.e., 23 3 voxels) was considered since voxel-scale calculations are performed. Consequently, edge effects dominate and the power -3/2 expected from Equation ( 12) is not found. Rather, a power -1 is observed suggesting that face boundaries are the main source of uncertainty.

The deformation of the considered texture is now analyzed. As previously, the local Young's modulus E is assumed to be a deterministic function of the gray level (an affine variation is chosen with a maximum contrast of 3 over the gray levels spanned by the reference volume, namely, the maximum value is equal to three times the minimum value of E). An elastic computation is first performed when a 1%-extension along the

x-direction is prescribed as in the previous case. The lateral faces are stress-free, and the displacement field is used to create an artificial deformed volume. The change of the dimensionless correlation residual as a function of the regularization length scale is shown in Figure 6 for the two different DVC techniques. All the trends observed in the previous case are also found for the present texture. In particular, a significant reduction of correlation residuals is achieved with the regularization, especially for small regularization lengths. A factor 4 gain is again observed. As in the previous case, a saturation of the correlation residuals occurs for small regularization lengths, ξ, signaling that the results become independent of that parameter. The main difference is related to the absolute levels themselves that are higher in the present case. This is due to the more rapidly varying displacement field induced by the more complex texture.

Further, it is shown that as the cut-off wavelength ξ b decreases, the correlation residuals decrease as well. This result allows to regularize the measurement problem with very small wavelengths for the boundary conditions. The fact that a plateau region when ξ < 4 voxels is observed for any tested values of ξ b allows for a choice of ξ independent of the value chosen for ξ b . Due to exogenous nature of the regularization used for the boundary conditions, ξ b has a slightly different role, namely, in the generic case (when the displacement close to the boundaries has to be considered) it should be as low as possible, but high enough to maintain the global convergence properties.

From the study of these two artificial cases, it appears that the proposed regularization fulfills the expected objective. It allows for a rapid convergence towards a displacement field that has a much lower residual level than what could be achieved traditionally. Moreover, it has been shown that for low ξ values, the residuals appeared to become independent of ξ. This means that the regularization is neutral and only provides a way of interpolating the displacement at small scales when another (correlation-based) information is missing, thus preserving the well-posedness of the problem, without compromise to the image registration algorithm.

Practical Case

The material studied herein is a ferritic cast iron whose composition is 3.4 wt.% C, 2.6 wt.% Si, 0.05 wt.% Mg, 0.19 wt.% Mn, 0.005 wt.% S and 0.01 wt.% P. After casting and heat treatment (ferritization at 880 • C, followed by air cooling) the microstructure obtained is a ferritic matrix (98 wt.% ferrite and 2 wt.% perlite, average grain size 50 µm) with a 14% volume fraction of graphite nodules. Carbon (nodules) and iron (matrix) have atomic numbers that are different enough to give a strong X-ray attenuation contrast. Therefore, the spherical nodules, which are homogeneously distributed inside the matrix, are easily imaged by tomography and they can be used as natural markers for image correlation in the reconstructed 3-D images. The average nodule diameter of 45 µm is large enough to allow for a relatively low resolution voxel size to be used for tomography. This also enables one to observe millimetric samples, the only limitation for the sample size being the overall attenuation of the material that has to allow for, at least, 10% transmission of the incoming X-Ray beam.

The experiment was performed at the European Synchrotron Radiation Facility (ESRF) in Grenoble (France) on Beamline ID19. The synchrotron beam is parallel so that the voxel size only depends on the optics used. To obtain a 3D image of the specimen, six hundred radiographs (referred to as a scan) were recorded while the sample was rotating over 180 • along a vertical axis. A Fast Readout Low Noise (FReLoN) 14-bit CCD camera with a resolution of 2048 × 2048 pixels was used [START_REF] Labiche | The fast readout low noise camera as a versatile X-ray detector for time resolved dispersive extended X-ray absorption fine structure and diffraction studies of dynamic problems in materials science, chemistry and catalysis[END_REF].

The time required to acquire each image was set to 3 s, resulting in a total scan time of about 42 minutes.

A specially designed in situ fatigue testing machine that allows for loading and cycling (up to a frequency of 50 Hz) of the specimen [START_REF] Buffière | Three-dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray micro-tomography[END_REF] was used. The specimen was loaded to a maximum level of 273 N, and then unloaded. Five tomographic scans were acquired during which the load level was equal to 22, 150, 200, 250 and 273 N, and one scan upon unloading (load level: 21 N). In the following, only the scan corresponding to the maximum load level will be considered in addition to the reference scan.

Reconstruction of the tomographic data was performed with a conventional filtered back-projection algorithm [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. It provides a reconstructed volume with a 32-bit dynamic range that is proportional to the local attenuation coefficient. The 32-bit volume is subsequently re-encoded onto an 8-bit range to reduce data size and computation time. This truncation does not degrade significantly the raw data since reconstruction artifacts are of higher amplitude than the error induced by the former [START_REF] Davis | Artefacts in X-ray microtomography of materials[END_REF][START_REF] Ketcham | New algorithms for ring artefact removal[END_REF]. The region of interest has dimensions of 200 × 340 × 512 voxels, i.e., with a voxel size of 5.06 µm in the reconstructed images, 1.01 × 1.72 × 2.59 mm 3 .

Numerical tests are analyzed on a ROI in the center of the volume of size 23 × 23 × 23 supervoxels (i.e., 13, 824 × 3 degrees of freedom for a V-DVC approach) after a pyramidal filter is applied (i.e., a recursive construction that consists in averaging groups of 2×2×2 voxels to form a "super-voxel"). The lengths ξ b and ξ are initially set to 16 supervoxels. In the present case, the physical size of each considered (super)voxel (or element since V-DVC is used) is equal to 10.12 µm. This coarsening was necessary since the nodule spacing was too big to get a sufficiently rich texture. A first step was applied to correct for rigid body translations (up to a level of 19.5 supervoxels along the loading direction). Then the V-DVC procedure itself is run. The convergence in terms of displacement corrections is shown in figure 7; a rather fast rate is observed.

The total CPU time is equal to 15 min on a Core i7 (2.9 GHz) computer.

The next check is given by the correlation residual that measures the registration quality. Figure 8 shows that there is a change only for the first three iterations. Afterwards, for ξ b and ξ both equal to 16 voxels, a saturation is observed at a level of 1.6 % of the dynamic range. This value is significantly lower than previously observed with C8-DVC [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF][START_REF] Hild | Three dimensional analysis of a compression test on stone wool[END_REF] and reaches values commonly observed in Q4-DIC [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]. The fact that the dimensionless residuals tend to such low values is an additional indication that the results are trustworthy. This value can be further decreased down to 1 % when ξ b and ξ are set to 6 voxels. In comparison, the correlation residual that was obtained with the C8-DVC approach on the same ROI is equal to 3 % when 9 elements are used, meaning that while improving the robustness, V-DVC allows for decrease of the correlation residual by a factor of 3. This result further validates the V-DVC approach.

The correlation residual map is shown in Figure 9 when ξ b and ξ are set to 6 voxels.

The maximum value is equal to 12 gray levels for a dynamic range of the scans of 256 gray levels. Except at very few locations, the residual levels are very low everywhere.

This may indicate deviations from the assumed behavior or a small texture gradient that does not allow for a good (local) registration. The scan g corresponding to the deformed state is given in the same figure for comparison purposes. The location of the nodules generally disappears in the correlation residual map.

Figure 10 shows the deformed volume at convergence of the V-DVC procedure.

A smooth displacement field is observed on the external surface. This result is made possible by the use of the wavelength ξ b . The displacement field is also regular inside the volume as can be seen on the cut shown in Figure 10

Conclusion

A new regularized Digital Volume Correlation method was introduced herein. It allows for the measurement of displacements at a voxel level, exploiting the assumption of an elastic behavior at a small scale, and a boundary regularization. When such a regularization is not used, even global DVC is not reliable for small element sizes, and eventually does not converge for very small element sizes since the problem becomes ill-posed. The regularization is performed using the Equilibrium Gap Method.

The correlation residuals are extremely small (much smaller than using a global DVC algorithm for any element size) showing the ability of the method to provide a faithful (and mechanically admissible) displacement field. To further lower the residuals, the measurement technique contains an "on the fly" identification procedure, which was not used herein. It consists in considering the elastic parameters as unknowns to be determined as was proposed for 2D applications [START_REF] Leclerc | Integrated Digital Image Correlation for the Identification of Mechanical Properties[END_REF]. A first step will consist in evaluating the sensitivity of the measurement technique to the sought mechanical parameters. Element size (voxels) Standard displacement uncertainty (voxel)

List of Figures

C8-DVC Interpolation

  for the displacement component along the tensile direction. It is worth emphasizing that the overall amplitude of the displacement field over the analyzed volume is equal to 1 voxel, or 10 µm. To see more clearly the impact of boundary and bulk regularization, Figure 10 also shows the component of the displacement along the tensile direction in the mid section of the volume. The corresponding texture as well as the correlation residual field of the same section are shown for comparison purposes in Figure 9. It is worth remembering that the displacement fluctuations are due to the heterogeneous nature of the studied material at the voxel scale.

from 1 to 3 . 27 5 28 6 29 7 30 8 31 9

 32728293031 The frame depicts the edges of the undeformed volume. . Standard displacement uncertainty as a function of the element size in C8-DVC, and the cut-off length ξ when ξ b = 16 voxels for V-DVC. . . . Dimensionless correlation residual, η, with two DVC methods as a function of finite element sizes for C8-DVC, and cut-off wavelengths for the proposed V-DVC method. The C8-DVC curve was stopped at the left because of lack of convergence. . . . . . . . . . . . . . . . . . . . . . . . Convergence of the present method in terms of RMS of the displacement increment between two iterations, as a function of the number of iterations for the chosen ferritic cast iron tomographies (load levels: 22 and 273 N). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dimensionless correlation residuals as a function of the number of iterations for the chosen ferritic cast iron tomographies (load levels: 22 and 273 N). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Correlation residual map (left) and deformed scan (right) when ξ b and ξ are set to 6 voxels. Note that the dynamic range of the residual map is significantly lower than that of the deformed scan (256 gray levels). The frame depicts the edges of the undeformed volume. . . . . . . . . . 32 10 Deformed volume with displacement amplitude (left) and along the tensile direction (right) expressed in voxels when ξ b and ξ are set to 6 voxels. The frame depicts the edges of the undeformed volume. . . . . . . . .

Fig. 1

 1 Fig. 1 Standard displacement uncertainty as a function of the element size in C8-DVC. The interpolation (solid line) is described by a power law B/ α , where α = 1.42 is good agreement with Equation (12).

Fig. 2

 2 Fig. 2 Deformed state of the artificial material used to generate a slowly varying texture. The gray level is a growing affine function of the local Young's modulus that varies from 1 to 100.The frame depicts the edges of the undeformed volume.

Fig. 3

 3 Fig.3Dimensionless correlation residual, η, with two DVC methods as a function of finite element sizes, , for C8-DVC, and of cut-off wavelengths, ξ, for the proposed V-DVC method.The C8-DVC curve was stopped at the left because of lack of convergence.

Fig. 5

 5 Fig. 5 Standard displacement uncertainty as a function of the element size in C8-DVC, and the cut-off length ξ when ξ b = 16 voxels for V-DVC.

Fig. 7

 7 Fig.7Convergence of the present method in terms of RMS of the displacement increment between two iterations, as a function of the number of iterations for the chosen ferritic cast iron tomographies (load levels: 22 and 273 N).

Fig. 8

 8 Fig. 8 correlation residuals as a function of the number of iterations for the chosen ferritic cast iron tomographies (load levels: 22 and 273 N).

Fig. 10

 10 Fig. 10 Deformed volume with displacement amplitude (left) and along the tensile direction (right) expressed in voxels when ξ b and ξ are set to 6 voxels. The frame depicts the edges of the undeformed volume.

  The frame depicts the edges of the undeformed volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3 Dimensionless correlation residual, η, with two DVC methods as a function of finite element sizes, , for C8-DVC, and of cut-off wavelengths, ξ, for the proposed V-DVC method. The C8-DVC curve was stopped at the left because of lack of convergence. . . . . . . . . . . . . . . . . . . 26 4 Deformed state of material used for the second test case. The gray level is a growing affine function of the local Young's modulus that ranges

	1	Standard displacement uncertainty as a function of the element size
		in C8-DVC. The interpolation (solid line) is described by a power law
		B/ α , where α = 1.42 is good agreement with Equation (12). . . . . . . 24
	2	Deformed state of the artificial material used to generate a slowly varying
		texture. The gray level is a growing affine function of the local Young's
		modulus that varies from 1 to 100.

Acknowledgments

This work was funded under the grant ANR-09-BLAN-0009-01 (RUPXCUBE Project).

It was also made possible by an ESRF grant for the experiment MA-501 on beamline ID19. The scans were obtained with the help of Drs. J.-Y. Buffière, A. Gravouil, N. Limodin, W. Ludwig, and J. Rannou.