The Maxwell radiation field is an essential physical characteristic of a galaxy. Here, an analytical model is built to simulate that field in an axisymmetric galaxy. This analytical model is based on an explicit representation for axisymmetric source-free Maxwell fields. In a previous work, the general applicability of this representation has been proved. The model is adjusted by fitting to it the sum of spherical radiations emitted by the composing “stars.” The huge ratio distance/wavelength needs to implement a numerical precision better than the quadruple precision. The model passes a validation test based on a spherically symmetric solution. The results for a set of “stars” representative of a disk galaxy indicate that the field is highest near the disk axis, and there the axial component of dominates over the radial one. This work will allow us in the future to check if the interaction energy predicted by an alternative theory of gravitation might be a component of dark matter.
Apart from pure magnetic fields, which are thought to be produced by a galactic dynamo action [1,2], the electromagnetic (EM) field in a galaxy is in the form of EM radiation, i.e., of propagating EM waves covering the whole spectrum, from radio to gamma. Of course, a lot of information can be found in the literature regarding the production of the EM radiation field by stars and by other astrophysical objects; about its interaction with dust, as well as other processes of radiative transfer; or about the EM wave spectrum and its dependence on the position in the galaxy, etc. [3,4]. However, it seems that little or nothing can be found about the description of the EM radiation field in a galaxy as an exact solution of the Maxwell equations. The aim of the present work is to propose and first check a method to obtain such relevant solutions.
Our main motivation for that work is to make a step toward testing the following prediction  of an alternative theory of gravity, which is a scalar theory with a preferred reference frame: in the presence of both a variable gravitational field and an EM field, the total energy(-momentum-stress) tensor is not the sum of the energy tensors of matter and the EM field – there must appear a specific interaction energy tensor, which should be distributed in space, and be gravitationally active. (Ref.  summarizes ref.  while also giving the main motivations for that alternative theory of gravitation.) Essentially, if instead one assumes the additivity of the energy tensors of matter and the EM field, then charge is not conserved in a variable gravitational field according to that theory. (See ref. , which contains also a summary of that theory.) Thus, one needs to introduce that interaction energy tensor, the form of which is determined by asking that it should be Lorentz invariant in the absence of a gravitational field [5,6]. That interaction energy could possibly contribute to the dark matter, because in addition to being gravitationally active and distributed in space, it has an “exotic” character, being different from both the energy of matter and that of the EM field. (It has a classical nature, however, as has the theory [5,7]. That is, the theory does not say anything about the possibility that the interaction energy might result from underlying quantum particles.) The interaction energy tensor is characterized by a scalar (field), which is determined by the gravitational and EM fields in the preferred frame of the theory; hence, it depends also on the velocity of the reference frame used, with respect to that preferred frame. More precisely, in order to calculate that scalar field in a weak and slowly varying gravitational field, one has to know the EM field and its first-order derivatives, as well as the time derivative and the spatial gradient of the latter . (Here is the Newtonian potential.)
The precise goal of the present work, therefore, is to build a representative analytical model for the EM radiation field in a galaxy, especially for the spatio-temporal variation of that EM field. This is in order to be later able to calculate the interaction energy predicted by the scalar theory , and to check if its distribution resembles a “dark halo.” However, the present goal: to describe the EM radiation field in a galaxy as an exact solution of the Maxwell equations, is important more generally, as that field is an essential physical characteristic of a galaxy. The work done in this article is independent of any theory of gravity. Section 2 describes the model that has been built in this work. Section 3 discusses its numerical implementation. The numerical results obtained so far are discussed in Section 4. Finally, we present our conclusions in Section 5.
Let us state and comment the two essential assumptions of the model, i.e., the ones which it would be difficult to change without going to a different model:
(i) The structure of the interstellar radiation field is determined by the sum of spherical scalar potentials emitted by the stars. The relevant potentials are defined below, equation (19), and the corresponding sum is given by equation (20). That sum defines the source input that we use to determine the parameters of a general solution of the Maxwell equations, see the end of this paragraph. Note that this assumption is different from assuming that the total interstellar radiation field is the sum of the EM fields emitted by the stars. The latter assumption would be inappropriate, because it would mean neglecting the effect of important physical processes like dust extinction (by absorption and scattering), and, more generally, radiative transfer. Nevertheless, according to Maciel , “The interstellar radiation field in the optical and ultraviolet comes essentially from integrated stellar radiation.” Note, moreover, that in directions which are roughly orthogonal to a galactic disk or at least are strongly inclined with respect to it, the light emitted by the stars of that galaxy will suffer much less alteration as compared with what happens close to the galactic plane. This is important in connection with the aim of checking if the “interaction energy” alluded to in Section 1 might significantly contribute to the dark matter halos, since, precisely, the most part of such a halo is outside the galactic disk. Even more important in this connection is an essential feature of our model: that it provides an EM field which is an exact solution of the Maxwell equations. In contrast, if one would try to account directly for absorption, emission, and scattering processes, this feature would be quite difficult to maintain. But on the other hand, in our model, the sum of the spherical potentials generated by the point sources (by the “stars”) is fitted by a general analytical solution, see the Theorem and see the statement of the model below. The sum just mentioned thus merely serves to determine the “shape” (the parameters) of that general solution. That analytical solution, in itself, is valid independently of whether the corresponding EM field has been generated directly or has undergone various radiative transfers such as absorption and scattering.
(ii) The distribution of the stars and (hence) the interstellar radiation field are axially symmetric. This is of course not exactly true (cf. e.g. the arms of a spiral galaxy), but it seems to be a reasonable simplification which should provide a correct first approximation. Except for peculiar cases, e.g., if there were a correlation between the intensity of the EM field emitted by a star and its angular position in the galaxy, and except for the vicinity of a star, the axisymmetry of the stars’ distribution should indeed imply that of the interstellar radiation field. In any case, we shall assume that both the stars’ distribution and the interstellar radiation field are axially symmetric.
Each star (or bright object) is schematized as a point determined by its cylindrical coordinates : distance to the symmetry axis, azimuth, altitude, thus .
A discrete set of such points, (that is, ), is got by random generation, as follows:
Also, an exponential distribution is assumed for : for any , we draw values ( being an even integer) by quasi-random generation with a probability law independent of :
For each value thus obtained, we introduce another value , i.e., we impose a perfect symmetry w.r.t. in the distribution of .
Finally, for any two and , we draw values , with a uniform distribution between 0 and (thus ensuring the axial symmetry of the distribution of the “stars,” as announced).
In ref. , two classes of time-harmonic axisymmetric solutions of the Maxwell equations were introduced. Although the aim of that paper  was to obtain in explicit form nonparaxial EM beams, these two classes are relevant for the present, very different work. The first class of solutions is got in the following way. One starts from a time-harmonic axisymmetric solution of the scalar wave equation, and one associates with it a vector potential by
The second class of solutions of the source-free Maxwell equations is deduced from the first one above by applying the EM duality to any solution of the first class, i.e., by setting :
In ref. , we showed that, by combining these two classes, one can define a method that allows one to get actually all totally propagating, time-harmonic, axisymmetric source-free Maxwell fields – and thus, by the appropriate summation on frequencies, all totally propagating axisymmetric source-free Maxwell fields. (The necessary restriction of the method to totally propagating fields turns out to be appropriate to describe the radiation field.) The main result that allows this is the following one:
Theorem . Let be any time-harmonic axisymmetric solution of the source-free Maxwell equations (whether totally propagating or not). There exist a unique solution of the first class (4)–(6) and a unique solution of the second class, both with the same frequency as has , and whose sum gives just that solution:
Of course, as usual, it is implicit that, in equations (4)–(6), , and are actually the real parts of the respective r.h.s. Therefore, when is totally propagating and hence has the form (7), we obtain (using the fact that ):
As we already mentioned, the case of a general time dependence is deduced from the case with harmonic time dependence by considering a frequency spectrum. We shall consider a discrete (and finite) spectrum for simplicity. A totally propagating solution of the wave equation with a finite frequency spectrum is got by summing solutions of the form (7):
If, in the axisymmetric time-harmonic solution (7), one puts
If we have a set of spherical sources situated at the points , all sources having the same amplitude and the same frequency spectrum, (19) becomes:
(2.2) Similarly, calculate the associated EM field of the second class, with , by using the EM duality (11).
Step 1 (especially) is delicate numerically, as we will see in Section 3. Therefore, we shall focus in this article on Steps 1 and (2.1). Thus, in the sequel of this paper, we shall omit step (2.2), i.e., we shall consider only solutions of the first class. This means that we shall obtain axisymmetric EM fields having . “Complete” axisymmetric EM fields, obtained by summing solutions of the two classes, will of course have to be considered in the future work. At the stage of the fitting (Step 1), one may think to consider two different frequency spectra for the two classes, e.g., “mutually interpenetrating” ones.
In order to determine the “spectra” in equation (16) by fitting the sum (20) to this equation, a second method is a priori conceivable, and has indeed been tried in this work: by inverse Fourier transform. Considering again one axisymmetric time-harmonic solution (7) of the scalar wave equation, and removing its dependence in , a formal inverse Fourier transform gives us
Instead, the method we finally used to determine is by the values
The ratio in equation (23), thus a number of the order of , gives the magnitude of the spatial argument in the complex exponential of equation (20) and the argument of the Bessel function in equation (22). However, the Bessel function, as well as the real and imaginary parts of , oscillate around 0 with a pseudo-period which is of the order of unity (exactly , for and ). So already to get only the correct sign, one needs to know their arguments to a precision better than . In view of the magnitude of the arguments: , it means that 25 significant digits are needed to know just the sign of and in equation (20) and the sign of in equation (22). Therefore, double precision (16 significant digits) is not enough: quadruple precision (32 significant digits) is needed – and even, it is not a luxury. Implementing quadruple precision, using the Matlab function (for “variable precision arithmetic”), increases drastically the computation time. But, fortunately, we could reduce significantly the computation time (by a factor of approximately 20 for our programs), by using the external toolbox “Multiprecision Computing Toolbox for Matlab,” of Advanpix. In that toolbox, we imposed the number of digits to be 41, in order to reach the same level of numerical precision as with the Matlab function with default precision, i.e., 32 digits plus 9 “guard digits.”
As with equation (21) and (22) for the potential: each of , , and is the sum of the time-harmonic components given by equations (13)–(15). And as we did with to obtain equations (34) and (35), we use the approximation (30) to calculate the integrals in equations (13)–(15). We thus get:
Suppose one starts from an exact solution, with a finite frequency spectrum , of the scalar wave equation: given by equations (21) and (22) with exact “wave-vector spectra” . Then equations (39)–(41) for the EM field give in general only an approximation of the associated EM field (defined by summing the contributions (13)–(15)): this is due to the discrete integration method (30), using the discrete values (28) of the exact functions .
However, the integration formula (30) is exact (i.e., the remainder is not only but exactly zero) if the integrand function is a polynomial of degree 3. This is because the remainder is proportional to for some . Moreover, given e.g. the first four arguments and the four corresponding function values ( ), there exists one and only one polynomial of degree 3, such that ( ). (This is the well-known “unisolvence theorem;” note that we are considering a fixed value of the frequency index .) By construction, is a multiple of 3, hence the whole integration interval is the union of adjacent subintervals, each covering three steps . Thus, due to the unisolvence theorem: in each of those subintervals, the four successive arguments and the four corresponding function values define a unique 3rd-degree polynomial, for which the integration (30) is exact. It follows that the integration (30) is exact for the piecewise polynomial function which continuously extends those polynomial functions to the whole interval . Therefore, equations (34) and (35) actually define an exact solution of the scalar wave equation, which corresponds with substituting the functions for in equations (21) and (22). Similarly, equations (39)–(41) for the EM field give the exact result of adding the contributions (13)–(15) for the different frequencies , when in these contributions one considers, for the frequency , the spectrum function . In other words, equations. (39)–(41) provide an exact solution of the source-free Maxwell equations, deduced from an exact solution (34) of the scalar wave equation – all corresponding with the spectrum functions , which are piecewise third-degree polynomials.
The formulas (34) and (39)–(41) were implemented numerically and that numerical implementation was tested for the case with spherical symmetry, as follows. We can define an exact solution of the source-free Maxwell equations by equations (4)–(6), with the spherically symmetric time-harmonic solution (18) of the scalar wave equation. This yields (after taking the real part):
Figures 1, 2 and 3 show, for the three different scales investigated, the relative average quadratic differences between the fields , , , , as calculated either “directly,” i.e., by equations (18) and (43)–(45), or “with the spectrum (46),” i.e., by equations (34) and (39)–(41) applied with the spectrum values (46). The different scales were as follows: ( ), with . For this test, the frequency was taken to be , thus . The average quadratic differences are in general evaluated on a regular three-dimensional spatio-temporal grid (36) for the variables . Thus, takes values between and , takes values between and , and takes values between and , with , , and . However, in the present case with harmonic time dependence, we took , thus one value of , which was fixed at , with . So here the grid is two-dimensional in fact. Also, we took here , , and , .
As expected, the errors (the relative delta’s on the different fields calculated either “directly” or “with the spectrum (46)”) decrease strongly as the discretization of the (here constant) spectrum function becomes finer, i.e., with increasing . (See equation (28).) This validates the correctness of our calculations. However, the errors increase quickly when the scale length is increased (even though it remains here enormously smaller than the galactic scale). This is because the integrals in equations (22) and (13)–(15) involve functions of that oscillate with a frequency or a pseudo-frequency which is proportional to the magnitude of the spatial variables and . As these integrals are approximated by the discrete sums (34) and (39)–(41), one would have to increase the discretization number in proportion of the scale length .
In addition to this validation test, which does test our programs but is based on the “unphysical” solution (18), we will also compare the outcome of equations (39)–(41) with the exact solution that corresponds to a set of spherical sources situated at the points . That is, one defines an exact solution of the source-free Maxwell equations by summing solutions obtained by equations (4)–(6) applied to the outgoing spherical wave (19). That sum is similar to the sum given by equation (20). (Thus, all sources have the same amplitude and the same frequency spectrum, but this can easily be changed.) The fields produced by the individual source at are denoted by , and . In the case of a single frequency , one first applies equations (4)–(6) with as given by equation (19):
A file of some randomly generated “stars” (as described in Section 2.2) has been used here, more precisely one with triplets .
The angular frequencies are the same in the expression to be fitted, equation (20), and in the analytical expression used to fit it, equation (21). They are regularly spaced and symmetric around a central frequency , thus
A few different spacetime domains (variables ) of galactic dimensions have been used. The adopted sizes of the domain for the calculations discussed here were as follows:
We have discretized that domain using grids (36) with , , – in short (4,8,7) –, or (5,12,11), or (7,14,13). The symbols and in equations (61)–(63) are meant to indicate that, e.g., for equation (61), the discrete variation of begins with and ends with the largest multiple of the time step that is smaller than . The time step is , so – and similarly for and in equations (62)–(63). (This is just the same kind of variation as for the validation test of Section 3.4, but here has a value that is relevant to a galaxy.) Thus, while we browse a very small total interval of time using a very small time step, we do scan a large spatial scale, representative of a disk galaxy. The reason for imposing this difference is that the variation of the fields has a quasi-periodic character in time. This has been checked for the potential by calculating the sum (20) with a very large time step, close to . However, as we will see below, the fields have a definite spatial variation. On the other hand, the following values were tried for the number in equations (28)–(34): .
We will compare the field components , , , as calculated either “directly,” i.e., from equations (56)–(58), or “from the model,” i.e., from equations (39)–(41), using in the latter case in equation (42) the spectrum values obtained from fitting the sum (20) by equation (34), equation (38). Thus, in both cases, the field derives from exact solutions of the scalar wave equation by equations (3) and (4)–(6): either the spherical functions entering the sum (20), or the function (34), which is obtained precisely from fitting the sum (20).
Figures 4, 5, 6, 7, 8 and 9 show, for the spatiotemporal grid, the contour levels of the field components , , , as calculated either “directly” or “from the model” – in the latter case, on the same spatiotemporal grid used for the fitting, and with . In order to save place, we selected somewhat arbitrarily, and independently for each component, 3 values of the time among the available 5 values. Also, very similar figures are obtained if one uses another spatiotemporal grid, like or .
If the increase of the error with scale, found for the smaller scales investigated in the validation test of Section 3.4 (up to ), would continue up to the scale relevant to a typical disk galaxy ( ), then the relative quadratic errors between the field components calculated either directly or from the model (e.g., ) would reach huge values for the galactic scale that we are investigating in this section. For the validation test, the increase of the error with scale was due to the fact that the discretization number was not increased in proportion of the scale length. The exact spectrum (17) was available, and the discretized spectrum values were taken from it, equation (46). In contrast, in this section, the discretized spectrum values are now obtained from fitting the sum (20) by equation (34) over a spatio-temporal grid with galactic-like spatial dimensions, equations (62)–(63). For the present calculation, the relative quadratic errors are rather close to unity. The qualitative features of the fields are the same for the fields calculated directly or with the model, i.e., from the spectra obtained by fitting the potential:
The fields are more intense close to the axis. This feature is true for both the direct calculation and the model, but it appears more clearly with the model. Moreover, using that model, based on equations (39)–(41), one can calculate the three components , and also for – which is not the case for the and components when one uses the direct calculation based on equations (56)–(58), see after those equations. The model gives for , because . In contrast, for , that model predicts very high values for , of the order with the spectrum values obtained by fitting the sum (20) by equation (34) on the spacetime domain (61)–(63), using e.g. the parameters described in the paragraph following equations (61)–(63) – but using these spectrum values to calculate the fields on a shifted grid, with starting at instead of .
The maximum intensities (positive and negative), calculated either directly or from the model, have quite similar values – but the positions of the maxima are generally different between the direct calculation and the model, except for the fact mentioned, that they are close to the axis.
Remind that the schematization which leads to the direct calculation (56)–(58) is a rather simple one: an axisymmetric disk-like distribution of point sources, each of which emitting an EM field deriving via equation (3) from a spherically symmetric solution of the scalar wave equation. As we saw at the end of Section 3.3, also the field got from the model is an exact axisymmetric solution of the source-free Maxwell equations. The set of point-like sources, that may represent a disk galaxy, and that leads directly to the calculation (56)–(58), is also used to adjust the model. The latter, however, is based on the nonsingular (“continuous”) equations (34) and (39)–(41) – in contrast with equations (20) and (56)–(58), that are singular at each source (i.e., for ). Therefore, the field obtained from the model is at least as representative of the EM radiation field in a disk galaxy as the field calculated directly can be.
In this work, an analytical model has been built for the Maxwell field in an axisymmetric galaxy, in particular for that field which results from stellar radiation. This model is based on a representation of any totally propagating axisymmetric source-free Maxwell field as the sum of two fields given explicitly: in the case of a time-harmonic field, the first field is given by equations (13)–(15), and the second one is deduced by the duality (11) from a field of the same form. In a previous work, the general applicability of this representation has been proved.
The model is adjusted by fitting to it the sum of spherical radiations emitted by a set of point-like “stars.” The distribution of these point-like objects is axisymmetric. It builds a flat disk, symmetrical with respect to a plane perpendicular to the symmetry axis, and the dimensions have been chosen to represent a disk galaxy similar to the Milky Way. The model provides an exact solution of the source-free Maxwell equations, also after the discretization that is used to calculate the relevant integrals.
The huge ratio distance/wavelength needs to implement a numerical precision better than the quadruple precision. The model and the corresponding software have passed a validation test based on an exact solution with spherical symmetry. The results for a disk galaxy indicate that the field is highest near to the axis, and there the component dominates over . In a further stage, it will be possible to adjust the model so as to very accurately describe the measured local EM spectrum. The model presented will then provide a theoretical prediction for the spatial variation of the EM spectrum in our Galaxy, which it will be possible to compare to other predictions, based on very different astrophysical models [3,4].
The author is grateful to Christian Boily and to Garrelt Mellema for useful remarks at the Computational Astrophysics Conference in Saint Petersburg, September 2019.
Conflict of interest: The author states that he has no conflict of interest.
 Beck R , Wielebinski R . Magnetic fields in the Milky Way and in galaxies. in: TD Oswalt , G Gilmore (eds), Planets, stars and stellar systems. vol. 5, Dordrecht: Springer, 2013, p. 641–723. Search in Google Scholar
 Chamandy L , Subramanian K , Shukurov A . Galactic spiral patterns and dynamo action I: a new twist on magnetic arms. Mon Not R Astron Soc. 2013;428:3569–89. Search in Google Scholar
 Porter TA , Strong AW . A new estimate of the galactic interstellar radiation field between 0.1μm and 1000μm . In: Proc. 29th International Cosmic Ray Conference, Pune. Mumbai:Tata Institute of Fundamental Research; Vol. 4; 2005. p. 77–80. Search in Google Scholar
 Maciel WJ . The interstellar radiation field. In: Astrophysics of the interstellar medium. Maciel WJ , editor. New York: Springer; Chapter 2, 2013. p. 17–31. Search in Google Scholar
 Arminjon M. On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy. Open Physics. 2018;16:488–98. Search in Google Scholar
 Arminjon M . Interaction energy of a charged medium and its EM field in a curved spacetime. In: Geometry, integrability and quantization XX. Mladenov IM , Pulov V , Yoshioka A , editors, Sofia: Avangard Prima; 2019. p. 88–98. Search in Google Scholar
 Arminjon M . Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation. Open Phys. 2016;14:395–409. Search in Google Scholar
 Kent SM , Dame TM , Fazio G. Galactic structure from the Spacelab infrared telescope. II. Luminosity models of the Milky Way. Astroph J. 1991;378:131–8. Search in Google Scholar
 Robin AC , Crézé M , Mohan V . The radial structure of the Galactic disc. Astron Astrophys. 1992;265:32–9. Search in Google Scholar
 Porcel C , Garzón F , Jiménez-Vicente J . The radial scale length of the Milky Way. Astron Astrophys. 1998;330:136–8. Search in Google Scholar
 Schneider P . Extragalactic astronomy and cosmology: an introduction. Berlin: Springer; 2016. p. 55. Search in Google Scholar
 Zamboni-Rached M , Recami E , Hernández-Figueroa HE. Structure of nondiffracting waves and some interesting applications. In: Localized waves. Hernández-Figueroa HE , Zamboni-Rached , Recami E , editors, Hoboken: John Wiley & Sons; 2008. p. 43–77. Search in Google Scholar
 Garay-Avendaño RL , Zamboni-Rached M. Exact analytic solutions of Maxwell’s equations describing propagating nonparaxial electromagnetic beams. Appl Opt. 2014;53:4524–31. Search in Google Scholar
 Arminjon M. An explicit representation for the axisymmetric solutions of the free Maxwell’s equations. Open Phys. 2020;18:255–63. Search in Google Scholar
 Jackson JD . Classical electrodynamics. 3rd ed. Hoboken: John Wiley & Sons; 1998. p. 360. Search in Google Scholar
 Durnin J. Exact solutions for nondiffracting beams. I. The scalar theory. J Opt Soc Am A. 1987;4:651–4. Search in Google Scholar
 Gradshteyn IS , Ryzhik IM . Table of Integrals, Series, and Products. 7th English Edition. Burlington (Mass., USA): Academic Press; 2007. § 6.677, p. 722. Search in Google Scholar
 Ala G , Francomano E , Viola F. A wavelet operator on the interval in solving Maxwell’s equations. Prog Electromag Res Lett. 2011;27:133–40. Search in Google Scholar
 Atkinson KE . An introduction to numerical analysis. 2nd Edition. New York: John Wiley & Sons; 1989. p. 257–8. Search in Google Scholar
© 2021 Mayeul Arminjon, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.