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

1 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 and Wolfendale [1]. 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 et al. [2], Gordon et al. [3], Robitaille [4], and Popescu et al. [5], 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, and hence, even less one that would be a solution of the Maxwell equations.

In a recent work [6], 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 that, as discussed earlier, 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 [7] 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. [6], its application to predict the variation of the spectral energy density (SED) in our Galaxy has been subjected to a first test [8]. To this purpose, the model has been adjusted by asking that the SED predicted for our local position in the Galaxy coincides with the SED determined from spatial missions by Henry et al. [9], Arendt et al. [10], Finkbeiner et al. [11], and Porter and Strong [12]. It has been observed in that most recent work [8] 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. [5], 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 article 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). Section 2 provides a summary of the model. Section 3 describes the improvement of the numerical scheme. In Section 4, 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 Section 6, asymptotic expansions are used to interpret the findings of the foregoing section. We present conclusions in Section 7, that is followed by Appendix A, which discusses the relation between the discrete and continuous descriptions of the SED.

2 Short presentation of the model

This model has been presented in detail in ref. [6]. 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 pseudo-random generation of their cylindrical coordinates ρ , ϕ , z with specific probability laws, ensuring that the distribution of ρ and z is approximately valid for the star distribution in the galaxy considered, and that the set { x i } is approximately invariant under azimuthal rotations of any angle ϕ [6]. In the present work, as in refs [6,8], 16 × 16 × 36 triplets ( ρ , z , ϕ ) were thus generated, so that i max = 9,216 , and the distribution of ρ and z is approximately 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 [13]. 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 time-harmonic 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 [14]. Hence, the ISRF is made of “totally propagating” EM waves, i.e., ones without any “evanescent” component [6,16]. 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 earlier, are themselves totally propagating. In that case, both A j z and A j z ′ have the explicit form [15,16]:

where ω j is the angular frequency, K j ≔ ω j / c , J 0 is the first-kind Bessel function of order 0, and 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 ′ [13].

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 centre 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 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:

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

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 ′ , and 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. [6], 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), (20), or (21) 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) and (ii) (e.g., the one corresponding with two equal potentials A i j z = A i j z ′ = ψ x i ω j ).

In the previous works [6,8], this fitting was done for all frequencies at once. That is, the following least-squares problem was considered:

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

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 ′ , whence ψ ω j S j = A j z = A j z ′ on the right-hand side 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

where

with δ j ≔ 2 K j / N , is a regular discretization of the interval [ − K j , + K j ] for k in the integral (1). Calculating those integrals with the “Simpson 3 8 composite rule,” (4) becomes the computable least-squares problem

with

The S n j s are the solved-for parameters in the least-squares problem (8). In Eq. (8), N must be a multiple of 3, and in Eq. (9), we have

Part (i) of the decomposition of the model ISRF then obtains as follows [6]:

with R n = K 0 2 − k n 2 and

(Here, k n and a n ( 0 ≤ n ≤ N ) are as k n j and a n j 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:

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

3 Frequency-by-frequency fitting of the potentials

Eq. (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. (8). Naturally, also the weight w j may then be removed from the left-hand side, by entering the inverse 1 / w j into the unknown spectrum function S j on the right-hand side. Eq. (8), thus becomes

At this point, one notes that both ψ x i ω j [Eq. (3)] and f n j [Eq. (9)] 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 :

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

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 ∣ / λ ∼ 1 0 25 , which enforces us to use a precision better than quadruple precision in the computer programs, thus leading to slow calculations [6]. Note that the “separate fitting,” i.e., the least-squares problem (21), is not exactly equivalent to the “grouped fitting,” i.e., the least-squares problem (8) (this will be confirmed by the numerical results): the two are slightly different ways of adjusting the global potentials (1).^[2] However, Eqs. (13)–(19) apply with the separate fitting as well – although the relevant values S n j 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 (21), with N para = N + 1 unknown parameters and N data = N ρ × N z data (the “data” are the values of the left-hand side of (21) on the spatial grid G ′ ). However, with the formerly used grouped fitting, we had to solve just one least-squares problem (8) with N para = ( N + 1 ) × N ω unknown parameters and N data = N t × N ρ × N z data.

On the other hand, through 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. [3]). 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 [8], 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 [8]

where the complex numbers C j ( q ) ( x ) ( q = 1 , … , 6 ) are the coefficients in the expansion, in time-harmonic functions, of each among the six components of the EM field:

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 × 1 0 9 ) in SI units). For an axisymmetric EM field, it is enough to consider the plane ϕ = 0 , thus x = x ( ρ , z ) , and we have

Using in that case the decomposition (i) and (ii), the expressions of three among the C j ( q ) coefficients follow directly from the expressions (14)–(16) of the corresponding components of the EM field [8]. Moreover, in the special subcase (18) considered here, the other components are given by (19), whence in the same way the three remaining C j ( q ) coefficients.

Now note that, in the least-squares problem (21), that we use to determine the values S n j allowing to compute the EM fields (14)–(16) and (19), no data relative to the intensity of the fields emitted by the point-like “stars” have been used until now. Hence, we may multiply the left-hand side of (21) by some number ξ j > 0 , thus obtaining now new values S n j ′ = ξ j S n j ( 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. [12] (Appendix A), and we take ρ loc = 8 kpc and z loc = 0.02 kpc , see, e.g., ref. [19]. 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, and these models are very different from the present model.

4 Results: maximum energy density

In the foregoing work [8], 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 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

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 j max ≃ 1 0 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 on the settings is governed by the “amount of overfitting”: less overfitting increases the range of the high values [8]. 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.

4.1 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 [8], 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 = 0 N f n j S n j , 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, 2, 3, 4 show these effects. The most salient result is that the extremely high values of u j max are now found with all calculations and in the whole domain of λ – except that on some curves, abrupt oscillations towards 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. [6].) 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 (Figure 3), increasing N does not allow us to eliminate the abrupt oscillations towards lower values: it even happens then, that increasing N from 192 to 384 actually deteriorates the u j max = f ( λ j ) curve. We interpret this as due to the fact that, when using a rougher spatial grid G ′ for the fitting, less data are provided (the values taken on the grid G ′ by the left-hand side of Eq. (21)) to determine the unknowns S n j on the right-hand side 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 j max depend quite little on N ω i.e., on the fineness of the frequency mesh: compare the bottom figures of Figures 1–3, and compare the three figures in Figure 4. Also, the levels of u j max depend quite little on whether the rough or the fine spatial grid is being used (see, e.g., Figure 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 j max = f ( λ j ) (Figure 6). In summary, we now find that, for the relevant values of N , say N = 192 or N = 384 , u j max decreases smoothly from ≃ 1 0 27 to ≃ 1 0 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 j max calculated using the present “separate fitting” have the same (extremely high) magnitude as those calculated with the former “grouped fitting” [8]. These observations lead us to conclude that: (i) the extremely high values of u j max (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 j max obtained for the higher values of λ when the “grouped fitting” was used [8] that were a numerical artefact.

Figure 1

Effect of discretization number N for: N ω = 23 , rough grid.

Figure 2

Effect of discretization number N for: N ω = 46 , rough grid.

Figure 3

Effect of discretization number N for: N ω = 76 , rough grid.

Figure 4

Effects of discretization number N and number of frequencies N ω , fine grid.

Figure 5

Comparison of two different grids.

Figure 6

Comparison of two different draws of the set of “stars.”

4.2 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 (Note 4), starting from its value for the calculations shown on Figures 1–3, i.e., ρ max = 9 kpc , and decreasing it to 1 , 1 0 − 1 , … , 1 0 − 14 kpc , using the “rough grid” parameters (see above), i.e., in particular, N ρ = 10 , and using the S n j ’s obtained with this rough grid with ρ max = 9 kpc – so that, for ρ max ≠ 9 kpc , those calculations 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 in Figure 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 ρ = 1 0 − 15 kpc , we have a good approximation u 1 ( ρ , z = 0 ) ≃ B / ρ , where

(27) B = u 1 ( ρ = 1 kpc , z = 0 ) ≃ 1 0 − 0.45 ( eV/cm 3 ) kpc .

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, Figure 7 shows that there is nothing special to j = 1 : we have u j ≲ u 1 , moreover for ρ ≳ 1 0 − 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

(28) u j ( ρ , z ) ≃ B / ρ for 1 0 − 15 kpc ≤ ρ ≤ 1 kpc ,

(29) u j ( ρ , z ) ≃ u j ( ρ ) ≲ B / ρ for ρ ≤ 1 0 − 15 kpc ,

with u j ( ρ ) a decreasing function of ρ . According to Eq. (A3) of the Appendix, this implies that, for ∣ z ∣ ≤ z max = 1 kpc , we have also

(30) f x ( ρ , z ) ( λ ) ≃ B / ρ for 1 0 − 15 kpc ≤ ρ ≤ 1 kpc ,

(31) f x ( ρ , z ) ( λ ) ≲ B / ρ for 0 ≤ ρ ≤ 1 0 − 15 kpc ,

independently of λ in the band

(32) λ ( 1 ) ≔ 0.1 μ m ≤ λ ≤ λ ( 2 ) ≔ 830 μ m .

With this approximation, we can assess the total EM energy (A5) contained in some disk

(33) D ( ρ 1 ) : ( 0 ≤ ρ ≤ ρ 1 , 0 ≤ ϕ ≤ 2 π , ∣ z ∣ ≤ z max ) ,

with ρ 1 ≤ 1 kpc , and in the wavelength band [ λ ( 1 ) , λ ( 2 ) ] . This energy is bounded, owing to (30) and (31), by

(34) W 1 − 2 , D ( ρ 1 ) ≲ Log λ ( 2 ) λ ( 1 ) × ∫ D ( ρ 1 ) B ρ ( x ) d 3 x = Log λ ( 2 ) λ ( 1 ) ∫ D ( ρ 1 ) B ρ ρ d ρ d ϕ d z ,

i.e.,

(35) W 1 − 2 , D ( ρ 1 ) ≲ Log λ ( 2 ) λ ( 1 ) × B ρ 1 × 2 π × 2 z max .

But consider, instead of the disk D ( ρ 1 ) , the ring R ( ρ 0 , ρ 1 ) : ( ρ 0 ≤ ρ ≤ ρ 1 , 0 ≤ ϕ ≤ 2 π , ∣ z ∣ ≤ z max ) , with ρ 0 ≥ 1 0 − 15 kpc . (Thus, a ring with a very narrow aperture.) Using this time only (30), the same calculation gives

(36) W 1 − 2 , R ( ρ 0 , ρ 1 ) ≃ Log λ ( 2 ) λ ( 1 ) × B ( ρ 1 − ρ 0 ) × 2 π × 2 z max .

Taking ρ 0 = 1 0 − 15 kpc , the conjunction of (35) and (36) shows that the contribution of the domain with 0 ≤ ρ ≤ ρ 0 is totally negligible, and hence, we may write

(37) W 1 − 2 , D ( ρ 1 ) ≃ Log λ ( 2 ) λ ( 1 ) × B ρ 1 × 2 π × 2 z max .

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 :

(38) δ U ≔ W 1 − 2 , D ( ρ 1 ) V 2 ≃ 4 Log λ ( 2 ) λ ( 1 ) × B ρ 1 ρ 2 2 .

(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) gives us

(39) ⟨ U ⟩ D ( 1 kpc ) ≃ 51 eV/cm 3 .

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 ≥ 1 0 6 eV/cm 3 – i.e., from (27) and (28), ρ ≤ ρ vh , where

(40) ρ vh = 1 0 − 6.45 ≃ 3.55 × 1 0 − 7 kpc ≃ 1.1 × 1 0 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 × 1 0 − 6 eV/cm 3 . In summary, the “very high” values of the SED are confined to the close neighbourhood of the Galaxy’s axis and contribute negligibly to the average energy density (39) in the disk D ( 1 kpc ) .

Figure 7

Decrease of the SED in the neighbourhood of ρ = 0 . ρ is given in kpc. For − 0.7 ≲ log 10 λ ≲ + 0.7 , the u j max and u j ( ρ = 0 , z = 0 ) curves are hidden but nearly coincide. The u j max curve is on Figure 1(b).

5 Results: spatial variation of the SED and 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 Section 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. [5], who used a radiation transfer model built by them. (Their model also assumes axial symmetry.) Figures 8, 9, 10, 11 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. [5]. The main difference is that our calculations oscillate somewhat strongly with the wavelength. The comparison of Figures 8–11 of this study with Figs. 2–5 in ref. [8] 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 [8]: the difference in log 10 ( u j ) between the results of our model and ref. [5] is ≲ 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 .) 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.

Figure 8

SED at ( ρ = 1 kpc , z = 0 ).

Figure 9

SED at ( ρ = 1 kpc , z = 1 kpc ).

Figure 10

SED at ( ρ = 8 kpc , z = 0 ).

Figure 11

SED at ( ρ = 8 kpc , z = 1 kpc ).

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. [5] (“K, B , UV”). Figures 12 and 13 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. [5], both in the levels and in the rough shape of the profiles. The most important difference is seen on the vertical profiles of Figure 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 Figure 12, for the λ = 2.29 μ m radiation, either at z ≃ 0 or at z = 1.25 kpc . By 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 ), Figures 14 and 15. The radial profiles of Figures 12 and 14 are consistent, although, in contrast with Figure 12, Figure 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 Figure 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 .

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 N z ) ( m = 1 , … , 10 ) , and then at ( ρ m , z 4 ) ( m = 1 , … , 10 ) are plotted.

Figure 13

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

Figure 14

Radial profiles of radiation fields. Fitting done on the same logarithmic grid as for Figure 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

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

6 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. [6]):

where K ≔ ω / c – the other components being obtained by the duality (18) 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 [20]:

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

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 [21]

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) and (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 j ( q ) ( 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 (7). After the discretization, each among the C j ( q ) coefficients in the expansions (24) of the components of the EM field has the form [8]: