Spectral energy density in an axisymmetric galaxy as predicted by an analytical model for the Maxwell field

An analytical model for the Maxwell radiation field in an axisymmetric galaxy, proposed previously, is first checked for its predictions of the spatial variation of the spectral energy distributions (SEDs) in our Galaxy. First, the model is summarized. It is now shown how to compute the SED with this model. Then the model is adjusted by asking that the SED predicted at our local position in the Galaxy coincide with the available observations. Finally the first predictions of the model for the spatial variation of the SED in the Galaxy are compared with those of a radiation transfer model. We find that the two predictions do not differ too much. This indicates that, in a future work, it should be possible with the present model to check if the"interaction energy"predicted by an alternative, scalar theory of gravitation, contributes to the dark matter. Keywords: Disk galaxy; interstellar radiation field; axial symmetry; electromagnetic field; Milky Way.


Introduction
The interstellar radiation field (ISRF) is a very important physical characteristic of a galaxy, and there is indeed a rich astrophysical literature that studies the modeling of the ISRF; see e.g. [1,2,3,4,5,6] and references therein. This subject is of importance not only in itself, but also because the radiation field produced by the stars and other bright objects has a strong interaction with the cosmic rays, which are another important component of the interstellar medium. The necessity of a modeling is obvious, since we can directly measure the radiation field only in the solar system, in fact in the neighborhood of the Earth [7,8,9]. The ISRF is significantly affected by its interaction with interstellar dust (e.g. [10,11]). In order to account for this interaction, equations for radiative transfer are solved: either by using Monte-Carlo methods, that simulate in a probabilistic framework the propagation of photons, and their absorption and scattering by dust, in a discretized spatial domain (e.g. [1,4]); or, by using "ray-tracing" algorithms that calculate the variation, due to absorption and scattering by dust, of the radiation's specific intensity along a finite set of directions (e.g. [3,4,6]). Thus, until now, the electromagnetic (EM) field, made of coupled electric and magnetic fields, has not been considered in the studies of the ISRF in galaxies, even less an EM field that would be a solution of the Maxwell equations -as should actually be the case for the real ISRF, however.
Obtaining for the ISRF such an exact solution of the Maxwell equations is somewhat demanding from both the mathematical and the numerical points of view, and has been the main purpose of a recent work [12]. A decisive motivation for that recent work was to prepare for checking an exotic prediction of an alternative, scalar theory of gravitation: namely, that the total energy(-momentum-stress) tensor, which provides the source of the gravitational field, is not in general the sum of the energy tensor of the EM field and the energy tensor of the material medium producing it. Instead, each time that an EM field is present together with a variable gravitational field (which is the general case), there should appear an additional energy tensor: the "interaction" energy tensor, T inter [13]. This speculative energy tensor T inter , being thus present virtually everywhere according to that theory, and having clearly an exotic character, could contribute to the dark matter. In order to calculate T inter , and to see if its distribution could be similar to dark matter halos, it is necessary to be able to calculate the EM field and its first order derivatives [13], whence the work [12]. In the latter work, a model has been formulated. It is based on an analytical representation of the general solution of the source-free Maxwell equations in the axisymmetric case [14]. Moreover, the model has been numerically implemented and has passed a numerical validation test based on an exact solution [12].
The aim of the present paper is to make a first check of the predictions of that "Maxwell model of the interstellar radiation field" for the spectral energy distribution (SED) and its spatial variation in a galaxy. Such SEDs are a very important output of the mainstream models for the analysis of the ISRF and, as was recalled at the beginning of this paper, those models depend essentially of the radiation transfer (mainly through dust). Note that, in contrast, the present model [12] does not directly take into account any radiation transfer. Section 2 summarizes that model. Section 3 details the calculation of the SED in that model. Section 4 presents the results obtained with that model for a simulation of the Milky Way. To obtain these results, the model is adjusted on the locally observed SED. In this way, it does indirectly take into account the radiation transfer, because the observed SED is indeed affected by the latter.

The Maxwell model of the interstellar radiation field
In this section, we give a new presentation of the model which was built in Ref. [12], and that uses results of Ref. [14] and references therein. This presentation uses an easy extension of the result [14]. Such an extension was used implicitly in Ref. [12], and is explicitly formulated in Appendix A.

General assumptions for the Maxwell field
We make two assumptions regarding the form of the Maxwell radiation field (E, B) (with E the electric field and B the magnetic field): i) We assume that it is an axially symmetric field. This is a relevant first approximation for many galaxies, even though it neglects some aspects like the arm structure of a spiral galaxy. Axial symmetry is often assumed in the literature on modeling the ISRF in galaxies and/or its interaction with dust (e.g. [2,3,6,10]).
ii) We search for a source-free solution of the Maxwell equations. This is because we seek to describe the EM field at the galactic scale, not the local EM field in the neighborhood of the stars that are the primary sources of that field.

The explicit representation for axisymmetric free Maxwell fields
This is based on the following recent Theorem [14]: any time-harmonic axisymmetric solution (E, B) of the source-free Maxwell equations has a unique decomposition where (E 1 , B 1 ) has the form (in cylindrical coordinates ρ, φ, z corresponding with the symmetry axis Oz): and where (E ′ 2 , B ′ 2 ) is deduced from a solution (E 2 , B 2 ) having just that same form (2)-(4), by the EM duality -i.e., (Here c is the velocity of light.) Note that the dual solution (5) has the complementary zero and non-zero components as compared with the direct solution (2)-(4), i.e., are time-harmonic with the same angular frequency ω as has (E, B). In Eqs. (2)-(4), the potential Ψ(t, ρ, z) = e −iωtΨ (ρ, z) is a time-harmonic axisymmetric solution of the standard scalar wave equation, which also has the same frequency ω as has (E, B); in general, the two solutions (E 1 , B 1 ) and (E 2 , B 2 ) correspond with two different scalar potentials Ψ and Ψ ′ [14]. The solution (2)-(4) actually derives from a vector potential having just its axial component A z = Ψ nonzero, i.e., A = Ψe z [14,15]. The Theorem [14] states the existence of the two potentials Ψ and Ψ ′ , not their uniqueness. In Appendix A, we extend this result to EM fields having a finite spectrum, which is easy, and allows us to give a more rigorous presentation of the model in Subsect. 2.3.
If a time-harmonic, axisymmetric solution of the scalar wave equation is totally propagating, it has the following explicit form [15,16]: with ω the angular frequency, K := ω/c, J 0 the first-kind Bessel function of order 0, and where S is some (generally complex) function of k ∈ [−K, +K]. "Totally propagating" means, in short, that Ψ has no "evanescent" component [15]; a detailed explanation is given in Ref. [12]. In turn, we shall define simply that a time-harmonic, axisymmetric EM field F = (E, B) is totally propagating if the two potentials Ψ and Ψ ′ , whose existence is guaranteed by the Theorem recalled above, can themselves be chosen totally propagating, i.e., having the form (6).

Statement of the model
We model a disk galaxy as a disk-like set of point-like "stars". 1 Each of them is assumed to produce an EM field with a finite spectrum: that field derives, in the way defined in the previous section, from potentials that obey the scalar wave equation, more precisely ones with spherical symmetry: Thus, the sum of the potentials emitted by the different stars at the different frequencies is: 1 The corresponding points are got by pseudo-random generation of their cylindrical coordinates. Here, the same set of 16 × 16 × 36 triplets (ρ, z, φ) was used as in Ref. [12], thus 9216 points. The distribution of ρ and z is approximately that valid for the star distribution in the Milky Way, and the distribution of φ is uniform, thus ensuring approximate axial symmetry.
Here r i := |x − x i |, with x the spatial position at which the function is evaluated and x i the spatial position of the i-th star; the numbers w j > 0, with j w j = 1, are the weights affected to the different frequencies. 2 Since the disk-like distribution of the "stars" is built (approximately) axisymmetric [12], the function (8) is also (approximately) axisymmetric. Moreover, this is a solution of the scalar wave equation, and it is clearly a totally propagating one, as is each of the individual spherical scalar waves in (8). It should thus be possible to put the function (8) -at least approximately -in the form of a sum, over the frequencies, of functions having the form (6). This is done by the least-squares method applied on a finite set of events (t, x) (with x = x(ρ, φ = 0, z) restricted to the plane φ = 0, due to the axial symmetry), making a regular spatio-temporal grid {(t l , ρ m , z p )}. We thus want to determine complex functions where the sign ∼ = indicates that the equality is in the sense of the least squares (on the relevant spatio-temporal grid, just mentioned). To the unique function on the r.h.s. of (9): Ψ = j ψ ω j S j , we may then associate: • either the EM field obtained by summing the time-harmonic contributions given by Eqs.

Discretization
To calculate the integrals in Eq. (6) and in Eqs.
This gives us for the EM field (10) [12], as also for the field (12): Even though these formulae follow from a discretization and from using an integration rule that is only approximate for the starting integral, they provide (without the O 1 N 4 remainder) an exact solution of the source-free Maxwell equations -see §3.3 in Ref. [12].
The arguments of the Bessel functions J 0 and J 1 in Eqs. (19), (20)-(22), as well as the spatial part of the argument of the complex exponential in Eqs. (19) and (23), are of the order of ρ/λ or z/λ. For a galaxy, these are numbers of the order of 10 25 , therefore the numerical implementation of the model needs to use a precision better than quadruple precision -which leads to long computation times [12].

Calculation of the spectral energy density
The volumic energy density of an EM field, in SI units, is with ǫ 0 = 1/(4π×9×10 9 ). In the relevant case with a frequency spectrum, the fields corresponding to the different frequencies do add, hence the quadratic expression (24) is not additive. The time average of such a quadratic expression is yet known to be additive in the case of a Fourier expansion [17]. Let us check if this remains true in the present case of an arbitrary finite frequency spectrum (ω j ) (j = 1, ..., N ω ). Let F (t) be some component of the EM field at some given spatial position x. We have, with some complex coefficients C j = C j (x): whence Unless But ω j + ω k = 0 does not occur, because ω j > 0, and ω j − ω k = 0 occurs only when j = k. Hence, (26) gives indeed the additivity for the time average of the square of a component admitting an expansion (25): For the EM field (10), its components j , respectively, which are directly read on Eqs. (20)-(22): From (24) and (27), we thus get for the field (10) a (time-averaged) energy density with a discrete spectrum: where the coefficients C (q) j (q = 1, 2, 3) are given by Eqs. (28)-(30), α 1 = ǫ 0 c 2 , α 2 = α 3 = ǫ 0 . As to the field (11), it is the dual (5) of the field (10). Therefore, the energy density (24) for the field (11) is given by the same Eq. (32). Finally, the energy density for the field (12) is hence just the double of this. Now, in the next section, the model will be adjusted so as to exactly describe the measured local EM spectrum, i.e., the values u j (x loc ) will be imposed to be the ones measured at our local position x loc in the Galaxy. In this context, the precise choice of the model, i.e., Eq. (10), or (11), or (12), is therefore totally neutral.

Results and comparison with a radiation transfer model
We first explain the final adjustment of the model. The scalar wave Ψ (x i ) (ω j ) (w j ) given by Eq. (8) depends on the numbers w j > 0, with j w j = 1. Each of them represents the relative weight of the corresponding frequency ω jthough in some average sense, since w j is common to all "stars" at the points x i . There is no exact way to assign values to these relative weights w j , which, together with the set (x i ) of the positions of the stars, determine the functions S j on the r.h.s. of Eq. (9) -or more exactly, determine the discretized values S nj in Eqs. (14) and (18). In turn, the values S nj (n = 0, ..., N), for a given value of the frequency index j, determine the intensity of the ω j -component of the scalar wave obtained after the fitting (the r.h.s. of (18)), and thus determine the intensity of the ω j -component of the EM field (20)-(22). However, since we start only from relative weights w j , this determination can be only up to a multiplying factor, and since the very values of the w j 's are assigned somewhat arbitrarily, that determination of the S nj 's (n = 0, ..., N) is really up to a factor depending on j, ξ j > 0. 3 We determine the values ξ j by imposing that the values u j given by Eq. (32), calculated at our local position x loc in the Galaxy, be equal to the ones measured. We take ρ loc = 8 kpc and z loc = 0.02 kpc (e.g. [18]). We took the local values of the SED from Fig. 1 of Ref. [2], using the software CurveUnscan to digitalize the curves.
Those local values were determined from the Apollo 17 mission [7] and from the COBE mission: the DIRBE experiment [8] and the FIRAS experiment [9]. Figure 1 shows the corresponding SED. For the calculations, in order to save computer time, we extracted 76 values of λ from the complete list of values of λ in the curve; the corresponding points are also shown on Fig. 1. Then, at any other spatial position x(ρ, z) (in the axially symmetric case considered), the model gives a prediction for the values u j (x) of the SED corresponding to the frequencies ω j (or the wavelengths λ j ) considered in our calculation. The spatio-temporal grid used was N t = 5, N ρ = 10, N z = 21 (thus 1050 events (t, ρ, z)), with t varying from 0 to T 0 = λ 0 /c with λ 0 = 10 µm, ρ varying from 0 to 10 kpc, and z from −1 to +1 kpc. Other grids have also been tested, giving very similar results. Moreover, unless otherwise mentioned, the discretization number N was N = 48. (This is the number of subintervals in which any of the integration intervals [−K j , +K j ] is subdivided, see Eqs. (13)- (14).) Figures 2 to 5 show the comparison of the SEDs, obtained at four different places in our Galaxy, according to either the thus-adjusted present model, or the radiation transfer model of Popescu et al. [6]. (Again, the relevant curves, Fig. 9 in Ref. [6], have been digitalized using CurveUnscan.)  Clearly, the general agreement between these two calculations is far from perfect. In particular, our calculation oscillates markedly as function of λ: at all wavelengths, for the two positions that are closer to the galactic centre (ρ = 1 kpc, Figs. 4 and 5); and merely at short wavelengths, for the two positions that are farther from the galactic centre (ρ = 8 kpc, Figs. 2 and 3). Nevertheless, this agreement is actually surprisingly good if one remembers how different are the two models. The SED predicted by our model is higher  than the one predicted by the radiation transfer model of Ref. [6] at the short wavelengths (λ 0.5µm, or log 10 λ −0.3) -except for the position (ρ = 8 kpc, z = 0 kpc), that is, relatively speaking, very close to our local position (ρ loc = 8 kpc, z loc = 0.02 kpc), where the adjustment is done. The higher levels found with our model at the short wavelengths might have a relation with the fact that the absorption by dust, not taken into account by our model, is strong for these radiations. However, this relation is not obvious, since our model is adjusted on the measured SED at our local position, and that local SED is of course affected by the dust absorption to which the radiations directed towards here have been subjected along their paths in the Galaxy. Actually, one notes also that, in the long wavelength domain as well as in the very short wavelength domain, the SED predicted by this model is, in average, higher than that predicted by the radiation transfer model, at the two positions which have a higher altitude (z = 1 kpc): Figs. 3 and 5.  that the present model with its current numerical settings predicts a slower decrease of the energy density as z (or rather |z|) increases than does the radiation transfer model. Of course the two models are very different, also in their modelling of the Galaxy: a discrete star distribution (for the present model) vs. an analytical, continuous variation of the stellar emissivity as function of ρ and z (for the radiation transfer model of Popescu et al. [6]). However, the ρ and z distribution of the stars in our model has been built to correspond (at least approximately) with that for the Milky Way [12], as has been built the ρ and z dependence of the stellar emissivity function for the radiation transfer model [6]. Of course also, the SEDs of a radiation transfer model (like that of Ref. [6]) are a prediction, too, not a measurement. Nevertheless, a further investigation is necessary to quantify the slower decrease of the SED with increasing z for the present model and to understand the reason for it.  A surprising prediction of our model is that for the maximum values u j max of the SED in our model of the Galaxy -the SED being calculated at each point on a regular (ρ, z) grid, as described above. Figure 6 shows these maximum values for three different spatial grids, with the same set of 23 wavelengths or frequencies. It is seen that the u j max values are quite similar for the different grids, and these values are extremely high at the short wavelengths, i.e., at high energies of the radiation field. They are found on the ρ = 0 axis, thus on the central axis of the galaxy. (Actually, the values u j (ρ = 0, z) are fairly independent of the value of z in the domain investigated, i.e., −1 kpc ≤ z ≤ 1 kpc.) However, if one refines the frequency mesh, going from N ω = 23 to N ω = 46, and then to N ω = 76, one observes that the very high values of u j max are obtained only at shorter and shorter wavelengths (Figs. 7 and 8). A similar phenomenon is observed when one  refines the discretization parameter N: for a given frequency mesh, increasing N reduces the domain of the high values of u j max towards smaller values of λ j , or, equivalently, decreasing N expands that domain towards larger values of λ j (Figs. 9-11). The likely common reason for the two effects is that, when one refines either the frequency mesh (N ω increased) or the discretization of the integration intervals (N increased), then the number of the solved-for parameters S nj (n = 0, ..., N; j = 1, ..., N ω ) increases. In most of the settings investigated in the present work, this number: N para = (N + 1) × N ω , is actually larger than the number N data of "data", that is the number of points in the spatio-temporal grid, 4 thus N data = N t × N ρ × N z . We have N data = 1050 for the grid used in this work, except for the two grids used for the first two curves in Fig. 6, which give respectively N data = 660 and N data = 550.
(We have N t = 5 also for those two grids.) Hence, we are often in a situation of "overfitting". 5 Higher values of the ratio R = N para /N data correspond with a more overfitted situation, in which one expects the predictions to be less accurate. Now, on each among Figs. 7-8 and 9-11, we saw that the domain of the low values of u j max increases when N para increases, the number N data being fixed at N data = 1050 for these five figures -thus we saw that the domain of the low values of u j max increases when R increases, i.e., when more overfitting occurs. Therefore, it is the low values of u j max which are likely to be a numerical artefact, not the extremely high values. Note that this effect can be seen also on Fig. 6: N and N ω , hence N para , are the same for the three curves. So the third and finest grid (N ρ = 10) × (N z = 21) leads to a lower R than do the other two grids (N ρ = 12)×(N z = 11) and (N ρ = 10)×(N z = 11) -and accordingly the domain of the high values of u j max is larger for the third curve. On the other hand, the effect of overfitting seems less important at a distance from the axis (where the values of the SED are not high at all): e.g., if one changes N ω from 76 to 23, the SEDs considered in Figs. 2-5 stay similar: they show a smaller number of oscillations due to the smaller number of points, but the amplitude of the oscillations is comparable. Tentatively, we may understand this by finding it natural that the very high numbers are more sensitive to overfitting. In summary, the status of the very high values of the maxima of u j (ρ, z): physical prediction of the model or numerical artefact, is not really clear at the present stage of this work -but the foregoing discussion, with different examples that all show an effect which can be attributed to the "amount of overfitting", seems to indicate that those high values might be a true prediction of the model.

Conclusion
In this work, the Maxwell model of the interstellar radiation field [12] has been first checked for its predictions of the spatial variation of the spectral energy distributions (SEDs) in our Galaxy. The model has been adjusted by asking that the SED predicted at our local position in the Galaxy, Eq. (32) with ρ = ρ loc := 8 kpc and z = z loc := 0.02 kpc, coincide with the observa- tions collected through different spatial missions [7,8,9]. Then the predictions of the model for other positions in the Galaxy (for which, of course, we have no real observation) have been compared with the predictions obtained by using a recent radiation transfer model [6]. The two predictions do not differ too much in magnitude, even though the predictions of the present model oscillate rather strongly as function of the wavelength, especially at short wavelengths and not very far (1 kpc) from the axis of the Galaxy. Also, the SED decreases more slowly with increasing altitude for the present model. The values of the SED reach their maximum on the axis itself of the Galaxy; that maximum is extremely high at short wavelengths. These surprisingly high values occur in a larger or smaller domain of wavelengths, depending on the settings of the calculation. However, this dependence is governed by the "amount of overfitting": less overfitting increases the range of the high values. This makes it plausible that the high values might be a true prediction of the model. The next work will aim at improving the (iii) to quantify the slower decrease of the SED with increasing altitude for the present model and to find the reason for it. Nevertheless, the present results: (1) show that the prediction, by the Maxwell model of the ISRF, of the spatial variation of the SED in the Galaxy, is quite comparable in magnitude to that provided by a radiation transfer model (except perhaps for high energies when the point is very close to the galaxy's axis); and therefore (2) indicate that an order-of-magnitude estimate of the "interaction energy" predicted [13] by an alternative, scalar theory of gravitation, should be indeed possible with the present model in a future work. A Appendix: Extension of the explicit representation to a finite spectrum In order to present this extension of the result [14] (the latter being recalled in Subsect. 2.2), it is necessary to write this result in terms of mappings. Consider the vector space S ω (respectively F ω ) made of those solutions Ψ of the scalar wave equation (respectively, made of those source-free Maxwell fields F = (E, B)) that are axisymmetric and time-harmonic with frequency ω. 6 We define two mappings Z 1 and Z 2 from S ω into F ω : given Ψ in S ω , Z 1 (Ψ) is the EM field given by Eqs.
(2)-(4), followed by Eq. (5). Obviously, Z 1 and Z 2 are linear. We can see that Z 1 (Ψ) belongs to F ω 1 , the subspace of F ω with E φ = B ρ = B z = 0, whereas Z 2 (Ψ) belongs to F ω 2 , the subspace of F ω with B φ = E ρ = E z = 0. The uniqueness of the decomposition (1) is just that of the EM fields Z 1 (Ψ) = F 1 ∈ F ω 1 and Z 2 (Ψ) = F 2 ∈ F ω 2 (if they exist), 6 As the domain of definition of the fields in S ω and F ω , we take the whole spacetime. More exactly, since we assume axial symmetry, we assume that they are defined for the whole range of values of the spacetime variables t, ρ, z, thus in the domain Ω: (t, z ∈ R, ρ ∈ R + ). such that the starting EM field F ∈ F ω is F = F 1 + F 2 . This results immediately from the complementarity of the non-zero components, noted in Subsect. 2.2, and inherent in the definition of F ω 1 and F ω 2 : the non-zero components of F 1 and F 2 must be just the corresponding components of F , hence are known from the data of F . Henceforth, for any EM field F , we will note F 1 (respectively F 2 ) the field whose components E φ , B ρ , B z are zero (respectively the field whose components B φ , E ρ , E z are zero), the other components being those of F . In contrast with the uniqueness, the existence of the decomposition (1) is a rather strong result, which amounts to state that the mapping Z from S ω × S ω to F ω defined by is surjective. This also amounts to state that, given any EM field F ∈ F ω , there exist Ψ and Ψ ′ in S ω , such that F 1 = Z 1 (Ψ) and F 2 = Z 2 (Ψ ′ ).
The general axisymmetric solutions of either the scalar wave equation or the source-free Maxwell equations are got by summing time-harmonic axisymmetric solutions of the respective equations. We restrict ourselves to such solutions having a finite set of frequencies, say (ω j ) (j = 1, ..., N ω ), for simplicity [12,14]. The axisymmetric solutions of the wave equation with this set of frequencies form the vector space S := S ω 1 ⊕ ... ⊕ S ω Nω . The axisymmetric solutions of the source-free Maxwell equations with this set of frequencies form the vector space F := F ω 1 ⊕ ... ⊕ F ω Nω . The sums are direct because the decomposition in a finite sum of time-harmonic fields is of course unique. Therefore, the linear mappings Z 1 and Z 2 extend immediately to ones from S to F 1 := F ω 1 1 ⊕ ... ⊕ F ω Nω 1 or respectively F 2 := F ω 1 2 ⊕ ... ⊕ F ω Nω 2 , by linearity, and so does the theorem. That is: for any EM field F = F (1) + ... + F (Nω ) in the direct sum F , the result valid for the timeharmonic case ensures that there exist 2N ω scalar potentials: Ψ (1) , Ψ ′(1) ∈ S ω 1 , ..., Ψ (Nω) , Ψ ′(Nω) ∈ S ω Nω , such that = Z 2 (Ψ ′(Nω ) ).