Mayeul Arminjon
email: mayeul.arminjon@3sr-grenoble.fr

Interstellar radiation as a Maxwell field: improved numerical scheme and application to the spectral energy density

The existing models of the interstellar radiation field (ISRF) do not produce a Maxwell field. Here, the recent model of the ISRF as a Maxwell field is improved by considering separately the different frequencies at the stage of the fitting. Using this improved procedure: (i) It is checked in detail that the model does predict extremely high values of the spectral energy density (SED) on the axis of a galaxy, that however decrease very rapidly when ρ, the distance to the axis, is increased from zero. (ii) The difference between the SED values (with ρ = 1 kpc or 8 kpc), as predicted either by this model or by a recent radiation transfer model, is reduced significantly. (iii) The slower decrease of the SED with increasing altitude z, as compared with the radiation transfer model, is confirmed. We also calculate the evolutions of the SED at large ρ. We interpret these evolutions by determining asymptotic expansions of the SED at large z, and also ones at large ρ.

Introduction

The interstellar radiation field (ISRF) in a galaxy is an electromagnetic (EM) field in a very high vacuum, hence it should be a solution of the Maxwell equations. However, the existing models for the ISRF do not take into account the full EM field with its six components coupled through the Maxwell equations. Consider, for example, the model of Chi & Wolfendale [START_REF] Chi | The interstellar radiation field: a datum for cosmic ray physics[END_REF]. It assumes an axisymmetric distribution of the volume emissivities j i (λ, ρ, z) of four stellar components (i) (i = 1, ..., 4): j i decreases exponentially with both the distance ρ to the galactic axis and the altitude z over the galactic central disk. The contribution of component (i) to the energy density of the ISRF at some position (ρ , z ) and wavelength λ is obtained by integrating j i (λ, ρ, z)g/l 2 over the whole galactic volume. Here l is the distance between the studied position and the running point in the galactic volume; g describes the dust absorption and is obtained by integrating the visual extinction per unit path length over the linear path joining the studied position and the running point in the galactic volume. Other models, e.g. by Mathis, Mezger and Panagia [START_REF] Mathis | Interstellar radiation field and dust temperatures in the diffuse interstellar matter and in giant molecular clouds[END_REF], Gordon et al. [START_REF] Gordon | The DIRTY model. I. Monte Carlo radiative transfer through dust[END_REF], Robitaille [START_REF] Robitaille | HYPERION: an open-source parallelized three-dimensional dust continuum radiative transfer code[END_REF], Popescu et al. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF], are based on similar principles: all of these models consider quantities such as the stellar emissivity and luminosity, and the dust opacity, and they evolve the light intensity emitted by the stars by taking into account (in addition to the distance) the radiative transfer, in particular by dust absorption/ reemission. Clearly, those models do not produce an EM field, hence even less one that would be a solution of the Maxwell equations.

In a recent work [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF], we proposed a model applicable to the relevant ideal case of an axisymmetric galaxy, and that provides for the ISRF such an exact solution of the Maxwell equations -a feature which, as discussed above, and to the best of our knowledge, appears to be fully new. This is indeed needed to study the relevance of a possible candidate for dark matter that emerges [START_REF] Arminjon | On the equations of electrodynamics in a flat or curved spacetime and a possible interaction energy[END_REF] from an alternative, scalar theory of gravity. However, it is also of astrophysical interest independently of the latter, since, as we noted, the ISRF must be an exact Maxwell field and this condition is not fulfilled by the existing models. As a step in checking the model proposed in Ref. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF], its application to predict the variation of the spectral energy density (SED) in our Galaxy has been subjected to a first test [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. To this purpose, the model has been adjusted by asking that the SED predicted for our local position in the Galaxy coincide with the SED determined from spatial missions by Henry, Anderson & Fastie [START_REF] Henry | Far-ultraviolet studies. vii. The spectrum and latitude dependence of the local interstellar radiation field[END_REF], Arendt et al. [START_REF] Arendt | The COBE Diffuse Infrared Background Experiment Search for the Cosmic Infrared Background. III. Separation of Galactic Emission from the Infrared Sky Brightness[END_REF], Finkbeiner et al. [START_REF] Finkbeiner | Extrapolation of Galactic Dust Emission at 100 Microns to Cosmic Microwave Background Radiation Frequencies Using FIRAS[END_REF], and Porter & Strong [START_REF] Porter | A new estimate of the galactic interstellar radiation field between 0.1µm and 1000µm[END_REF]. It has been found in that most recent work [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF] that the spatial variation of the SED thus obtained with our model does not differ too much in magnitude from that predicted by the recent radiation transfer model of Ref. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF], but that the SED predicted by our model: (i) is extremely high on the axis of the Galaxy -i.e., on the axis of the axial symmetry that is assumed for the model of the Galaxy; (ii) has rather marked oscillations as function of the wavelength; and (iii) seems to decrease more slowly when the altitude z increases (or rather when |z| increases), as compared with the radiation transfer model.

The aim of this paper is to present an improved numerical scheme to operate that "Maxwell model of the ISRF", and to apply this improved scheme to check the findings (i)-(iii) above. Section 2 provides a summary of the model. Section 3 describes the improvement of the numerical scheme. In Sect. [START_REF] Robitaille | HYPERION: an open-source parallelized three-dimensional dust continuum radiative transfer code[END_REF], we check whether the model really predicts extremely high values of the SED on the axis of the Galaxy. Section 5 studies the spatial variation of the SED and compares it with results of the literature. In Sect. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF], asymptotic expansions are used to interpret the findings of the foregoing section. The Conclusion section 7 is followed by Appendix A, which discusses the relation between the discrete and continuous descriptions of the SED.

Short presentation of the model

This model has been presented in detail in Ref. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF]. An axisymmetric galaxy is modelled as a finite set of point-like "stars", the azimuthal distribution of which is uniform. Those points x i (i = 1, ..., i max ) are obtained by pseudorandom generation of their cylindrical coordinates ρ, φ, z with specific probability laws, ensuring that the distribution of ρ and z is approximately that valid for the star distribution in the galaxy considered, and that the set {x i } is approximately invariant under azimuthal rotations of any angle φ [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF]. In the present work, as in Refs. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF][START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF], 16 × 16 × 36 triplets (ρ, z, φ) were thus generated, so that i max = 9216, and the distribution of ρ and z is approximately that valid for the star distribution in the Milky Way.

The ISRF is also assumed axisymmetric, and thus depends only on ρ and z. Since we want to describe, not the field inside the sources and in their vicinity, but instead the smoothed-out field at the intragalactic scale, we search for a solution of the source-free Maxwell equations. In the axisymmetric case, any time-harmonic source-free Maxwell field is the sum of two Maxwell fields: (i) one deriving from a vector potential having just the axial component A z non-zero, with A z obeying the standard wave equation, and (ii) one deduced from a solution of the form (i) by EM duality [START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF]. We consider for simplicity a model ISRF that has a finite frequency spectrum (ω j ) j=1,...,Nω , hence we may apply the foregoing result to each among its timeharmonic components (j), and then just sum these components. Moreover, we envisage the ISRF as being indeed an EM radiation field, thus excluding from consideration the purely magnetic part of the interstellar EM field [START_REF] Beck | Magnetic fields in the Milky Way and in galaxies[END_REF]. Hence the ISRF is made of "totally propagating" EM waves, i.e., ones without any "evanescent" component [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF][START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]. Specifically, we assume that the two scalar potentials A j z and A j z that define the decomposition (i)-(ii) of each time-harmonic component (j), mentioned above, are themselves totally propagating. In that case, both A j z and A j z have the explicit form [START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]:

ψ ω j S j (t, ρ, z) = e -iω j t +K j -K j J 0 ρ K 2 j -k 2 e ik z S j (k) dk, (1)
with ω j the angular frequency, K j := ω j /c, J 0 the first-kind Bessel function of order 0, and where S j is some (generally complex) function of k ∈ [-K j , +K j ]. For a totally propagating, axisymmetric EM field, but otherwise fully general, the two potentials A j z and A j z may be different, i.e., may correspond with different "spectra" in Eq. ( 1), say S j and S j [START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF].

To determine these potentials, that is, to determine the spectrum functions S j , we use a sum of potentials emitted by the "stars". We assume that every "star", each at some point x i , contributes to the global potential A j z of a given frequency ω j (j = 1, ..., N ω ) by a spherically symmetric scalar wave of the same frequency ω j , whose emission center is its spatial position x i -in order that all the directions starting from the star be equivalent. Thus, consider time-harmonic spherically symmetric solutions of the wave equation that have a given angular frequency ω. It is easy to check by direct calculation that they can be either an outgoing wave, an ingoing wave, or the sum of an ingoing wave and an outgoing one, and that, up to an amplitude factor, the following is the only outgoing wave:

ψ ω (t, x) = e i(Kr-ωt) Kr , K := ω c , r := |x| . (2)
Clearly, only that outgoing solution is relevant here, given that the point-like "stars" must be indeed sources of radiation. 1 Thus, the contributions of the i-th star to the potentials A j z and A j z can differ only in amplitude, since both must be a multiple of

ψ x i ω j (t, x) := ψ ω j (t, x -x i ) = e i(K j r i -ω j t) K j r i , (3)
where

K j := ω j c , r i := |x -x i |.
But there is no apparent physical reason to affect different amplitudes to the contribution of the i-th star to A j z and to A j z , hence we assume both of them to be equal to ψ x i ω j . To determine the global potentials A j z and A j z (j = 1, ..., N ω ), that generate the axisymmetric model ISRF with a finite frequency spectrum (ω j ), the sum of the spherical potentials (3) emanating from the point stars is fitted to the form (1). As noted in Ref. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF], this is not assuming that the ISRF is indeed the sum of the radiation fields emitted by the different stars (which is not correct, due to the radiation transfers) -because (a) the equalities ( 4), [START_REF] Dieudonné | Calcul infinitésimal[END_REF] or [START_REF]Méthode de la phase stationnaire[END_REF] below are not exact equalities but ones in the sense of the least squares, and (b) nothing is really assumed regarding the EM field of the "star" itself, in particular we actually do not need to assume that it has the form (i)-(ii) above (e.g. the one corresponding with two equal potentials A i j z = A i j z = ψ x i ω j ).

In the previous works [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF][START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF], this fitting was done for all frequencies at once. That is, the following least-squares problem was considered:

Nω j=1 imax i=1 w j ψ x i ω j ∼ = Nω j=1 ψ ω j S j on G, (4)
where the sign ∼ = indicates that the equality is in the sense of the least squares (the arguments of the functions varying on some spatio-temporal grid G), and where the numbers w j > 0 are the weights affected to the different frequencies. In view of the axial symmetry, the spatial position x is restricted to the plane φ = 0, so x = x(ρ, z) and

G = {(t l , ρ m , z p ), 1 ≤ l ≤ N t , 1 ≤ m ≤ N ρ , 1 ≤ p ≤ N z }. ( 5
)
spherically symmetric solution of the wave equation that satisfies the Sommerfeld radiation condition [START_REF] Sommerfeld | Die Greensche Funktion der Schwingungsgleichung[END_REF]. However, the Sommerfeld condition aims precisely at selecting a boundary condition in order to find only "physical", i.e., outgoing solutions for the Helmholtz equation. The latter equation applies to general time-harmonic solutions of the wave equation.

In the spherically symmetric case, the time-harmonic solutions are easy to find and the outgoing solutions are immediate to recognize.

Since the contributions of the i-th star to A j z and to A j z have both been assumed to be equal to ψ x i ω j , there is no possibility to distinguish between A j z and A j z , either -whence ψ ω j S j = A j z = A j z on the r.h.s. of (4). The unknowns of the problem are the spectrum functions S j , j = 1, ..., N ω . We determine S j by the (generally complex) values

S nj := S j (k nj ) (n = 0, ..., N ), (6)
where

k nj = -K j + nδ j (n = 0, ..., N ), (7)
with δ j := 2K j /N , is a regular discretization of the interval [-K j , +K j ] for k in the integral [START_REF] Chi | The interstellar radiation field: a datum for cosmic ray physics[END_REF]. Calculating those integrals with the "Simpson

w j ψ x i ω j ∼ = Nω j=1 N n=0 f nj S nj on G, (8)
with

f nj (t, ρ, z) = a nj J 0 ρ K 2 j -k 2 nj exp [i (k nj z -ω j t)] . (9)
The S nj 's are the solved-for parameters in the least-squares problem [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. In Eq. ( 8), N must be a multiple of 3, and in Eq. ( 9) we have

a nj = (3/8) δ j (n = 0 or n = N ), (10)
a nj = 2 × (3/8) δ j (mod(n, 3) = 0 and n = 0 and n = N ), [START_REF] Finkbeiner | Extrapolation of Galactic Dust Emission at 100 Microns to Cosmic Microwave Background Radiation Frequencies Using FIRAS[END_REF]

a nj = 3 × (3/8) δ j otherwise. (12)
Part (i) of the decomposition of the model ISRF then obtains as follows [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF]:

E φ = B ρ = B z = 0, (13)
B φ (t, ρ, z) = N n=0 Nω j=1 R n J 1 ρ ω j ω 0 R n Re [F nj (t, z)] + O 1 N 4 , (14)
E ρ (t, ρ, z) = N n=0 Nω j=1 c 2 ω 0 k n R n J 1 ρ ω j ω 0 R n Re [F nj (t, z)] + O 1 N 4 , (15)
E z (t, ρ, z) = N n=0 Nω j=1 c 2 ω 0 k 2 n -ω 0 J 0 ρ ω j ω 0 R n Im [F nj (t, z)] + O 1 N 4 , ( 16
) with R n = K 2 0 -k 2 n and F nj (t, z) = ω j ω 0 2 a n exp i ω j ω 0 k n z -ω j t S nj . (17)
(Here k n and a n (0 ≤ n ≤ N ) are as k nj and a nj in Eqs. ( 7) and ( 10), replacing K j by K 0 = ω 0 c , with ω 0 some (arbitrary) reference frequency.) Since we assume A j z = A j z for the global potentials generating the model ISRF, part (ii) of its decomposition is deduced from the first part by the EM duality:

E = cB, B = -E/c. (18)
It follows from this and from (13) that the model ISRF, sum of these two parts, has the components ( 14)-( 16), and that the other components are just

E φ = cB φ , B ρ = -E ρ /c, B z = -E z /c. ( 19
)
3 Frequency-by-frequency fitting of the potentials Equation (4) may be split into the different frequencies (marked by the index j), simply by removing the sum on j from both sides of either equation. The same is true for Eq. [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. Naturally, also the weight w j may then be removed from the l.h.s., by entering the inverse 1/w j into the unknown spectrum function S j on the r.h.s. Equation ( 8) thus becomes imax i=1

ψ x i ω j ∼ = N n=0 f nj S nj on G (j = 1, ..., N ω ). ( 20
)
At this point, one notes that both ψ x i ω j [Eq. ( 3)] and f nj [Eq. [START_REF] Henry | Far-ultraviolet studies. vii. The spectrum and latitude dependence of the local interstellar radiation field[END_REF]] have the same dependence on time, exp(-iω j t), which we can hence remove also, to obtain a least-squares problem with merely the spatial variables ρ and z:

imax i=1 e iK j r i K j r i ∼ = N n=0 g nj S nj on G (j = 1, ..., N ω ), (21)
where

G = {(ρ m , z p ), 1 ≤ m ≤ N ρ , 1 ≤ p ≤ N z }
is the spatial grid, and

g nj (ρ, z) = a nj J 0 ρ K 2 j -k 2 nj exp (ik nj z) . (22)
The separation, into the different frequencies, of the fitting of the sum of the potentials emitted by the "stars", is consistent with the linearity of the wave equation and the Maxwell equations. Moreover, the elimination of the time variable from the fitting represents an appreciable gain in computing time.

We recall that, for the EM field in a galaxy, the arguments of the Bessel function J 0 and the angular exponential, e.g. in Eq. ( 22), have the huge magnitude |x| /λ ∼ 10 25 , which enforces us to use a precision better than quadruple precision in the computer programs, thus leading to slow calculations [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF]. Note that the "separate fitting", i.e. the least-squares problem [START_REF]Méthode de la phase stationnaire[END_REF], is not exactly equivalent to the "grouped fitting", i.e. the least-squares problem (8) (this will be confirmed by the numerical results below): the two are slightly different ways of adjusting the global potentials (1). 2 However, equations ( 13)-( 19) apply with the separate fitting as well -although the relevant values S nj are different. The separate fitting is more appropriate, because solutions corresponding with different frequencies behave independently in the Maxwell equations, and each frequency can be treated with more precision by considering it alone. Indeed a very important point is that, by switching to the separate fitting, we improve the situation regarding the "overfitting", i.e., we decrease the ratio R of the number of parameters N para to the number of data N data : now, for each value of the frequency index j, we have to solve the least-squares problem [START_REF]Méthode de la phase stationnaire[END_REF], with N para = N + 1 unknown parameters and N data = N ρ × N z data (the "data" are the values of the l.h.s. of ( 21) on the spatial grid G ). Whereas, with the formerly used grouped fitting, we had to solve just one least-squares problem [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF] with

N para = (N + 1) × N ω unknown parameters and N data = N t × N ρ × N z data.
On the other hand, through the processes of radiative transfer, there are indeed transfers of radiation intensity from some frequency domains to other ones, e.g. the interaction with dust leads to a transfer from higher to lower frequencies (see e.g. Fig. 3 in Ref. [START_REF] Gordon | The DIRTY model. I. Monte Carlo radiative transfer through dust[END_REF]). But these processes are not directly taken into account by the present model: not any more with the grouped fitting than with the separate fitting. They are indirectly taken into account through the adjustment of the energy density [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF], which we briefly recall now.

The time-averaged volumic energy density of an EM field having a finite set of frequencies, (ω j ) j=1,...,Nω , is given by [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]

U (x) := δW δV (x) = Nω j=1 u j (x), u j (x) := 1 4 6 q=1 α q C (q) j (x) 2 , ( 23
)
where the complex numbers C (q) j (x) (q = 1, ..., 6) are the coefficients in the expansion, in time-harmonic functions, of each among the six components of the EM field:

F (q) (t, x) = Re Nω j=1 C (q) j (x)e -iω j t (q = 1, ..., 6); (24)
and where α q = 0 for an electric field component, whereas α q = 0 c 2 for a magnetic field component (here 0 is the vacuum permittivity, with 0 = 1/(4π × 9 × 10 9 ) in SI units). For an axisymmetric EM field, it is enough to consider the plane φ = 0, thus x = x(ρ, z), and we have

C (q) j = C (q) j (ρ, z), u j = u j (ρ, z). (25)
Using in that case the decomposition (i)-(ii), the expressions of three among the C (q) j coefficients follow directly from the expressions ( 14)-( 16) of the corresponding components of the EM field [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. Moreover, in the special subcase [START_REF] Sommerfeld | Die Greensche Funktion der Schwingungsgleichung[END_REF] considered here, the other components are given by [START_REF] Majaess | Characteristics of the Galaxy according to Cepheids[END_REF], whence in the same way the three remaining C (q) j coefficients. Now note that, in the least-squares problem [START_REF]Méthode de la phase stationnaire[END_REF], that we use to determine the values S nj allowing to compute the EM field ( 14)-( 16) and [START_REF] Majaess | Characteristics of the Galaxy according to Cepheids[END_REF], no data relative to the intensity of the fields emitted by the point-like "stars" has been used until now. Hence, we may multiply the l.h.s. of ( 21) by some number ξ j > 0, thus obtaining now new values S nj = ξ j S nj (n = 0, ..., N ) as the solution of ( 21). 3 Therefore, to adjust the model, we determine the numbers ξ j (j = 1, ..., N ω ) so that the values u j (x loc ) of the SED for our local position x loc in the Galaxy and for the frequencies ω j , as given by Eq. ( 23), coincide with the measured values, as determined from spatial missions. We take the measured local values f x loc (λ j ) as plotted in Ref. [START_REF] Porter | A new estimate of the galactic interstellar radiation field between 0.1µm and 1000µm[END_REF] (see Appendix A), and we take ρ loc = 8 kpc and z loc = 0.02 kpc, see e.g. Ref. [START_REF] Majaess | Characteristics of the Galaxy according to Cepheids[END_REF]. The model thus adjusted then allows us to make predictions: in particular, predictions of the spatial variation of the SED in the Galaxy. Such predictions may then be compared with predictions of the mainstream models of the ISRF, which models are very different from the present model.

Results: maximum energy density

In the foregoing work [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF], the same adjustment just described was used in the framework of the "grouped fitting" (i.e. the least-squares problem ( 8)). A surprising result was that found for the values of the maximum of the energy density u j (x) in the Galaxy -thus, owing to the axial symmetry (25), for the values of

u jmax = Max{u j (ρ m , z p ); m = 1, ..., N ρ , p = 1, ..., N z }, (26)
found for the different spatial grids investigated, all having ρ varying regularly from ρ 0 = 0 to ρ max 10 kpc and z varying regularly from z 0 = 0 or z 0 = -z max to z max ≤ 1 kpc. 4 These maximum values, which are always found at ρ = 0, thus on the axis of symmetry, are extremely high for lower wavelengths λ j , with u jmax 10 27 eV/cm 3 . Moreover, the value of u j (ρ = 0, z) depends little on z in the domain investigated. These surprisingly high values occur in a larger or smaller range of wavelengths, depending on the settings of the calculation. Therefore, the question arises whether these extremely high values are a physical effect or a numerical artefact. However, the dependence 3 We might even affect different weights ξ ij to the radiations of frequency ω j emitted by the different stars (i), in order to account for different luminosities. However, given that our aim is to determine the spectra S j , each of which characterizes according to Eq. (1) the axisymmetric radiation of frequency ω j in the galaxy, we feel that this would not likely change the results very significantly. 4 Precisely: ρ 0 := ρ m=1 ; ρ max := ρ m=Nρ with, in this subsection, ρ max = 10 kpc×

Nρ-1

Nρ ; z 0 := z p=1 and z max := z p=Nz with, in this paper, z max = 1 kpc and z 0 = -z max .

on the settings is governed by the "amount of overfitting": less overfitting increases the range of the high values [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. This makes it plausible that the high values might be a true prediction of the model. We will now try to check whether this is indeed the case.

Robustness of the high values on the axis

In the present work based on the separate fitting (which, we argued, is more appropriate), we investigated rather systematically the question which we just asked. Since the influence of the spatial grid was found weak in the foregoing work [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF], only two grids were tried: an (N ρ = 10, ρ 0 = 0) × (N z = 21, z 0 = -1 kpc, z max = 1 kpc) grid (hereafter "rough grid"), and an (N ρ = 20, ρ 0 = 0)×(N z = 23, z 0 = -1 kpc, z max = 1 kpc) grid (hereafter "fine grid"). However, we investigated the influence of the fineness of the frequency mesh (N ω ) and the influence of the discretization number N quite in detail. [That integer N is used to compute the integrals over the wavenumber k, e.g. the integral (1) approximated to N n=0 f nj S nj , see Eq. ( 8).] The effect of choosing N ω = 23, N ω = 46, or N ω = 76, was studied simultaneously with the effect of choosing N = 12, or N = 24, 48, 96, 192, 384, and this was done for the two different grids.

Figures 1 to 4 show these effects. The most salient result is that the extremely high values of u jmax are now found with all calculations and in the whole domain of λ -except that on some curves, abrupt oscillations toward lower values of the energy density are present for some wavelengths. By looking at the set of these figures, it is manifest that such abrupt oscillations occur when an inappropriate choice of parameters is done: essentially, the discretization number N has to be large enough. (This is certainly expected, and this expectation is confirmed by the validation test in Ref. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF].) Indeed, for a given value of N ω , those oscillations are maximum for the lowest value of N in the trial (N = 12) and progressively disappear when N is increased. What is a "large enough" value of N is not strongly dependent of the fineness of the spatial grid (i.e., of whether the "rough" one or the "fine" one is used) and that of the frequency mesh (N ω ). However, when using the finest frequency mesh (N ω = 76) for the "rough" spatial grid (Fig. 3), increasing N does not allow us to eliminate the abrupt oscillations toward lower values: it even happens then, that increasing N from 192 to 384 actually deteriorates the u jmax = f (λ j ) curve. We interpret this as due to the fact that, when us-ing a rougher spatial grid G for the fitting, less data are provided (the values taken on the grid G by the l.h.s. of Eq. ( 21)) to determine the unknowns S nj on the r.h.s. of (21) -while, of course, increasing N increases the number of unknowns and thus asks for more data. On the other hand, it is seen that (for the relevant values of N , say N = 192 or N = 384, so that little or no oscillations are present), the levels of u jmax depend quite little on N ω i.e. on the fineness of the frequency mesh: compare the bottom figures between Figs. 1, 2, and 3, and compare the three figures in Fig. 4. Also, the levels of u jmax depend quite little on whether the rough or the fine spatial grid is being used (see e.g. Fig. 5). We also checked that the results are little dependent of the pseudo-random "draw" of the set of point-like "stars": another draw of 16 × 16 × 36 triplets (ρ, z, φ) gives very similar curves u jmax = f (λ j ) (Fig. 6). In summary, we now find that, for the relevant values of N , say N = 192 or N = 384, u jmax decreases smoothly from 10 27 to 10 21 eV/cm 3 when λ j varies in the domain considered, i.e., from λ 0.11µm to 830µm. We note moreover that, for the low values of λ j , the values of u jmax calculated using the present "separate fitting" have the same (extremely high) magnitude as those calculated with the former "grouped fitting" [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. These observations lead us to conclude that: (i) the extremely high values of u jmax (in the whole domain of λ considered) are really what the "Maxwell model of the ISRF" predicts for this model of the Galaxy. (ii) Somewhat surprisingly, it is the low values of u jmax obtained for the higher values of λ when the "grouped fitting" was used [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF] that were a numerical artefact.

Decrease of the energy density away from the axis

Recall that the maxima of the u j 's, which are extremely high, are always obtained for ρ = 0, i.e. on the axis of the (model of the) Galaxy, and that the energy density for ρ = 0 depends little on z in the domain investigated. The next questions are therefore: which is the extension around ρ = 0 of the domain of the very high values? Do such values imply that "too much EM energy" would be contained there? To answer these questions, we calculated the SED with successive lower and lower values of ρ max (see Note 4), starting from its value for the calculations shown on Figs. 1 to 3, i.e., ρ max = 9 kpc, and decreasing it to 1, 10 -1 , ..., 10 -14 kpc, using the "rough grid" parameters (see above), i.e., in particular, N ρ = 10, and using the S nj 's obtained with this rough grid with ρ max = 9 kpc -so that, for ρ max = 9 kpc, those cal-culations are not ones on the fitting grid. We looked in detail to the values u j (ρ m=2 , z p=1 ) = u j (ρ max /9, z = 0). The main results are shown on Fig. 7: even for very small values of ρ = 0, the values of u j are much smaller than u j max . That is, u j (ρ, z) decreases very rapidly when ρ is increased from 0. Actually, we found on the example of the smallest wavelength, corresponding with j = 1, that, from ρ = 1 kpc down to ρ = 10 -15 kpc, we have to a good approximation u 1 (ρ, z = 0) B/ρ, with

B = u 1 (ρ = 1 kpc, z = 0) 10 -0.45 (eV/cm 3 ).kpc. ( 27
)
This behaviour is not valid until ρ = 0, because for ρ → 0, u 1 (ρ, z = 0) tends towards u 1 (0, 0) < ∞, so we may assume u 1 (ρ, 0) B/ρ. On the other hand, Fig. 7 shows that there is nothing special to j = 1: we have u j u 1 , moreover for ρ 10 -15 kpc, u j depends only moderately on λ j . We observed in our calculations that, for ρ ≤ 1 kpc, u j (ρ, z) depends quite little on z with |z| ≤ z max = 1 kpc. Thus we may give the following approximation (which is likely an overestimate) to u j : for all j, and for |z| ≤ z max = 1 kpc, we have

u j (ρ, z) B/ρ for 10 -15 kpc ≤ ρ ≤ 1 kpc, (28)
u j (ρ, z) u j (ρ) B/ρ for ρ ≤ 10 -15 kpc, (29)
with u j (ρ) a decreasing function of ρ. According to Eq. (51) of the Appendix, this implies that, for |z| ≤ z max = 1 kpc, we have also

f x(ρ,z) (λ) B/ρ for 10 -15 kpc ≤ ρ ≤ 1 kpc, (30)
f x(ρ,z) (λ) B/ρ for 0 ≤ ρ ≤ 10 -15 kpc, (31)
independently of λ in the band

λ (1) := 0.1µm ≤ λ ≤ λ (2) := 830µm. ( 32
)
With this approximation, we can assess the total EM energy (53) contained in some disk

D(ρ 1 ) : (0 ≤ ρ ≤ ρ 1 , 0 ≤ φ ≤ 2π, |z| ≤ z max ), (33)
with ρ 1 ≤ 1 kpc, and in the wavelength band [λ (1) , λ (2) ]. This energy is bounded, owing to (30)-(31), by

W 1-2, D(ρ 1 ) Log λ (2) λ (1) × D(ρ 1 ) B ρ(x) d 3 x = Log λ (2) λ (1) × D(ρ 1 ) B ρ ρ dρ dφ dz, (34) i.e. W 1-2, D(ρ 1 ) Log λ (2) λ (1) × B ρ 1 × 2π × 2z max . (35)
But consider, instead of the disk D(ρ 1 ), the ring R(ρ 0 , ρ 1 ) :

(ρ 0 ≤ ρ ≤ ρ 1 , 0 ≤ φ ≤ 2π, |z| ≤ z max )
, with ρ 0 ≥ 10 -15 kpc. (Thus a ring with a very narrow aperture.) Using this time only (30), the same calculation gives

W 1-2, R(ρ 0 ,ρ 1 ) Log λ (2) λ (1) × B (ρ 1 -ρ 0 ) × 2π × 2z max . (36)
Taking ρ 0 = 10 -15 kpc, the conjunction of ( 35) and (36) shows that the contribution of the domain with 0 ≤ ρ ≤ ρ 0 is totally negligible, hence we may write

W 1-2, D(ρ 1 ) Log λ (2) λ (1) × B ρ 1 × 2π × 2z max . (37)
We can calculate the contribution δU that it gives to the average density of the EM energy in some disk D(ρ 2 ) of the Galaxy, with ρ 2 ≥ ρ 1 , making a volume V 2 = πρ 2 2 z max :

δU := W 1-2, D(ρ 1 ) V 2 4 Log λ (2) λ (1) × B ρ 1 ρ 2 2 . (38)
(Note that we may leave B in (eV/cm 3 ).kpc and ρ 1 and ρ 2 in kpc.) To give figures, let us first take ρ 1 = ρ 2 = 1 kpc, so that the corresponding value of δU is just the average volumic energy density U D(1 kpc) in the disk D(ρ 1 = 1 kpc). Then (38) with (27) give us

U D(1 kpc) 51 eV/cm 3 . ( 39
)
Note that this value is not very high. Another interesting application of Eq. ( 38) is to assess the effect, on that average value in the same domain D(1 kpc), of the domain of the "very high" values of the SED, say the domain for which u j ≥ 10 6 eV/cm 3 -i.e., from ( 27) and (28), ρ ≤ ρ vh , with ρ vh = 10 -6.45 3.55 × 10 -7 kpc 1.1 × 10 10 km,

which is almost twice the average distance Sun-Pluto, but still very small on a galactic scale. Taking to this effect ρ 1 = ρ vh and ρ 2 = 1 kpc in Eq. ( 38), the numerical values ( 27), ( 32) and (40) give us δU 4.54 × 10 -6 eV/cm 3 . In summary, the "very high" values of the SED are confined to the close neighborhood of the Galaxy's axis and contribute negligibly to the average energy density (39) in the disk D(1 kpc).

Results: spatial variation of the SED & comparison with the literature

This model's prediction for the spatial variation of the SED in the Galaxy was investigated, using again the separate fitting and the adjustment of the local SED on the measured values (both being described in Sect. 3). It was shown by using two different types of representations.

First, we plotted the SED at four different points in the Galaxy, and we compared the results with those obtained by Popescu et al. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF], who used a radiation transfer model built by them. (Their model also assumes axisymmetry.) show this comparison, our model being used here with N = 192 and N ω = 76. (Other choices of parameters that we tried gave similar figures.) It can be seen that the predictions of the present model do not differ very significantly from those of the radiation transfer model of Ref. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF]. The main difference is that our calculations oscillate somewhat strongly with the wavelength. The comparison of Figs. 8-11 here with Figs. 2345in Ref. [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF] shows that the difference between the results of the two models is significantly smaller now than it was in our previous work, in which the calculations were based on the grouped fitting [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]: the difference in log 10 (u j ) between the results of our model and Ref. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF] is here 1, whereas it went beyond 3 and even 4 in the previous calculations. However, the new calculations oscillate with the wavelength also at higher wavelengths. Whereas, when the grouped fitting was used, there was virtually no oscillation for λ 10 µm at the two positions at ρ = 8 kpc. (There were oscillations in the whole range of λ for the two positions at ρ = 1 kpc.) In order to check if those calculations inside the spatial domain of small values of the SED could be "polluted" by the extremely high values on the galaxy's axis, we investigated the effect of doing the fitting on a "shifted" grid with ρ ≥ 1 kpc. This did not lead to less oscillations. The general reason for these oscillations may be simply that this model takes fully into account the wave nature of the EM field.

Second, we plotted the radial and vertical profiles of the radiation fields at three wavelengths close to the ones considered in Fig. 7 of Popescu et al. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF] ("K, B , UV"). Figures 12 and13 show these profiles as they are calculated at points (ρ, z) belonging to the "logarithmic" grid on which the fitting was done for this calculation (see the legend). Those profiles of the radiation fields obtained with the present model on the fitting grid are relatively similar to those that are predicted with the very different model of Ref. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF], both in the levels and in the rough shape of the profiles. The most important difference is seen on the vertical profiles of Fig. 13: according to the Maxwell model of the ISRF, the energy density level decreases more slowly when the altitude z increases -or even, for the λ = 2.29µm radiation at ρ = 0.059 kpc or the λ = 0.113µm radiation at ρ = 7.5 kpc, the level of the SED does not decrease in the range considered for z. A similar lack of decrease is found on the radial profiles of Fig. 12, for the λ = 2.29µm radiation, either at z 0 or at z = 1.25 kpc. Using that same fitting done on a logarithmic grid, we also calculated and plotted the radial and vertical profiles of the same radiations, but this time for regularly spaced values of ρ (or respectively z), and in a larger range for ρ (or respectively z), Figs. 14 and15. The radial profiles of Figs. 12 and 14 are consistent, although, in contrast with Fig. 12, Fig. 14 plots the SED at points (ρ, z) which were not involved in the fitting, and which moreover involve an extrapolation to a larger range for ρ as compared with the fitting. The vertical profiles of Fig. 15, which also correspond with points which were not involved in the fitting, and also involve an extrapolation to a larger range for z as compared with the fitting, show an oscillating behaviour without any tendency to a decrease at large z.

Asymptotic behaviour at large ρ and at large z

To help understanding the behaviours just noted, in this section we study the asymptotic behaviour of the expressions of the components of the EM field and of the SED, as they are given by the Maxwell model of the radiation field.

The expressions ( 14)-( 16) that are implemented in the numerical model, are deduced from the exact integral expressions of the EM field for a given angular frequency ω, after the summation over the frequencies, and after the discretization ( 7) is done. Hence, we begin with the exact integral expressions of the EM field for a given angular frequency. These expressions, which are valid for any totally propagating, axisymmetric, time-harmonic EM field, are (Eqs. ( 13)-( 15) in Ref. [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF]):

B φ ω S = Re e -iωt +K -K √ K 2 -k 2 J 1 ρ √ K 2 -k 2 S(k) e ikz dk , (41)
E ρ ω S = Re -i c 2 ω e -iωt +K -K √ K 2 -k 2 J 1 ρ √ K 2 -k 2 ik S(k) e ikz dk , (42)
E z ω S = Re ie -iωt +K -K J 0 ρ √ K 2 -k 2 ω - c 2 ω k 2 S(k) e ikz dk , (43)
where K := ω/c -the other components being obtained by the duality [START_REF] Sommerfeld | Die Greensche Funktion der Schwingungsgleichung[END_REF] from the components (41)-( 43), with in the most general case an other spectrum function S (k).

The dependence in ρ of the components (41)-( 43) is determined by that of the Bessel functions J 0 and J 1 , and by the form of the integrals which involve them. At large x we have the asymptotic expansion [START_REF] Dieudonné | Calcul infinitésimal[END_REF]

J α (x) = 2 πx cos x - απ 2 - π 4 + O x -3 2 . ( 44
)
However, the argument of the Bessel functions in Eqs. ( 41)-( 43)

is x = ρ √ K 2 -k 2 .
Hence, as ρ → ∞, x does not tend towards ∞ uniformly, depending on the integration variable k: we even have x ≡ 0 independently of ρ, for k = ±K. Therefore, it is not obvious to see if the integrals (41)-(43) do have an expansion at fixed z as ρ → ∞.

As to the behaviour at fixed ρ and at large z: up to the real part, and for a fixed value of ρ, the components (41)-( 43) are expressions of the form e -iωt I(z), with

I(z) = b a f (k)e izg(k) dk, ( 45
)
and where, specifically, a = -K, b = +K, and the phase function is simply g(k) ≡ k, which has no stationary point. (The regular function k → f (k) depends on the component being considered, and also on ρ as a parameter.)

In that case, we have the asymptotic expansion [START_REF]Méthode de la phase stationnaire[END_REF]

I(z) = f (K) iz e izK - f (-K) iz e -izK + O 1 z 2 . ( 46
)
So at large z and for a fixed value of ρ, all components of any totally propagating, axisymmetric, time-harmonic EM field are order 1 z (unless the coefficient of 1 z in this expansion is zero, which is true only in particular cases -then the relevant component is higher order in 1 z ). This applies indeed to the part (i) of the decomposition (i)-(ii), that is given by Eqs. ( 41)-( 43), but also to the part (ii), since it is obtained from ( 41)-( 43) by applying the EM duality (18) (with, in the most general case, a different spectrum function S (k)). Hence the SED (23) is order 1 z 2 at large z, for any fixed value of ρ -when the C (q) j (x) coefficients correspond with the exact expressions (41)-( 43). [The explicit expression of the coefficient of 1 z 2 , depending on ρ , K, S(K), S(-K) (and, in the most general case, of the values S (K), S (-K) of the spectrum function S corresponding to the part (ii) of the decomposition (i)-(ii)) might easily be obtained from (23), (41)-(43), and (46).] The foregoing result applies to a general spectrum function S(k) (and S (k)). By summation on the frequency index j, it extends to an EM field having a finite set of frequencies.

Let us now investigate the asymptotic behaviour of the EM field and the SED, still in the totally propagating case with axial symmetry, but now after the summation over the frequencies and the discretization [START_REF] Arminjon | On the equations of electrodynamics in a flat or curved spacetime and a possible interaction energy[END_REF]. After the discretization, each among the C (q) j coefficients in the expansions (24) of the components of the EM field has the form [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]:

C (q) j = C (q) j (ρ, z) = N n=0 R (q) n J α ρ ω j ω 0 R n G nj (z) (α = 0 or α = 1), (47) where R
(q) n > 0 (except for R (q) 0 and R (q)
N , both of which turn out to be zero) are constant numbers, and where G nj (z) = exp (iω j t) F nj (t, z) is just the function F nj in Eq. ( 17) hereabove, deprived of its periodic timedependence (and thus is a periodic function of z). Together with (44), Eq. (47) shows that, at a given value of z, we have

C (q) j = O(1/ √ ρ) as ρ → ∞.
The SED for an EM field having a finite set of frequencies is given by Eq. ( 23).

For any given frequency (j), u j is a quadratic form of the C (q) j coefficients, hence

u j (ρ, z) = O 1 ρ (ρ → ∞). (48)
This is compatible with the curves shown on Fig. 14.

Passing to the behaviour at large z: in Eq. ( 47), the dependence in z is entirely contained in the functions G nj (z) which, we noted, are periodic. Hence, the coefficients C (q) j (ρ, z), each of which involves a linear combination of these functions (with coefficients that depend on ρ), are almost-periodic functions of z [START_REF] Ameriosi | Almost-periodic functions and functional equations[END_REF], and the same is true for the components (24) of the EM field. Moreover, for any given value of ρ, each u j in Eq. ( 23) is hence the sum of the square moduli of periodic complex functions of z. Therefore [START_REF] Ameriosi | Almost-periodic functions and functional equations[END_REF], the SED is an almost-periodic function of z, too. This result allows us to understand the lack of a decrease with z, observed on the vertical profiles of Fig. 15, which involve an extrapolation to a larger range for z as compared with the domain used for the fitting: an almost-periodic function f does not tend towards zero at infinity, unless f ≡ 0.5 As to Figs. 12 and 13, they involve no extrapolation, thus the relevant coefficients result from the fitting done on the very domain to which the curves belong. Hence the asymptotic behaviour of u j (whether at large z or at large ρ) is not relevant to them.

Discussion and conclusion

In this paper, we developed an improved numerical scheme to adjust the Maxwell model of the ISRF in a galaxy, which was proposed in a foregoing work [START_REF] Arminjon | An analytical model for the Maxwell radiation field in an axially symmetric galaxy[END_REF]. Namely, at the stage of fitting the radiations emitted by the many different point-like "stars" which make the model galaxy, we are now considering each time-harmonic component separately, which is more precise. This allows us as a bonus to eliminate the time variable at this stage, Eq. ( 21)thus reducing the computer time.

We used that "separate fitting" procedure, first, to check if the extremely high values of the spectral energy density (SED), which were predicted by this model on the axis of our Galaxy with the former "grouped fitting" [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF], are a physical prediction or a numerical artefact. A rather detailed investigation led us to conclude that these extremely high values are indeed what the model predicts -see Sect. 4.1. However, we find also that the SED decreases very rapidly when one departs from the galaxy's axis, see Fig. 7. Moreover, the average energy density of the EM field in, for example, a disk of diameter 1 kpc and thickness 2 kpc, is not very high, Eq. ( 39). The extremely high values of the SED on the axis of our Galaxy (and likely also in many other galaxies) are a new and surprising prediction for the ISRF. Recall that our model is adjusted so that the SED predicted for our local position in the Galaxy coincide with the SED determined from spatial missions, and thus is fully compatible with what we see of the ISRF from where we are. The prediction of the present model may be interpreted as a kind of self-focusing of the EM field in an axisymmetric galaxy. On the other hand, as we mentioned in the Introduction, the existing (radiation-transfer) models for the ISRF do not consider the EM field with its six components coupled through the Maxwell equations. These models consider paths of photons or rays and do not take into account the nature of the EM radiation as a field over space and time, subjected to specific PDE's. So a self-focusing cannot be seen with those models. It is difficult to assess the degree to which this prediction depends on the specific assumptions of the model, in particular the axial symmetry. If this prediction is at least partly confirmed, this will have important implications in the study of galaxies.

Second, we studied the spatial variation of the SED predicted by our model with the new procedure, and compared it with the predictions of a recent radiation transfer model [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF]. The difference between the results of the two models is much smaller now than it was [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF] with the older procedure. However, the SED predicted by our model still oscillates as function of the wavelength (or the frequency) also with the new, "separate fitting" procedure, although the different frequencies are then fully uncoupled. We also plotted the radial and vertical profiles of the radiation fields at three wavelengths. We confirm the slower decrease at increasing altitude z as compared with the radiation transfer model of Ref. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF], indicated by the previous work [START_REF] Arminjon | Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field[END_REF]. Actually, when the vertical profiles of the radiation fields are calculated and plotted in a domain that involves an extrapolation to a (three times) larger domain of z, a slightly oscillating behaviour without a decrease at large z is observed. This is explained by our study of the asymptotic behaviour of the analytical expressions of the EM field and the corresponding SED: we show that the SED calculated by the implemented model, that involves a discretization of the wave number, is a quasi-periodic function of z -although the exact SED obtained from the integral expressions (41)-( 43) is order 1/z 2 at large z. Thus, extrapolation on the domain of z should be used parsimoniously with the current numerical implementation based on a discretization of the wave number.

A Appendix: Discrete vs. continuous descriptions of the spectral energy density

The SED, u or rather u x (it depends on the spatial position x), is normally a continuous density with respect to the wavelength or the frequency: the time-averaged volumic energy density of the EM field at some point x and in some wavelength band [λ (1) , λ (2) ] is given by

U 1 2 (x) := δW 1 2 δV (x) = λ (2)
λ (1) u x (λ) dλ.

However, in many instances, including the present work, one is led to consider a discrete spectrum, thus a finite set of frequencies, (ω j ) (j = 1, ..., N ω ), hence a finite set of wavelengths. It leads also to a discrete energy density, Eq. ( 23). This raises the question of how to relate together these discrete and continuous descriptions of the SED. To answer this question, we note first that the u j 's in Eq. ( 23) are indeed volumic energy densities. Whereas, u x in Eq. (49) has physical dimension [U ]/[L], i.e., it is

f x (λ) := λu x (λ) (50)
which is a volumic energy density. And it is indeed f x that is being considered by Popescu et al. [START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF], when plotting the SEDs at different places in the Galaxy, or when plotting the radial or axial profiles of the radiation field at some selected wavelengths.

As is more apparent with the "separate fitting" used now (see Sect. 3), the discrete set of frequencies ω j , considered in the Maxwell model of the ISRF, represents just a finite sampling of the actual continuous distribution. The link between the two descriptions is hence given simply by the following relation:

u j (x) = f x (λ j ) = λ j u x (λ j ). ( 51
)
Consider a bounded spatial domain D. The total EM energy contained in the domain D and in the wavelength band [λ (1) , λ (2) ] is given, according to Eq. ( 49), by

W 1-2, D := D δW 1 2 δV d 3 x = D λ (2)
λ (1) u x (λ) dλ d 3 x,

i.e., using (50):

W 1-2, D = D λ (2)
λ (1)

f x (λ) d (Logλ) d 3 x. (53)
If we are using a model considering a fine-enough finite set of wavelengths (λ j ), j = 1, ..., N ω , with λ 1 = λ (1) and λ Nω = λ (2) , we may use (51) to estimate the integral over λ in Eq. (53) as a finite sum, e.g. a Riemann sum:

λ (2) λ (1) f x (λ) d (Logλ) Nω-1 j=1 f x (λ j )(Logλ j+1 -Logλ j ) Nω-1 j=1 u j (x)(Logλ j+1 -Logλ j ) (54
) or a better approximation (trapezoidal, Simpson,...). Maximum of SED on (N ρ =10,ρ 0 =0)×(N z =21,z 0 =-1kpc,z max =1kpc) grid Maximum of SED on (N ρ =20,ρ 0 =0)×(N z =23,z 0 =-1kpc,z max =1kpc) grid (N ρ =10,ρ 0 =0)×(N z =21,z 0 =-1kpc,z max =1kpc) grid, N=192 (N ρ =20,ρ 0 =0)×(N z =23,z 0 =-1kpc,z max =1kpc) grid, N=192 (N ρ =20,ρ 0 =0)×(N z =23,z 0 =-1kpc,z max =1kpc) grid, N=384

N ω =76, N=192 N ω =76, N=384
N ω =76, N=96 N ω =76, N=192 N ω =76, N=384

  Maximum of SED on (N ρ =10,ρ 0 =0)×(N z =21,z 0 =-1kpc,z max =1kpc) grid Maximum of SED on (N ρ =10,ρ 0 =0)×(N z =21,z 0 =-1kpc,z max =1kpc) grid

Figure 1 : 1 -

 11 Figure 1: Effect of discretization number N for: N ω = 23, rough grid

Figure 2 : 1 -

 21 Figure 2: Effect of discretization number N for: N ω = 46, rough grid

Figure 3 : 1 -

 31 Figure 3: Effect of discretization number N for: N ω = 76, rough grid

Figure 4 :

 4 Figure 4: Effects of discretization number N and number of frequencies N ω , fine grid

Figure 5 :

 5 Figure 5: Comparison of two different grids

Figure 6 :Figure 7 :

 67 Figure 6: Comparison of two different draws of the set of "stars"

  Predicted energy density at ρ = 1 kpc, z = 0 kpc(N ρ =10,ρ 0 =0)×(N z =21,z 0 =-1kpc,z max =1kpc) grid, N ω =76, N=192[START_REF] Popescu | A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics[END_REF]

Figure 8 :

 8 Figure 8: SED at (ρ = 1 kpc, z = 0)

Figure 9 :

 9 Figure 9: SED at (ρ = 1 kpc, z = 1 kpc)

Figure 10 :

 10 Figure 10: SED at (ρ = 8 kpc, z = 0)

Figure 11 :

 11 Figure 11: SED at (ρ = 8 kpc, z = 1 kpc)

Figure 12 :

 12 Figure 12: Radial profiles of radiation fields. Fitting done on a logarithmic grid:ρ 1 = ρ max , ρ m = ρ m-1 × q (m = 2, ..., N ρ ); z 1 = z max , z k = z k-1 × q (k= 2, ..., N z ); q = 0.5. SED values at (ρ m , z Nz ) (m = 1, ..., 10), then at (ρ m , z 4 ) (m = 1, ..., 10), are plotted.

Figure 13 :

 13 Figure 13: Vertical profiles of radiation fields. Fitting done on the same logarithmic grid as for Fig. 12. SED values at (ρ 10 , z k ) (k = 1, ..., N z ), then at (ρ 3 , z k ) (k = 1, ..., N z ), are plotted.

Figure 14 :

 14 Figure 14: Radial profiles of radiation fields. Fitting done on the same logarithmic grid as for Fig. 12. SED values at regularly spaced values of ρ, starting at 0.1 kpc, and for z = 0, then z = 1.25 kpc, are plotted.

Figure 15 :

 15 Figure 15: Vertical profiles of radiation fields. Fitting done on the same logarithmic grid as for Fig. 12. SED values at regularly spaced values of z, starting at 0, and for ρ = 0.1, then ρ = 7.5 kpc, are plotted.

The outgoing wave (2) can also be characterized[START_REF]Sommerfeld radiation condition. Wikipedia, The Free Encyclopedia[END_REF] as being the only time-harmonic

If Eq.[START_REF] Dieudonné | Calcul infinitésimal[END_REF], or equivalently Eq. (21), were an exact equality instead of being an equality in the sense of the least squares, then of course it would imply Eq. (8) (with w j ≡ 1) as an exact equality.

This results from the most common definition of an almost-periodic function f[START_REF] Ameriosi | Almost-periodic functions and functional equations[END_REF]: the existence, for any > 0, of a relatively dense set of almost-periods. I.e., for any > 0, there exists a length l such that, for any x ∈ R, there is at least one numberT ∈ [x, x + l [ such that for any t ∈ R, |f (t + T ) -f (t)| ≤ . If f is not identically zero, let a ∈ R, such that f (a) = α = 0. Taking = α2 in the definition above, we thus have|f (a + T ) -f (a)| ≤ α 2 , hence |f (a + T )| ≥ α 2 .Since x can be taken arbitrarily large and since T ≥ x, this proves that f does not tend towards zero at infinity.


Matthieu Labat

Monika Woloszyn

Géraldine Garnier

Jean Jacques Roux

Assessment of the air change rate of airtight buildings under natural conditions using the tracer gas technique. Comparison with numerical modelling

Keywords: air change rate, single-room building, wind and thermal effects, tracer gas measurements, numerical modelling

Under natural conditions, air change rates are very sensitive to specific building elements and to climate conditions, even more so in the case of airtight buildings. Consequently, applying general correlations to such cases may lead to inaccurate predicted air change rates. Still, this approach remains valuable because of its simplicity compared to other methods such as wind tunnelling and CFD simulations. In this paper, the tracer gas concentration decay technique was selected to contribute additional information to classical air-tightness measurements. The measurements were used to fit the coefficients of a general single-node pressure model. Simulation results were found to be consistent with tracer gas measurements most of the time. However, the closeness of the fit is strongly related to the average pressure coefficients from the literature, which were estimated more precisely using the other techniques mentioned above. From a general point of view however, it would seem promising to extend this method to other buildings.

Introduction

As saving energy is so important in buildings nowadays, building envelopes need to be increasingly efficient in order to meet current regulation codes. The use of efficient thermal insulation materials has made it possible to reduce energy consumption by reducing heat conduction through the envelope. As this type of heat loss has progressively decreased over the years, the focus has shifted to other heat transfer modes, in particular, convective heat transfers. This type of transfer can be reduced by making the envelope air tight in order to limit air transfers between indoor and outdoor areas. While energy losses are reduced however, other considerations arise such as moisture accumulation and indoor air quality [START_REF] Blondel | Screening of formaldehyde indoor sources and quantification of their emission using a passive sampler[END_REF]Plaisance 2011, Karlson and[START_REF] Karlsson | A comprehensive investigation of a low energy building in Sweden[END_REF]. One way to deal with these issues is to improve ventilation systems [START_REF] Woloszyn | The effect of combining a relative-humidity-sensitive ventilation system with the moisture buffering capacity of materials on indoor climate and energy efficiency of buildings[END_REF][START_REF] Koffi | Analyse multicritère des stratégies de ventilation en maisons individuelles[END_REF]. On the other hand, natural ventilation is also useful given that it uses no energy at all. Natural ventilation can be easily achieved by considering openings [START_REF] Santamouris | Experimental investigation of the air flow and indoor carbon dioxide concentration in classrooms with intermittent natural ventilation[END_REF] or dedicated building elements such as a solar chimney [START_REF] Afonso | Solar Chimneys: simulation and experiment[END_REF].

Under natural conditions, air change rates vary greatly because of shifting climates. Indeed, the driving forces are wind and thermal effects, which are both time-dependent. Using constant air change rates would lead to a high level of inaccuracy and detailed approaches are therefore preferable. Numerous studies have been conducted on such cases and have shown that natural driven air change rates are difficult to assess, even when considering simple buildings [START_REF] Larsen | Single-sided natural ventilation driven by wind pressure and temperature difference[END_REF]Heiselberg 2008, Li et al. 2001). As a consequence, it would seem difficult to extend the results reported in the literature to existing buildings.

The input of additional information on air transfers would make it possible to select appropriate results from the literature and to enhance modelling. The main objective of this paper is to propose and to test a simple experimental technique which will help improve a classic single-node pressure model. The results from both the measurements and the simulation will be presented and discussed.

This study is based on an existing airtight building with a single opening. The air tightness of existing buildings is commonly measured using depressurizing methods, which provide preliminary information. However, these methods do not provide information on the driving forces of natural ventilation. Wind tunnelling measurements provide useful additional information about wind effect, which is highly dependent on building geometry, envelope permeability and the surroundings [START_REF] Koinakis | The effect of the use of openings on interzonal air flows in buildings: an experimental and simulation approach[END_REF]). Nevertheless, this method is expensive and time-consuming and in many cases, therefore, other solutions are preferable, such as computational fluid dynamics (CFD) modelling [START_REF] Van Hoof | On the effect of wind direction and urban surroundings on natural ventilation of a large semi-enclosed stadium[END_REF]. Still, natural ventilation is time-dependent because it is strictly dependent on climate conditions, which can vary greatly within a single day. Numerous climate conditions should be simulated so that accurate results can be assessed. Another experimental method called the tracer gas technique is described by [START_REF] Cheong | Airflow measurements for balancing of air distribution systemtracer gas technique as an alternative?[END_REF]. This is an interesting alternative, as it works properly under natural pressure differences. It has already been applied to measure the air change rate in existing buildings by [START_REF] Santamouris | Experimental investigation of the air flow and indoor carbon dioxide concentration in classrooms with intermittent natural ventilation[END_REF] and Karlson and Moshfegh (2007).

The tracer gas technique

Various tracer gas techniques have been examined. SF6 is the most commonly used tracer gas as it is detectable at very low concentrations [START_REF] Afonso | Solar Chimneys: simulation and experiment[END_REF][START_REF] Chao | Ventilation performance measurement using constant concentration dosing strategy[END_REF][START_REF] Janssens | Effects of wind on the transmission heat loss in duo-pitched insulated roofs: a field study[END_REF], Laporthe et al. 2001, Santamouris et al. 2008). However, the technique is not limited to a single gas and other gases such as CO2 [START_REF] Baker | A method for estimating carbon dioxide leakage rates in controlled-environment chamber using nitrous oxide[END_REF][START_REF] Blondel | Screening of formaldehyde indoor sources and quantification of their emission using a passive sampler[END_REF][START_REF] Collignan | Use of metabolic-related carbon dioxide as tracer gas for assessing air renewal in dwellings[END_REF] and N2O [START_REF] Baker | A method for estimating carbon dioxide leakage rates in controlled-environment chamber using nitrous oxide[END_REF], Karlson and Moshfegh 2007[START_REF] Koinakis | The effect of the use of openings on interzonal air flows in buildings: an experimental and simulation approach[END_REF], Laporthe et al. 2001) have also been found suitable in many cases.

The criteria defined by [START_REF] Cheong | Airflow measurements for balancing of air distribution systemtracer gas technique as an alternative?[END_REF] are used to help compare the advantages and drawbacks of these three gases. The results are presented in Table 1 below. Laporthe et al. (2001) experimentally compared SF6 and N2O in order to investigate the effects of the higher density of SF6. However, this did not seem to affect the results significantly. In fact, one conclusion drawn was that SF6 may be preferable because it was more detectable. [START_REF] Baker | A method for estimating carbon dioxide leakage rates in controlled-environment chamber using nitrous oxide[END_REF] experimentally compared N2O and CO2, concluding that both were as accurate when experimental solutions were set to limit the drawbacks of using a tracer gas which is already present in the atmosphere. In the end, CO2 was selected here to obtain the tracer gas measurements as it has less impact on the environment (greenhouse effect) and it is also the least expensive solution.

Concerning the gas injection, several techniques can be found in the literature. [START_REF] Cheong | Airflow measurements for balancing of air distribution systemtracer gas technique as an alternative?[END_REF] defines four of them precisely:

• Constant injection: since the flow release rate is known, it is possible to achieve a tracer gas balance so that air-leakage rates are shown. It is a high gas consumption method, but it is well adapted to long-term observations [START_REF] Afonso | Solar Chimneys: simulation and experiment[END_REF]. • Constant concentration: the principle is very similar, but the system is far more complex (Koinakis 2005, Janssens and[START_REF] Janssens | Effects of wind on the transmission heat loss in duo-pitched insulated roofs: a field study[END_REF]. When used correctly, such systems can lead to a high level of accuracy. It is also the most expensive of the methods presented here. • Pulse injection: a small quantity of gas is injected into a duct and its concentration is measured further along [START_REF] Tanner | Use of carbon dioxide as a tracer gas for determining in package airflow distribution[END_REF]. Consequently, this method is of no interest here. • Concentration decay: a high quantity of gas is released and mixed into the system. Since the tracer gas concentration is measured along with the decrease, air leakages can be calculated [START_REF] Blondel | Screening of formaldehyde indoor sources and quantification of their emission using a passive sampler[END_REF].

Concentration decay is the most widespread method used to study the air change rate in a building because it is the easiest to set up and gives accurate results. This method was compared experimentally with the constant concentration method by [START_REF] Chao | Ventilation performance measurement using constant concentration dosing strategy[END_REF]. Both were found to achieve the aim, but the author warned about the importance of achieving a widespread sampling in the building, particularly if the measurements take place in a multizonal building. Similar conclusions were also reached by [START_REF] Van Buggenhout | Influence of sampling positions on accuracy of tracer gas measurements in ventilated spaces[END_REF]. In order to measure the air change rate in our experimental house, it was finally decided to use the concentration decay method.

The air change rate is deduced from the general tracer gas mass balance equation, given in [START_REF] Cheong | Airflow measurements for balancing of air distribution systemtracer gas technique as an alternative?[END_REF] and presented in (1). The calculation of the air change rate will be discussed further in 3.2.

( ) ( ) ( ) ( ) ( ) ( ) t C t C t Q t S dt t dC V In Out V In -  + =  (1)
where CIn and COut are the CO2 concentration of indoor and outdoor air, respectively (ppm), V is the indoor air volume (m 3 ), QV the air leakage rate (m 3 .s -1 ) and S is the amount of gas released indoors.

Experimental set-up 3.1 A brief description of the experimental set-up

The experimental house is situated in Grenoble, France (latitude: 45.2°E, longitude: 5.77°N) and is exposed to a natural climate. It consists of a single-room building, the dimensions of which are designed to be representative of a living room (4.56×4.55×2.36 m 3 interior dimensions). The indoor air volume is approximately equal to 49 m 3 (this value will be used further to compute the air change rate of the whole test house, as presented in (1)).

The structure is made of spruce studs (section: 0.07×0.165 m 2 ), positioned every 0.60 m. The walls are composed of gypsum board on the indoor side, cellulose wadding used as an insulation material between the spruce studs, particle board on the exterior side, a rain screen, a ventilated air gap and wooden cladding. The ceiling and the floor are heavily insulated and vapour transfer is prevented by two layers of vapour barrier. Further details about the test house can be found in [START_REF] Piot | Experimental wooden frame house for the validation of whole building heat and moisture transfer numerical models[END_REF].

The roof is a typical tiled French roof with two 30° slopes facing north and south. The floor has been elevated to 0.57 m above the ground in order to simplify the boundary conditions for the numerical models. The door is located in the middle of the northern face and has been insulated with polystyrene in order to achieve a homogeneous level of insulation.

The test house is also equipped with a ventilation system that is located in the attic: the indoor opening is located in the ceiling while the outdoor opening is located close to the eaves on the western face, 3.43 m up the wall of the experimental house, as shown in Fig. 1(a andb). Another opening is located on the door but has been sealed to increase the air tightness of the building. The indoor temperature and relative humidity are also measured, and a weather station, located close to the test house, records the climatic conditions every 10 min. In particular, wind speed and direction are measured 8 m above the ground.

Tracer gas equipment and setting the initial and final indoor concentrations

The outdoor CO2 concentration varies throughout the day, even though the experimental house is not set too close to any human activity area. In order to monitor this, a first gas analyser was set up outside at the weather station. Preliminary measurements showed that the outdoor CO2 concentration remained between 400 and 600 ppm for several months, as shown in Fig. 2. A second analyser was set indoors, while the outdoor CO2 concentration remained monitored. Both devices measure the infrared absorption of gases and are calibrated to detect CO2. The indoor analyser is equipped with a pump which continuously samples indoor air. To enhance the sampling efficiency, four ventilators were installed to mix indoor air and four sampling tubes were added to the analyser. Carbon dioxide is released from a bottle set inside the experimental house.

In order to minimize the uncertainty related to the outdoor CO2 concentration, the indoor concentration should be set far higher than the outdoor concentration. At the same time, it is better to reach the gas analyser's saturation rate in order to take advantage of its sensitivity. These two limits are used to define the initial and final indoor concentration, respectively CIn(t=0) and CIn(t=tStop). To first order, they can be used as boundary conditions to solve (1), considering the average outdoor concentration Out C , a constant air change rate and no gas release. As a consequence, the elapsed time tElapsed while the indoor concentration decreased from CIn(t=0) to CIn(t=tStop) depends on these two limits, as presented in (2).

( ) ( )         - = - =  = Out Stop In Out In V Elapsed C t t C C t C Q V t 0 ln (2)
It is preferable to conduct a single experiment over several hours so that the outdoor conditions variations have a detectable influence and are easier to identify. Based on preliminary experiments, the average air change rate and outdoor CO2 concentrations were estimated equal to 6 m 3 .h -1 and 500 ppm, respectively. The indoor gas analyser range is [0; 10,000 ppm]. The initial concentration was set at 9,500 ppm (this value is preferred to 10,000 ppm to avoid extreme scale values). The final CO2 concentration was set at 2,000 ppm, so the indoor concentration is at least four times higher than outside. As a consequence, the average elapsed time for a single experiment should reach 15 h. Similarly, the average air change rate can be deduced from (2), measuring the elapsed time only. This is a very convenient way to discuss the general influence of outdoor conditions, yet a detailed analysis would still contribute more information. This approach will be discussed further in 5.1.

Before carrying out air change rate measurements, several preliminary tests were conducted and aimed at checking:

• The efficiency of the sampling system (ventilators and sampling tubes);

• The low impact of indoor air mixing on the global air change rate;

• A good match with estimated extracted rates when the ventilation system was active (three air changes per hour); • Low reactivity between carbon dioxide and indoor sides and with water vapour.

Air change rate modelling

General equations

Air can move into and out of the building by several means: through openings (here, the ventilation opening located on the western face), through walls because of their air permeability (here, the four vertical walls) and through cracks or at junctures (here, at the door frame).

Once the N air paths have been identified, it is necessary to consider the pressure field all around the building and to compute a pressure balance in order to correctly model air transfers throughout the whole building [START_REF] Koffi | Analyse multicritère des stratégies de ventilation en maisons individuelles[END_REF]. For each single i air path (with i from 1 to N), both wind and thermal effects have to be estimated. Next, the pressure difference between indoor and outdoor air dPi (Pa) can be evaluated using (3). As wind and thermal effects vary quickly within a day, each pressure value is assumed to be time-dependent.

( ) ( ) ( ) ( )

t
+ = + (3)
where PInt is the indoor air pressure (Pa), PDyn is the dynamic pressure induced by wind (Pa) and dPTh is the pressure difference induced by the thermal effect (Pa).

Next, related airflow rates can be estimated using the general relationship (4).

( ) (

)

i n i i i V t dP C Q  = , (4)
where QV,i is the air flow (m 3 .h -1 ), C is the discharge coefficient (m 3 .s -1 .Pa -n ), n is a non-dimensional coefficient which characterizes the flow (1 is for laminar flows, 0.5 for turbulent flows). Every value is related to the i air path being considered.

Finally, the N airflows are used to compute a mass balance which leads to the calculation of the indoor pressure PInt used in (3). To avoid numerical problems, (3) is first computed using the indoor pressure from the previous time step.

As already mentioned, wind tunnel measurements would provide accurate and useful information, but these take too long to set up. General results from the literature have therefore been used to describe the wind effect. This can be estimated very simply by computing the dynamic wind pressure as presented in (5):

( ) ( ) ( ) ( ) 2 , , , 2 1 z t V C t t P Z i P Out i Dyn    =   (5)
where ρOut is the density of outdoor air (kg.m -3 ), CP is a non-dimensional coefficient depending on wind incidence (θ) and VZ (m.s -1 ) is the wind speed at wall altitude z (m) away from the test house.

In order to evaluate the pressure difference induced by thermal effects (the so-called stack effect), the general relationship (6) was used. This requires the definition of a single parameter; the altitude zi (m) of the openings and leakages through which air can move.

( ) ( ) ( ) ( ) t t z g t dP In Out i i Th   -   - = , ( 6
)
where g is gravity (= 9.81 m.s -2 )

Using air tightness measurements to estimate air flow rates

Air tightness tests (referred to below) rely on creating significant pressure gradients between indoor and outdoor climates. Here, this has been achieved using a device called a Permeascope ® , developed by the French company ALDES. This device was selected over the well-known Blowerdoor ® technique because its design is much better suited to small airtight rooms. Indoor air is extracted by means of the ventilation duct (see Fig. 3(a)) and the extracted air flow rate is measured. Even if these kinds of tests take place under greater pressure differences (from 20 to 100 Pa) than those occurring under natural conditions (from 2 to 20), [START_REF] Walker | A comparison of the power law to quadratic formulations for air infiltration calculations[END_REF] have shown that ( 4) is still valid.

As pressure and air flow rates are both measured, Ci and ni values such as the ones used in (4) can be estimated by performing linear regressions on the measurements. This is not immediate however, because a single test includes the air flows coming from several air paths (here, from the four vertical walls and from the leakages at the door frame).

A vapour barrier is another material commonly used in wood-framed houses to prevent moisture damage, which also strongly reduces air transfers and contributes to achieving airtight buildings. Here, a vapour barrier which has an impact on moisture transfers equivalent to that of an 18-m-thick air layer (Sd = 18 m), was placed in the vertical walls. A second air tightness test was carried out and the results are presented in Table 2. To allow quick comparison, the I4 value (from the French standard) was included. This is the ratio of the air change rate occurring under 4 Pa to the building surface (21 m 2 ). Indeed, 4 Pa is roughly representative of pressure differences occurring under natural conditions. Even if lower I4 values can now be achieved with passive buildings, the envelope air tightness level is high compared to most existing buildings, as presented in [START_REF] Sfakianaki | Air tightness measurements of residential houses in Athens, Greece[END_REF]. When assuming that the leakages at the walls are reduced to almost naught when the vapour barrier is in place, it becomes apparent that half the leakages are located at the door frame. A similar conclusion was also reached by [START_REF] Piot | Hygrothermique du bâtiment : expérimentation sur une maison à ossature bois en conditions climatiques naturelles et modélisation numérique[END_REF] considering other measurements. When assuming that the remaining leakages (the second half) are equally distributed over the four vertical walls, (4) can be rewritten for the door as (7a) and for a single vertical wall as (7b) (the pressure difference is bound to be different for every single vertical wall).

( ) ( ) (

)

n Door Door V t dP C t Q  = 2 , ( 7a
) ( ) ( ) ( ) n Wall Wall V t dP C t Q  = 8 , ( 7b
)
The air flow rate which occurs through the ventilation opening remains to be defined. A third test (see Fig. 3(b) referred to later as the global airflow was therefore carried out, during which the air was extracted through the opening located at the door (the seal was removed for this test only). Leakages at the door frame and through the walls were also taken into account, as shown in Fig. 3(a). The results are presented in Table 3. When assuming that the leakages at the wall and the door frame are the same in both experiments, the airflow rate crossing the ventilation opening can be deduced by subtracting the air tightness measurements from the global airflow measurements. So can be rewritten as ( 8) in order to model the air flow through the ventilation opening.

Uncertainties on wind and thermal effects

Modelling wind and thermal effects may be difficult to achieve for a real building. Firstly, the thermal effect can be easily computed at the ventilation opening given that its height hopening is clearly identified. This is not the case when evaluating the thermal effect for leakages spread over the entire height of the wall, as illustrated in Fig. 4(a).

Secondly, the pressure field all around the building is non-homogeneous and depends on wind direction: this is mathematically represented by the use of the CP values in (5). However, these values are also strongly dependent on the building geometry, as shown by [START_REF] Uematsu | Wind pressure acting on low rise buildings[END_REF]. The dependency is greatest when looking at specific locations, such as the roof and its edges. Unfortunately, this is where the ventilation opening is located. Moreover, the wind effect is very significant on an opening, which results in higher uncertainties on the airflow rate modelled at the ventilation opening, as illustrated in Fig. 4(b). At the same time, the vertical walls are covered with an open jointed cladding, which may also affect the pressure field acting on air transfers. An alternative way of modelling thermal effects is to set two equivalent heights hWall and hDoor which would represent the global thermal effect at the walls and at the door, respectively. Concerning the wind effect, it was decided to use general CP values (see Table 4), which were already proposed in the HAM-Tools model for vertical walls. It is well known that pressure coefficients vary greatly with the building geometry, as presented in [START_REF] Uematsu | Wind pressure acting on low rise buildings[END_REF] and [START_REF] Larsen | Single-sided natural ventilation driven by wind pressure and temperature difference[END_REF]. (Costolla et al. 2010) advise conducting detailed studies to assess these pressure coefficients accurately, even when the building geometry studied is simple, as here. Indeed, the authors showed that using average pressure coefficient values can lead to substantial errors. Moreover, errors were found to be the greatest when the openings were the smallest; the ratio exceeded 5 in their study. However, the authors also indicated that few methods are able to achieve this goal (wind tunnel experiments and detailed CFD studies, for example). Most of the time, these means are both expensive and time-consuming. As a consequence, the results reported in the literature could still be used and should fit with average trends.

Local wind speed in ( 5) is simply estimated, as presented in [START_REF] Hagentoft | Introduction to Building Physics[END_REF]: it is based on formula (9), which considers a logarithmic wind speed profile. It requires the use of a wind speed measurement V10 obtained 10 m above the ground. In order to model the influence of the surroundings, two constants are introduced and called "terrain coefficients" (k; a). The values given by the author are reported in Table 5.

a Z z k V V   = 10 (9)

Results and discussion

Tracer gas measurements and general analysis

Twelve experiments, conducted from October 2010 to January 2011, are listed in Table 6. The last two columns indicate average values so as to compare experiments, yet further calculations use timedependent values. The average airflow rate (Qv) calculation is based on (2) and uses the elapsed time measured.The very last column indicates the temperature difference between indoors and outdoors, representative of the thermal effect. The moisture content difference is also influential and will be taken into account in further calculations. In all cases the electric heater was active and maintained indoor air at a constant temperature. Finally, two distinct regimes were observed during experiment no. 1, separated 3 hours 15 minutes after the beginning of the experiment, so results were split in two lines in Table 6. Results for this specific case will be discussed further. 1) Based on the indoor CO2 decrease from 9,500 to 5,720 ppm (2) Based on the indoor CO2 decrease from 5720 to 2,000 ppm (3) This experiment was stopped when the indoor CO2 concentration reached 3,100 ppm

The results clearly show that the elapsed time can be very different from one experiment to another: it is more than 30 h long for experiment no. 2 and can diminish to almost 5 h for experiment no. 10. As the average airflow rate is computed directly from the decay of the carbon dioxide concentration, it closely follows elapsed time values and varies from 3 to 18 m 3 .h -1 (average air change rates ranging from 0.06 to 0.37 h -1 ). Given the rather low indoor air volume and low air permeability of the test house to air, results will be discussed further by considering the average airflow rather than the air change rate.

It should be noted that considering average weather conditions may lead to poor interpretation. Indeed, the temperature difference can strongly vary during a single experiment. For example, it increased from almost 0 to 15°C during experiment no. 2 (see Fig. 5(b)). This results in an increase of the air change rate, accelerating the CO2 concentration decrease, as it can be seen in Fig. 5(a). When the temperature difference remains the same and under low wind speed conditions, the air change rate does not vary much (see experiments no. 6 and no. 8). This results in a constant slope on the CO2 concentration over time (when plotted with half-logarithmic scales). The highest average wind speed value was recorded during experiment no. 10, which corresponds to the highest measured airflow rate. However, since wind speed did not remain constant during a single experiment, but indeed could vary substantially, more information was obtained by subdividing some experimental runs into short period. Experiment no. 1 was split in two parts to resolve one period lasting 3 hours where the mean wind speed was recorded at 6.6 m.s -1 and a longer second period during which the wind speed had decreased to 1.2 m.s -1 (see Fig. 6(b)). This resulted in a significantly slower decrease rate of the indoor concentration as presented in Fig. 6(a).

The air change rate for experiment no. 10 was observed to be greater than for experiment no. 1a, despite comparable and high wind speeds for both periods. The main difference was noted to be wind direction: for no. 1a the wind blew mainly towards the east whereas for no. 10 it blew towards the west. From a general point of view, it is very convenient to consider average weather conditions and the total elapsed time to give a quick overlook of the experimental conditions and results. To explain the air change rate variations observed however, one should consider weather conditions in great detail. In the following sections, instantaneous weather conditions will be used to model air flows.

Air transfer model

Equations presented in part 4 were implemented in a Matlab / Simulink environment and a variable time-step solver was used to solve the mass balance at the building scale. Because of the uncertainties related to wind and thermal effects, the model was tested repeatedly with different values for hWall, hDoor and {k; a}. In order to compare each solution, computed air change rates were used to recalculate the tracer gas decrease and computed elapsed times were set against measurements. Somehow, this represents the model's mean accuracy.

A more detailed approach was used to complement the average performance. Indeed, in some cases comparisons based only on the elapsed time are not relevant. For example, an important error occurring at a single time step may result in a permanent bias on the elapsed time, as illustrated in Fig. 7(a). Equation ( 2) was changed to (10) to compute the hourly averaged air change rate from tracer gas measurements. The results are plotted in Fig. 7(b) as "Measurements".

( ) ( )         -  + -   = Out In Out In V C t t C C t C t V Q ln ( 10
)
where ∆t (h) is the time step, here equal to 1 h.

Preliminary simulations showed that the model was very sensitive to the equivalent door height hDoor, poorly sensitive to hWall and behaved differently from terrain coefficient values. This could be a consequence of the rather complex interaction between wind and thermal effects, as shown in [START_REF] Li | Some examples of solution multiplicity in natural ventilation[END_REF]. When considering both the elapsed time and the averaged error on the air change rate, the best solution used a "scattered" terrain coefficient; hDoor = 2 m; hWall = 0.75 m. Detailed errors are presented in Table 7 for the best coefficient combination. A negative error on elapsed time means that the model overestimates the airflow rate. The average error on elapsed time is 15 %, 23 % when looking at the air change rate. Experiment no. 7 shows the necessity of studying errors on both the elapsed time and the airflow rates: the computed airflow rates are 25 % accurate and can be either overestimated or underestimated. The errors at different time steps compensate each other, so that the computed elapsed time almost matches the measured value.

Still, these results are not homogeneous, meaning that the model has possibly encountered limits. Indeed, the airflow rate error exceeds 2.2 m 3 .h -1 for experiment no. 9 and may have been overestimated by 82 % for experiment no. 1.a when the measured airflow rate was 8.4 m 3 .h -1 .

When looking at weather conditions, it becomes apparent that the wind speed V10 is greater than 5 m.s -1 when the highest errors are computed. However, computed airflow rates were correctly estimated for experiment no. 10, which also occurred under high wind speed. This could be a consequence of the use of CP values (see Table 4) to model the influence of wind incidence: this is consistent with earlier observations (see 4.1). In order to illustrate this statement, CP values from [START_REF] Hagentoft | Introduction to Building Physics[END_REF] were used to replace the previous values. Both references are compared in Fig. 8 and the elapsed times are presented in Table 8. The influence of CP values is not homogeneous, which attests to how sensitive the model is. However, the trend is not very clear: better results can be achieved using (Hagentoft 2001) CP values (experiments no. 1.a, 1.b, 6, 11) as well as less accurate results (experiments no. 3, 9, 10). Moreover, the influence can be sizeable with lower wind speed conditions if the thermal effect is less dominating (experiments no. 3 and no. 6). With high wind speed conditions, even if an improvement is observed for experiment no. 1, the results are less accurate for experiments no. 9 and no. 10. From a general point of view, better results were obtained with CP values coming from (Sharag-Eldin 2007).

Outlook

Depending on the desired level of accuracy, one could investigate the numerical fitting further by setting CP values in great detail. This could be achieved by studying the air flow around the tested building, which seems to be reasonable with respect to recent numerical improvements as presented in [START_REF] Obrecht | Towards urban scale flow simulations using the Lattice Boltzmann method[END_REF].

In this case however, the air transfer model developed will be implemented in a more global model called HAM-Tools [START_REF] Kalagasidis | The international building physics toolbox in Simulink[END_REF] to study the interactions between heat, air and moisture transfers. Weather data was analysed over a 6-month period. This has revealed that 95 % of the measured V10 values are below the 5 m.s -1 limit. This would indicate that the 12 experiments were carried out under a representative climate. The model is thought to be accurate enough to further investigate the coupling with such transfers, estimating the potential of natural ventilation to limit indoor vapour excess, for example. In the same vein, this approach could be extended to other pollutants.

Using tracer gas measurements has been found to be effective when fitting a general model, and its applicability to other buildings seems promising. However, a few limitations must be mentioned. Firstly, using tracer gas measurements with multizonal buildings needs to be validated further. Secondly, it is important to underline the importance of achieving tracer gas measurements under representative climates; another experimental campaign would need to be carried out if different average weather conditions were to be considered. So, the legitimate domain is the second limit of this approach.

Conclusions

The tracer gas method was used to obtain air change rate measurements in a 50 m 3 natural ventilated building. Twelve experiments were conducted from October to January under varying wind and thermal conditions, with the air change rate ranging from 3 to 18 m 3 .h -1 . As these rates are rather low, the experiment time lasted from 5 to 30 h. This is long enough to ignore short-term weather variations and allow us to consider general weather and wind data.

A classical single-node pressure model was used to estimate the global air change rate. General correlations were used to represent both wind and thermal effects. Experimental results were used to fit numerically obtained wind and thermal coefficient values. The results are 15% accurate when considering the elapsed time and 23 % when considering the global air change rate. As the uncertainties are related to the use of general pressure coefficients from the literature, better results could probably be obtained using wind-tunnelling or CFD simulations. Still, this method should remain valuable because of its simplicity and cost.

Fig. 1 :

 1 Fig. 1: Scheme of the experimental house with the door on the northern side and ventilation opening on the western side (a) and picture of the ventilation opening (b)

Fig. 2 :

 2 Fig. 2: Outdoor carbon dioxide concentration measured from 22/09/2010 to 20/10/2010

Fig. 3 :

 3 Fig. 3: Scheme of the air tightness test principle (a) and the global airflow test principle (b)

Fig. 4 :

 4 Fig. 4: Difficulties assessing the thermal effect on a widespread area (a) and wind effect at a specific location (b)

Table 4 :

 4 CP values depending on the wind incidence θ on the walls according to (Sharag-Eldin 2007) 69 0.54 0.35 0.08 -0.23 -0.61 -0.88 -0.92 -0.77 -0.58 -0.50 -0.42

Fig. 5 :

 5 Fig. 5: Indoor concentration decay for a few experiments with half-logarithmic scales (a) and temperature difference (b) occurring during experiments no. 2, 4, 6 and 8.

Fig. 6 :

 6 Fig. 6: Indoor concentration decay for a few experiments with half-logarithmic scales (a) and wind speed at 10 m (b) occurring during experiments no. 1, 5, 10 and 12.

Fig. 7 :

 7 Fig. 7: Comparison of measured tracer gas decay (a) and airflow rate (b) ("scattered" terrain coefficient; hDoor = 0.01 m; hWall = 0.75 m) for experiments no. 1 and no. 2.

Fig. 8 :

 8 Fig. 8: Comparison of CP values according to (Sharag-Eldin 2007) and (Hagentoft 2001)

Table 1 : Comparison of gases used as tracer gas based on the criteria defined by Cheong (2001)

 1

	Criteria	SF6	N2O	CO2
	No hazard	Greenhouse effect	< 300 ppm	< 30,000 ppm
	No chemical reactivity	✓	✓	✓
	Density (kg.m -3 )	5.11	1.53	1.87
	Distinctive from atmospheric gases	✓	✓	no -[300; 600] ppm
	Affordable and available	✓	✓	Cheapest solution

Table 2 : Air tightness measurements obtained with and without the vapour barrier

 2

	Test type	Vapour barrier	C (m 3 .s -1 .Pa -n ) n (-)	I4 (m 3 .h -1 .m -2 )
	Air tightness	No	5.64	0.69	0.62
	Air tightness	Yes	2.92	0.73	0.27

Table 3 : Comparison of the air tightness measurements with the global airflow measurement

 3

	Test type	Vapour barrier	C (m 3 .s -1 .Pa -n ) n (-)	I4 (m 3 .h -1 .m -2 )
	Air tightness	Yes	2.92	0.73	0.38
	Global airflow	Yes	13.03	0.60	1.43

Table 5 : Terrain coefficients used in (4) (Hagentoft 2001)

 5

	Terrain coefficient

Table 6 : Tracer gas measurements from October 2010 to January 2011, elapsed time while the indoor concentration decreased from 9,500 ppm to 2,000 ppm and average weather conditions

 6

				Average	Average	Average
	no.	Date	Elapsed time (h)	airflow rate	wind speed	temp. difference
				(m 3 .h -1 )	(m.s -1 )	(°C)
	1a 1b	04/10/2010	3h15 (1) 24h04 (2)	8.4 2.6	6.6 1.2	2.9 7.4
	2	06/10/2010	30h22	3.0	1.0	7.5
	3	08/10/2010	26h39	3.4	1.3	10.5
	4	11/10/2010	26h25	3.4	1.2	10.8
	5	12/10/2010	17h52	5.0	1.3	18.3
	6	15/10/2010	17h13	5.2	1.5	21.2
	7	03/01/2011	14h59	6.0	2.2	28.2
	8	04/01/2011	15h54	5.6	1.7	26.4
	9	06/01/2011	13h48 (3)	4.6	2.9	19.8
	10	07/01/2011	4h58	18.0	5.1	10.5
	11	10/01/2011	19h14	4.7	1.6	21.5
	12	25/01/2011	13h39	6.6	1.8	27.4
	(

Table 7 : Error computed on the time duration and the air flow rate for each experimental sequence with the best parameter combination

 7

	N°	1.a	1.b	2	3	4	5	6	7	8	9	10	11	12
	Measured Qv (m 3 .h -1 )	8.4	2.6	3.0	3.4	3.4	5.0	5.2	6.0	5.6	4.6	18.0 4.7	6.6
	Error Elapsed time	-54 % -22 % -14 % -7 % -18 % -7 % 1 % 2 % 8 % -31 % 3 % 5 % 1 %

Table 8 : Error on the elapsed time using CP values according to (Sharag-Eldin 2007) and to (Hagentoft 2001)

 8

	N°	1.a	1.b	2	3	4	5	6	7	8	9	10	11	12
	Measured Qv (m 3 .h -1 )	8.4	2.6	3.0	3.4	3.4	5.0	5.2	6.0	5.6	4.6	18.0 4.7	6.6

Sharag-Eldin 2007 82 % 27 % 24 % 11 % 24 % 14 % 12 % 25 % 19 % 47 % 15 % 20 % 18 % Hagentoft 2001 63 % 12 % 19 % 32 % 26 % 19 % 6 % 26 % 16 % 99 % 25 % 15 % 16 %

Acknowledgements

The authors would like to thank ADEME (the French Environment and Energy Management Agency, OPTI-MOB project, n°0704C0099) for its financial support.