Alexandre Gallenne
email: alexandre.gallenne@obspm.fr

Antoine Mérand

Pierre Kervella

Julien H V Girard

Spatially extended emission around the Cepheid RS Puppis in near-infrared hydrogen lines Adaptive optics imaging with VLT/NACO ⋆

Keywords: circumstellar matter, Stars: variables: Cepheids, Stars: individual: RS Pup, Instrumentation: adaptive optics, Infrared: stars, Atmospheric effects

come

Spatially extended emission around the Cepheid RS Puppis in near-infrared hydrogen lines. Adaptive optics imaging with VLT/NACO

Introduction

For almost a century, Cepheid stars have been used as distance indicators thanks to the well-known period-luminosity relation (P-L), also called Leavitt's law. A good calibration of this relation is necessary to obtain an unbiased estimate of the distance. After the detection of circumstellar envelopes (CSEs) around many Cepheids in the near-and mid-infrared (Kervella et al. 2006;Mérand et al. 2006Mérand et al. , 2007;;Kervella et al. 2009a), it was shown that these CSEs could bias the angular diameter determination using near-infrared interferometry (Mérand et al. 2006). A circumstellar envelope could have an impact on the distance estimate if we do not take into account its flux contribution. It could lead to a bias both in interferometric and photometric measurements and therefore influence the accuracy of the estimated distance. Because Cepheids are used as standard candles to estimate distances in the universe, it is important to understand the properties of the CSEs around this kind of stars.

We report new observations with VLT/NACO in the IB 2.18 and NB 1.64 filters using a fast imaging mode in order to characterize the extended emission around RS Puppis (HD 68860) at sub-arcsecond angular resolution (≈0.01 pc at RS Pup distance). Data and reduction methods are described in Sect. 2. In Sects. 3 and 4 we present two different data processing analyses: the shift-and-add method, and a method based on a statistical analy-sis of the speckle noise to enable the detection of the CSE in the speckle cloud surrounding the star's image.

This paper is focused on RS Pup's envelope, which has been known for many years (Westerlund 1961;Havlen 1972) for it is very extended: it has a spatial extension of around 2 ′ in V band, which corresponds to about one parsec at RS Pup's distance (about 2 kiloparsecs; Kervella et al. 2008). Its origin is not precisely known; it is especially unclear whether it is made of material ejected by the star or of a material that predates the formation of the star. Its actual shape has been widely discussed: discrete nebular knots have been invoked by Kervella et al. (2008) to estimate the distance using light echoes, whereas Feast (2008) argued that a disk-like CSE will bias this kind of estimate. Alternately, a bipolar nebula has been considered by Bond & Sparks (2009), making a ratio between two images of different epochs. We were able to spatially resolve the CSE close to the star. The morphology of the CSE is discussed in Sect. 4.3.

Observations and data reduction

RS Pup was observed with the NACO instrument installed at the Nasmyth B focus of UT4 of ESO VLT. NACO is an adaptive optics system (NAOS, [START_REF] Rousset | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF]) and a high-resolution near IR camera (CONICA, [START_REF] Lenzen | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF], working as imager or as spectrograph in the range 1-5 µm. We used the S13 camera for our RS Pup observations (FoV of 13.5 ′′ ×13.5 ′′ and 13.26 ± 0.03 mas/pixel, [START_REF] Marengo | The Evolving ISM in the Milky Way and Nearby Galaxies[END_REF] with the narrow band filter NB 1.64 (λ 0 = 1.644 µm with a width ∆λ = 0.018 µm) and the intermediate band filter IB 2.18 (λ 0 = 2.180 µm with a width ∆λ = 0.060 µm). We chose the cube mode in order to apply the shift-and-add technique. We chose windows with 512×514 pixels and 256×258 pixels respectively for the narrow and the intermediate band filters. Each cube (10 for RS Pup and 4 for the reference star) contains 460 images at 1.64 µm and 2 000 at 2.18 µm. Each frame has the minimum integration time and readout type for this mode, i.e 109 ms and 39 ms respectively for the 1.64 and 2.18 µm cubes. With these DITs and filters we are in a quasi-monochromatic regime with a very good speckle contrast (Roddier 1981). The large number of exposures allows us to select the frames with the brightest image quality.

Data were obtained on 2009 January 7, with the atmospheric conditions presented in Table 1. The PSF (point spread function) calibrator stars were observed just before or immediately after the scientific targets with the same instrumental configurations. Each raw image was processed in a standard way using the Yorick 1 language: bias subtraction, flat-field and bad pixel corrections. The negligible sky background was not subtracted.

To compare the AO system behavior and artefacts on RS Pup, we used the star Achernar and its PSF reference star observed the same night. These data were also acquired in cube mode with the same camera in the NB 1.64 and NB 2.17 (λ 0 = 2.166 µm with a width ∆λ = 0.023 µm) filters. Only one cube was acquired per filter, containing 20 000 frames for Achernar and 8 000 for the calibrator. The window detector size is 64×66 pixels in both filters with an integration time of 7.2 ms. The slightly different second filter and the shorter exposure time compared with RS Pup is not critical for our analysis because we are interested in flux ratios. The smaller window size cuts the AO PSF wings, but we will show in Sect. 3 that this does not affect our conclusion significantly.

Once the pre-processing was over, we carried out a precentering as well as a sorting according to the maximum intensity of the central peak to reject the images for which the correction was the least efficient. On the one hand, we applied the shift-and-add technique to obtain the best angular resolution possible. On the other hand, we applied a statistical study aimed at extracting a possible diffuse (i.e. completely resolved) component inside of the speckle cloud of the central core (the star). These methods are explained in Sects. 3 and 4.

Shift-and-add process

In this section we use the shift-and-add method and then compare the RS Pup images with another star observed on the same night for which we do not expect an envelope.

RS Pup

Both stars, RS Pup and HD 74417, are unresolved by the telescope: diffraction limits for an 8.2m UT are 55 mas and 41 mas in filters IB 2.18 and NB 1.64 respectively, while angular diameters are approximately one milli-arcsecond (1 mas) for both stars.

The shift-and-add method was originally proposed by Bates & Cady (1980) as a technique for speckle imaging. Applied to NACO's cube mode, it enhances the Strehl ratio by selecting the frames that are the least altered by the atmospheric turbulence. It also reduces the halo contribution when used with adaptive optics by selecting the best AO-corrected frames. The efficiency of the cube mode vs. the standard long exposure mode is discussed in Kervella et al. (2009b). [START_REF] Girard | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF] showed for the 1 http://www.maumae.net/yorick/doc/index.php RS Pup's reference star at 2.18 µm that by using this method we have an absolute gain of 9 % in Strehl ratio.

Our processing steps are as follows: we select the 10 % best frames according to the brightest pixel (as a tracer of the Strehl ratio) in our 10 data cubes (4 for the calibrator). We then spatially resample them by a factor 4 using a cubic spline interpolation and co-align them using a Gaussian fitting on the central core, at a precision level of a few milliarcseconds (this method is described in detail in Kervella et al. 2009a). Each cube is then averaged to obtain ten final average images (4 for the calibrator). We finally compute the mean of the 10 average images of RS Pup (or the 4 for the calibrator).

The result of this method is shown in Fig. 1 for the two filters. The two upper frames are HD 74417 and RS Pup as seen through the IB 2.18 filter, while the two lower are in the NB 1.64 filter. The spatial scale in all images is 3.3 mas/pix and we chose a logarithmic intensity contrast with adjusted lower and higher cut-off to show details in the halo. On the 1.64 µm image of RS Pup compared with the HD 74417 image we can see that there is a resolved circumstellar emission up to ∼ 1 ′′ from the Cepheid. This feature is fainter in the 2.18 µm filter. In both filters, because of this emission, diffraction rings are less visible around RS Pup than around HD 74417.

For each final co-added frame (i.e the image shown in Fig. 1) we computed a ring median and normalized to unity. The resulting curves are presented in Fig. 2 for two filters. The error bars are the standard deviation of values inside rings and are plotted only for few radii for clarity. We notice a significant difference between the two curves, more pronounced in the NB 1.64 filter.

Achernar

We carried out the same processing for Achernar and its reference star δ Phe (HD 9362). These two stars are unresolved by the UTs (angular diameters are approximately 2 mas). With this window size (64×66 pixels) the field of view is four times smaller than in the RS Pup's frames, and even if we can also see some Airy rings in all the images, the background is dominated by the AO halo.

We plot in Fig. 2 the radial profiles of the stacked images, normalized to 10 (for clarity in the graph).

Halo fitting

We now investigate the presence of a circumstellar envelope by studying the radial profile of the stacked images we obtained for both RS Pup and its PSF calibrator and compare them with the observations of Achernar. The main idea behind this study is to compare the variations of encircled energy (the ratio of the core's flux to the total flux) to the expected variations caused by atmospheric condition changes. If RS Pup has indeed a circumstellar envelope, its encircled energy will be significantly lower than that of our reference stars.

We use an analytical function to fit the residual light out of the coherent core (i.e. the halo). We used a turbulence degraded profile from Roddier (1981)

I halo (r) = f 0.488 ρ 2        1 + 11 6 r ρ 2        -11 6 , (1)
where the fitting parameters are the FWHM ρ and f proportional to the flux of the envelope. In this case, the parameter ρ corresponds to the equivalent seeing disk for long exposure. NACO's 2009-01-07 05:26:57-05:31: Table 1. Log of the NACO observations. seeing DIMM are the seeing measurements in V from the Differential Image Motion Monitor (DIMM) station, while seeing λ are interpolated to our wavelength. The uncertainties denote the standard deviations over the number of cubes. The start-end UTC indicated denotes the time when the first and the last cube was recorded. The phase φ of the Cepheid was computed using (P, T 0 ) from Kervella et al. (2008). AM denotes the airmass, while Nbr is the number of cubes recorded. r 0 and t 0 are the Fried parameter and the coherence time given by the AO real time calculator. resulting core, even for an unresolved source, is fairly complex and suffers from residual aberrations, this is why we decided to fit the halo at a large radial distance of the core. We will later define the core as the residual of the fit of the halo function to our data. We show the resulting fits in Fig. 4. They were applied for r > 0.26 ′′ at 2.18 µm and 2.17 µm while the range was for r > 0.2 ′′ at 1.64 µm for RS Pup (and its calibrator) and r > 0.26 ′′ for Achernar (and its calibrator). The black lines are our data (ringmedian), while the gray lines represent the adjusted curves. The vertical scale is logarithmic and normalized to the total flux. We also evaluated the normalized encircled energy (E enc ), i.e. the ratio of the flux in the core (i.e within 1.22λ/D) to the total flux. We can already see from these curves that the flux in the halo is different in the images of RS Pup and Achernar. These differences can be linked to the emission of the CSE around RS Pup or a modification of the seeing between the two observations. This point will be discussed bellow with a statistical study of the encircled energy.

For RS Pup in the NB 1.64 filter we found from our fit comparing the two FWHMs, i.e (ρ RSPup -ρ HD74417 )/ρ HD74417 , a change of halo size of 11% and encircled energies E enc (RSPup) = 0.327 and E enc (HD74417) = 0.452 (see Table 2). This yields an encircled energy ratio ε ∼ 72 %. The local DIMM station registered a seeing variation of around 26 %. Different values can be explained by the station location, which is not inside the telescope dome, but a few tens of meters below, rendering it more susceptible to ground layer effects that are known to affect the seeing [START_REF] Sarazin | [END_REF]. As a result the UT seeing is consistently better than the one measured by the DIMM. For IB 2.18 we show the fitting curves in Fig. 4 and the resulting fit parameters in Table 2. Our fit gives a seeing variation of around 9 %, while the DIMM measurements give a 29 % variation, leading to the same conclusion as before about the accuracy of the DIMM station seeing estimates. The encircled energies measured in IB 2.18 are E enc (RSPup) = 0.417 and E enc (HD74417) = 0.517, giving an encircled energy ratio ε ∼ 80 %.

We evaluated the uncertainties according to the variations of the encircled energy in the cubes. Because we already selected the 10 % best frames, we kept these samples and plotted the histogram of the encircled energy to check the distribution function (Fig. 3). At the first order, the E enc distribution can be approximated by a Gaussian distribution. We then computed the relative standard deviation σ(E enc )/E enc for both stars. We assessed the uncertainty of the encircled energy ratio ε with

σ(ε) ε = σ(E enc ) E enc 2 RSPup + σ(E enc ) E enc 2 HD74417 , (2)
We assume that the relative standard deviation of the encircled energy is a good approximation of the seeing variations. We measured a relative standard deviation [σ(E enc )/E enc ] RSPup ∼ 8 % at 2.18 µm and ∼ 10 % at 1.64 µm while [σ(E enc )/E enc ] HD74417 ∼ 4 % and ∼ 6 % respectively. We also evaluated the statistical uncertainties with a bootstrap method that gave σ ∼ 0.5 % for both stars. These statistical uncertainties do not take into account the seeing variations. We therefore chose Eq. 2 as a reliable estimate of the final uncertainty. This gives σ(ε)/ε ∼ 9 % and ∼ 12 % respectively in IB 2.18 and NB 1.64.

For Achernar, the fit gives a seeing variation ∼ 3% both in NB 1.64 and IB 2.17. The encircled energies measured in IB 2.18 are E enc (Achernar) = 0.495 and E enc (δ Phe) = 0.505, giving an encircled energy ratio ε ∼ 98 % (see Table 2). In NB 1.64 we have 0.374 and 0.376 respectively, which gives ε ∼ 99 %. The encircled energies are nearly identical, giving ratios close to 100 %. Using Eq. 2 as previously, we found a relative standard deviation ∼ 15 % in IB 2.17 and ∼ 11 % in NB 1.64. Because Achernar does not have an extended circumstellar component in NB 1.64 and IB 2.17, we attribute these variations to the Strehl variations and the AO correction quality. The encircled energy ratio seems stable between the two stars at the 2 % level in IB 2.18 (see Table 2).

We can notice that the normalized curves of Achernar and its calibrator almost overlap, and this is what is expected for unresolved stars without envelope taken in same atmospheric and instrumental conditions. It is not the case for RS Pup (see Fig. 2 and4). The additional emission coming from the envelope of the Cepheid adds a component that is clearly not negligible. In the same atmospheric conditions and AO correction between the star and its calibrator, the encircled energy ratio should be close to 100 %. That is what we obtain from Achernar's data. Scaling RS Pup's images to the same window size (i.e. r 0.42 ′′ ) slightly changes E enc because we normalized to the total flux, but we still have a discrepancy of 100 % for the encircled energy ratio (ε ∼ 79 % and ε ∼ 87 % respectively in NB 1.64 and IB 2.18). A smaller window on RS Pup makes the total flux lower because we do not integrate all the envelope flux, consequently the normalized encircled energy is larger and therefore the E enc ratio as well.

The drop in E enc of RS Pup compared with the statistics over Achernar's data and the fact that this drop is larger than the one directly estimated from the seeing halo variations, motivate the following section where this drop is analyzed under the assumption that it is caused by CSE emission.

Encircled energy variation

In this part we evaluate the flux of the envelope using the fact that the theoretical encircled energy ratio should be 100% and that an additional component to the flux should decrease this ratio. The encircled energy of RS Pup can be written as

E enc = F coh + αF env F + F env ,
where F is the total flux of the star, without the extended flux, F coh is the flux in the coherent core, F env is the emission of the environment and α is the fraction of the envelope flux that lays in the core. For simplicity, we will assume α << 1, which means that we assume that the CSE is significantly larger than the core. Accordingly we can write

E enc ≈ F coh /F 1 + F env /F ,
Assuming that the encircled energy is stable in NACO for our observing conditions, we can calibrate the previous formula to obtain the extended emission flux:

ε = E enc (Cepheid) E enc (reference) ≈ 1 1 + F(env)/F(Cepheid)
The relative emission F(env)/F(Cepheid) can be estimated with the encircled energy ratio from our measured values (Table 2). We can also take the uncertainty values assessed in the previous section, defined as being the relative standard deviation on 10% of best frames in the cubes. We found a relative extended emission around RS Pup contributing to 38 ± 17 % of the star's flux in the NB 1.64 band and 24 ± 11 % in IB 2.18.

Morphological analysis

We now present a qualitative analysis of the envelope morphology. In this analysis, based on the statistical properties of the noise owing to speckle boiling, we presented the ratio σ(F)/F as a function of the factor (g). This parameter is intrinsic to speckle temporal and spatial structures, is not expected to depend on azimuth, and is proportional to the envelope flux. We mapped this factor, and any azimuthal asymmetry should be attributed to the morphology of the envelope.

Analysis method

The principle of the method is to recover the flux of the envelope in the short exposures. The difference between a simple unresolved star and a star surrounded by an envelope is that the envelope will prevent the flux between the speckles from reaching zero.

Statistical descriptions of stellar speckles have been proposed by many authors (Racine et al. 1999;Fitzgerald & Graham 2006), who showed that speckle noise is not negligible. The distribution of speckle statistics in AO were also studied by Canales & Cagigal (1999). For example, Racine et al. (1999) (hereafter Ra99) showed that the speckle noise dominates other noise sources within the halo of bright stars corrected with AO systems. We know that the total variance in the cube short exposures after the basic corrections (i.e corrected for bias, bad pixels, flat-field and recentered) is

σ 2 (r) = σ 2 s (r) + σ 2 ph (r) + σ 2 RON + σ 2 Strehl (r), (3)
where σ 2 s is the variance of the speckle noise, σ 2 ph the variance of the photon noise, σ 2 RON the variance of the readout noise and σ 2

Strehl the noise owing to Strehl ratio variations. We neglected the variance of the sky photon noise, which is negligible in our wavelength for bright targets. Because the Strehl ratio variations occur in the coherent core and because we are particularly interested in the halo part, we omit these variations in the rest of this section.

We present for one calibrating cube the variance observed in the top panels of Fig. 5 for both filters. We plotted the radial profile computed using median rings on all azimuths. The dashed curve denotes the variance of the readout noise. We can estimate some of the components in Eq. 3 from our data. For example, from the stacked image that we obtained by summing the cube, we obtain the core and halo profile, as we did in Sect. 3. We then obtain the photon noise variance using its Poissonian property, knowing that the variance is equal to the mean. We define for the stacked image the radial intensity profile I(r) and we estimate the term σ 2 ph of Eq. 3 with

σ 2 ph (r) = γ I(r) = γ I(r) N ,
where N is the number of frames in the cube (1999 frames in IB 2.18 and 457 in NB 1.64) and γ the gain of the detector in order to convert noises in electron unit. The readout noise is estimated from the outer part of the images (see in Fig. 5 the top panels). We estimate the speckle noise from the total variance of the cubes as

σ 2 s (r) = σ 2 (r) -σ 2 ph (r) -σ 2 RON ,
In the bottom panels of Fig. 5 we present the relative importance of the different noise sources. We see that the different variances show a different radial behavior, and we see in particular that the photon noise of the core never dominates. More precisely, we can define two regimes:

σ(r) ≃ σ s (r) for r r 1 , in this regime, the speckle noise dominates the other variances (σ s > 3 σ RON )

-σ(r) ≃ σ 2 s (r) + σ 2 RON for r r 1
with r 1 ≃ 0.22 ′′ in NB 1.64 band and r 1 ≃ 0.24 ′′ in IB 2.18. The fairly weak (∝ λ -12/5 ) wavelength dependency of the speckle noise may explain the same order of magnitude between the limit radius in IB 2.18 and NB 1.64. Indeed there is another regime in the range 0 r < r 1 where the Strehl variations occur. The PSF core is mainly affected by these variations. As said previously, we are interested in the PSF halo to detect the envelope between the speckles, which is why we neglect this term.

Speckle noise invariant

We use an interesting property of the speckle noise, namely that its variance is proportional to the squared flux in the cube (Ra99), F/σ is therefore an invariant. More precisely:

σ s = g(r, λ, ∆t, S , ...) F ⋆ (λ),
where the function g is a function depending on the wavelength λ, the integration time ∆t, the Strehl ratio S and atmospheric parameters.

We checked this property by plotting in Fig. 6 the mean g function for each star, i.e we evaluated a g function for each cube that we then averaged (the error bars correspond to the standard deviation between the cubes). We took care to remove the flux contribution from the circumstellar envelope with the values found in Sect. 3. Obviously there is no clear difference between the two stars except in the core. This can be explained by the Strehl ratio variations we omitted before in the core. We can assess these variations by evaluating the encircled energy for each frame of each calibrating cube (here we take 100 % of frames, without selection) and by using it as an estimator of the Strehl ratio in function of time. By fitting a Gaussian distribution to E enc for each star and both filters, we estimate the relative standard deviation for RS Pup to be σ(E enc )/E enc ∼ 15 % in IB 2.18 and ∼ 17 % in NB 1.64. For HD 74417 we found ∼ 8 % and ∼ 9 % respectively at 2.18 µm and 1.64 µm.

The total invariance of the g function, i.e for the complete radius, depends on the Strehl ratio variations. Nevertheless, because these variations occur in the coherent core, we can neglect We can see a larger difference in the core owing to the Strehl variations. From r > r 1 the functions are similar at a 9 % level in NB 1.64 and at 6 % in IB 2.18. them if we go far away from the central part. For r > r 1 we found a mean absolute relative difference (see Fig. 6) of 5 % in IB 2.18 and 8 % in NB 1.64. So for r > r 1 the invariance of the g function is verified with a good accuracy.

We can measure g using our reference star. Assuming it is stable in time, we can define the observable Γ as the average of the flux of a given pixel in the cube, divided by the standard deviation. Defining F cal = αF RS Pup = αF ⋆ and making the difference between Γ sci and Γ cal (where the index sci and cal refer to RS Pup and HD 74417 respectively) we have

ξ = Γ sci -Γ cal = F ⋆ + F env g 2 ⋆ F 2 ⋆ + σ 2 RON - αF ⋆ g 2 ⋆ α 2 F 2 ⋆ + σ 2 RON , (4)
All variables in Eq. 4 are functions of the pixel coordinates (x, y), but we omit this notation for clarity.

The parameter α can be estimated with the results from Sect. 3. The flux ratio between the two final images of RS Pup and HD 74417 is a function of α and F(env)/F ⋆ :

F(RS Pup) F(HD 74417) = F ⋆ + F env α F ⋆ = 1 + F env /F ⋆ α , (5)
By measuring the total flux in our final images (Sect. 3) and using F env /F ⋆ estimated in Sect. 3.4 we get α = 1.06 ± 0.13 at 1.64 µm and α = 0.96 ± 0.08 at 2.18 µm. We also used the total flux variations in our 10 (for RS Pup) and 4 (for the reference) averaged images to assess the uncertainties, but they turn out to The white circle represents the r 1 limit of Eq. 7 separating the two regimes. The scale intensity is logarithmic. be small (< 2 %) compared to F env /F ⋆ uncertainties. As a first approximation we can simplify Eq. 4 with α ∼ 1 to obtain

ξ = F env g 2 ⋆ F 2 ⋆ + σ 2 RON . ( 6
)
ξ is thus directly proportional to the flux of the CSE. We expect this function to be null if there is no envelope. The denominator is a radial parameter, therefore any spatial asymmetry in ξ will be linked to an asymmetry of the CSE.

From our different data cubes we computed a mean ξ by subtracting to each Γ sci a mean Γ cal , i.e:

ξ(x, y) = 1 n n j=1 (Γ sci,j -Γ cal ) and Γ cal = 1 m m i=1 Γ cal,i ,
where n and m are the number of cube for RS Pup and HD 74417 respectively. This process provides us with an average image proportional to the flux of the CSE in each filter.

Symmetry of the CSE

In Fig. 7 we show the ξ function for the NB 1.64 and IB 2.18 band. For radii larger than r 1 , we expect the function to be null if there is no CSE (see Eq. 6), whereas we clearly obtain a positive image, indicative of a likely detection of a CSE. The central part of the image with r < r 1 (white circle for r = r 1 ) denotes the part where the previously defined function g is not invariant because of Strehl ratio fluctuations. It is difficult to estimate and therefore does not supply reliable informations close to the star. This method is useful for studying the morphology of the CSE, in particular whether or not it has a central symmetry, like one expects for a uniform shell structure or a face-on disk; or an asymmetrical shape, like in an inclined disk-like structure for instance. Despite adaptive optics artefacts, the images in Fig. 7 seem to present a uniform intensity distribution in agreement with a shell or a face-on disk structure. More precisely, we estimated a symmetry level using a residual map obtained by subtracting 90 degree rotated versions of themselves from the images presented in Fig. 7. The result of this subtraction is shown in Fig. 8. In the initial image we assessed the average value in a ring, while in the subtracted image (0 • -90 • ) we compute an average value in a small window (see Fig. 8). Making the ratio of these two averages (i.e. the average value in the small window over the average value in the ring), we exclude a departure from symmetry larger than 7 % in NB 1.64 and 6 % in IB 2.18. Top panels: on the left is plotted the ξ parameter at 2.18 µm with a ring of radius r 1 < r < r 2 = 0.92 ′′ where we computed a mean value. On the right we plot the difference of the initial (0 • ) and 90 • rotation images where we cut a small window to assess a mean residual value. The ratio between the mean residual value and the mean value give us a level of symmetry. Bottom panels: same process at 1.64 µm with a ring of radius r 1 < r < r 2 = 1.1 ′′ . The scale intensity is logarithmic.

Conclusion

We present a study of the close environment (1000-2000 AU) of the Cepheid RS Pup using near-infrared AO-assisted imaging. We used two different techniques to 1) measure the photometric emission due to the CSE and 2) obtain an image proportional to the CSE through a statistical analysis of the large number of short exposures in which speckles are still visible.

Using a simple image modeling method, based on the unresolved PSF calibrator, we estimated (Sect. 3) a likely photometric emission of 38 ± 17 % and 24 ± 11 % of the star flux in the NB 1.64 and IB 2.18 bands respectively for RS Pup's CSE. Apart from the actual interest to better understand the mass loss in Cepheids, the presence of such an extended emission may have an impact on the distance estimate with the Baade-Wesselink method.

Using an original statistical study of our short exposure data cubes, we qualitatively showed (Sect. 4) that the envelope in the two bands (NB 1.64 and IB 2.18) has a centro-symmetrical shape (excluding a departure from symmetry larger than 7 %), suggesting that the CSE is either a face-on disk or a uniform shell.

Our observations were obtained in two filters isolating the Brackett 12-4 (1.641 µm, NB 1.64 filter) and Brackett γ 7-4 (2.166 µm, IB 2.18 filter) recombination lines of hydrogen. We therefore propose that our observations show the presence of hydrogen in the CSE of RS Pup. The star is naturally relatively faint in these spectral regions, because of deep hydrogen absorption lines in its spectrum. This appears as a natural explanation to the high relative contribution from the CSE compared to the stellar flux. Another contribution to the CSE flux may also come from free-free continuum emission, or Rayleigh scattering by circumstellar dust. However, this contribution is probably limited to a few percents of the stellar flux. The high relative brightness of the detected CSE therefore implies that the spectral energy distribution of the detected CSE is essentially made of emission in hydrogen spectral lines instead of in the form of a continuum. RS Pup's infrared excess in the broad NB 1.64 and IB 2.18 bands is limited to a few percents (Kervella et al. 2009a), which excludes a continuum CSE emission of 40% in the NB 1.64 filter and 24% in the IB 2.18 filter as measured in our NACO images.

Our discovery is also remarkably consistent with the nonpulsating hydrogen absorption component detected in the visible Hα line of RS Pup by Nardetto et al. (2008), together with a phase-dependent component in emission. The extended emission we observe in the infrared is therefore probably caused by the same hydrogen envelope as the one that creates the Hα emission observed in the visible. It should also be mentioned that our result is consistent with the proposal by Kervella et al. (2009a) that the envelope of RS Pup is made of two distinct components: a compact gaseous envelope, which is probably at the origin of the observed NACO emission, and a very large (≈ 2 ′ ) and cold (≈ 45 K) dusty envelope (Kervella et al. 2008) that is mostly transparent at near-infrared wavelengths owing to the reduced scattering efficiency. These two envelopes appear to be largely unrelated because of their very different spatial scales. However, it should be noted that [START_REF] Marengo | The Evolving ISM in the Milky Way and Nearby Galaxies[END_REF] recently detected large envelopes around several Cepheids using Spitzer. The presence of hydrogen close to the Cepheids of their sample could be investigated at high angular resolution using the observing technique presented here.

  43 -NB 1.64 HD 74417 0.54 ± 0.01 ′′ 0.43 ± 0.01 ′′ 1.04 4 20.3 ± 0.1 5.1 ± 0.1

Fig. 1 .

 1 Fig. 1. Final frames from shift-and-add process. In the two upper panels the PSF calibrator star HD 74417 and RS Pup are shown in the IB 2.18 filter. The bottom panels show the same stars in the NB 1.64 filter. The scale intensity is logarithmic.

Fig. 2 .Fig. 3 .

 23 Fig. 2.On the left are presented the radial profiles at 2.18 µm for RS Pup, Achernar and their reference stars. On the right the same at 1.64 µm. Radial profiles were computed using the median value on one-pixel thick circular rings. The gray color denotes Achernar and δ Phe profiles, shifted by a factor 10 for clarity. The scale is logarithmic for both axes.

Fig. 4 .

 4 Fig. 4.Halo fit for all filters with logarithmic ordinates. The fluxes are normalized to the total flux. Left panel: halo fit from the radial profile of images in the two filters for RS Pup. Right panel: halo fit from the radial profile of images in the two filters for Achernar. The fit in IB 2.18 was for r > 0.26 ′′ and in NB 1.64 for r > 0.2 ′′ for RS Pup and for r > 0.26 ′′ for Achernar. The dashed lines denote the reference stars, while the thick lines denote the Cepheid and Achernar. The fitted curves are plotted in gray.

Fig. 5 .

 5 Fig. 5.Left panel: at the top we plot the variance (circles) and the average value (triangle) of one data cube of HD 74417 at 2.18 µm. Bottom: contribution of the different noise sources. Right panel: the same for the 1.64 µm filter. The scale for all plots is logarithmic in ordinates. The dotted curves show the 3σ limits of the readout noise, while the gray curves denote the halo fit.

Fig. 6 .

 6 Fig.6. Mean g function for each star. The error bars correspond to the RMS. The bottom plot represents the absolute relative difference. We can see a larger difference in the core owing to the Strehl variations. From r > r 1 the functions are similar at a 9 % level in NB 1.64 and at 6 % in IB 2.18.

Fig. 7 .

 7 Fig.7. ξ parameter (Eq. 6) proportional to the envelope flux at 2.18 µm (left) and at 1.64 µm (right). The white circle represents the r 1 limit of Eq. 7 separating the two regimes. The scale intensity is logarithmic.

Fig. 8 .

 8 Fig.8. Estimate of the symmetry level. Top panels: on the left is plotted the ξ parameter at 2.18 µm with a ring of radius r 1 < r < r 2 = 0.92 ′′ where we computed a mean value. On the right we plot the difference of the initial (0 • ) and 90 • rotation images where we cut a small window to assess a mean residual value. The ratio between the mean residual value and the mean value give us a level of symmetry. Bottom panels: same process at 1.64 µm with a ring of radius r 1 < r < r 2 = 1.1 ′′ . The scale intensity is logarithmic.

Table 2 .

 2 Best

	filter	name	ρ ( ′′ )	f	E enc
	NB 1.64	RS Pup	0.46 ± 0.01 0.58 ± 0.01 0.327
		HD 74417 0.43 ± 0.01 0.43 ± 0.01 0.452
		Achernar	0.49 ± 0.01 1.00 ± 0.01 0.374
		δ Phe	0.46 ± 0.01 0.97 ± 0.01 0.376
	IB 2.18	RS Pup	0.56 ± 0.01 0.45 ± 0.01 0.417
		HD 74417 0.50 ± 0.01 0.34 ± 0.01 0.517
	NB 2.17	Achernar	0.52 ± 0.01 0.59 ± 0.01 0.495
		δ Phe	0.54 ± 0.01 0.61 ± 0.01 0.505

-fit results for the residual halo where. ρ is the FWHM and f is a parameter proportional to the flux of the halo normalized to the total flux. E enc represents the flux in the central core (i.e. within 1.22λ/D) to the total flux. Note: The fitted I halo analytical function is not normalized to f , i.e I halo rdr f .

Acknowledgements. We thank Drs. Jean-Baptiste Le Bouquin and Guillaume Montagnier for helpful discussions. We also thank the ESO Paranal staff for their work in visitor mode at UT4. We received the support of PHASE, the high angular resolution partnership between ONERA, Observatoire de Paris, CNRS, and University Denis Diderot Paris 7. This work made use of the SIMBAD and VIZIER astrophysical database from CDS, Strasbourg, France and the bibliographic informations from the NASA Astrophysics Data System. Data processing for this work have been done using the Yorick language, which is freely available at http://yorick.sourceforge.net/.


Eric Darrigrand
email: eric.darrigrand-lacarrieu@univ-rennes1.fr

Peter Monk
email: monk@math.udel.edu

Combining the Ultra-Weak Variational Formulation and the Multilevel Fast Multipole Method *

Keywords:

Because of its practical significance, many different methods have been developed for the solution of the time-harmonic Maxwell equations in an exterior domain at higher frequency. Often methods with complimentary strengths can be combined to obtain an even better method. In this paper we provide a numerical study of a method for coupling of the Ultra-Weak Variational Formulation (UWVF) of Maxwell's equations, a volume based method using plane wave basis functions, and an overlapping integral representation of the unknown field to obtain an exact artificial boundary condition on an auxiliary surface that can be very close to the scatterer. Combining the new algorithm with a multilevel fast multipole method we obtain an efficient volume based solver with an exact auxiliary boundary condition, but without the need for singular integrals.

Introduction

At progressively higher frequencies, the problem of approximating the time-harmonic Maxwell equations in an exterior domain is increasingly challenging because of the need for an increasingly large number of unknowns per wavelength to maintain accuracy. Although many numerical schemes exist for this problem we will focus on hybridizing two approaches with complementary advantages: a volume based method using finite element grids and plane wave basis functions and an overlapping integral equation formulation. Volume methods are derived from a three dimensional discretization of a portion of the exterior domain close to the scatterer and limited by an artificial boundary. If a convenient local boundary condition is used the distance between the scatterer and the artificial boundary must be taken sufficiently large to ensure accuracy. In addition the artificial boundary must be convex. This leads to the expensive discretization of free space. However volume based methods can easily handle non-constant material coefficients, and give rise to sparse matrices. In contrast an integral formulations enables the restriction of the problem to the surface of the scatterer. However, integral techniques require a homogeneous external media, and the evaluation of singular integrals. Both methods are often coupled leading to a hybrid method that inherits the strengths of both component methods.

In this paper, we focus on coupling the Ultra-Weak Variational Formulation (UWVF) and a direct boundary integral representation using a multilevel fast multipole method. The UWVF is a volume based numerical method for solving the time-harmonic Maxwell system on a bounded domain developed by Després and Cessenat [START_REF] Cessenat | Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation[END_REF][START_REF] Cessenat | Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniquesproblèmes de Helmholtz 2D et de Maxwell 3D[END_REF][START_REF] Huttunen | Solving Maxwell's Equations Using the Ultra Weak Variational Formulation[END_REF]. It uses local plane wave solutions on a finite element mesh to approximate impedance traces of the field on the skeleton of the mesh. By varying the number of plane wave basis functions from element to element the UWVF can discretize the electromagnetic field with a coarser volume mesh in comparison to more classical methods low order finite elements or finite differences. It can handle general piecewise constant material properties. However, to approximate scattering on an unbounded domain, the UWVF requires an artificial boundary denoted Γ ext sufficiently far from the scatterer. The original UWVF uses a low order absorbing boundary condition on Γ ext (but see [START_REF] Huttunen | Solving Maxwell's Equations Using the Ultra Weak Variational Formulation[END_REF] for use of the perfectly matched layer). To overcome the drawback of needing a large distance between the scatterer and a convex artificial boundary, we suggested using an integral representation of the unknown field on Γ ext .

In particular we shall use an overlapping scheme. A standard Galerkin volume equation combined with an overlapping integral equation system was first proposed for Maxwell's equations by Hazard and Lenoir [START_REF] Hazard | On the Solution of Time-Harmonic Scattering Problems for Maxwell's Equations[END_REF]. A finite element discretization of the method using edge elements was then analyzed by Hsiao, Monk and Nigam [START_REF] Hsiao | Error analysis of a finite element-integral equation scheme for approximating the time-harmonic Maxwell system[END_REF]. In particular, Hazard and Lenoir proposed to use an integral representation of the field near the artificial boundary Γ ext computed from unknown field values on a second artificial boundary Σ taken closer to the boundary of the scatterer Γ int (see Fig. 1 for a graphic of the relationship of various boundaries). This couples the Cauchy data on Σ to the solution on Γ ext . The main constraint is that the domain between Σ and Γ ext be homogeneous (i.e. the background medium). However the domain within Σ can be inhomogeneous. Upon discretization, the well-known integral operators involved in this representation lead to expensive evaluations by direct means (but do not involve singular integrals, or near interaction terms). Indeed the expense and memory needed to compute the linear system resulted in no numerical results being presented in [START_REF] Hsiao | Error analysis of a finite element-integral equation scheme for approximating the time-harmonic Maxwell system[END_REF]. Clearly the problem can potentially be solved using a fast multipole implementation of the integral equations. Note that for standard edge finite elements this was the approach used by Liu and Jin [START_REF] Liu | A Novel Hybridization of Higher Order Finite Element and Boundary Integral Methods for Electromagnetic Scattering and Radiation Problems[END_REF]. However, the UWVF can be more efficient than finite elements at higher wave number, and this motivates our investigations here.

The main aim of coupling the UWVF and a Fast Multipole Method (FMM) is to reduce the overall numerical cost, and we have already introduced the coupled UWVF and FMM in several papers (see for example [START_REF] Coifman | The Fast Multipole Method for the Wave Equation: A Pedestrian Prescription[END_REF], [START_REF] Song | Multilevel Fast Multipole Algorithm for Solving Combined Field Integral Equations of Electromagnetic Scattering[END_REF], [START_REF] Darve | The Fast Multipole Method: Numerical Implementation[END_REF]). In particular in [START_REF] Darrigrand | Coupling of the Ultra-Weak Variational Formulation and an Integral Representation using a Fast Multipole Method in Electromagnetism[END_REF], we used a one-level FMM. The numerical results show that the one-level FMM leads to an algorithm with a cost comparable to the UWVF used in its original form. In [START_REF] Darrigrand | Ultra-Weak Variational Formulation and Efficient Integral Representation in Electromagnetism: a Thorough Study of the Algorithm Complexity[END_REF], we gave a thorough study of the complexity of that algorithm using either the one-level or multilevel fast multipole method in different configurations. This paper is devoted to numerical tests of the multilevel version of the algorithm. We shall show that the method can effectively solve several model problems, and give numerical evidence to support the complexity estimates.

The plan of this paper is as follows. The next section gives a brief presentation of the UWVF. In Section 3, we describe the use of the integral representation within the hybrid UWVF and integral representation (UWVF-IR) algorithm, and recall the results on the complexity that were developed in [START_REF] Darrigrand | Ultra-Weak Variational Formulation and Efficient Integral Representation in Electromagnetism: a Thorough Study of the Algorithm Complexity[END_REF]. The last section, which is the focus of this paper, presents numerical results obtained using a multilevel FMM (MLFMM) and demonstrates agreement with the study of complexity. These results show the improved efficiency of the novel algorithm compared to the original UWVF.

Ultra-Weak Variational Formulation

Suppose we wish to solve the time harmonic Maxwell system in the unbounded exterior of a bounded scatterer (assumed to have connected complement and polyhedral surface Γ int ). Introducing a polyhedral artificial boundary Γ ext containing the scatterer in its interior, we denote by Ω the bounded annular domain between Γ int and Γ ext (see Fig. 1). We denote by ν the outward normal to Ω.

To solve the time-harmonic Maxwell equations in a domain Ω we need to find the electric field E and magnetic field H such that the following equations hold:

∇ ∧ E -ıωµH = m , ∇ ∧ H + ıωεE = j , in Ω, (1)
where m and j are given data vector functions specifying the volume sources, ε and µ are positive piecewise constant functions of position and ω > 0 is the angular frequency of the field. Although not required by the method, we usually assume a source free region and select m = j = 0 (our numerical test will conform to this).

For the UWVF, it is convenient to specify the boundary condition on ∂Ω = Γ int ∪ Γ ext in the following non standard form ( [START_REF] Cessenat | Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation[END_REF])

Ω Γ int Γ ext Σ Figure 1:
The computational domain Ω is the region outside the boundary of the impenetrable scatterer Γ int and inside the artificial boundary Γ ext . A second artificial boundary Σ inside Γ ext and containing the scatterer (both impenetrable and penetrable components) in its interior is used in the integral representation of the field. If the exterior domain is entirely homogeneous and isotropic, we can take Σ = Γ ext .

where Q = 0 gives the standard low order absorbing boundary condition on Γ ext and g is computed from the incident wave. Since we model the total field, we choose g = 0 on Γ int and use Q to set the boundary condition. For example choosing Q = 1 gives the perfectly conducting boundary condition, while |Q| < 1 gives an impedance condition. The choice Q = -1 gives a magnetic wall condition that is useful for symmetric structures. We choose Q = 0 on Γ ext and g will be given by an integral operator.

The UWVF is based on the decomposition of the domain Ω into tetrahedra {Ω k } k=1,...,K and it computes the impedance trace of the solution on the boundaries of all these tetrahedra (which can then be post processed to give the solution in the entire element). This variational formulation is defined on the Hilbert space V = K k=1 L 2 t (∂Ω k ) where L 2 t (∂Ω k ) is the space of square integrable tangential fields on ∂Ω k the boundary of Ω k . For this space the scalar product is given by

(X , Y) V = k ∂Ω k X /∂Ω k Y /∂Ω k .
Under the assumption that ε and µ are positive constants on each element Ω k , (E, H) is found through the restriction of the field (E k , H k ) to ∂Ω k , where (E k , H k ) = (E, H) /Ω k . The method then solves for an unknown function X ∈ V , defined element by element by the impedance trace X /∂Ω k ∈ L 2 t (∂Ω k ) on ∂Ω k and

X /∂Ω k = ε /∂Ω k (E k ∧ ν k ) + µ /∂Ω k ((H k ∧ ν k ) ∧ ν k ) . (3)
where ε /∂Ω k and µ /∂Ω k are quantities defined by the values of ε and µ on each side of ∂Ω k (see [START_REF] Darrigrand | Coupling of the Ultra-Weak Variational Formulation and an Integral Representation using a Fast Multipole Method in Electromagnetism[END_REF] for details), and ν k is the exterior normal to ∂Ω k . The UWVF involves two operators Π and F defined in [START_REF] Cessenat | Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation[END_REF][START_REF] Cessenat | Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniquesproblèmes de Helmholtz 2D et de Maxwell 3D[END_REF]. The operator Π : V → V switches boundary traces across faces shared by two tetrahedra and involves the boundary condition (2) through the function Q on faces which contribute to the boundaries Γ ext and Γ int . More precisely if elements Ω j and Ω k meet at a face

f then ΠX j | f j,k = X k | f j,k
and on a boundary face of an element

Ω k , ΠX k | ∂Ω = QX k | ∂Ω .
The operator F : V → V is a local operator that links the outgoing and incoming impedance traces of the degrees of freedom on the interfaces between tetrahedra. If (u, v) satisfy the Maxwell system on Ω K and

X k = ε /∂Ω k (u ∧ ν k ) + µ /∂Ω k ((v ∧ ν k ) ∧ ν k ) . then F (X k ) = ε /∂Ω k (u ∧ ν k ) -µ /∂Ω k ((v ∧ ν k ) ∧ ν k ) .
The UWVF of Maxwell's equations is as follows [START_REF] Cessenat | Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation[END_REF][START_REF] Cessenat | Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniquesproblèmes de Helmholtz 2D et de Maxwell 3D[END_REF]:

Find X ∈ V such that (X , Y) V -(ΠX , F Y) V = ( b, Y) V for all Y ∈ V, (4)
for all Y ∈ V where g is the extension by zero of g to a function in V .

Thus by taking a finite dimensional subspace V h ⊂ V and using basis functions Z i , i ∈ J for V h , a Galerkin discretization of the formulation (4) leads to problem of finding

X h = i∈J X i Z i ∈ V h such that (X h , Y h ) V - (ΠX h , F Y h ) V = ( g, Y h ) V for all Y h ∈ V h .
Equivalently, in matrix/vector form, we seek to compute

X = [X 1 , • • • , X card(J) ] T such that (D -C)X = b , (5)
where D is the matrix with (i, j)th entry (Z j , Z i ) V and C has (i, j)th entry given by (ΠZ j , F Z i ) V . The data vector b is derived from the right hand side above in the same way.

As usual for the UWVF, to facilitate calculating the action of F , on each element Ω k we use a basis generated by taking the impedance trace of p k plane waves satisfying the adjoint Maxwell system on Ω k (p k /2 directions with two polarizations for each direction). In particular to discretize the problem, we follow [START_REF] Cessenat | Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation[END_REF] and use boundary functions given by

Y /∂Ω k = ε /∂Ω k (E ′ k ∧ ν k ) + µ /∂Ω k ((H ′ k ∧ ν k ) ∧ ν k )
where the fields (E ′ k , H ′ k ) are taken from the span of a set of p k plane waves that satisfy the adjoint Maxwell problem

∇ ∧ E ′ k -ıωµ Ω k H ′ k = 0 in Ω k , ∇ ∧ H ′ k + ıωε Ω k E ′ k = 0 in Ω k .
In ( 4), b ∈ V is derived from the right hand side of (1) and from g given in [START_REF] Cessenat | Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniquesproblèmes de Helmholtz 2D et de Maxwell 3D[END_REF]. At least six plane waves (and usually more) are used per element, counting polarizations.

The UWVF then leads to a sparse square system of size (

K k=1 p k ).
The number of plane waves p k is chosen depending on the local wavelength and diameter of the element (see [START_REF] Huttunen | Solving Maxwell's Equations Using the Ultra Weak Variational Formulation[END_REF]).

Compared to more classical volume methods, the UWVF enables one to reduce the number of elements in the mesh. The complexity of the method is then linked to the number of elements in the mesh and the number of basis functions per element. For concreteness, suppose the electromagnetic parameters of the domain are constant and define the wave number κ = ω √ εµ. Before a closer consideration of the complexity, we need to introduce another parameter: K 0 denotes the average number of tetrahedra taken in one dimension so that K ∼ K 3 0 . As a volume method, the UWVF method leads to a sparse system: The number of degrees of freedom is of order K 3 0 p and the complexity of the algorithm is O(K 3 0 p 2 ) where p denotes the average number of basis functions per tetrahedra which typically satisfies K 0 p ∼ κ [START_REF] Huttunen | Solving Maxwell's Equations Using the Ultra Weak Variational Formulation[END_REF].

Use of an Integral Representation within the UWVF

To simplify the presentation we shall suppose that ε = µ = 1 so that the exterior domain is entirely homogeneous and we use the perfectly conducting boundary condition on Γ int (i.e. the scatterer is not penetrable and the exterior medium is homogeneous). In this case we may take Σ = Γ int (extensions to more general domains in which Σ = Γ int is then easy to understand). In this case the artificial boundary Γ ext can then be taken very close to the boundary of the obstacle. Following Hazard and Lenoir [START_REF] Hazard | On the Solution of Time-Harmonic Scattering Problems for Maxwell's Equations[END_REF], the hybrid scheme consists in replacing the low order absorbing boundary condition

-E ∧ ν + (H ∧ ν) ∧ ν = -E 0 ∧ ν + (H 0 ∧ ν) ∧ ν on Γ ext by the boundary condition -E ∧ ν + (H ∧ ν) ∧ ν = -E s ∧ ν +(H s ∧ ν) ∧ ν -E 0 ∧ ν +(H 0 ∧ ν) ∧ ν ,
where (E s , H s ) are given by the Stratton-Chu formula ( [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]) in terms of ν × H and ν × E on Σ (i.e. under our assumptions on Γ int ) via

E s (x) = ∇ x ∧ Γ int G(x, y) ν Σ (y) ∧ E(y) dγ(y) - 1 ı ω ∇ x ∧ ∇ x ∧ Γ int G(x, y) ν Σ (y) ∧ H(y) dγ(y) , (6)
H s (x) = ∇ x ∧ Γ int G(x, y) ν Σ (y) ∧ H(y) dγ(y) + 1 ı ω ∇ x ∧ ∇ x ∧ Γ int G(x, y) ν Σ (y) ∧ E(y) dγ(y) , (7)
where ν Σ is the exterior normal to the surface Σ = Γ int and G(x, y) = exp(ıκ|x -y|)/(4π|x -y|) is the fundamental solution for the Helmholtz equation. Thanks to the structure of the unknowns of the UWVF, as shown in [START_REF] Darrigrand | Coupling of the Ultra-Weak Variational Formulation and an Integral Representation using a Fast Multipole Method in Electromagnetism[END_REF], the fields in the integrands above can be computed directly from the degrees of freedom of the UWVF (3) taking into account the convention for the direction of normals, and the boundary condition on Γ int [START_REF] Cessenat | Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniquesproblèmes de Helmholtz 2D et de Maxwell 3D[END_REF] where g = 0.

The system (5) becomes (D -C -C)X = b where C couples the degrees of freedom on Γ int (more generally Σ) and Γ ext . The matrix C can be split into different discrete integral operators C i , i = 1, ..., 4 of the form

( C i X h ) kl = Σ ext kk c k S i (X h ) • F Y kl dγ ext , where • Σ ext
kk is the face on Γ ext of a tetrahedron intersects the exterior boundary, • c k depends only on ε and µ on Σ ext kk , • F is the local operator introduced in (4),

• S i is a global operator which comes from the right hand side of ( 6)-( 7), for instance

(S 1 (X ))(x) = - Γ int f Q (y) ∇ y G(x, y) ∧ X (y) dγ(y) ∧ ν(x) ,
where f Q is a function involving Q and ε. These integral operators can be evaluated by the FMM.

In this paper, the solution of the new system (D -C -C)X = b is obtained by the same method (BiCGStab) as used for the classical UWVF system (D -C)X = b, considering C + C as a small perturbation of C. It is important to note that we never compute C explicitly, but evaluate its action as needed by the FMM.

The integral representation aims to reduce the distance of the absorbing boundary from the scatterer to a number of elements independent of κ. We then have a number of elements in the mesh of order K 2 0 . This reduces the complexity related to the volume calculation. The FMM is used to control the cost of the integral calculation related to the integral operators which give rise to large dense blocks in the discrete system.

A rigorous expression of the algorithm complexity involves several parameters: the size of the mesh K ∼ K 3 0 or K 2 0 , the average number of basis functions per tetrahedron p, the wave number κ, with the correlation K 0 p ∼ κ. Moreover, the integral representation involves oscillating functions such as the Green function, and the exterior normal to the boundaries. Evaluation of these operators requires an accurate model of the boundary which is assured by our assumption of a polyhedral domain. For more general surfaces it a possible modification of our algorithm would be to use a coarse volume mesh for the UWVF and a finer surface mesh for the integral calculations (the mesh used for the resolution of classical integral equations); this follows the double-mesh concept introduced by Zhou et al. [START_REF] Abboud | Improvement of the Integral Equation Method for High Frequency Problems[END_REF] and used in [START_REF] Darrigrand | Coupling of Fast Multipole Method and Microlocal Discretization for the 3-D Helmholtz Equation[END_REF]. Furthermore, note that the FMM used here does not involve very close interactions. Indeed, the integral operators link points on one boundary to points on another boundary so the source points are not in the same domain as the target points. All these considerations lead to different algorithms with different complexities. In Tab. 1, we give the results developed in [START_REF] Darrigrand | Ultra-Weak Variational Formulation and Efficient Integral Representation in Electromagnetism: a Thorough Study of the Algorithm Complexity[END_REF] which suggest that the UWVF with MLFMM and no close interactions should have attractive work estimates.

Numerical Results

In this section, we give numerical results with the new algorithm which couples the UWVF and the integral representation using a MLFMM. Our first results are for the problem of approximating scattering by a perfectly conducting unit sphere for Γ int ) with the wave numbers κ = 3m -1 , 4m -1 , 6m -1 , 10m -1 , 15m -1 . The wavelength is then λ = 2π/3 ≈ 2.09m, π/2 ≈ 1.6m, π/3 ≈ 1.05m, π/5 ≈ 0.63m, 2π/15 ≈ 0.42m. This very simple example has the advantage that a Mie series solution is available for comparison.

We also note that the method easily extends to disconnected scatterers. If the scatterers are close together (closer than approximately 4 elements in the mesh), we need to surround the scatterers by a single auxilary boundary. In this case multiple scattering is taken into account by the UWVF. If the scatterers are sufficiently separated, we can use a thin layer of elements around each scatterer. In this case multiple scattering is taken into account via the fast multipole method and hence dispersion error is reduced. We have not yet found it necessary to precondition in a special way in these cases.

K 3 0 ∼ (κ/p) 3 K 3 0 p 2 ∼ K 0 κ 2 UWVF + IR K 2 0 ∼ κ 2 K 2 0 p 2 + K 4 0 p 2 ∼ κ 2 + κ 4 UWVF + IR + 1-level FMM K 2 0 ∼ κ 2 K 2 0 p 2 + K 3 0 p ∼ κ 2 + κ 3 UWVF + IR + multilevel FMM K 2 0 ∼ κ 2 K 2 0 p 2 + K 2 0 ln 2 (K 0 ) p ∼ κ 2 + κ 2 ln 2 κ UWVF + IR + multilevel FMM with double mesh K 2 0 ∼ (κ/p) 2 K 2 0 p 2 + κ 2 ln 2 (κ) p ∼ κ 2 + κ 2 ln 2 (κ)p UWVF + IR + 1-level FMM without close interaction K 5/2 0 ∼ κ 5/2 K 5/2 0 p 2 + K 3 0 p ∼ κ 2 + κ 3 UWVF + IR + multilevel FMM without close interaction K 2 0 ∼ κ 2 K 2 0 p 2 + K 2 0 ln 2 (K 0 ) p ∼ κ 2 + κ 2 ln 2 κ
Results from three codes are presented: the classical UWVF with a Silver-Müller type low order auxilliary boundary condition of order 0 (i.e the standard UWVF with the boundary condition described earlier setting Q = 0), the code UWVF+FMM using a 1-level FMM and the code UWVF+MLFMM using a multi level FMM. All the resuts in this paper were obtained on an Apple Mac Pro with 2 × 3GHz quad-core Xeon processors using 16 Gb RAM.

The unit sphere

For approximating scattering by the unit sphere, the exterior boundary Γ ext is taken to be a concentric sphere. Obviously these are not polyhedral domains, but can be approximated using a sufficiently fine mesh near the boundary. We have experimented with several exterior boundaries giving rise to different meshes as defined in [START_REF] Darrigrand | Coupling of the Ultra-Weak Variational Formulation and an Integral Representation using a Fast Multipole Method in Electromagnetism[END_REF]. Table 2 describes those used in this paper. The names "Sxxx" denote the different meshes, where "xxx" denotes one of the numbers 400, 200, 100, 075, 050, 025, 040, 017. This number gives the distance between Γ int and Γ ext in centimeters. For κ = 4 the distance from the perfect conductor to the artificial boundary ranges from 0.16λ (S025) to 2.5λ (S400) where λ = 2π/κ is the wavelength.

All the meshes have been generated using COMSOL Multiphysics version 3.5a. The meshes S400 and S200 are appropriate for a classical use of the UWVF. The meshes S025, S040, S017, S010 and S007 have been generated optimizing the ratio between the average edge-length h and the wavelength λ to h ≈ λ/5. These meshes are quite uniform (the S400 and S200 meshes are graded to give larger elements away from the scatterer). The large number of tetrahedra in S025 might appear to be a disadvantage for the UWVF+MLFMM method (the number of tetrahedra is comparable to S400). However the number of Degrees of Freedom (DoF) is much less than for S400 because fewer plane waves are used per element due to the smaller size of elements as shown in the bottom row of Table 2.

We use the radar cross section (RCS) as a measure of the efficiency of the method since the RCS is often We start by comparing the classical UWVF, the UWVF-FMM and UWVF-MLFMM at constant wave number κ = 4. Figure 2 compares the results of our codes as a function of angle with the exact Mie series solution for TE polarization, and Fig. 3 shows results for TM polarization, but only for UWVF-MLFMM. As is well known, the classical UWVF code requires a large distance between Γ int and Γ ext to obtain an accurate solution due to its low order absorbing boundary condition. Figure 2 clearly shows the impact of the integral representation on the accuracy of the results for the case κ = 4. The code UWVF+FMM gives results which fit with the Mie series solution, even with the thin mesh S025. The results obtained with the code UWVF+MLFMM are comparable to the those of the code UWVF+FMM, indicating that use of the multilevel scheme does not degrade the FMM accuracy. CPU-time, memory requirements and error measures are given in the case of the TE polarization in Table 3. The results are given for the meshes S400 and S200 using the classical UWVF code and for the mesh S025 using the codes UWVF+FMM and UWVF+MLFMM. In the table we use the following notation (units are seconds and Giga-bytes): • Nit = Number of iterations for the bi-conjuguate gradient.

• mem = Memory required by the code.

• Rel Err 2 = Relative quadratic error on RCS. Table 3: Computational costs comparisons between the UWVF (S400; S200), UWVF+FMM (S025 -F) and UWVF+MLFMM ( S025 -M) for scattering by a unit sphere at constant κ. For more details of the meshes see • Rel Err ∞ = Relative infinity error on RCS.

• RMS error = Root mean square error on the scattering amplitude. Like in [START_REF] Cui | Efficient MLFMA, RPFMA, and FAFFA algorithms for EM scattering by very large structures[END_REF],

RMS error = 1 n n i=1 (RCS(i) -Exact-RCS(i)) 2
is used as a classical RCS error ( [START_REF] Cui | Efficient MLFMA, RPFMA, and FAFFA algorithms for EM scattering by very large structures[END_REF]).

The RCS curve obtained with S200 using the UWVF code could be acceptable, however its accuracy is quite poor in comparison with the one obtained with S025 using the codes UWVF+FMM and UWVF+MLFMM hence we should compare results for S400 and S025 in Table 3. That case shows that the UWVF-MLFMM code reduces computer time by approximately 50% and the memory needed by almost 75%. for roughly the same RMS error.

Next we investigate the κ dependence of the method to determine how well Table 1 reflects the practical work needed by the method. Of course, as κ increases the method needs to add degrees of freedom to maintain accuracy. Moreover our crude geometry approximation requires a finer mesh as κ increases. Hence we use meshes S040, S035, S017, S010 and S007. Figures 3, 4 and5 show the angular dependence of the RCS from our UWVF-MLFMM code compared to the exact Mie series.

Memory, CPU and time requirements, together with quantitative error results are shown in Table 4. We have tried to keep the RMS error roughly constant to enable comparison. We then plot the resulting number of degrees of freedom, CPU and memory requirements against κ in Fig. 8. We also consider reference lines to enable comparison with results predicted in Table 4. These curves support the theoretical results on the complexity.

Multiple scatterers

In this subsection, we show results obtained on two different multiple-scatterers examples with the wave number κ = 4. The first result (Figure 10) is obtained with UWVF+MLFMM on a scatterer consisting of two closely separated spheres using a mesh named 2S025. This scatterer consists of two unit spheres centered at (±1.125, 0, 0), and compared to UWVF with 2S200 and 2S400 where the name of the meshes were chosen as before. For the different cases, the artificial boundary is taken around the two spheres such that the domain is as in the two leftmost panels of Figure 9. Because the scatterers are close together in comparison to our meshing constraints (i.e. closer than 4 elements apart), we need to use a single mesh surrounding both scatterers. For the meshes 2S200 and 2S400 that surround two balls, the error is given compared to the results obtained with UWVF+MLFMM on 2S025. Indeed, in the case of scattering from two spheres, we consider our solution on 2S025 as the reference solution. The second example (results in Figure 11) is a set of four unit spheres centered at (±1.8, 0, ±1.8), and the mesh is called 4S025. In this case, the spheres are well separated and so we can take the artificial boundary to consist of four spheres of radius 1.5 m and centered at the same points as the spheres which contribute to the scatterer. In this situation, the domain is a four-components domain, and multiple scattering is taken into account by the integral equation. Figure 9 shows the domain. No error estimate is available for the case of scattering by four balls (mesh 4S025).

Conclusion

We have described how to couple the UWVF and multi-level fast multipole scheme via an overlapping strategy. This allows us to use non-convex artificial boundaries close to the scatterer, and hence decrease the memory and CPU requirements for solving an exterior scattering problem. The main content of the paper is a series of numerical tests that 1) demonstrate that the new method improves on the original UWVF when applied to simple scatterers, 2) supports the predicted wave number dependence of the work estimates and 3) shows how the method can handle multiple scatterers.

Clearly much more work is needed to make this method generally applicable. Error estimates should be proved and a better preconditioner for the UWVF derived. In addition a subgrid boundary model needs to be implemented to allow for a coarser volume grid near the boundary, and to allow a better approximation of curved boundaries.

Figure 2 :

 2 Figure2: The TE-polarized RCS as a function of polar coordinate θ computed using the three codes in the study, with κ = 4. Left: results for the classical UWVF code with meshes S400 and S200. Right: results for the codes UWVF+FMM and UWVF+MLFMM with the mesh S025.

•

  CPU time: T00 = precalculation of the matrices C and D. T0c (resp. T0f) = precalculation of the close (resp. far) interactions related to C. TD (resp. TC) = one multiplication by D -1 (resp. C) (average). TCc (resp. TCf) = one multiplication by C close (resp. far) (average). Tcg = Total CPU time for the solution of the system. Ttot = Total CPU time required by the code.

Figure 3 :

 3 Figure3: The TM-polarized RCS against polar angle θ computed using the code UWVF+MLFMM (dashed line) with the mesh S025 and the exact solution (solid line). Here κ = 4.

Figure 4 :

 4 Figure 4: The RCS computed using UWVF+MLFMM with the mesh S040, with κ = 3. Left: TE-polarized RCS. Right: TM-polarized RCS.

Figure 5 :Figure 6 :

 56 Figure 5: The RCS computed using UWVF+MLFMM with the mesh S017, with κ = 6. Left: TE-polarized RCS. Right: TM-polarized RCS.

Figure 7 :

 7 Figure 7: The RCS computed using UWVF+MLFMM with the mesh S007, with κ = 15. Left: TE-polarized RCS. Right: TM-polarized RCS.

Figure 8 :

 8 Figure 8: Number of degrees of freedom (left), CPU-time (middle) and memory requirements (right) with respect to the wave number κ, using a log-log scale. We show the data from Table 4 marked with a symbol in each case. The best fit line is shown as a solid line, and lines corresponding to O(κ), O(κ 2 ) and O(κ 3 ) for comparison.

Figure 9 :

 9 Figure 9: Meshes 2S025, 2S200 and 4S025.

Figure 10 :

 10 Figure 10: The RCS computed using UWVF+MLFMM with the mesh 2S025 and using UWVF with the meshes 2S200 and 2S400, with κ = 4. Left: TE-polarized RCS. Right: TM-polarized RCS.

Figure 11 :

 11 Figure 11: The RCS computed using UWVF+MLFMM with the mesh 4S025, with κ = 4. Left: TE-polarized RCS. Right: TM-polarized RCS.

Table 1 :

 1 Complexity estimates for the various versions of the UWVF + IR

	Method	Number of elements	Cost of the solution
	UWVF

Table 2 :

 2 The meshes for the unit sphere.

	Name	S400	S200	S040	S025	S017	S010	S007
	κ	4	4	3	4	6	10	15
	Radius of Γ ext in m	5	3	1.40	1.25	1.17	1.10	1.07
	Distance between Γ int and Γ ext Number of tetrahedra	≈ 2.5λ 16179	≈ 1.3λ ≈ λ/7.5 ≈ λ/6.3 ≈ λ/6.4 ≈ λ/6.3 ≈ λ/6.6 14526 5822 11008 8449 22630 44459
	Number of basis functions
	per tetrahedron	8 to 128 8 to 72 10 to 22 10 to 24 16 to 30 16 to 28 14 to 30
	Number of DoF	880200 508450 101814	178146	183108	480226	989238

Table 2 .

 2 The RCS error is computed via the Mie series.

	Case	κ T00	T0c	T0f	TD	TC TCc TCf Tcg Nit
	S400	4	326	-	-	1.2	7.6	-	-	1399 156
	S200	4 73.7	-	-	0.4 1.99	-	-	388 155
	S025 -F	4	34	138.6 4.3	0.1	0.2	1	3.5	656 135
	S025 -M 4 10.4	52	2.98 0.11 0.3 0.36 4.5	733 139
	Case S400	κ Ttot mem Rel Err 2 Rel Err ∞ RMS error 4 1725 6.2 1.5 10 -2 1.1 10 -2 9.5 10 -2
	S200	4	462	1.8	3.2 10 -2	2.5 10 -2		2 10 -1
	S025 -F	4	838	2.8	4.6 10 -3	5.1 10 -3		2.8 10 -2
	S025 -M 4	798	1.6	1.7 10 -2	1.1 10 -2		1.1 10 -1

Table 4 :

 4 Computational costs comparison for UWVF+MLFMM (S040 ; S025 -M ; S017 ; S010 ; S007) as κ changes. For more details of the meshes see Table2. All results are for the unit sphere and the RCS error is computed via the Mie series.

	Case		κ	T00 T0c T0f	TD	TC TCc TCf Tcg Nit
	S040		3	6.81	80			2.7	0.06 0.18 0.5		3		469 125
	S025 -M	4	10.4	52		2.98 0.11 0.3 0.36 4.5		733 139
	S017		6	13.5 141		5.7	0.1	0.4 0.95 7.9		894	95
	S010		10		35	467		15	0.3	1	3			19.6 2252 93
	S007		15		76	921		66	1.3	4.5 16.7 61.5 8880 105
	Case S040		κ 3	Ttot mem Rel Err 2 Rel Err ∞ RMS error 558 1.5 1.2 10 -2 1.2 10 -2 7.3 10 -2
	S025 -M	4		798	1.6		1.7 10 -2	1.1 10 -2			1.1 10 -1
	S017		6	1054 2.99		6.4 10 -3	4.5 10 -3			4.1 10 -2
	S010		10 2770 8.38		5.5 10 -3	4.1 10 -3			3.6 10 -2
	S007		15 9944	17		1.3 10 -2	1.1 10 -2			8.5 10 -2
	10 6	Number of DoF with UWVF+IR+MLFMM ; Spheres κ = 3, 4, 6, 10, 15			10 4	CPU requirements with UWVF+IR+MLFMM ; Spheres κ = 3, 4, 6, 10, 15				Memory requirements with UWVF+IR+MLFMM ; Spheres κ = 3, 4, 6, 10, 15
	Number of DoF								CPU										Memory
	slope(nb DoF)								slope(CPU)									slope(MEM)
	κ									κ										κ
	κ 2									κ 2										κ 2
	κ 3									κ 3									10 1	κ 3
	Number of DoF								CPU (s)									Memory (Gb)
										10 3
	10 5	4	6	κ	8	10	12	14	16		4	6	κ	8	10	12	14	16	0 10	4	6	κ	8	10	12	14	16

Table 5 :

 5 Computational costs comparison for the basic UWVF (2S400 ; 2S200) and the UWVF+MLFMM (2S025; 4S025) for scattering from multiple objects.

	Case	T00 T0c	T0f	TD	TC TCc TCf	Tcg Nit
	2S400	250	-	-	1.12 7.4	-	-	1172 136
	2S200	93.5	-	-	0.5 2.68	-	-	424 130
	2S025 17.29 187	5.73	0.17 0.52 1.28 7.25 1087 118
	4S025	28	361 12.25 0.25 0.79 2.45 16.03 1675 86
	Case 2S400 1422 6.12 Ttot mem Rel Err 2 Rel Err ∞ RMS error 4.3 10 -2 5.7 10 -2 3.6 10 -1
	2S200	518	2.28	7.8 10 -2	9.9 10 -2		6.4 10 -1
	2S025 1297 3.34	-		-		-
	4S025 2077 6.73	-		-		-

-| √ ε | E ∧ ν + (| √ µ | H ∧ ν) ∧ ν = Q(| √ ε | E ∧ ν + (| √ µ | H ∧ ν) ∧ ν) + g ,(2)

The research of PM is partially supported by grant number FA-9550-08-1-0138 from the US AFOSR. This paper is dedicated to Professor George Hsiao on the occasion of his 65 th birthday. Professor Hsiao has inspired us by his love for integral equations.


Dr Laurence Galet
email: laurence.galet@mines-albi.fr

Dr Yamina Ouabbas

Alain Chamayou

Dr Philippe Grosseau

Michel Baron

Prof. Dr Gérard Thomas

Dr A Chamayou

Prof M Baron

Surface analysis of silica gel particles after mechanical dry coating with magnesium stearate

Keywords: AFM phase imaging, adhesion force, dry coating, silica, magnesium stearate

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Dry particle coating to change the surface properties of powders is very important in many industries. Typical applications include, but are not limited, to flowability, wettability (hydrophobic/hydrophilic properties), solubility, dispersibility, flavour, shape, electrostatic, optical, electric, magnetic, etc. In dry particle coating processes, materials with relatively large particle size (host particles; 1-500 µm) are mechanically coated with fine particles (guest particles; 0.1-50 µm) in order to create new functionality or to improve their initial characteristics [START_REF] Yoshihara | Hybridization -Technology for surface modification of powders without binders[END_REF]. Since the size of the guest particles are so small, van der Waals interactions are strong enough to keep them firmly attached to the host particles. Thus, either a discrete or continuous coating of guest particles can be achieved depending on a variety of operating conditions including processing time, rotation speed, weight fraction of guest to host particles and particle properties [START_REF] Pfeffer | Synthesis of engineered particulates with tailored properties using dry particle coating[END_REF]. Figure 1 below is a simple schematic illustrating the process of dry particle coating.

Over the last few years, we have performed experimental investigations of applications of the dry coating technique to study the effect of mechanical dry coating on the surface properties of powders. Our recent work is related to the modification of the flowability and the wettability properties of silica gel particles coated with different weight ratios of magnesium stearate by using a high energy impact coater Hybridizer from Nara Machinery [START_REF] Ouabbas | Surface modification of silica particles by dry-coating: Characterization and powder ageing[END_REF] and a high shear mixer Cyclomix from Hosokawa [START_REF] Ouabbas | Particle-particle coating in a cyclomix impact mixer[END_REF]. For example, the flowability of the silica gel powder was significantly decreased after treatment in the Cyclomix mixer with 15% of MgSt 2 . Moreover, it has been found that the coating by hydrophobic MgSt 2 in the Cyclomix can reduce the high affinity between silica gel and water after treatment with 5% and 15% of MgSt. In those papers, it has been demonstrated that a dry particle coating technique can be used to modify the properties of silica gel powder by coating with small quantities of hydrophobic magnesium stearate in both the Hybridizer and the Cyclomix. The more uniform coating has been obtained after treatment in the Hybridizer device. In this paper, we report an investigation of the surface analysis of the coated silica gel particles using atomic force microscopy in tapping mode (TM-AFM) and contact mode (CM-AFM). The TM-AFM method, based on the measurements of the oscillations of a cantilever probe in contact on the sample, has been used for several years to characterize the surface topography of particles. In particular, phase contrast images are related to characterize the attractive and repulsive probe sample interactions [START_REF] Chen | Phase contrast and attraction-repulsion transition in tapping mode atomic force microscopy[END_REF], the sample surface components [6] and the morphology of composite particles [START_REF] Gun | Composite powders with titania grafted onto modified fumed silica[END_REF]. The CM-AFM method consists in using a colloidal probe fixed on the cantilever to quantify the adhesion force between the particles. As an example, Kani et al report an interesting study on the silicamica adhesion force as a function on the presence of impurities on the silica particles [START_REF] Kani | Influence on the surfaceadhered nanoparticles and nanoporous structure on particle-particle interaction of silica[END_REF].

Meincken et al. report a relationship between contact angle analysis and adhesive forces measurements by the AFM technique in a study of the surface hydrophobicity of polyurethane coatings [START_REF] Meincken | Influence of the viscosity and the substrate on the surface hydrophobicity of polyurethane coatings[END_REF]. This technique has been also used to characterize the adhesion properties of pharmaceutical carriers [START_REF] Louey | Characterization of adhesional properties of lactose carriers using atomic force microscopy[END_REF]. In our case, we performed AFM phase-contrast and adhesion forces measurements to characterize the surface of silica particles as a function of the magnesium stearate mass ratio (1 to 30%).

Experimental

Powders

Silica gel powder supplied by Merck and usually used for filling chromatography columns has been chosen as host particles for dry coating (Silica gel 60, 0.040-0.063 mm Table1: Properties of host and guest particles.

Coating process

The hybridization system (HB, type NHS-0; Nara Machinery Co.), a high-speed dry impact blending coater, has been used to make the silica gel-MgSt 2 composite particles. The reactor consists in a very high-speed rotating rotor with six blades, a stator and a powder recirculation circuit. The coating chamber is surrounded with a jacket in which coolant is circulated [START_REF] Yoshihara | Hybridization -Technology for surface modification of powders without binders[END_REF]. The coating process can be summarized as follows: the powder mixture (host and guest particles) is subjected to high impaction and dispersion due to collisions with blades and the walls of the device and continuously re-circulates in the machine through the cycle tube. Particle coating is achieved due to the embedding or filming of the guest particles onto the host particles by high impaction forces and friction heat. Since the rotor of the hybridizer can rotate from 5000 to 16000 rpm, very short processing time is required to achieve coating. The operating conditions used in experiments are 4800 rpm for 5 min. The preparation of the coated particles is described in detail in previous papers [START_REF] Ouabbas | Surface modification of silica particles by dry-coating: Characterization and powder ageing[END_REF][START_REF] Ouabbas | Particle-particle coating in a cyclomix impact mixer[END_REF]. Coating experiments have been carried out with 1%, 5%, 15% and 30% of mass fraction of guest MgSt 2 particles.

Particles characterization

The uncoated and coated silica gel particles were examined by environmental scanning electron microscopy (ESEM) to study the surface morphology and the covering by the MgSt 2 particles (XL 30 Philips). The ESEM images reveal that greater MgSt 2 coverage is observed on the surface of silica gel particles as the percentage of MgSt 2 is increased (Figure 2). Additional investigations were done to analyse the surface atomic composition of the coated particles by X-ray spectroscopy confirming the presence of magnesium atoms when analysis is performed on a black spot [START_REF] Ouabbas | Procédés mécaniques d'élaboration à sec de particules composites à propriétés d'usage contrölées -Caractérisation et stabilité d'un gel de silice[END_REF]. At 15%, the surface coverage of the MgSt 2 is intense but seems discrete. TG-DSC experiments have been performed on silica gel-MgSt 2 mixtures made by manual mixing (1%, 3%, 15%, 30% and 50% mass ratios). The loss of mass obtained between 250°C and 600°C as a function of the MgSt 2 ratio is reported in the Figure 5. This result is the mean of two measurements. The results show a very good relation between the AGT-DSC mass loss (%) as a function of the MgSt 2 ratio (%). The linear relation obtained with a high correlation coefficient has been used to calculate the real composition of the coated particles after processing in the reactor with a precision of less than ±1%. AFM measurements were made with a Multimode Nanoscope IIIA from Veeco Metrology Group. In tapping mode, the cantilever is driven to oscillate. Height, amplitude and phase images were simultaneously recorded for each sample. Beyond the standard height imaging, the phase contrast mode is quite efficient to exhibit the main surface property differences between de silica gel and MgSt 2 .

The amplitude and the frequency of the oscillations change when the tip scans on a surface.

TM-AFM images are produced by imaging the oscillating contacts of the tip with the sample surface: height and amplitude images give information on the sample surface topography and phase images give information on the atomic surface composition (interaction tip-surface). The TM-AFM measurements were performed with a MPP 11100

Silcon phosphorus doped silicium tip with a spring constant of 40 N/m. Large surface scale of surface sample has been scanned: from 250x250 nm 2 to 5x5 mm 2 , for silica gel and silica gel-MgSt 2 coated particles (1%, 5% and 15% MgSt 2 mass ratio).

In contact mode AFM, the tip scans the sample in close contact with the surface. The force between the tip and the surface sample is measured as the distance, by maintaining a constant deflection. To examine the mechanism of interaction between silica gel and MgSt 2 particles, the CM-AFM measurements were performed with a NP silicon nitride tip with a spring constant of 0.32 N/m on which a MgSt 2 particle was glued under optical microscope.

The adhesion forces were measured between the MgSt 2 particles attached on the end of the cantilever and the different samples of uncoated and coated particles with 1% to 30% MgSt 2 mass ratio. For each sample 450 force curves were obtained. The presence of the magnesium stearate particle on the tip was checked after each series of measurements.

Influence of the coating on the surface property of the silica gel particles was also studied.

The sessile drop method was used to study the wettability of the different samples. We report the contact angle value of a water drop (10 µl) deposited on different particulate systems: pure silica gel, pure MgSt 2 and coated particles. The method performed is described in details in previous papers [START_REF] Ouabbas | Surface modification of silica particles by dry-coating: Characterization and powder ageing[END_REF][START_REF] Ouabbas | Particle-particle coating in a cyclomix impact mixer[END_REF].

Results and Discussion

The real concentration of the MgSt 2 after coating is calculated from the TG-DSC measurements and the calibration straight line shows in Figure 5. The results are reported in the Figure 6. This loss of MgSt 2 , due to a deposit of MgSt 2 on the coating chamber surface, varies as a function of the introduced MgSt 2 mass ratio. That is from 30% for the sample containing initially 30% of MgSt 2 in the mixture, to 40% for the sample containing initially 5% MgSt 2 in the mixture. The phase imaging TM-AFM characterisations of the samples are reported in Figure 7.

Concerning the MgSt 2 , the amplitude image reveals an irregular surface topography. The phase contrast image shows various levels of phase contrast (black and white areas), which can reach 50°. The silica gel particles have a more regular surface topography, formed by spherical micro beads (about 100 to 600 nm) and a low phase contrast only due to the topography (from 0° to 10°). The analyses of the silica-MgSt 2 coated particles clearly reveal a modification of the surface due to the presence of the magnesium stearate, even at 1% of MgSt 2 introduced. With a MgSt 2 initial mass ratio of 15%, both the amplitude and the phase contrast images evolve into the pure MgSt 2 images. For the coated particles, the part of the image characterized by high phase angles can be attributed to the presence of MgSt 2 layers on the surface of silica gel. The contact mode analyses bring some pertinent information about the silica gel covering by the magnesium stearate. Figure 8 gives an adhesion force-displacement curve for the pure products. The force measured between a particle of MgSt 2 fixed onto the AFM tip is about 10 nN for the pure silica gel particle. That is near 120 nN for the MgSt 2 -MgSt 2 , revealing a high interaction. The adhesion force distributions for the different particulate systems are reported in the Figure 9. These distributions were obtained from 450 measurements of each particulate system. This force varies proportionally to the quantity of MgSt 2 introduced and confirms the increasing presence of the MgSt 2 on the silica gel surface when the mass ratio increases.

For the coated particles, the number of high adhesion forces increased as the mass ratio of MgSt 2 increases. The mean adhesion force between the MgSt 2 particles attached on the end of the cantilever and the MgSt 2 sample (68 nN) was stronger than for the silica gel sample ( 8nN). The large distribution of the adhesion forces confirm the discrete deposition of the MgSt 2 particles onto the surface silica gel particles as observed by the ESEM images.

The mean adhesion force calculated from the 450 measurements for each system is reported in the Table 2. The results show an increasing of the mean value of the force as a function of the increasing of the MgSt 2 concentration. The mean adhesion force seems to evolve progressively to the mean value of the pure MgSt 2 (Figure 10). This no linear evolution could be due to a random distribution of the MgSt 2 particles on the silica gel surface. The affinity of water is evaluated by the behaviour of a water drop deposited on the particulate systems. The initial contact angle values are reported in Table 3. The Figure 11 shows the evolution of the contact angle as a function of the MgSt 2 ratio. We observe a increasing of the contact angle of the water drop as the presence of the MgSt 2 on the silica gel surface. This result is of course in agreement with the hydrophobicity of this compound.

However, even with a low content of MgSt 2 , the contact angle measured is significatively modified: 70° and 110° for 0.1% and 1.1% MgSt 2 respectively. The surface property of the silica gel particles is modified by dry coating even with a very low ratio of MgSt 2 . This result reveal that at the water drop scale, a very few presence of MgSt 2 is enough to modify the hydrophilic property of the silica gel.

Conclusion

Dry particle coating can be used to create new generation materials with altered properties or new functionalities. This treatment creates a discrete or a continuous coating of fine particles on the surface of the large particle surface. In this study a Nara Hybridizer has been used as a dry coating device for the coating of silica gel particles with different mass ratios

Legends of figures

Figure 1 .

 1 Figure 1. Schematic principle of dry coating.

Figure 2 .

 2 Figure 2. ESEM images of silica gel and MgSt 2 -silica gel coated particles.

Figure 3 .

 3 Figure 3. TG-DSC diagrams for silica gel.

Figure 4 .

 4 Figure 4. TG-DSC diagrams for MgSt 2 .

Figure 5 :

 5 Figure 5: Experimental relation between the TG-DSC mass loss (from 250°C to 600°C) and

Figure 6 :

 6 Figure 6: Real mass fraction of MgSt 2 after the coating process calculates from the mass

Figure 7 .

 7 Figure 7. TM-AFM images of MgSt 2 , silica gel and silica coated with 1% and 15 % initial mass ratios of MgSt 2 (2x2 µm 2 ).

Figure 8 :

 8 Figure 8: Force curves F (N/m) as a function of the distance z (µm) measured with MgSt 2

Figure 9 :

 9 Figure 9: CM-AFM adhesion force distribution obtained with 450 analyses for different

Figure 10 :

 10 Figure 10: Mean adhesion force (nN) as a function of the sample MgSt 2 concentration.

Figure 11 :

 11 Figure 11: Contact angle measurements (°) as a function of the sample MgSt 2 concentration.

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig. 1: Schematic principle of dry coating.

Fig. 4 :Fig. 5 :Fig. 6 :

 456 Fig. 4: TG-DSC diagrams for MgSt 2 .

Fig. 7 :

 7 Fig. 7: TM-AFM images of MgSt 2 , silica gel and silica coated with 1% and 15 % initial mass ratios of MgSt 2 (2x2 µm 2 ).

Fig. 8 :Fig. 9 :Fig. 10 :

 8910 Fig. 8: Force curves F (N/m) as a function of the distance z (µm) measured with MgSt 2

Figure 11 :

 11 Figure 11: Contact angle measurements (°) as a function of the sample MgSt 2 concentration.

Fig. 1 :

 1 Fig. 1: Schematic principle of dry coating.

Fig. 2 :Fig. 3 :

 23 Fig. 2: ESEM images of silica gel and MgSt 2 -silica gel coated particles.

Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 :

 4567 Fig. 4: TG-DSC diagrams for MgSt 2 .

Fig. 8 :Fig. 9 :

 89 Fig. 8: Force curves F (N/m) as a function of the distance z (µm) measured with MgSt 2

Table 2 :

 2 Mean adhesion force (nN) for different particulate systems.

Table 3 :

 3 Contact angle measurements (°) for different particulate systems.

Table 2 :

 2 of MgSt 2 particles. ESEM observations show that MgSt 2 particles form a discrete coating onto the surface silica gel particles. The MgSt 2 surface coverage increases and the surface becomes more and more uniform as the mass ratio of MgSt 2 is increased. In addition to ESEM observations, AFM images and distributions curves of adhesion forces confirm the distribution of MgSt 2 in discrete layers onto the silica gel surface. The mean adhesion force obtained for the different particulate systems reveal a no linear evolution to the adhesion force measured for the pure MgSt 2 . From these results, we are currently working on the development of a distribution model of MgSt 2 particles onto the surface of silica gel based on the adhesion force results. The sessile drop method reveals that a discrete deposition of the MgSt 2 allows a real modification of the hydrophilicity of the silica gel. Less than 3% of MgSt 2 is enough to increase the contact angle near the value of the pure MgSt 2 . This additional approach confirms the efficiency of the dry coating technique to design some new composite systems with controlled end use properties. Mean adhesion force (nN) for different particulate systems.

	5

Table 3 :

 3 Contact angle measurements (°) for different particulate systems.

Table 2 :

 2 Mean adhesion force (nN) for different particulate systems.

	Particles	Size (d 50 ) (µm)	Solid density (g/cm 3 )	Specific surface area	Pore Volume
		(Mastersizer 2000)	(Helium Pycnometer)	(SBET (m 2 /g)	(cm 3 / g)
	Silica gel	55 (D host )	2.07 (ρ host )	475	0.69
	MgSt	5 (d guest )	1.04 (ρ guest )	7.7	0.02

Acknowledgements

The authors would like to acknowledge Sylvie Delconfetto for the TG-DSC measurements and Christine Rolland for the ESEM images.