T Dudok De Wit
email: ddwit@cnrs-orleans.fr

SOME REGRESSION PROBLEMS IN SOLAR-TERRESTRIAL SCIENCES: LEARNING FROM MISTAKES

We address three timely regression analysis problems in solarterrestrial observations: the identification of trends in observations that exhibit a high level of internal variability, the choice of explanatory variables in the multilinear regression of climate data, and the identifi-cation of power laws in power spectral densities. In all three of them we focus on some common mistakes, and on how these may help facilitate critical reading of research in the field.

What not to do in regression analysis

.

Here, we concentrate on three timely issues that have been hotly debated in recent years. The first one is about identifying a drift in ionospheric observations, which is important for assessing the existence of long-term changes that may be related to global climate change. Technically, the problem is about the choice of the explanatory variables. The second issue is about the identification of a solar signature in climate records. Here the problem deals with the impact of collinearity on the regression analysis. The third and last example is about the identification of power laws in power spectral densities, with in addition the location of cutoff

frequencies. This is again a classical regression problem, with too many cases of erroneous conclusions. These simple illustrations are meant to facilitate critical reading of research in the field, and to help avoid being led astray by mistakes we all make.

2 Trend determination: Is the sky falling down?

The increasing concentration of greenhouse gases in the atmosphere causes enhanced cooling of the upper atmosphere, which should result in changes of atmospheric parameters. One of these is a trend in the ionosphere. This at-mospheric layer, however, is highly variable and driven by solar and geomagnetic activity, whose levels also changed throughout the 20 th century. If a trend exists at all, then it is most likely hidden in the natural variability. Disentangling these various signatures is an interesting but challenging regression problem (Laˇstoviˇcka et al. 2011). One the most representative ionospheric parameters is the critical frequency, called foF2, of the F2 layer. This layer has the densest electron concen-tration in the ionosphere, and extends from about 200 km to 500 km.

There has been a long quest for trends in foF2 (Laˇstoviˇcka et al. 2006). Ideas have progressively evolved as the pitfalls in the analyses have came to light. Figure 1 shows one of the longest records available, which is from the Julius-ruh station in Northern Germany. The trend in foF2 is heavily dominated by an 11-year solar cycle modulation, by an annual and semi-annual modulation that is related to the inclination of the Earth, and by more impulsive events that are due to geomagnetic activity. The earliest studies started by fitting a simple line to the Fig. 2. Same plot as Figure 1, showing in addition four trends (dashed lines) that were estimated by considering intervals running between different solar minima. time series, with the model

foF2(t) = a + bt + (t)
and then by interpreting the value of b; here, (t) stands for the residuals. However, it rapidly became evident that the slope b would heavily depend on the phasing with the solar cycle modulation. Both positive and negative trends were reported, depending on the endpoints of the time interval.

A first improvement consisted in applying trend analyses only to full periods of the solar cycle, as illustrated in Figure 2. However, it then became evident that the results would be biased by the non stationary properties of the solar cycle. Later models thus included the sunspot number S(t) as a proxy for solar activity foF2(t) = a + b t + c S(t) + (t) with variants. For example, should the slope b be considered as the trend, or should we fit instead the simpler model

foF2(t) = a + c S(t) + (t)
and search for a trend in the residuals (t)? It took some time to realise that the last two models lead to different conclusions because missing influential variables bias the results. Successive improvements followed unabated, and several hun-dreds of models were tested, including more advanced ones, involving for example multiscale decompositions.

Today, there is a global consensus for the trend in foF2 at Juliusruh to be relatively weak and negative, of the order of -0.02 to -0.015 MHz/year. However, these results are not conclusive yet, neither from a statistical, nor from a physical point of view. Some of the statistically important questions that remain to be addressed in order to move from a mere description of observations to actual statistical inference (with hypothesis tests), are:

• Measurement errors, and their propagation to the model parameters, have been largely ignored so far. The errors (and the natural variability) in foF2 are heteroscedastic, reaching a maximum at the time when solar activity also peaks.

• There is no single good proxy for solar activity, but rather an ensemble of them.

The most frequently used proxies are the sunspot number and the solar radio flux at 10.7 cm, which have quite different statistical properties. One emanates from a counting process, with Poisson-like noise, whereas the other has normal noise. Not surprisingly, the two lead to different results.

• The most critical issue probably is the proper choice of the influential vari-ables.

Missing influential variables are known to inflate type II errors, which are the failure to reject a false null hypothesis. aerosols, and the internal variability of the climate system, whose dominant mode is the El Ni˜no Southern Oscillation (ENSO). On top of these comes the solar forcing, which is represented by the sunspot number. This simple approach has received considerable attention in the literature. Many cases of misuse have been reported too [START_REF] Benestad | [END_REF]. A typical case is illustrated in Figure 3, in which the monthly-averaged global sur-face temperature anomaly ∆T is expressed as a combination of four contributions

∆T(t) = b0 + b1GHG(t) + b2AEROSOL(t) + b3ENSO(t) + b4S(t) + (t)
which are GHG, aerosols, ENSO and the Sun; stands for the residual. For a com-plete description of the various observables, see for example (Lean & Rind 2008). Fortunately, all regressors are almost totally uncorrelated here, so that collineary is not an issue. However, examples abound wherein the number of regressors, or collinearity, become serious issues (Crooks & Gray 2005). Well-established ap-proaches, such as partial least squares then apply [START_REF] Montgomery | Introduction to Linear Regression Analysis[END_REF]), but are often ignored in practice.

There are three frequent problems which we wish to concentrate on. The first one deals with the identification of the solar signal through bandpass filtering. On the time scales of interest (months to decades), the solar contribution mostly consists of a monochromatic signal, in contrast to all other explanatory variables.

The solar forcing indeed has a period of about 11 years, and a slowly varying amplitude. A simple strategy for emphasising its signature in the highly dynamic climate signal then consists in applying a 11-year bandpass filter, see for example (Gleisner & Thejll 2003).

Let ∆T * (t) then denote the bandpass filtered version of the dependent variable (in our case the global surface temperature anomaly) and S * (t) ≈ S(t) that of the sunspot number. In that case, we should have approximately ∆T * (t) = b4S * (t) + * (t).

The major advantage of this filtering is the elimination of all other regressors, assuming that they have little spectral power content around 11 years. Unfortu-nately, such a filtering artificially enhances the correlation between ∆T * (t) and the solar signal, thus giving the false impression that there truly is a solar signature hidden in the climate noise. Examples abound, wherein false correlations were reported because of such filtering [START_REF] Coughlin | [END_REF]. One solution would be to consider the regression coefficient b4, which is less impacted.

A second common problem arises from adding to the regressors all known contributions to the temperature anomaly, based on the widespread belief that this should eventually reduce the residual error, and consequently improve the description of the solar contribution. This is yet another example wherein model bias is reduced at the expense of a larger variance. Typical regressors, in addition to the four ones we just presented, are various internal modes of the climate system, such as the North Atlantic Oscillation (NAO). Unfortunately, these modes are not fully independent, and actually occasionally synchronise. Therefore, the model rapidly tends to become ill-conditioned. Orthogonalisation may cure here the numerical problem, but will not solve the physical one. What we need is a robust strategy for identifying the most significant regressors. This is once again a common problem in linear regression, but certainly not a trivial one.

Various interesting solutions have been developed for that purpose. One of them is based on the error reduction ratio (ERR), which is a criterion for pruning the set of regressors (Korenberg et al. 1988). Using matrix notation, let y be the vector containing the dependent variable, X the matrix of regressors, b the model coefficients, and e the residuals:

y = Xb + e.
The matrix X = WA can be decomposed into the product of an upper triangular matrix A and a matrix with orthogonal columns W. Then The strategy then consists in determining the ERR for each regressor, and select-ing those that offer the highest ERRs until their sum reaches a predefined level. This criterion is widely used for trimming linear and nonlinear parametric models [START_REF] Billings | Nonlinear System Identification: NARMAX Methods in the Time, Frequency and Spatio-Temporal Domains[END_REF], but dissemination outside of that field has been slow.

A third interesting problem is the frequent association between regression models and measures of correlation. On many occasions, when the influential variables are not properly known, or too numerous, then the correlation between one of them and the dependent variable is used as an alternative for their link, see for example (Sfˆıc˘a & Voiculescu 2014). The usual gauge of correlation is Pearson's parametric correlation coefficient. Spearman's non parametric rank correlation coefficient has received surprisingly little attention, although it is more appropriate for handling variables that have a monotonic but nonlinear relationship (Feigelson & Babu 2012).

Unfortunately, correlation coefficients are frequently interpreted as a measure of the influence of the explanatory variables on the dependent variable, which is incorrect (King 1986). Indeed, nothing attributes causal or sensitivity assumptions to the correlation coefficient. [START_REF] Cook | Applied Regression Including Computing and Graph-ics[END_REF] express this idea in a more elaborate way when stating that regression analysis aims at understanding "as as far as possible with the available data how the conditional distribution of the response y varies across subpopulations determined by the possible values of the predictor or predictors".

There are few cases, however, in which the correlation coefficient does bring some value to other statistical tools. One of them is the comparison of two equa-tions with the same dependent variable but not the same set of explanatory vari-ables. Fortunately, this is precisely what many of us are after when seeking to identify the most influential variables.

4 Power laws: Are there lines everywhere?

Power laws are ubiquitous in nature, and are often considered as a signature of selfsimilarity, possibly caused by critical processes [START_REF] Sornette | Critical Phenomena in Natural Sciences: Chaos, Fractals Self-Organization and Disorder: Concepts and Tools[END_REF]Aschwanden 2013). A quantity x follows a power law if it is drawn from a probability distribution p(x) ∝ x -α .

In astrophysics, power laws often show up in distributions, for which the Poisson statistics apply (Maschberger & Kroupa 2009;Andreon & Hurn 2013). Here, we consider instead empirical data for which p(x) is not necessarily a distribution. The noise properties differ, and so do the solutions for estimating the slope, or scaling parameter α. Typical examples p(x) are power spectral densities, for which p(ω) ∝ ω -β .

In most power law regression problems there are two objectives: one is to estimate the slope α, such as the -5/3 spectral index in Kolmogorov's model of fully developed turbulence. A second, and often even more important objective is about locating the transition from this power law to another functional expression, which then tells us about the characteristic scale of the system. There should be a third objective, which is to determine whether the power law model is actually the best one, but this aspect is often left out.

The quest for power laws is particularly frequent in solar wind turbulence. Here, the power spectral densities of the electric or magnetic fields are the key to the understanding of the underlying microphysical processes. Of particular inter-est are their extension down to small (so-called kinetic) scales, where dissipation processes set in (Sahraoui et al. 2009;[START_REF] Alexandrova | [END_REF].

Figure 4 illustrates a typical power spectral density from the solar wind. Fol-lowing common practice, we identify power laws by the approximately straight line fit in a loglog representation. The human eye is indeed remarkably good in detecting such straight lines. It can also be easily led astray by other types of distributions, such as log-normal or exponential. Examples abound, in which so-called evidence for power laws could be equally well, or even better explained by other distributions. As a rule of thumb, straight line fits should not be applied to ranges that cover less than a decade. Even a decade, however, may be too optimistic, not to mention the excessive number of digits used to express the value of the slopes. Let {xi}i=1,...,n be the observed values of a random variable, which is known to follow a power law in the interval xmax ≥ xi ≥ xmin. Under the condition where xmax xmin, the maximum likelihood estimator of the slope is

αˆ = 1 + n n ln xmin i -1 i=1 .
x

Unfortunately, this estimator is very sensitive to the initial guess of xmin, whose underestimation can severely impact the value of αˆ. [START_REF] Clauset | [END_REF] provide useful guidance for properly estimating both the interval and the slope by max-imum likelihood. But this is only part of the challenge. The second and more difficult challenge is their validation, which still remains a largely unsolved prob-lem when considering the bounds of the interval. Bootstrapping is one of the ap-proaches that may help us estimate confidence intervals (Feigelson & Babu 2012). This problem of estimating the power law is equally challenging when p(x) is not a distribution. However, a new and even more challenging problem arises: the error-invariable now cannot be neglected anymore. That is, the error in the independent variable x can substantially affect the conclusions, in addition to the error in the dependent variable p(x).

Various solutions have been proposed for dealing with this error-in-variable case. For least squares regression, there is the natural generalisation to total least squares. Let us, however, take one step back. With standard maximum likelihood approaches, we are answering the question: "How do the data I observed match my model?". Actually, we should rather ask "What do I know about the model parameters even before collecting the data?". The difference is not just philosophical, because the two questions lead to substantially different approaches, a frequentist one in the first case, and a Bayesian one in the second case [START_REF] Gelman | Bayesian Data Analysis[END_REF]. A discussion on their differences is definitely beyond the scope of this text. The main asset of the Bayesian approach lies in its ability to explicitly add prior information that may help further constrain the power law model. In cases where the statistics on the number of events is poor, such as with low photon counts, the Bayesian approach excels in extracting information from noisy data [START_REF] Van Dyk | [END_REF]). In our case, the Bayesian setting enables the estimation of the best range [xmin,xmax] and slope α simultaneously, rather than sequentially.

However, it would be fair to conclude that the numerical implementation of the solution can be challenging too.

Conclusion

Not surprisingly, regression analysis is key concept solar-terrestrial sciences, with implications that may sometimes have far-reaching societal consequences. In many applications, there is a widespread belief that issues such as collinearity, and lack of knowledge of what the influential variables should be, are some state of nature against which nothing can be done. The examples we have shown, however, suggest that there is much room left for improvementprovided we can learn from these mistakes and not repeat them. Therefore, in some sense, the mistakes that were risen in this short article, should be seen as a strong incentive for going back to basics and determining whether our model matches our true assumptions. There is no doubt here that the Bayesian framework will strongly help in a near future, as it has already in astrophysics.

Fig. 1 .

 1 Fig. 1. Raw critical frequency foF2, measured daily at noon at the Juliusruh station (54.6N, 13.4E). The straight line has been computed without any prior reduction of data to solar or geomagnetic activity.

Fig. 3 .

 3 Fig.3. The global surface temperature anomaly (the dependent variable) and its four re-gressors, from top to bottom: the concentration of GHG, the concentration of tropospheric aerosols, the strength of ENSO, and the level of solar activity. Monthly averages are used. All quantities are in arbitrary units, except for the temperature anomaly.

i

  The three terms respectively describe the variance of the dependent variable, of the explanatory variables, and the unexplained variance. From this, we define the ERR

Fig. 4 .

 4 Fig. 4. Power spectral density of a component of the electric field (Ey) and a component of the magnetic field (Bz) as measured in solar wind turbulence, respectively by the EFW and FGM/STAFF-SC instruments onboard the Cluster 4 spacecraft. The straight lines are power law fits of the spectra. Figure adapted from Sahraoui et al. (2009).

The first two questions suggest that Simpson's paradox(Pearl 2009) may also play a role here. In this paradox, a trend may change quite substantially depending on how the data are aggregated. Since both the statistical and the physical properties of the observations change with their phasing with respect to the solar cycle, the search for a better model, including a rethinking of the way the data are sampled, is definitely needed.3 Multilinear regression: Desperately searching for solar signatures One of the key issues in climate research is to scientifically ascertain the mecha-nisms that are responsible for recent climate change. The Sun is in the spotlight because the level of its contribution to global warming has far reaching political implications. Although there is a clear consensus today on the prevalent effect of anthropogenic greenhouse gases in global warming, the quantification of the solar contribution remains a difficult and ongoing challenge[START_REF] Stocker | Climate Change 2013 -The Physical Science Basis. Work-ing Group I Contribution to the Fifth Assessment Report of the IPCC[END_REF]. The reason for this is the extreme complexity of the dynamic response of the at-mosphere and the couplings of the mechanisms involved, which preclude simple sensitivity analyses.In this context, several simplified approaches have been developed for testing the sensitivity of the climate system to the solar forcing. A common one is a multilinear regression analysis wherein a climate signal, typically the global surface temperature anomaly, is modelled as a linear superposition of various contributions(Lean & Rind 2008;[START_REF] Gray | [END_REF]. Foremost among these is the concentration of greenhouse gases (GHG), volcanic cooling expressed in terms of tropospheric


Doru Constantin
email: constantin@lps.u-psud.fr

Brigitte Pansu

Marianne Impéror

Patrick Davidson

Franc ¸ois Ribot

Repulsion Between Inorganic Particles Inserted Within Surfactant Bilayers

Keywords: numbers: 87.16.dt, 61.05.cf, 82.70.Dd

We study by synchrotron small-angle X-ray scattering highly aligned lamellar phases of a zwitterionic surfactant, doped with monodisperse and spherical hydrophobic inorganic particles as a function of particle concentration. Analysis of the structure factor of the two-dimensional fluid formed by the particles in the plane of the bilayer gives access to their membrane-mediated interaction, which is repulsive, with a contact value of about 4 kBT and a range of 14 Å. Systematic application of this technique should lead to a better understanding of the interaction between membrane inclusions.

In the last decades, much effort was dedicated to the understanding of self-assembled membranes, and in particular of the interaction between membrane inclusions and the host bilayer. This is a challenging problem, since the membrane must be considered as a many-particle system, its properties being collectively determined by the assembly and not by the chemical properties of the individual molecules [1]. Notwithstanding the complexity of the system, the concepts developed in soft matter physics should be operative in this context and even 'simplified' models could yield valuable information. For this reason, a considerable body of work deals with the theoretical modelling and numerical simulation of such mixed systems, aiming in particular to determine the membranemediated interaction between inclusions [2]. Thorough understanding of these interactions could make a substantial contribution to topics as diverse as the formulation of new composite systems and understanding the activity of membrane proteins.

However, this theoretical work was not yet matched by enough experimental results. The first such data was obtained by directly measuring the radial distribution function of membrane inclusions using freeze-fracture electron microscopy (FFEM) [3]. These data were compared to liquid state models and could be described by a hard-core model with, in some cases, an additional repulsive or attractive interaction [4]. FFEM was not extensively used, undoubtedly due to the inherent experimental difficulties; moreover, it is not obvious that the distribution measured in the frozen sample is identical to that at thermal equilibrium.

Considering the typical length scales to be probed, X-ray and neutron scattering techniques [5] are uniquely adapted to the study of this problem. As an example, the pores formed by the antimicrobial peptide alamethicin in dimyristoylphosphatidylcholine bilayers were shown to repel each other [6]. These studies are however hindered by the low scattering contrast of the proteins (and, in many cases, by their scarcity), whence the practical (and also the fundamental) interest of finding out whether other -perhaps more adapted-particles can be inserted within membranes.

The purpose of this Letter is to show that self-assembled bilayers can be doped with significant amounts of (hydropho-bic and charge-neutral) hybrid nano-objects and that these probes can be used to determine accurately the membranemediated interaction. Their use presents substantial advantages: they are "rigid" (fixed atomic configuration) and perfectly monodisperse, therefore imposing a well-defined membrane deformation, whereas membrane proteins can assume various conformations; their scattering contrast is high (due to the presence of metal atoms); their surface properties can be tailored by changing the nature of the grafted ligands. Conceptually, these inclusions are also easier to model, since they do not "break" the surface of the monolayer and hence there is no contact line, where the choice of the boundary condition would be somewhat delicate [7].

This approach can help clarify long-standing questions, such as: Is a continuous model sufficient for an accurate description of the membrane and, if so, down to what length scale? What are the relevant parameters and how can they be measured? What are the specificities of mixed bilayers, and in particular of lipid membranes (microscale separation, raft formation etc.)? In the long run, systematic studies should yield a clearer picture of membranes as two-dimensional complex systems.

The nanoparticles used here are butyltin oxo clusters {(BuSn) 12 O 14 (OH) 6 } 2+ (4-CH 3 C 6 H 4 SO - 3 ) 2 , denoted by BuSn12 in the following. They were synthesized and characterized as described in reference [8] (for their structure, see Figure 1 in this reference). They were dissolved in ethanol at a concentration of 23.47 wt.%.

The zwitterionic surfactant, dimethyldodecylamine-Noxide (DDAO) was purchased from Sigma-Aldrich and dried in vacuum (using a liquid nitrogen-cooled solvent trap) for 20 h. No weight loss was observed after this step, so we conclude that the surfactant was dry as supplied (see [9] for a detailed discussion). DDAO was then dissolved in isopropanol at a concentration of 23.02 wt.%.

We mixed the BuSn12 and DDAO solutions to yield the desired particle concentrations, and then dried the mixtures in vacuum; the final (dry) mass was about 200 mg for each sample. We then added water at a concentration of 15-27 wt.% of the final hydrated mixtures, which are thus in the fluid lamellar L α phase (see the phase diagram of the undoped system in [9]). The molecular weight of DDAO is 229.40 (Sigma-Aldrich), its density is 0.84g/cm 3 and the thickness of the bilayer is 25 ± 1 Å [10], yielding an area per surfactant molecule A DDAO = 37.8 Å2 . For the BuSn12 particles, we take a molecular weight of 2866.7 and a density of 1.93g/cm 3 [8].

Using these values and neglecting the increase in bilayer surface due to the inserted particles yields the (two-dimensional) number density of particles in the plane of the membrane, n.

The samples were prepared in flat glass capillaries (Vitro-Com Inc., Mt. Lks, N.J.), 100 µm thick and 2 mm wide by gently sucking in the lamellar phase using a syringe. The capillaries were flame-sealed. Good homeotropic alignment (lamellae parallel to the flat faces of the capillary) was obtained by thermal treatment, using a Mettler FP52 heating stage. The samples were heated into the isotropic phase (at 130 • C) and then cooled down to the lamellar phase at a rate of 1 • C/min.

The small-angle x-ray scattering measurements were performed at the European Synchrotron Radiation Facility (ESRF, Grenoble, France) on the bending magnet beamline BM02 (D2AM), at a photon energy of 11 keV. See reference [11] for more details on the setup. The data was acquired using a CCD Peltier-cooled camera (SCX90-1300, from Princeton Instruments Inc., New Jersey, USA) with a resolution of 1340 × 1300 pixels. Data preprocessing (dark current subtraction, flat field correction, radial regrouping and normalization) was performed using the bm2img software developed at the beamline.

The incident beam was perpendicular to the flat face of the capillary (parallel to the smectic director, which we take along the z axis.) Thus, the scattering vector q is mostly contained in the (x, y) plane of the layers, and the measured scattered signal I(q) probes inhomogeneities of the electron density in this plane. Since the bilayers form a two-dimensional liquid, the scattering pattern exhibits azimuthal symmetry: I = I(q = |q|). We also measured the scattering intensity of two BuSn12/ethanol solutions in the same type of capillary. The accessible scattering range was 0.04 < q < 0.9 Å-1 . Since the electron density of the butyl chains is similar to that of the dodecyl chains within the bilayers and to that of ethanol, we expect the electron contrast of the particles to be due exclusively to their inorganic core, which is slightly ovoidal, with an average radius of 4.5 Å. Indeed, the intensity at higher scattering vectors (q > 0.5 Å-1 ) is well described for all samples by the form factor of a sphere |F f (R, q)| 2 , with a radius R = 4.5 ± 0.2 Å used as a free fitting parameter. The interaction between particles is described by the structure factor, defined as S(q) = I(q)/|F f (R, q)| 2 [START_REF] Chaikin | Principles of Condensed Matter Physics[END_REF]. The structure factors thus obtained are shown in Figure 1 for all in-plane concentrations, listed in the caption.

The additional interaction is viewed as a perturbation with respect to the hard core (disk or sphere) model, taken into account via the random phase approximation (RPA) [START_REF] Andersen | [END_REF]. In this approach, one obtains the direct correlation function of the perturbed system c(r) from that of the reference system FIG. 1: Structure factors S(q) of the 2D fluid formed by the BuSn12 particles in the plane of the membranes, for different concentrations, from top to bottom: n = 0.217, 0.231, 0.429, 0.434, 0.451, 0.584, 0.801, 1.102, 1.275, 1.976, 2.147, and 2.304 10 -3 Å-2 . The curves are shifted vertically in steps of 0.2. Gaps in the curves correspond to the presence of (weak) lamellar peaks, due to occasional alignment defects, mainly taking the form of oily streaks. c ref (r) as: c(r) = c ref (r) -βU (r) [START_REF] Hansen | Theory of Simple Liquids[END_REF] or, equivalently:

nβ U (q) = S -1 (q) -S -1 ref (q) (1)
with β = (k B T ) -1 .

In three dimensions, the reference structure factor S 3D ref for BuSn12 particles dissolved in ethanol is given by a hard sphere interaction (in the Percus-Yevick approximation [START_REF] Wertheim | [END_REF]16]) with a hard-core radius of 4.5 Å; the numerical particle density n 3D (in 3D) is determined from the mass concentration of the solutions.

In two dimensions, an analytical form for the structure factor S 2D ref of hard disks was given by Rosenfeld [17]; we use the same core radius of 4.5 Å as above. The Fourier transform of the interaction potential, U (q), obtained by applying relation (1) to the data in Fig. 1, is shown in Fig. 2 for all concentrations. For ease of calculation, both in two and three dimensions we model the interaction U (r) by a Gaussian, with amplitude U 0 and range ξ:

U (r) = U 0 exp -(r/ξ) 2 /2 (2)
FIG. 2: The Fourier transform of the interaction potential, U (q), (open dots) obtained from the structure factors in Figure 1 in the RPA approximation (1). The solid line is a Gaussian model. The shaded area is chosen to cover most experimental points and is also delimited by two Gaussians (see text for details). The data points far above and far below the shaded range correspond to n = 0.429 and 0.584 10 -3 Å-2 , respectively (third and sixth curve from the top in Figure 1).

with U (q) the Fourier transform of U (r), given by:

U (q) =    2π U 0 ξ 2 exp -(qξ) 2 /2 in 2D (2π) 3/2 U 0 ξ 3 exp -(qξ) 2 /2 in 3D (3)
With the exception of two curves (for n = 0.429 and 0.584 10 -3 Å-2 ) that are anomalously high or low, respectively, all data points are reasonably well covered by the shaded area in Figure 2, limited by two Gaussians (2D case in Eq. 3), with U 0 = 3.22 k B T and ξ = 15.71 Å for the lower bound and U 0 = 6.93 k B T , ξ = 12.86 Å for the upper bound. The solid line, giving approximately the midline of the shaded area, is described by U 0 = 4.75 k B T and ξ = 14.14 Å.

The apparent negative values of U (q) around q = 0.2 Å-1 in Figure 2 are probably due to the enhancement of the interaction peak of S(q) with respect to the reference potential by the presence of the repulsive interaction; this feature is not captured by the simple RPA treatment. The (real-space) interaction potentials U (r) corresponding to the aforementioned values of U 0 and ξ are plotted in Figure 3, using the same convention. Thus, the solid line is the estimate for the interaction and the shaded area the uncertainty due to the value spread in Figure 2. Of course, for r < 2R = 9 Å the repulsion is due to the hard core (solid vertical line).

To make sure that the interaction is induced by the membrane, we also show the repulsion measured for BuSn12 particles in ethanol (lower solid line in Figure 3); its amplitude and range are clearly lower than for the interaction in the membrane (for reference, U 0 = 4 k B T and ξ = 7 Å). It is probably due to the steric repulsion between the butyl chains grafted onto the inorganic core. The main result of this work is hence that the (twodimensional) interaction between BuSn12 particles inserted within DDAO bilayers can be described by a potential of the form (2), with U 0 = 5 ± 1 k B T and ξ = 14 ± 1 Å. Below, we discuss briefly the various theoretical predictions, but a direct comparison cannot be made since they generally concern cylindrical inclusions that traverse the membrane, without being covered by the monolayers.

The simplest approach is to model the membrane as a continuous medium and write the free energy of the bi-layer+inclusions system in terms of elastic deformation, described by the Helfrich Hamiltonian [18], where the molecular properties of the bilayer are abstracted into mesoscopic parameters, such as the bending and stretching moduli and the spontaneous curvature of the monolayer. A systematic study was performed in the group of Pincus [19][20][21]. Their results highlight the importance of the spontaneous curvature of the monolayers c 0 : the interaction is attractive when c 0 vanishes (the elastic energy is minimized by the aggregation of inclusions) but can become repulsive when the deformation induced by the inclusion is such that the curvature of the individual monolayer has the sign of c 0 . In our case, this would mean that the DDAO monolayers tend to have c 0 > 0 (convex head-groups, see the sketch in Figure 3) However, inclusions also perturb the structure of the membrane by restricting the conformation of the lipid chains in their vicinity [22], and this can lead to significant interaction, even in the absence of "large scale" bilayer deformation (no hydrophobic mismatch.) In particular, Lagüe et al. [23] find that "smooth" hard cylinders repel each other in some lipid bilayers; the amplitude and range of this interaction is in qualitative agreement with our results (for cylinders with a 5 Å radius in dioleoylphosphatidylcholine bilayers, the repulsive lipid-mediated interaction has a maximum value of 7 k B T and extends 20 Å from contact.) Nonetheless, the interaction is highly dependent on the chemical structure of the lipid, so a meaningful comparison is difficult to make using the available data.

Recently, some efforts were made [24,25] to account for both effects (elastic energy and restrictions on chain conformations) within an extended model, including as variables both the variation in membrane thickness and the local molecular tilt. Their coupling removes the symmetry between the positive and negative spontaneous curvature values, leading to repulsion only for c 0 > 0 combined with negative hydrophobic mismatch (bilayer "pinching".) At first sight this prediction seems to be at odds with our system, but the difference in geometry plays an important role: in our case, we expect the particles to be covered by positively curved monolayer caps (energetically favourable), while the calculation in [25] considers the inclusions as vertical cylinders. A more accurate comparison should be very interesting.

We have shown that nano-objects can be used to probe the properties of self-assembled bilayers; conversely, one can envision using ordered surfactant phases (combining a high degree of order with easy processability and very good wetting properties) to organize and align such objects in view of applications.

The ESRF is gratefully acknowledged for the award of beam time (experiment 02-01-732) and we thank C. Rochas for competent and enthusiastic support. A. Dessombz and A. Poulos are acknowledged for helping with the synchrotron experiments.

FIG. 3 :

 3 FIG.3: Interaction potential U (r) of the BuSn12 particles within the bilayers, obtained by taking the inverse Fourier transform of the solid line and the shaded area in Figure2(see text for the numerical values). The lower curve is the interaction potential of the particles in ethanol. The solid vertical line marks the hard core interaction with radius 4.5 Å


Keevin Béneut

Doru Constantin
email: constantin@lps.u-psud.fr

Patrick Davidson

Arnaud Dessombz

Corinne Chanéac

Magnetic nanorods confined in a lamellar lyotropic phase

The dilute lamellar phase of the nonionic surfactant C 12 EO 5 was doped with goethite (iron oxide) nanorods up to a fraction of 5 vol%. The interaction between the inclusions and the host phase was studied by polarized optical microscopy (with or without an applied magnetic field) and by small-angle x-ray scattering. We find that when the orientation of the nanorods is modified using the magnetic field, the texture of the lamellar phase changes accordingly; one can thus induce a homeotropic-planar reorientation transition. On the other hand, the lamellar phase induces an attractive interaction between the nanorods. In more concentrated lamellar phases (under stronger confinement) the particles form aggregates. This behaviour is not encountered for a similar system doped with spherical particles, emphasizing the role of particle shape in the interaction between doping particles and the host phase.

Introduction

In recent years, the field of soft matter physics has witnessed a surge of activity in the area of hybrid organic-inorganic materials. 1 This sustained interest was of course motivated by the manifold applications of these systems, 2 and also by novel fundamental issues related to the interaction between the two components. In many cases, the materials are obtained by dispersing solid nanoparticles in a "soft" continuous matrix, formed by surfactants, polymers, emulsions etc. By a careful choice of the components, one tries to combine the specific properties of the particles (catalytic, optical, magnetic etc.) and the processability of the host phase.

Obviously, such dispersions also represent a new field among soft matter systems; its novelty with respect to the "classical" colloidal solutions resides in the complexity of the matrix. Beyond the theoretical interest of this study there is a very practical one: for what parameter values (particle size and shape, structure and elastic moduli of the matrix etc.) is the dispersion stable? Can the confinement imposed by the host phase lead to ordering of the inclusions? What are the ensuing applications? None of these questions can be answered without a thorough understanding of the interaction between the host phase and the inclusions.

In the case of nanoparticles dispersed in a liquid crystalline matrix, one should naturally consider the effect of the elastic and anisotropic medium on the interparticle potential. For lamellar phases, the effect of smectic elasticity was modeled in detail, [3][4][5] but the experimental studies are still unsubstantial. Conversely, the inclusions can change the interaction potential between membranes, and thus its repeat distance 6 and elastic moduli. 7 Clearly, in a composite system one must study:

• The influence of the confinement (due to the host phase) on the inclusions.

• The changes induced by the particles in the structure of the matrix.

Both these aspects are illustrated in a system that we formulated recently: the host lamellar phase is the C 12 E 5 /hexanol/water system, with C 12 E 5 the nonionic surfactant penta(ethylene glycol) monododecyl ether, and the inclusions are iron oxide nanorods (goethite) with complex magnetic properties. 8 We demonstrate the attractive interaction between the particles induced by the lamellar matrix and show that, when the confinement becomes too strong (i.e. the lamellar repeat distance is too small) the particles aggregate, even when their concentration is very low. On the other hand, we show the action of the inclusions on the texture of the lamellar phase: when the nanorods are oriented using an applied magnetic field, the bilayers 'follow' and (at high field) they align perpendicular to the field. Perfect planar monodomains can thus be obtained, and the alignment persists after removal of the field. The magnetic field has no effect on the pure lamellar phase (without inclusions).

Lamellar lyotropic phases doped with small and spherical magnetic particles have already been formulated 9,10 and their structural 7,11 and magnetic 12,13 properties were studied in detail long ago. In contrast, our study deals with large, anisotropic particles and with their interaction due to the confinement in the lamellar phase. It is also noteworthy that gold nanorods can be confined in lamellar phases of block copolymers, as reported recently. 14 This hybrid system appears promising for the preparation of surface layers of magnetic nanoparticles with well-defined spacing and orientation (controlled by the host-induced interaction), with applications, for instance, in high-density storage media. 15,16 The viscoelastic properties of the lamellar phase are also interesting in view of fine-tuning the deposition process (e.g., by spin-coating).

Experimental

Goethite (α -FeOOH) is an iron oxyhydroxide, widely used as a pigment. The nanorods were synthesized according to well-established protocols. 17,18 Their dimensions are of the order of 150 × 25 × 10 nm 3 (length × width × height). 19 The surface of the particles is hydroxylated, with a surface charge of 0.2 C m -2 at pH = 3 and with an isoelectric point corresponding to pH = 9; see 19 for more details. Bulk goethite is antiferromagnetic, 20 but the nanorods bear a permanent magnetic dipole µ ∼ 1200 µ B along their long axis, probably due to uncompensated surface dipoles (µ B = 9.274 10 -24 J/T is the Bohr magneton). Furthermore, the easy magnetization axis is perpendicular to this direction so that, at high applied fields, the induced magnetic moment overtakes the permanent one and the orientation of the rods switches from parallel to perpendicular to the field at a critical value B ∼ 250 mT.

The surfactant, C 12 EO 5 , was acquired from Nikko and the hexanol from Fluka; they were used without further purification. The phase diagram of the C 12 EO 5 /H 2 O mixture was determined more than 20 years ago. 21 Ever since, it was extensively studied due to the presence of several mesophases, and especially of a lamellar phase that can be swollen up to a few percent of membrane fraction. This dilute lamellar phase appears at fairly high temperature, but it was shown that it can be brought down to room temperature by the addition of a co-surfactant such as hexanol. 22,23 We used a hexanol/C 12 EO 5 ratio of 0.35 by weight, corresponding to a molar ratio of 1.3 (hexanol molecules for each surfactant molecule). The temperature domain of the lamellar phase changes with dilution, but it extended at least between 17 and 32 • C for all our samples. The membrane thickness is δ ≈ 2.9 nm. 22 Silica particles with a nominal diameter of 27 nm were obtained from Sigma-Aldrich as concentrated colloidal suspensions (Ludox TMA 34) in deionized water (34 wt.%). We found a pH of 7 for the initial suspension.

Concentrated stock solutions of C 12 EO 5 /hexanol/H 2 O were mixed with colloidal suspensions (of goethite or silica) and deionized water to yield the desired volume fractions of membranes and doping particles. The samples were contained in flat glass capillaries, 50-100 µm thick (Vitrocom) and aligned in homeotropic anchoring (by thermal annealing).

The magnetic field was applied using a home-made setup based on permanent magnets with a variable gap. Fields of up to 0.9 T can be achieved. For polarized microscopy observation we used an Olympus BX51 microscope with a rotating stage and long working distance objectives (with a 5× or 10× magnification.)

Small angle x-ray scattering (SAXS) experiments were performed at the ID02 station of the European Synchrotron Radiation Facility synchrotron in Grenoble, France. The incident beam had a wavelength λ = 0.0995 nm, and the sample-detector distance was 5 m. The scattered x-rays were detected with a specially developed CCD camera. A detailed description of the experimental setup can be found in reference. 24 The q range over which the data could be reliably collected was 0.018 < q < 0.6 nm -1 . The flat faces of the capillaries were set perpendicular to the x-ray beam.

Results and Discussion

Phase diagram of the doped system

The first step of the study was determining the phase diagram of the system, more specifically the range of confinement (controlled by the membrane volume fraction) for which the particles can be added to the phase without demixing and their maximum concentration. We started by preparing mixtures with a volume fraction of goethite φ g = 0.5, 1, 1.5 and 2 % and a membrane volume fraction φ m = 4.4, 7.2, 10.2 and 14.4 %. φ g is defined as the ratio between the volume of goethite particles and the total volume. φ m is the sum of the hexanol and surfactant volumes divided by the total volume (we assume that mixing volumes are negligible). The particles remained well dispersed in the dilute lamellar phase (φ m = 4.4 and 7.2 vol. %) for all values of φ g investigated. In the concentrated phases, on the other hand, particle aggregation was discernable after a few hours and was very clear after a few days, even at the lowest particle concentration (see Figure 1). We then prepared samples with φ m = 7.2 vol. % and φ g up to 5 vol. %. All these samples have been stable for months. 25 We conclude that a significant amount of goethite can be dispersed in the lamellar phase as long as the membrane fraction remains below a threshold in the interval 7.2 < φ m < 10.2 vol %, corresponding to a lamellar repeat distance 28 < d < 40 nm. The upper transition temperature of the lamellar phase (towards the sponge phase) exhibits no significant variation as a function of the doping fraction φ g . The mixing of the nanorods and the lamellar phase presumably leads to an energy gain owing to the formation of hydrogen bonds between the hydrated surface of the particles and the surfactant heads. 26 At this point, we have no convincing explanation for the threshold value of the repeat distance. The most plausible connection is that, as the lamellar phase becomes more concentrated, its elastic moduli increase and so does the interaction between particles, 3 to the point of inducing aggregation.

It should also be noted that the threshold value is of the order of the particle width. An alternative explanation would therefore be that aggregation occurs when rotation about the long axis of the particles is hindered (and the particle loses a degree of freedom). It is however not clear whether this explanation is compatible with a strong interaction between the particles and the surfactant heads.

Magnetic field effect

Magnetic field measurements were performed on flat glass capillaries, 50 µm thick and 1 mm wide. The field was applied in the plane of the capillary, perpendicular to its long axis.

We started by applying an increasing field (from 0 to about 0.8 T) to a sample oriented in very good homeotropic anchoring (obtained by annealing overnight close to the transition temperature to the sponge phase). The membrane volume fraction was φ m = 7.2 % and the goethite volume fraction φ g = 1.5 %. A few very thin oily streaks persisted. The succession of images is shown in Figure 2. At low field, the transmitted intensity increased with the field up to about 0.15 T (nanorods aligned along the field); it then decreased to 0 at 0.25 T and increased again at higher field, as the rods aligned perpendicular to the field. Starting from the initial homeotropic anchoring, above 0.3 T the existing oily streaks became more pronounced and new ones nucleated; the texture gradually switched to planar anchoring, with the smectic director along the field. This crossover corresponds to the value at which the particle orientation changes from parallel to perpendicular to the field in aqueous solution. 8 Consequently, we infer that, as the particles turn, the lamellae follow, presumably due to the strong association between the goethite nanorods and the surfactant heads.

Very good planar anchoring can thus be obtained, as shown in Figure 3. The field was applied overnight; during this time, the sample was kept at a temperature about 1 • C below the transition temperature to the sponge phase, in an oven. The temperature was then slowly decreased to its ambient value. Figure 3a shows the sample, in natural light and between crossed polarizers, under high field. It contains two domains separated by a wall. By rotating the sample between the polarizers (which remain parallel to the the sides of the photo) each domain can be extinguished (Fig. 3b). The disorientation between the domains can thus be estimated at 3.6 • . The field was then progressively decreased (in steps of 0.1 T every 10 minutes). Some focal conic textures developed during the process, but they annealed after a couple of hours. The resulting texture at zero field is shown in Fig. 3c; it was stable for days.

The effect of the magnetic field is similar to that observed in ferrosmectic phases obtained by doping dilute lamellar phases with small ferromagnetic particles. 9,12 Indeed, these authors also observe a reorientation transition, signaled by the appearance of focal conic defects in homeotropic samples under the influence of a magnetic field applied along their director axis (normal to the layers). However, in their system the layers tend to align along the field, while in our case they prefer to be perpendicular to it above the critical value.

Interaction induced by the lamellar phase

We used x-ray scattering to study the interaction between colloidal particles in the lamellar phase and in solution. It is well-known 27 that the intensity scattered by a collection of identical particles can usually be written as the product I(q) = |F (q)| 2 × S(q) of a form factor, |F (q)| 2 , dependent only on the size and shape of the individual particle, and a structure factor S(q) quantifying the interactions between particles (S(q) = 1 in the absence of interactions).

As form factor we used the intensity scattered by a very dilute dispersion (φ g = 0.066%) in the lamellar phase at the same membrane concentration as for the curves shown in Figure 4, namely φ m = 7.2%. After background subtraction, the scattering curves for (L α + goethite) systems presented in the following were divided by this signal and normalized to one at large q vectors.

The first observation is that all structure factors shown in Figure 4 exhibit a marked increase at small angles (below 0.1 nm -1 ), clear sign of a strong Goethite in the lamellar phase phase. The membrane volume fraction is φ m = 7.2% and the volume fraction of goethite φ g is indicated for each curve. The vertical dashed line q = 0.018 nm -1 is the limit below which the data exhibits marked uncertainty (due mainly to background subtraction). Solid lines are exponential fits (see text).

f g = 0.5 % f g = 1.5 % f g = 3.0 % f g = 4.0 % f g = 5.0 %
attractive interaction. A quick estimate of the interaction range ξ can be obtained by fitting the data to an exponential decrease:

S(q) = 1 + A exp(-qξ) , (1)
yielding 30 nm < ξ < 50 nm. The attractive range is similar to both the width of the nanorods and the lamellar repeat distance. A more detailed study for different dilutions of the lamellar phase is needed to assess the nature of the interaction. The most interesting feature of this interaction is that it only appears under confinement (in the lamellar phase), and for anisotropic particles, as discussed below. As a reference system, we studied water dispersions of goethite particles at similar concentrations, see Figure 5. The form factor used was obtained from a very dilute aqueous solution, φ g = 0.066%. At lower particle concentrations the structure factor is negligible; at φ g = 7.3% a typical shape for hard-core systems, with a well-defined peak and an oscillation at higher q starts to appear, but its amplitude is moderate and the shape very different from that measured in the lamellar phase. We conclude that the effects described above (Fig. 4) are due to the presence of the confining lamellar phase.

The X-ray scattering data for the doped lamellar phases indicate that the lamellar phase induces an attractive interaction between the goethite nanorods. In order to assess the effect of shape we studied the same lamellar phase doped with a comparable concentration of silica beads, a system already described in the literature. [START_REF] Salamat | [END_REF] The membrane fraction was φ m = 7 vol% and the volume fraction of beads was φ s = 0.5, 1, 2 and 3 vol%. The samples were stable and homogeneous for months, although in reference, [START_REF] Salamat | [END_REF] the doped phase was found to be stable only for volume fractions up to φ s = 0.8 vol%. This discrepancy might be related to the difference in the presentation of the silica beads: they used Ludox TM solutions, at pH 9 and with relatively high salt concentrations.

Homeotropic doped samples were studied using the same procedure as for the goethite-containing phases. The resulting structure factors are shown in Figure 6. The scattering intensity I(q) for the dilute sample (φ s = 0.5 vol%) was used as form factor; it is well described by a polydisperse sphere model, with a radius R ∼ 13 nm and polydispersity p = σ/R ∼ 0.1. Similar values are obtained for the aqueous dispersions using an in-house rotating anode setup (data not shown).

The first observation is that for the silica beads dispersed in the lamellar phase the structure factors only exhibits a slight increase at small angles. fraction of goethite φ g is indicated for each curve. The vertical dashed line q = 0.018 nm -1 is the limit below which the data exhibits marked uncertainty (due mainly to background subtraction). If present, the induced attraction is thus much weaker for spheres than for rods. For wave vectors q > 0.06 nm -1 , the structure factors for φ s = 2 and 3 vol% are well described by a three-dimensional hard sphere interaction (in the Percus-Yevick approximation 29,30 ) with an effective hard-core radius of 19 nm, see Figure 6. In conclusion, the presence of the lamellar phase has no discernible effect on the interaction between silica spheres. Rigorously speaking, for spheres confined between rigid planes a two-dimensional (harddisk) interaction would be a more adequate description. We performed such an analysis using the analytical form for the structure factor given by Y. Rosenfeld 31 and obtain similar results, with a hard disk radius of 17 nm. More concentrated systems would be needed to discriminate between the 2D and 3D cases.

Conclusion

In conclusion, we formulated a nonionic lamellar phase doped with large magnetic nanorods (in comparison with the interlayer distance). The inclusions experience an attractive interaction under confinement, a feature absent in simple aqueous solutions of similar concentration or in systems of confined silica spheres. Under even higher confinement (membrane concentration), the nanorods aggregate. The interaction between the particles and the host phase is also apparent in the orienting effect of the inclusions on the lamellar phase.

Figure 1 :

 1 Figure 1: Lamellar phase doped with a goethite concentration φ g = 0.5 %, for a membrane concentration φ m = 7.2 % (a) and 10.2 % (b), two weeks after preparation. Left: in natural light. Right: between crossed polarizers. Aggregation of the nanorods is observed in the more confined system (b).

Figure 2 :

 2 Figure 2: A magnetic field was applied across a capillary of lamellar phase (membrane volume fraction 7.2 vol. %) doped with 1.5 vol. % nanorods. The field was increased from 0 to about 0.8 T (the field direction is shown in the first image and its values are given below each image). The images are taken between crossed polarizers parallel to the sides of the image.

Figure 3 :

 3 Figure 3: Two-domain area in a sample with φ m = 4.4 % and φ g = 1 %. a) Under strong field (834 mT), in natural light (left) and between crossed polarizers parallel to the sides of the image (right). b) By rotating the sample, total extinction is obtained for each of the domains. c) After field removal, the planar orientation persisted.

Figure 4 :

 4 Figure 4: Structure factors for goethite particles dispersed in the lamellar L α

Figure 5 :

 5 Figure 5: Structure factors for water dispersions of goethite particles. The volume

  lamellar phasef s = 1 vol % f s = 2 vol % f s = 3 vol %

Figure 6 :

 6 Figure 6: Structure factors for silica spheres dispersed in the lamellar phase. The volume fraction of silica particles φ s is indicated for each curve. The data for φ s = 2 and 3 vol% is fitted with a Percus-Yevick hard sphere model with an interaction radius of 19 nm.

Acknowledgements P. Panine, E. Belamie and A. Poulos are acknowledged for helping with the synchrotron experiments.

For Table of Contents Use Only


Salvatore Capozziello
email: capozziello@unina.it

Tomislav Prokopec
email: t.prokopec@uu.nl

Alessandro D A M Spallicci
email: spallicci@cnrs-orleans.fr

HAL is

From the infinitesimal scale of particle physics to the cosmic scale of the universe, there have been spectacular advances in physics during the past century. The understanding of the universe is largely based on electromagnetic observations. As photons are the messengers, it is very important for fundamental physics to test the foundations of electromagnetism. The same stands for gravitation, the dominant interaction at large scale, for which we have just witnessed the opening of the high frequency gravitational wave window.

The basic properties of the Universe are often believed to be well understood, with only minor adjustments left in front of us. The observations which may repre-sent a problem are interpreted to make them concordant or treated as peculiarities to be later dealt within the framework of an otherwise correct model. Further, assum-ing our current description of gravity and electromagnetism, we sometimes run into extremely complex and sophisticated cosmological models, partly based on purposely built "dark" or unknown ingredients. Ironically, we are awkwardly confronted with a largely unknown universe, and yet we adopt the term "precision cosmology".

Indeed, cosmology faces puzzling observations. First, galactic rotation curves, lensing and clustering do not appear to be generated by conventional baryonic matter.

Second, it is commonly, although not universally, accepted that cosmic expansion is accelerating. We are clearly still facing aspects of the universe that we do not under-stand.

While the standard approach might lead to the ultimate answers as promised, con-currently the scientific community has to explore other avenues. Comparing alternative and standard views, it is the sole way that science can advance.

Thus, a debate on the foundations and on the theories of physics with a bearing on cosmology, or conversely a debate on cosmology (theoretical or observational) with a bearing on the foundations and theories of physics are the two tracks along which we collected these four contributions, while others may follow. In this issue, the papers deal with massive photons, the microwave background, the cosmological red shift and inflation.


Florian Blachère
email: florian.blachere@utt.fr

Christophe Chalons
email: christophe.chalons@uvsq.fr

Rodolphe Turpault
email: rodolphe.turpault@u-bordeaux.fr

Very high-order asymptotic-preserving schemes for hyperbolic systems of conservation laws with parabolic degeneracy on unstructured meshes

Keywords: asymptotic-preserving schemes, diusion limit, nonlinear hyperbolic systems, high order nite volumes schemes. AMS : 35L50, 65M08

HAL is

Introduction

We are interested in the numerical approximation of the solutions to nonlinear hyperbolic systems of conservation laws with possibly sti relaxation source terms. Such systems are ubiquitous in many physical applications, like for instance in the modelling of gas dynamics with gravity and friction terms, multiphase ows or radiative transfer. In the case of sti relaxation source terms and under compatibility conditions with the convective terms, see [START_REF] Berthon | Late-time/sti-relaxation asymptotic-preserving approximations of hyperbolic equations[END_REF] for more details, such systems may have some typical asymptotic (long-time) behaviors and degenerate into parabolic type equations. The concept of asymptotic-preserving schemes has been introduced in Jin [START_REF] Jin | Ecient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF],

Gosse and Toscani [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] and the aim is to preserve this asymptotic behavior at the numerical level. This means in particular that the consistency and stability of the discretization should be uniform with respect to the sti parameters, or equivalently they should be preserved also in the asymptotic limit. There is an extensive literature available on this topic and the design of rst-order asymptotic-preserving schemes in one space dimension is now well understood, see the recent book [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF], the references therein, and without any attempt to be exhaustive and for instance [3], [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]. . . However, turning now to the multi-dimensional case with high-order accuracy, there is a large gap and only a few numerical techniques have been developed up to our knowledge. As far as the multi-dimensional case is concerned, the 1D algorithms can be easily extended on cartesian grids and admissible meshes, but it is much more dicult on unstructured meshes since most schemes loose consistency in the asymptotic limit in this case. Regarding the high-order extension, the main challenges are to preserve the set of admissible states and the high-order accuracy in the asymptotic limit at the same time. We refer for instance the reader to [START_REF] Blachère | An admissibility and asymptoticpreserving scheme for systems of conservation laws with source term on 2D unstructured meshes[END_REF]2] and the references therein for rst-order schemes on unstructured meshes, and to [START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2D unstructured meshes with high-order MOOD reconstruction[END_REF][START_REF] Chalons | High-order asymptotic-preserving schemes for linear systems: Application to the goldstein-taylor equations[END_REF] and the references therein for high-order extensions. Note that the literature is now extremely large on the topic and that the proposed numerical strategies may also depend on the type of equations under consideration. Regarding kinetic equations for instance, we can also quote [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with sti sources[END_REF][START_REF] Latte | A High-Order Asymptotic-Preserving Scheme for Kinetic Equations Using Projective Integration[END_REF][START_REF] Boscarino | High-order asymptoticpreserving methods for fully nonlinear relaxation problems[END_REF][START_REF] Boscarino | Implicit-explicit Runge-Kutta Reference HLL-θ-P3 HLL-DLP-AP-P3 (d) Fine mesh with P3 schemes Figure 6: Density results for the blast at time t = 0.1 on the line (0.5, 0.4) -(0.5, 1.0) schemes for hyperbolic systems and kinetic equations in the diusion limit[END_REF]. . . It is therefore the purpose of this paper to design a class of high-order and asymptotic-preserving schemes on unstructured meshes for the numerical approximation of nonlinear systems with sti source terms. The present work extends in some sense the recent contribution [START_REF] Chalons | High-order asymptotic-preserving schemes for linear systems: Application to the goldstein-taylor equations[END_REF] devoted to highorder and asymptotic-preserving schemes in one space dimension and for linear systems. The proposed strategy is very simple and consists of a mild modication of the usual HLL scheme originally proposed in [START_REF] Harten | On upstream dierencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] and such that the consistency error of the scheme stays uniform with respect to the sti parameters. As we will see, this approach can be understood as a numerical diusion reduction technique and is especially well-adapted to the high-order extensions on unstructured meshes.

Throughout this paper, we will consider as a typical example the 2D isentropic Euler model with friction given by

∂ t U + div(F (U )) = σ(U )S(U ), (1)
where t usually denotes the time variable, x the location in space, and where we have set

U =   ρ ρu ρv   , F (U ) =   ρu ρv ρu 2 + p(ρ) ρuv ρuv ρv 2 + p(ρ)   , S(U ) =   0 -ρu -ρv   . (2)
In the following, we will also use the notation V = (u, v) for the velocity vector. The friction coecient σ is such that σ(U ) = κ(ρ) > 0 and the pressure law is assumed to satisfy p (ρ) > 0, so that the convective system is hyperbolic with eigenvalues given by V. nc, V. n and V. n + c, for any vector n ∈ R 2 and where the sound speed c is given by c = p (ρ). Recall that the characteristic eld associated with the two extreme eigenvalues are genuinely nonlinear, while the last one is linearly degenerate. The set of admissible states of this model is

A = {U = (ρ, ρV) ∈ R 3 , ρ > 0}.
It was proved in [4] that when κt → +∞, this system degenerates into the following diusion equation

∂ t ρ -div p (ρ) κ(ρ) ∇ρ = 0, (3)
with convergence speeds given for p ≥ 2 by

∂ β x ρ L p = O(1 + σt) -(1-1/p)-β/2 ∂ β x ∂ t ρ L p = O(1 + σt) -(1-1/p)-β/2-1/2 ∂ β x ρu L p = O(1 + σt) -(1-1/p)-β/2-1/2 ∂ β x ∂ t ρu L p = O(1 + σt) -(1-1/p)-β/2-1 (4)
The outline of the paper is as follows. In the next Section, we rst explain why the classical HLL scheme fails in preserving the asymptitc limit (3) by clearly showing the misleading term, and we propose a simple correction in multiple space dimensions. Section 3 discusses the extension to high-order accuracy and Section 4 illustrates the behaviour of our scheme by considering several numerical experiments.

2. First-order numerical scheme

Basic notations

Let us rst introduce some classical notations. We suppose that the computational domain Ω ⊂ R 2 is covered by N polygonal cells K. We consider e, a face of the cell K, and we suppose that the following admissibility assumptions are satised:

either there exists a single polygonal cell L such that e = K ∩ L = ∅.

In this case, e can be either a vertex or a single face of the mesh, and in the case of a single face, we note e = K ∩ L, either e ⊂ ∂Ω and in this case we will use the same notation e = K ∩ L where L is assumed to be a ghost cell to impose boundary conditions.

We note |K| the area of cell K, E K the set of interfaces of K, |e| the length of the interface of K and L (e = K ∩ L, L being the neighbour of K by the edge e), and n K,e the unit normal vector to e = K ∩ L pointing out of K.

At last, the perimeter p K of a cell K is dened by p K = e=K∩L∈E K |e|, and

δ K = |K| p K
is the so-called space step of the cell K.

The classical HLL scheme

Before going into the details, we rst recall the general framework of rst-order explicit nite volume scheme for (1) in several space dimensions.

We consider that for all n and K, U n K represents an approximation of the average value of the exact solution at time t n and on the volume K. Invoking the rotational invariance of the Euler equations, the standard nite volume approach to update the solution from time t n to time t n+1 writes

U n+1 K = U n K - ∆t |K| e=E K |e|F e • n K,e + ∆tσ K S K . (5)
In this scheme, σ K = σ(U n K ) and S K = S(U n K ) approximate the average of σ and S on K, and F e is a numerical ux which approximates the time integral of the exact ux at the interface e. An instance of such a numerical ux is given by classical HLL approximate Riemann solvers [START_REF] Harten | On upstream dierencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] such as the Rusanov ux leading to

F e • n K,e = 1 2 F (U K ) + F (U L ) • n K,e - b e 2 U L -U K , (6)
where b e > 0 denotes an upper bound of the absolute value of all wave speeds at interface e. We refer the reader to [START_REF] Harten | On upstream dierencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] for more details.

Using now the divergence formula

e∈E K |e|n K,e = 0 ⇐⇒ e∈E K |e|F K • n K,e = 0, (7)
with F K = F(U K , U K ) and to be plugged into (5), we get

U n+1 K = U n K - ∆t |K| e∈E K |e|(F e -F K ) • n K,e + ∆tσ K S K .
Finally, introducing the parameter w K,e = |e| p K > 0 such that

e∈E K w K,e = 1,
the scheme can be equivalently recast as

U n+1 K = e∈E K w K,e U n K - ∆t δ K (F e -F K ) • n K,e + ∆tσ K S K ,
which can thus be understood as a convex combination of quasi-1D schemes.

As far as the denition of the time step ∆t, we impose the following natural CFL condition,

∆t ≤ min K,e δ K 4b e + σ K δ K , (8)
See for instance [START_REF] Blachère | An admissibility and asymptoticpreserving scheme for systems of conservation laws with source term on 2D unstructured meshes[END_REF] in the homogeneous case (σ K = 0) and [START_REF] Chalons | High-order asymptotic-preserving schemes for linear systems: Application to the goldstein-taylor equations[END_REF] for a rigorous proof of stability in the linear non homogeneous case.

Failure of the classical HLL scheme

In this section, we aim at proving that the classical HLL scheme fails in preserving the asymptotic behavior of the solutions of (1), which is given by (3) when κt → +∞. A nice way to do that is to reproduce the limit κt → +∞ by rst considering the diusion scaling given by ∆t ← ∆t/ε and κ ← κ/ε in (5) and then letting ε tend to zero. In this section, we assume for the sake of simplicity that κ is constant. After easy manipulations, the scheme (5) rst writes:

ρ n+1 K = ρ n K - ∆t ε|K| e∈E K |e| 2 (ρu) n K + (ρu) n L (ρv) n K + (ρv) n L • n K,e -b e ρ n L -ρ n K , (ρu) n+1 K = 1 - κ∆t ε 2 (ρu) n K - ∆t ε|K| e∈E K |e| 2 (ρu 2 + p) n K + (ρu 2 + p) n L (ρuv) n K + (ρuv) n L • n K,e -b e (ρu) n L -(ρu) n K , (ρv) n+1 K = 1 - κ∆t ε 2 (ρv) n K - ∆t ε|K| e∈E K |e| 2 (ρuv) n K + (ρuv) n L (ρv 2 + p) n K + (ρv 2 + p) n L • n K,e -b e (ρv) n L -(ρv) n K .
Then, considering the Chapmann-Enskog expansions

(ρV) n K = (ρV) n,0 K + ε(ρV) n,1 K + O(ε 2 ), a rst identication of the terms in ε -2 yields (ρu) n,0 K = (ρv) n,0 K = 0. (9)
Using these equalities, an identication of the terms in ε -1 gives

(ρV) n,1 K = - 1 κ (∇ h p) K , (10)
where (∇ h p) K is an approximation of the pressure gradient in cell K and given by

(∇ h p) K = 1 |K| e∈E K |e| p n K + p n L 2 n K,e . (11)
Considering now the evolution equation on the density ρ, we easily get

ρ n+1 K = ρ n K - ∆t |K| e∈E K |e| 2 (ρu) n,1 K + (ρu) n,1 L (ρv) n,1 K + (ρv) n,1 L • n K,e + O(ε)+ + ∆t |K| e∈E K |e| 2 b e ρ n L -ρ n K ε which, using (10) becomes ρ n+1 K = ρ n K + ∆t |K| e∈E K |e| κ 1 2 (∇ h p) K + 1 2 (∇ h p) L • n K,e + O(ε)+ + ∆t |K| e∈E K |e| 2 b e ρ n L -ρ n K ε .
It is thus clear that the last term of the right-hand side prevents the limit scheme from being a consistent approximation of the limit equation (3), except of course if the mesh size is small compared to ε in order to compensate the 1/ε factor. Therefore, the classical HLL scheme ( 5)-( 6) is not asymptotic-preserving. However, it is important to observe that this wrong asymptotic behaviour comes from the diusion term of the HLL numerical ux. This remark motivates the denition of very simple scheme satisfying the asymptotic preserving and given in the next section.

A very simple asymptotically consistent correction

The results of the last section allow us to see that the problem comes from the numerical diusion in the density equation. The convergence of this term towards the diusion regime is not fast enough to be negligible in comparison with the terms that drive the behavior of the system of equations ( 1)-(2). In order to x this problem and to obtain the asymptotic-preserving property, we follow [START_REF] Chalons | High-order asymptotic-preserving schemes for linear systems: Application to the goldstein-taylor equations[END_REF] and [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] and propose to replace the classical HLL ux (6) at each interface by

F e • n K,e = 1 2 F (U K ) + F (U L ) • n K,e - b e 2 θ e U L -U K , (12)
where θ e = (θ e , 1, 1) T , which means that only the rst equation on the density ρ is concerned by the modication, and θ e = O(ε 1+α ) with α > 0 in the limit ε → 0 in order to recover the asymptotic-preserving property by forcing the numerical diusion to be negligible in the rst equation close to the diusion regime. Indeed, the limit scheme writes in this case (still assuming that κ is constant)

ρ n+1 K = ρ n K + ∆t |K| e∈E K |e| κ 1 2 (∇ h p) K + 1 2 (∇ h p) L • n K,e ,
which is a consistent approximation of the limit equation (3), without any specic requirement on the mesh. Note that the asymptotic-preserving property is also valid when κ is not constant but the pressure gradient in the limit scheme are of course weighted in this case. Therefore, the adjustment proposed here provides a numerical scheme that naturally preserves the asymptotic in the diusion regime.

At last, notice that at this stage the relation θ e = O(ε 1+α ) does not provide us with an explicit denition of θ e , especially because of the use the dimensionless form of the equations and therefore the presence of ε. This motivates to check the asymptotic-preserving property by comparing the continuous convergence speeds given in (4) with their discrete counterparts.

With this in mind, we rewrite the proposed scheme in the following form

ρ n+1 K -ρ n K ∆t = - 1 |K| e∈E K |e| 2 (ρu) n K + (ρu) n L (ρv) n K + (ρv) n L • n K,e + 1 |K| e∈E K |e| 2 b e θ e ρ n L -ρ n K , (13)
(ρu

) n+1 K -(ρu) n K ∆t = -κρu n K - 1 |K| e∈E K |e| 2 (ρu 2 ) n K + (ρu 2 ) n L (ρuv) n K + (ρuv) n L • n K,e - 1 2 (∇ h p) K + (∇ h p) L x + 1 |K| e∈E K |e| 2 b e (ρu) n L -(ρu) n K , ( 14
) (ρv) n+1 K -(ρv) n K ∆t = -κρv n K - 1 |K| e∈E K |e| 2 (ρuv) n K + (ρuv) n L (ρv 2 ) n K + (ρv 2 ) n L • n K,e - 1 2 (∇ h p) K + (∇ h p) L y + 1 |K| e∈E K |e| 2 b e (ρv) n L -(ρv) n K , (15)
where we have used the notations

(∇ h p) K + (∇ h p) L =   (∇ h p) K + (∇ h p) L x (∇ h p) K + (∇ h p) L y   .
Using standard Taylor expansions combined with the convergence speeds (4) in the L ∞ norm (see again [4]), we have

(ρu) n+1 K -(ρu) n K ∆t ≈ ∂ t (ρu)(t n , x K ) = O(1 + κt) -2 , ρu n K = O(1 + κt) -3/2 , 1 |K| e∈E K |e| 2 (ρu 2 ) n K + (ρu 2 ) n L (ρuv) n K + (ρuv) n L • n K,e ≈ div(ρV ⊗ V ) x (t n , x K ) = O(1 + κt) -2 , 1 2 (∇ h p) K + (∇ h p) L x ≈ ∂ x p(t n , x K ) = O(1 + κt) -3/2 , 1 |K| e∈E K |e| 2 b e (ρu) n L -(ρu) n K ≈ div( ∇ρu) = O(1 + κt) -2 ,
where ∇ is a rst-order operator which corresponds to a gradient twisted in the direction KL instead of n K,e on each interface.

Once these developments are completed, the dominant term can be identied

to obtain ρu = - 1 κ ∂ x p. ( 16
)
The same strategy applies on ρv to get

ρv = - 1 κ ∂ y p. (17)
Considering now the mass conservation equation ( 13), we get

ρ n+1 K -ρ n K ∆t ≈ ∂ t ρ(t n , x K ) = O(1 + κt) -3/2 , 1 |K| e∈E K |e| 2 (ρu) n K + (ρu) n L (ρv) n K + (ρv) n L • n K,e ≈ div(ρV )(t n , x K ) = O(1 + κt) -3/2 , 1 |K| e∈E K |e| 2 b e θ e ρ n L -ρ n K ≈ div(θ e ∇ρ),
so that a denition of θ e such that θ e = O(1 + κt) -k for a given k > 0 and the leading terms would give

∂ t ρ = -div(ρV),
and hence, using ( 16)-( 17) the expected diusive limit, namely

∂ t ρ = div 1 κ ∇p .
In practice, we will set

θ e = 1 1 + κt ,
in order to keep as much numerical diusion as possible without modifying the convergence speed. Indeed, choosing k ∈ (0, 1) would still give an asymptotic preserving scheme but with a wrong speed of convergence towards the limit.

Note that even if the limit scheme is self-imposed and therefore any target scheme cannot be chosen, it is remarkable that it is directly consistent even on unstructured meshes. We will also see that it is very easy to extend to obtain a uniformly high order scheme.

However, it is worth noting two minor disadvantages in practice. The rst one is that the CFL condition ( 8) is in general more restrictive CFL in general compared to the unmodied scheme. The second one is that due to the denitions of (∇ h p) K , the stencil of the limit scheme is wider than usual and can cause checkerboard eects in a regime very close to diusion.

Indeed, it is easily seen in 1D where the diusion operator is approximated by

ρ i+2 -2ρ i + ρ i-2 ∆x 2 .
Hopefully, the checkerboard eect only happens on specic meshes when κ is so large that the diusion limit is reached instantly. Therefore, it is seldom a real issue.

High-order numerical scheme

This section is devoted to the high-order extension of the numerical scheme. The particular form of correction makes it possible to consider doing so using most conventional techniques. We will focus here on the MOOD method but keeping an arbitrary order. In this context, the high order is achieved through a polynomial reconstruction in each cell. The reconstructed polynomial can then be expressed as follows, considering a generic polynomial degree d on cell K

ŨK (x; d) = U K + d |α|=1 R α (x -c) α - 1 |K| K (x -c) α dx ( 18
)
where c is the centroid of the cell K, R α are the polynomial coecients, α = (α 1 , α 2 ) ∈ N 2 is a multi-index with |α| = α 1 + α 2 , and (xc) α = (x 1c 1 ) α 1 + (x 2c 2 ) α 2 , where x = (x 1 , x 2 ) and c = (c 1 , c 2 ). Notice that this form of polynomial insures that the following property holds by construction

K ŨK (x; d)dx = U K .
This property is indeed mandatory to obtain eective orders greater than 3.

In order to determine the coecients R α , a standard least-square approximation is used on a given compact stencil associated with cell K, and the reader is referred for instance to [START_REF] Diot | La méthode MOOD Multi-dimensional Optimal Order Detection: la première approche a posteriori aux méthodes volumes -nis d'ordre très élevé[END_REF] for more details.

Once the reconstructed polynomials are dened on each cell, the following generic high-order (in space) nite volume formula is considered

U n+1 K = U n K - ∆t |K| e∈E K |e| R r=1 ξ r F e,r • n K,e + ∆tσ K S K , (19)
with

F e,r • n K,e = 1 2 F (U e,r K ) + F (U e,r L ) • n K,e - b e,r 2 θ e,r U e,r L -U e,r K . (20)
In these formulas, U e,r K and U e,r L represent the high-order approximations of U obtained by evaluating the in-cell reconstructed polynomials at quadratures points q e,r LR located on the edge e. The ξ r are the R quadrature weights.

Obviously, the quadrature formula has to be exact for polynomials up to degree d. As far as σ K and S K are concerned, it is dened by the average of σ and S on K using the appropriate quadrature formula. At last, the highorder accuracy in time, is obtained by using a standard Strong Stability Preserving Runge-Kutta method (see [START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF] for instance).

This results in a high order scheme that can easily be implemented generically by choosing appropriate quadrature rules and suitable reconstruction stencils. It remains to be decided which criteria will be used to correct the numerical results, particularly in the vicinity of discontinuities, in order to obtain a satisfactory numerical approximation. The concept of the MOOD method is to test the prediction of the high order numerical scheme at each time step with respect to these criteria and to a posteriori correct the cells that do not validate them by locally recalculating with a lower order scheme (possibly the rst order parachute scheme which must therefore verify the criteria). The validation criteria will be recalled for each numerical simulation, but it should be noted that it is always necessary to use at least the PAD (Physical Admissibility Detector) criterion. As its name suggests, it checks whether the prediction belongs to the set of physically admissible states and is therefore essential, otherwise the code may crash.

It now remains to check the asymptotic-preserving character of this highorder extension. In order to do so, one can reproduce exactly the same arguments as in Sections 2.3 and 2.4. More precisely, using the diusive scaling and focusing for instance on the rst-order accuracy in time (but high-order in space) for the sake of simplicity, ( 19)-( 20) can be recast as

ρ n+1 K = ρ n K - ∆t ε|K| e∈E K R r=1 |e| 2 ξ r (
) n+1 K = 1 - κ∆t ε 2 (ρv) n K - ∆t ε|K| e∈E K R r=1 |e| 2 ξ r (ρuv) n,e,r K + (ρuv) n,e,r L (ρv 2 + p) n,e,r K + (ρv 2 + p)
n,e,r L

• n K,e -b e θ e,r (ρv) n,e,r L -(ρv) n,e,r K .

Then, considering again the Chapmann-Enskog expansions

(ρV) n K = (ρV) n,0 K + ε(ρV) n,1 K + O(ε 2 ),
by an identication of the terms in ε -2 and ε -1 we get (ρu) n,0 K = (ρv) n,0 K = 0 and (ρV) n,1 K = -1 κ (∇ h p) K , where (∇ h p) K still approximates the pressure gradient in cell K but is now given by

(∇ h p) K = 1 |K| e∈E K R r=1 |e|ξ r p n,e,r K + p n,e,r L 2 n K , e. (21)
Finally, considering the evolution equation on the density ρ, the limit scheme is again given by (assuming that κ is a constant)

ρ n+1 K = ρ n K + ∆t |K| e∈E K |e| κ 1 2 (∇ h p) K + 1 2 (∇ h p) L • n K,e ,
which is a consistent approximation of the limit equation (3), again without any specic requirement on the mesh. Furthermore, this approximation is naturally an high-order one due to the properties of the reconstruction polynomials. Indeed, the MOOD nature of the extension is not required to obtain this property and therefore, as mentioned earlier, other techniques based on polynomial reconstructions (e.g. (W)ENO) may be used to extend the scheme to high-order.

A uniformly high order numerical scheme was thus obtained, which preserves both the diusion asymptotic and the set of admissible states. This scheme is relatively simple to obtain through a correction that does not require any particular implementing eort. The only price to pay is a CFL condition which is slightly more restrictive than the one of the scheme described in [START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2D unstructured meshes with high-order MOOD reconstruction[END_REF]. However, this drawback may be overcame by using a slightly coarser mesh due to its uniformly high-order nature (the scheme designed in [START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2D unstructured meshes with high-order MOOD reconstruction[END_REF] loses its high-order nature when the friction is large).

Numerical results

In this section we numerically check the behaviour of our scheme in various situations. Besides, we also compare the new HLL-θ scheme ((5) with ux ( 12)) to the classical HLL (scheme ( 5) with ux (6), i.e. HLL-θ with θ e = 1 everywhere) and to the HLL-DLP-AP schemes (from [START_REF] Blachère | An admissibility and asymptoticpreserving scheme for systems of conservation laws with source term on 2D unstructured meshes[END_REF][START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2D unstructured meshes with high-order MOOD reconstruction[END_REF]).

For the high-order schemes the proper time integrators are used, more precisely SSPRK(1, 1), SSPRK(2, 2) and SSPRK(5, 4) using the notations of [START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF] for the P0, P1 and P3 schemes. Where the Pn notation denotes a reconstruction of degree n in (18) using the appropriate quadrature formulas [START_REF] Diot | La méthode MOOD Multi-dimensional Optimal Order Detection: la première approche a posteriori aux méthodes volumes -nis d'ordre très élevé[END_REF].

Let us notice, that to have a proper convergence with the HLL-θ-P3 scheme, we impose θ e = 1 in the ux of the boundary interfaces.

In the MOOD loop, the PAD (Physical Admissibility Detector) criterion is always enabled to check if the solution lies in the set of admissible state A. To avoid spurious oscillations and to keep smooth extrema the DMP (Discrete Maximum Principle) and u2 criteria may be enabled with the P 3 reconstruction. For more details about criteria, the reader is referred to [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Clain | The MOOD method for the non-conservative shallow-water system[END_REF].

Convergence to a manufactured solution

The aim of the rst numerical experiment is to check the convergence with respect to the mesh size. In order to dene an exact solution of our system (2), we consider a manufactured solution of the form ρ(x, y, t) = exp(-(x + y)),

(ρu)(x, y, t) = (ρv)(x, y, t) = exp(-t) k , p(x, y, t) = ρ(x, y, t) γ , (22)
which corresponds to the following source term

σ(x, y, t) = -1 (ρu) ∂ t (ρu) + ∂ x p -2 (ρu) 2 ∂ x (ρ) ρ 2 . ( 23
)
With this exact solution, the tests are run with a nal time t f = 1.

Dirichlet boundary conditions are imposed with the exact solution on each sides of the unit square. The rst mesh is composed of triangles created with Gmsh (see [START_REF] Geuzaine | Gmsh: A 3-D nite element mesh generator with built-in pre-and post-processing facilities[END_REF]), while all the ner meshes used are renements of the rst one.

In Figure 1, the L 2 -norm of the error on ρ is drawn regarding to the number of cells in the mesh. The left plot corresponds to the choice k = 20 in ( 22) and the right one to k = 640. The new HLL-θ scheme [START_REF] Diot | La méthode MOOD Multi-dimensional Optimal Order Detection: la première approche a posteriori aux méthodes volumes -nis d'ordre très élevé[END_REF] with various polynomial reconstructions (P0, P1, P3) is compared to the classical HLL scheme [START_REF] Blachère | An admissibility and asymptotic preserving scheme for systems of conservation laws with source term on 2D unstructured meshes with high-order MOOD reconstruction[END_REF] and the HLL-DLP-AP scheme developed in [START_REF] Blachère | An admissibility and asymptoticpreserving scheme for systems of conservation laws with source term on 2D unstructured meshes[END_REF].

In the MOOD loop, only the PAD criterion is used as the solution is very smooth and no spurious oscillations are created with the schemes. Figure 1: Rate of convergence regarding to the mesh with a manufactured solution [START_REF] Toro | Riemann solvers and numerical methods for uid dynamics[END_REF] Figure 1 shows a good agreement with the theoretical convergence rates for all the schemes. Besides, the new schemes oer a better convergence rate especially with large source terms.

Convergence to the diusive limit

The next test case illustrates the behaviour of the numerical schemes with respect to the long time behaviour and the corresponding convergence speeds (4). With this in mind, we plot on Figure 2 the evolution of ρ(t) L p and (ρu)(t) L p for p = 2 and p = ∞ with respect to 1 + σt.

The initial condition is a Gaussian curve, namely

ρ 0 (x, y) = exp -(x -50) 2 -(y -50) 2 + 1 V 0 (x, y) = 0. (24)
An unstructured mesh of [0; 100] 2 composed of 6.6 • 10 3 triangles and a space step ∆x = min K δ K = 1.7 • 10 -1 are used. Neumann boundary conditions The results are obtained with various schemes (HLL-θ, HLL and HLL-DLP-AP). The high-order schemes only use the PAD criterion in the MOOD loop. We observe that the proposed asymptotic-preserving schemes fully respect the convergence rates whereas the classical HLL scheme gives the correct convergence rate when the high-order reconstructions are used.

Blast with friction

The last numerical test shows the behaviour of our scheme in a more complex conguration. In the unit square two zones are set: one inside a circle of radius 0.3 centered at (0.5, 0.4) with the following initial condition ρ 0 = 5, V 0 = 0, and one outside the circle with:

ρ 0 = 1, V 0 = 0.
All sides of the square are considered as walls and the source term is controlled by: σ(U ) = 5 ρ 5

3

The mesh created with Gmsh is composed of 4 • 10 4 triangles for a space step of 8 • 10 -4 . Finally, the MOOD loop used for the high-order schemes use the PAD, DMP and u2 criteria [START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF].

The results are presented in Figures 4 and5 In order to make the comparison, reference solutions obtained on ne meshes of 6 • 10 5 cells (∆x = 2 • 10 -4 ) and 2.7 • 10 6 cells (∆x = 1 • 10 -4 ) with the HLL-θ-P0 and HLL-θ-P3 schemes, are presented in Figure 3.

For the P0 schemes in Figure 4 the precision is a bit better with the HLL-DLP-AP scheme. Whereas, the HLL-θ-P3 scheme gives a better description of the central structure in Figure 5. This behaviour is expected as HLL-θ is a true high-order scheme while HLL-DLP-AP only use high-order polynomial reconstruction and get down to a rst-order scheme when targeting the limit.

The total run to the nal time t f = 1 took 2.2•10 4 iterations and 1h with the HLL-θ-P0 scheme with an average time step of ∆t = 4.5 • 10 -5 whereas the HLL-DLP-AP-P0 took 5 • 10 3 iterations and 2h with ∆t = 2 • 10 -4 . This dierence can be explained as the HLL-DLP-AP scheme use the classical hyperbolic CFL but needs to compute nonlinear coecients (ν in [START_REF] Blachère | An admissibility and asymptoticpreserving scheme for systems of conservation laws with source term on 2D unstructured meshes[END_REF]) at each time step.

For the P3 schemes the dierence between the two computational times is greater (10h and 30h) as there are more coecients to compute for the HLL-DLP-AP scheme.

In Figure 6, a comparison is made with a reference solution obtained as in [START_REF] Toro | Riemann solvers and numerical methods for uid dynamics[END_REF]. This reference solution is computed with a 1D code using the HLL-θ scheme from [START_REF] Chalons | High-order asymptotic-preserving schemes for linear systems: Application to the goldstein-taylor equations[END_REF] adapted to the isentropic Euler model with friction (2). In order to take account of the cylindrical symmetry of the test case, we added the geometrical source term : S g (U ) = -1 r ρ ρu .

with r the distance from the center.

  ) n,e,r L • n K,e -b e θ e,r (ρu) n,e,r L -(ρu) n,e,r K , (ρv

  k = 640, σ ∈ [64; 2600]

Figure 2 :

 2 Figure 2: Decays in time of ρ L p and ρu L p

  at two nal times (0.35 on the left and 0.75 on the right). The HLL-θ-P0 and HLL-DLP-AP-P0 are used in Figure 4 while the Figure 5 shows the results of the HLL-θ-P3 and HLL-DLP-AP-P3. The density values lie in [0.5; 5.1] at time 0.35 and [1; 3.5] at time 0.75.

10 5 Figure 3 :Figure 4 :

 534 Figure 3: Density results on ne meshes for the blast at various times