Open Access Published by De Gruyter Open Access March 16, 2021

# An analytical model for the Maxwell radiation field in an axially symmetric galaxy

Mayeul Arminjon
From the journal Open Physics

## Abstract

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 E 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.

### 1 Introduction

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 [5] 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. [6] summarizes ref. [5] 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. [7], 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 T U and the spatial gradient of the latter [5]. (Here U 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 [5], 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.

### 2 Description of the model

#### 2.1 Main assumptions

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 [4], “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.

#### 2.2 Distribution of the stars

• Each star (or bright object) is schematized as a point x determined by its cylindrical coordinates ( ρ , ϕ , z ) : distance to the symmetry axis, azimuth, altitude, thus x = x ( ρ , ϕ , z ) .

• A discrete set of such points, S = { x i } = { x ( ρ m , ϕ m p q , z m p ) } (that is, i = i ( m , p , q ) ), is got by random generation, as follows:

1. An exponential distribution is assumed for ρ > 0 , i.e., n ρ values ρ m ( m = 1 , , n ρ ) are got by quasi-random generation with a probability:

(1) P ( a < ρ m < b ) = 1 h a b e ρ h d ρ ,

where h is the scale length, with h = 3 kpc in the numerical computations. Such a distribution, with this value of h , is approximately correct for the Milky Way [8,9,10].

2. Also, an exponential distribution is assumed for z > 0 : for any m = 1 , , n ρ , we draw n z / 2 values z m p ( n z being an even integer) by quasi-random generation with a probability law independent of m :

(2) P ( a < z m p < b ) = 1 h z a b e z h z d z ,

with h z = 0.2 kpc in the numerical computations. This also roughly corresponds to the Milky Way [11].

3. For each value z m p > 0 thus obtained, we introduce another value z m p , i.e., we impose a perfect symmetry w.r.t. z = 0 in the distribution of z .

4. Finally, for any two m = 1 , , n ρ and p = 1 , , n z , we draw n ϕ values ϕ m p q , with a uniform distribution between 0 and 2 π (thus ensuring the axial symmetry of the distribution of the “stars,” as announced).

#### 2.3 Explicit representation for axisymmetric EM radiation fields

In ref. [13], two classes of time-harmonic axisymmetric solutions of the Maxwell equations were introduced. Although the aim of that paper [13] 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 Ψ ( t , ρ , z ) = e i ω t Ψ ˆ ( ρ , z ) of the scalar wave equation, and one associates with it a vector potential A by

(3) A Ψ e z , or A z Ψ , A ρ = A ϕ = 0 .

(We shall denote by ( e ρ , e ϕ , e z ) the standard, point-dependent, direct orthonormal basis associated with the cylindrical coordinates ρ , ϕ , z .) In the time-harmonic case considered for the moment, such a vector potential defines uniquely – as it was shown in ref. [13] and, in more detail, in ref. [14] – the following exact solution of the source-free Maxwell equations (i.e., the Maxwell equations in a domain where the source (the 4-current) is zero[1]):

(4) B ϕ = A z ρ , E ϕ = 0 ,

(5) E ρ = i c 2 ω 2 A z ρ z , B ρ = 0 ,

(6) E z = i c 2 ω 2 A z z 2 + i ω A z , B z = 0 .

In ref. [13], this class was defined only when the scalar wave Ψ has the following form:

(7) ψ ω S ( t , ρ , z ) = e i ω t K + K J 0 ρ K 2 k 2 e i k z S ( k ) d k ,

with ω the angular frequency, K ω / c ,[2] and J 0 the first-kind Bessel function of order 0. ( c is the velocity of light.) Any totally propagating,[3] time-harmonic, axisymmetric solution of the scalar wave equation can be set in the form (7) [12,13]. However, as noted in ref. [14], this form is not necessary at this stage and the solution (4)–(6) applies whether A z is totally propagating or not. On the other hand, given any integrable function S ( k ) ( S L 1 ( [ K , + K ] ) , equation (7) defines an axisymmetric solution of the standard wave equation for the whole range of values of the spacetime variables t , ρ , z , thus in the domain Ω : ( t , z R , ρ R + ) (i.e., ρ 0 ) [12,14]. Therefore, given any frequency ω > 0 and any integrable function S ( k ) , equations (4)–(6), applied with A z = ψ ω S given by equation (7), define an axisymmetric solution of the source-free Maxwell equations in the whole spatiotemporal domain, Ω .

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 [13]:

(11) E = c B , B = E / c .

In ref. [14], 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 [14]. Let ( A , E , B ) be any time-harmonic axisymmetric solution of the source-free Maxwell equations (whether totally propagating or not). There exist a unique solution ( E 1 , B 1 ) of the first class (4)–(6) and a unique solution ( E 2 , B 2 ) of the second class, both with the same frequency as has ( A , E , B ) , and whose sum gives just that solution:

(12) E = E 1 + E 2 , B = B 1 + B 2 .

Of course, as usual, it is implicit that, in equations (4)–(6), B ϕ , E ρ and E z are actually the real parts of the respective r.h.s. Therefore, when A z is totally propagating and hence has the form (7), we obtain (using the fact that d J 0 / d x = J 1 ( x ) ):

(13) B ϕ ω S = ℛe e i ω t K + K K 2 k 2 J 1 ρ K 2 k 2 S ( k ) e i k z d k ,

(14) E ρ ω S = ℛe i c 2 ω e i ω t K + K K 2 k 2 J 1 ρ K 2 k 2 i k S ( k ) e i k z d k ,

(15) E z ω S = ℛe i e i ω t K + K J 0 ρ K 2 k 2 ω c 2 ω k 2 S ( k ) e i k z d k ,

where K ω / c .

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 ( ω j ) ( j = 1 , , N ω ) is got by summing solutions of the form (7):

(16) Ψ ( ω j ) ( S j ) ( t , ρ , z ) = j = 1 N ω ψ ω j S j ( t , ρ , z ) ,

where each ψ ω j S j is a time-harmonic solution having the form (7). Note that the different frequencies do not necessarily have the same weight, since any given “wave vector spectrum” S j can be multiplied by a factor w j . In other words, the weights are contained in the functions S j .

#### 2.4 The case of spherical waves

If, in the axisymmetric time-harmonic solution (7), one puts

(17) S ( k ) c 2 ω K < k < K ω c ,

then one gets a spherical time-harmonic solution of the scalar wave equation: this yields indeed [13,17]

(18) ψ ω S c 2 ω ( t , ρ , z ) = e i ω t sinc K r , K ω c ,

where sinc θ sin θ θ , r x = ρ 2 + z 2 . However, because sin K r = ( e i K r e i K r ) / ( 2 i ) , this solution contains both an outgoing wave and an ingoing wave, and therefore it does not satisfy the Sommerfeld radiation condition. Up to a multiplying coefficient (amplitude), there is only one outgoing spherical solution of the scalar wave equation that is time-harmonic with a given frequency ω , and of course it is well known:

(19) ψ ω ( t , ρ , z ) = e i ( K r ω t ) K r , K ω c .

Obviously, this is a totally propagating solution. Nevertheless, unfortunately, it does not seem to be amenable to the integral form (7) with an analytical function S ( k ) – in contrast to equation (18). However, for a spherical source (as opposed to a sink), we must use equation (19).

If we have a set of spherical sources situated at the points x i , all sources having the same amplitude and the same frequency spectrum, (19) becomes:

(20) Ψ ( x i ) ( ω j ) ( S j ) ( t , ρ , z ) = i = 1 i max j = 1 N ω S j e i ( K j r i ω j t ) K j r i ,

where r i x x i , K j ω j / c , and setting again the initial phases to zero for simplicity. Of course one might also make the amplitude of the source, as well as the weights affected to the different frequencies ω j , depend on the source, i.e., on the index i , by giving a dependence on i to the positive numbers S j , which would thus become S i j .

#### 2.5 The model

• Step 1. Determine a relevant axisymmetric solution Ψ of the scalar wave equation, by fitting to the form (16) the sum (20) of the spherical radiations emitted by the “stars” that make the “galaxy.”

• Step 2.

1. (2.1) Calculate the associated EM field of the first class, ( E ρ , E z , B ϕ ) with B ρ = B z = E ϕ = 0 , by summing its time-harmonic contributions given by equations (4)–(6), with successively A z ψ ω j S j ( j = 1 , , N ω ), where the functions ψ ω j S j ( j = 1 , , N ω ) are the result of step 1.

2. (2.2) Similarly, calculate the associated EM field of the second class, ( B ρ , B z , E ϕ ) with E ρ = E z = B ϕ = 0 , 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 B ρ = B z = E ϕ = 0 . “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.

### 3 Numerical implementation

#### 3.1 Precise object of the fitting

The fitting of the sum (20) is done to determine the “wave vector spectra” S j in equation (16):

(21) Ψ ( ω j ) ( S j ) ( t , ρ , z ) = j ψ ω j S j ( t , ρ , z ) ,

where [equation (7) with ω = ω j and K j = ω j / c ]

(22) ψ ω j S j ( t , ρ , z ) = e i ω j t K j + K j J 0 ρ K j 2 k 2 e i k z S j ( k ) d k .

Thus, we have one spectrum S j for each value of the index j = 1 , , N ω , the latter specifying the frequency ω j . To determine these spectra, several methods could be a priori envisaged. However, a difficulty comes from the huge ratio

(23) galactic distances wavelength kpc μ m 3 × 1 0 25 ,

which is the order of magnitude of the arguments of the Bessel function J 0 and the complex exponential in equation (22). This huge number discards several possibilities. First, it turns out to be not tractable at all here to determine each S j by its Fourier coefficients C j n – as proposed (for a very different problem) by Garay-Avendaño and Zamboni-Rached [13]. Indeed, considering (to begin with) one axisymmetric time-harmonic solution (7) of the scalar wave equation, this method leads to the following expansion equation (8) in ref. [13] :

(24) ψ ω S ( t , ρ , z ) = 2 K e i ω t n = C n sinc ( h ω n ( ρ , z ) ) ,

with K = ω / c and

(25) h ω n ( ρ , z ) = K 2 ρ 2 + ( z K + π n ) 2 .

Introducing the wavelength λ = 2 π c / ω = 2 π / K , we have from (25):

(26) h ω n ( ρ , z ) 2 = π 2 2 ρ λ 2 + 2 z λ + n 2 .

On the r.h.s. of equation (26), ρ / λ and z / λ have the huge magnitude (23). Therefore, the functions h ω n , hence also the functions sinc ( h ω n ) which form the basis in the expansion (24), are practically independent of n for relevant values of the spatial variables ρ and z , unless n would take similarly huge values. In that case, presumably, an integer of an akin value should give the number of the different values n to be taken in order to have an accurate expansion (24) – which of course is not tractable. Anyway, this method has been tried in this work and has not allowed us to get an accurate fitting of the sum (20).

In order to determine the “spectra” S j 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 t , a formal inverse Fourier transform gives us

(27) S ( k ) = 1 2 π J 0 ρ K 2 k 2 + ψ ( ρ , z ) e i k z d z .

This should thus be independent of ρ when e i ω t ψ ( ρ , z ) is indeed an axisymmetric time-harmonic solution, with frequency ω , of the wave equation. What we found numerically using the Matlab software is that, for the solution (18) corresponding to a spherical wave, with the frequency ω 0 being defined in Section 4: (i) the spectrum S obtained by equation (27) depends on ρ , which it should not. (ii) For a given value of ρ , S ( k ) varies with k K , + K instead of being the constant S ( k ) c 2 ω = 1 2 K . (iii) For values of ρ in the investigated range ( 1 0 3 kpc , , 1 0 2 kpc ), S ( k ) is much smaller than that value 1 2 K – by a factor of at least 1 0 5 to 1 0 10 at given ρ , and depending on ρ . We tried also to determine S ( k ) using an inverse Fourier transform starting from the EM field instead of the EM potential A z : e.g., starting from the component E z [equation (15) in the time-harmonic case]. This also was not successful in the numerical application with the relevant numbers. Another possibility could be to investigate an approach based on wavelets (e.g., [18]). We did not study the feasibility of such an approach for the present problem.

Instead, the method we finally used to determine S j is by the values

(28) S n j S j ( k n j ) ( n = 0 , , N , j = 1 , , N ω ) ,

where

(29) k n j = K j + n δ j ( n = 0 , , N )

is a regular discretization of the integration interval [ K j , + K j ] for k , corresponding with the frequency ω j , equation (22). (We remind that K j ω j / c ; moreover, δ j 2 K j / N is the size of the discretization interval.) Using the so-called “Simpson 3 8 composite rule,” integrals like the one in equation (22) are approximated by discrete sums, as follows:

(30) K j + K j f ( k ) d k = n = 0 N a n j f ( k n j ) + O 1 N 4 ,

where N must be a multiple of 3, and

(31) a n j = ( 3 / 8 ) δ j ( n = 0 or n = N ) ,

(32) a n j = 2 × ( 3 / 8 ) δ j ( mod ( n , 3 ) = 0 and n 0 and n N ) ,

(33) a n j = 3 × ( 3 / 8 ) δ j otherwise .

Using the approximation (30) to calculate Ψ [equations (21) and (22)], we get:

(34) Ψ ( t , ρ , z ) = n = 0 N j = 1 N ω f n j ( t , ρ , z ) S n j + O 1 N 4 ,

with

(35) f n j ( t , ρ , z ) = ω j ω 0 a n J 0 ρ ω j ω 0 K 0 2 k n 2 exp i ω j ω 0 k n z ω j t ,

where the real numbers a n 0 and k n ( 0 n N ) are as a n j and k n j in equations (31) and (29), replacing K j by K 0 = ω 0 c , so that

a n j = ω j ω 0 a n , k n j = ω j ω 0 k n .

To determine the spectra S j , i.e., the unknown complex numbers (28), we fit the sum of the scalar potentials of the individual “stars,” equation (20), by equation (34). To do that, we evaluate the sum (20) at a discrete set G of values of t , ρ and z , that makes a regular three-dimensional grid of points of spacetime:

(36) G = { ( t l , ρ m , z p ) , 1 l N t , 1 m N ρ , 1 p N z } .

We group the indices l , m , p as a single index J = J ( l , m , p ) ( 1 J J max = N t × N ρ × N z ) . We denote the corresponding values of the sum (20) by D J :

(37) D J = i = 1 i max j = 1 N ω S j exp ( i ( K j r i m p ω j t l ) ) K j r i m p , r i m p x ( ρ m , ϕ = 0 , z p ) x i .

(Recall that x ( ρ , ϕ , z ) is the spatial point with cylindrical coordinates ( ρ , ϕ , z ) .) Similarly, we denote A J n j = f n j ( t l , ρ m , z p ) . The fitting of the sum (20) by equation (34) on the spatiotemporal grid G amounts to solving the linear system

(38) n = 0 N j = 1 N ω A J n j S n j = D J ( J = 1 , , J max )

in the sense of the least squares, hence getting the complex numbers S n j as the output. These calculations are implemented on a PC, using the Matlab language and software.

#### 3.2 Quadruple precision is needed

The ratio in equation (23), thus a number of the order of 1 0 25 , gives the magnitude of the spatial argument K j r i in the complex exponential e i θ of equation (20) and the argument of the Bessel function J 0 in equation (22). However, the Bessel J 0 function, as well as the real and imaginary parts of e i θ , oscillate around 0 with a pseudo-period which is of the order of unity (exactly 2 π , for cos θ and sin θ ). So already to get only the correct sign, one needs to know their arguments to a precision better than O ( 1 ) . In view of the magnitude of the arguments: O ( 1 0 25 ) , it means that 25 significant digits are needed to know just the sign of cos K j r i and sin K j r i in equation (20) and the sign of J 0 ρ ω j 2 c 2 k 2 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 vpa (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 vpa with default precision, i.e., 32 digits plus 9 “guard digits.”[4]

#### 3.3 Calculation of the EM field and its exact character

As with equation (21) and (22) for the A z potential: each of B ϕ , E ρ , and E z is the sum of the time-harmonic components given by equations (13)–(15). And as we did with A z to obtain equations (34) and (35), we use the approximation (30) to calculate the integrals in equations (13)–(15). We thus get:

(39) B ϕ ( t , ρ , z ) = n = 0 N j = 1 N ω R n J 1 ρ ω j ω 0 R n ℛe F n j ( t , z ) + O 1 N 4 ,

(40) E ρ ( t , ρ , z ) = n = 0 N j = 1 N ω c 2 ω 0 k n R n J 1 ρ ω j ω 0 R n ℛe F n j ( t , z ) + O 1 N 4 ,

(41) E z ( t , ρ , z ) = n = 0 N j = 1 N ω c 2 ω 0 k n 2 ω 0 J 0 ρ ω j ω 0 R n ℐm F n j ( t , z ) + O 1 N 4 ,

with R n = K 0 2 k n 2 and

(42) F n j ( t , z ) = ω j ω 0 2 a n exp i ω j ω 0 k n z ω j t S n j .

Suppose one starts from an exact solution, with a finite frequency spectrum ( ω j ) , of the scalar wave equation: A z = Ψ ( ω j ) ( S j ) given by equations (21) and (22) with exact “wave-vector spectra” S j . 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 S j .

However, the integration formula (30) is exact (i.e., the remainder is not only O ( 1 / N 4 ) but exactly zero) if the integrand function f is a polynomial of degree 3. This is because the remainder is proportional to f ( 4 ) ( ξ ) for some ξ K j , + K j [19]. Moreover, given e.g. the first four arguments k n j and the four corresponding function values S n j ( n = 0 , , 3 ), there exists one and only one polynomial P of degree 3, such that P ( k n j ) = S n j ( n = 0 , , 3 ). (This is the well-known “unisolvence theorem;” note that we are considering a fixed value of the frequency index j .) By construction, N is a multiple of 3, hence the whole integration interval [ K j , + K j ] is the union of N / 3 adjacent subintervals, each covering three steps δ j . Thus, due to the unisolvence theorem: in each of those subintervals, the four successive arguments k n j and the four corresponding function values S n j 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 S ˜ j which continuously extends those N / 3 polynomial functions to the whole interval [ K j , + K j ] .[5] Therefore, equations (34) and (35) actually define an exact solution Ψ of the scalar wave equation, which corresponds with substituting the functions S ˜ j for S j 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 ω j , when in these contributions one considers, for the frequency ω = ω j , the spectrum function S = S ˜ j . 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 S ˜ j , which are piecewise third-degree polynomials.

#### 3.4 Validation test

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 A z = ψ ω S c 2 ω the spherically symmetric time-harmonic solution (18) of the scalar wave equation. This yields (after taking the real part):

(43) B ϕ ( t , ρ , z ) = ρ cos ω t r 2 sin K r K r cos K r ,

(44) E ρ ( t , ρ , z ) = c 2 ω ρ z r 3 sin ω t 3 cos K r r + 3 K r 2 K sin K r ,

(45) E z ( t , ρ , z ) = c 2 ω sin ω t r 1 3 z 2 r 2 cos K r r + 3 z 2 K r 3 1 + K 2 z 2 K r + K r sin K r r ,

with K ω / c and r ρ 2 + z 2 . The outputs of the (exact) equations (43)–(45) were compared to those of equations (39)–(41) applied with the relevant constant spectrum S ( k ) c 2 ω , thus in equation (42):

(46) S n S ( k n ) = c 2 ω ( n = 0 , , N ) .

(We are considering the single angular frequency case: N ω = 1 , hence the index j is omitted since it takes only one value: j = 1 ; moreover, we set ω 1 ω 0 ω in equations (39)–(42).) Note that, with the exact spectrum S ( k ) c 2 ω , equations (13)–(15) provide just the same exact fields as do equations (43)–(45). However, even with the exact spectrum values (46), equations (39)–(41) provide only an approximation of those exact fields, due to the discrete integration (30).

Figures 1, 2 and 3 show, for the three different scales investigated, the relative average quadratic differences between the fields Ψ = A z , B ϕ , E ρ , E z , 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: scale = 1 0 n λ ( n = 1 , 2 , 3 ), with λ c / ν . For this test, the frequency was taken to be ν ω / ( 2 π ) = 100 MHz , thus λ = 3 m . The average quadratic differences are in general evaluated on a regular three-dimensional spatio-temporal grid (36) for the variables t , ρ , z . Thus, t takes N t values between t 0 and t 0 + ( N t 1 ) δ t , ρ takes N ρ values between ρ 0 and ρ 0 + ( N ρ 1 ) δ ρ , and z takes N z values between z 0 and z 0 + ( N z 1 ) δ z , with δ t = T / N t , δ ρ = scale / N ρ , and δ z = scale / ( 10 N z ) . However, in the present case with harmonic time dependence, we took N t = 1 , thus one value of t = t 0 , which was fixed at T / 8 , with T = 1 / ν = 1 0 8 s . So here the grid is two-dimensional in fact. Also, we took here N ρ = 14 , N z = 13 , and ρ 0 = δ ρ , z 0 = 5 δ z .

### Figure 1

Average quadratic differences on a ( ρ , z ) grid. Scale: 10 λ .

### Figure 2

Average quadratic differences on a ( ρ , z ) grid. Scale: 1 0 2 λ .

### Figure 3

Average quadratic differences on a ( ρ , z ) grid. Scale: 1 0 3 λ .

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 S ( k ) becomes finer, i.e., with increasing N . (See equation (28).) This validates the correctness of our calculations. However, the errors increase quickly when the scale length scale 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 k that oscillate with a frequency or a pseudo-frequency which is proportional to the magnitude of the spatial variables ρ and z . As these integrals are approximated by the discrete sums (34) and (39)–(41), one would have to increase the discretization number N in proportion of the scale length scale .

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 x i . 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 Ψ ( x i ) ( ω j ) ( S j ) 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 x i are denoted by B ϕ i , E ρ i and E z i . In the case of a single frequency ω , one first applies equations (4)–(6) with A z = ψ ω as given by equation (19):

(47) B ϕ ω ( t , ρ , z ) = ρ r 2 sin φ + cos φ K r , φ K r ω t ,

(48) E ρ ω ( t , ρ , z ) = c 2 ω ρ z r 3 3 cos φ r + K 3 K r 2 sin φ ,

(49) E z ω ( t , ρ , z ) = c 2 ω r 3 z 2 r 2 1 cos φ r + 3 z 2 K r 3 + 1 + K 2 z 2 K r K r sin φ r ,

but with ρ , z and r being replaced by

(50) ρ i ( x x i ) 2 + ( y y i ) 2 ,

(51) Z i z z i ,

(52) r i x x i = ( x x i ) 2 + ( y y i ) 2 + ( z z i ) 2 = ρ i 2 + Z i 2 .

The definition of B ϕ i , E ρ i and E z i is that the exact fields produced at x by the source at x i are decomposed on the orthonormal direct basis made of

(53) e ρ i = e ρ ( x ; x i ) = ( x x i ) e x + ( y y i ) e y / ρ i ,

(54) e ϕ i = e ϕ ( x ; x i ) e z e ρ ( x ; x i ) = ( x x i ) e y ( y y i ) e x / ρ i ,

and e z . When each spherical source has the same finite frequency spectrum ( ω j , S j ) , as in equation (20), it is understood that the components B ϕ i , E ρ i and E z i involve the corresponding weighted sum, e.g.,

(55) B ϕ i = j = 1 N ω S j B ϕ i ω j ,

where B ϕ i ω j is given by equation (47) with ρ i , Z i , r i , ω j , K j in the place of ρ , z , r , ω , K . The total exact fields, sums of these different fields, are (reminding that B ρ i = B z i = E ϕ i = 0 ):

(56) B ϕ B.e ϕ = i B ϕ i e ϕ i .e ϕ ,

(57) E ρ E.e ρ = i E ρ i e ρ i .e ρ ,

(58) E z E.e z = i E z i .

In general, the other components of the total fields may be non-zero also, but we assume that the distribution of the identical spherical sources is axisymmetric (see Section 2.2). In that case, we expect that E ϕ = B ρ = B z = 0 . Note also that e ρ and e ϕ , hence also the scalar products in equations (56)–(57), and therefore also the B ϕ and E ρ components, are not defined if ρ = 0 .

### 4 Results and discussion

A file of some 1 0 4 randomly generated “stars” (as described in Section 2.2) has been used here, more precisely one with 16 × 16 × 36 triplets ( ρ , z , ϕ ) .

The angular frequencies ω j 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 ω 0 , thus

(59) ω j = ω 0 Δ ω + ( j 1 ) Δ ω N inter , ( j = 1 , , N ω = 2 N inter + 1 ) ,

where ω 0 = 2 π c λ 0 , Δ ω < ω 0 . In the calculations, we took N inter = 5 (hence N ω = 11 ), λ 0 0.5 μ m , Δ ω = ω 0 / 2 . The weights S j affected to the different frequencies in equation (20) have the form

(60) S j ω j exp 1 2 σ 2 ω j ω 0 2 ( j = 1 , , N ω )

with σ = Δ ω , and are normed so that j S j = 1 .

A few different spacetime domains (variables t , ρ , z ) of galactic dimensions have been used. The adopted sizes of the domain for the calculations discussed here were as follows:

(61) 0 t < T 0 λ 0 c ,

(62) ρ 0 ρ < scale ,

(63) scale / 20 < z < scale / 20 ,

with scale = 3.086 × 1 0 20 m 10 kpc and ρ 0 = scale / 1 0 6 in these calculations.

We have discretized that domain using grids (36) with N t = 4 , N ρ = 8 , N z = 7 – 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 t begins with t = 0 and ends with the largest multiple of the time step that is smaller than T 0 . The time step is δ t = T 0 / N t , so t = ( i t 1 ) δ t ( i t = 1 , , N t ) – and similarly for ρ and z in equations (62)–(63). (This is just the same kind of variation as for the validation test of Section 3.4, but here scale 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 A z potential by calculating the sum (20) with a very large time step, close to δ ρ / c . 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 N in equations (28)–(34): N = 6 , 12 , 24 , 48 , 96 , 192 .

We will compare the field components B ϕ , E ρ , E z , 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 S n j 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 ( 5 , 12 , 11 ) spatiotemporal grid, the contour levels of the field components B ϕ , E ρ , E z , as calculated either “directly” or “from the model” – in the latter case, on the same spatiotemporal grid ( 5 , 12 , 11 ) used for the fitting, and with N = 48 . 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 ( 4 , 8 , 7 ) or ( 7 , 14 , 13 ) .

### Figure 4

B ϕ calculated directly or from the model; t = 2 T 0 / 5 .

### Figure 5

B ϕ calculated directly or from the model; t = 4 T 0 / 5 .

### Figure 6

E ρ calculated directly or from the model; t = 2 T 0 / 5 .

### Figure 7

E ρ calculated directly or from the model; t = 3 T 0 / 5 .

### Figure 8

E z calculated directly or from the model; t = T 0 / 5 .

### Figure 9

E z calculated directly or from the model; t = 3 T 0 / 5 .

If the increase of the error with scale, found for the smaller scales investigated in the validation test of Section 3.4 (up to 1 0 3 λ ), would continue up to the scale relevant to a typical disk galaxy ( 3 × 1 0 25 λ ), then the relative quadratic errors between the field components calculated either directly or from the model (e.g., δ B ϕ / B ϕ ) 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 N 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 S n j 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 A z potential:

1. The fields are more intense close to the z 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 B ϕ , E ρ and E z also for ρ = 0 – which is not the case for the B ϕ and E ρ components when one uses the direct calculation based on equations (56)–(58), see after those equations. The model gives B ϕ = E ρ = 0 for ρ = 0 , because J 1 ( 0 ) = 0 . In contrast, for E z , that model predicts very high values for ρ = 0 , of the order E z = O ( 10 ) 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 ρ 0 = 0 instead of ρ 0 = scale / 1 0 6 .

2. 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 z axis.

An important point has to be noted in this connection. As we saw, the fields have a definite spatial variation at the galactic scale, in contrast to their quasi-periodic time variation with a very small time period T 0 = λ 0 / c . However, to that large-scale spatial variation, is superimposed an oscillatory variation at the very small scale of the main wavelength λ 0 = 0.5 μ m . This results again from the analytical expressions of the fields, e.g., it is easy to see in equations (39)–(41) for the model.[6] The oscillatory spatial variation at the wavelength scale implies a high sensitivity of the details of the calculations on a big spatiotemporal grid like (61)–(63) to small variations of the parameters. This is seen, for example, when the spectrum functions S j obtained from a fitting are used to calculate the fields on a different spatiotemporal grid than the one used for the fitting; or, when two different orders N are used for the same grid. However, the main features of the calculations, as described at points (i) and (ii) above, are robust, and the relative quadratic differences between different calculations on the same grid usually remain of the order of unity. An exception is if a too low order N is used (e.g., N = 6 used to fit and calculate with the model on the ( 5 , 12 , 11 ) grid), in which case larger differences can exist.

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 x = x i ). 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.

### 5 Conclusion

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 z axis, and there the E z component dominates over E ρ . 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].

## Acknowledgements

The author is grateful to Christian Boily and to Garrelt Mellema for useful remarks at the Computational Astrophysics Conference in Saint Petersburg, September 2019.

1. Conflict of interest: The author states that he has no conflict of interest.

#### References

[1] 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. 10.1007/978-94-007-5612-0_13Search in Google Scholar

[2] 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. 10.1093/mnras/sts297Search in Google Scholar

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

[4] Maciel WJ . The interstellar radiation field. In: Astrophysics of the interstellar medium. Maciel WJ , editor. New York: Springer; Chapter 2, 2013. p. 17–31. 10.1007/978-1-4614-3767-3_2Search in Google Scholar

[5] 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. 10.1515/phys-2018-0065Search in Google Scholar

[6] 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. 10.7546/giq-20-2019-88-98Search in Google Scholar

[7] Arminjon M . Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation. Open Phys. 2016;14:395–409. 10.1515/phys-2016-0045Search in Google Scholar

[8] 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. 10.1086/170413Search in Google Scholar

[9] Robin AC , Crézé M , Mohan V . The radial structure of the Galactic disc. Astron Astrophys. 1992;265:32–9. 10.1063/1.43976Search in Google Scholar

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

[11] Schneider P . Extragalactic astronomy and cosmology: an introduction. Berlin: Springer; 2016. p. 55. 10.1007/978-3-642-54083-7Search in Google Scholar

[12] 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. 10.1002/9780470168981.ch2Search in Google Scholar

[13] 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. 10.1364/AO.53.004524Search in Google Scholar PubMed

[14] Arminjon M. An explicit representation for the axisymmetric solutions of the free Maxwell’s equations. Open Phys. 2020;18:255–63. 10.1515/phys-2020-0117Search in Google Scholar

[15] Jackson JD . Classical electrodynamics. 3rd ed. Hoboken: John Wiley & Sons; 1998. p. 360. Search in Google Scholar

[16] Durnin J. Exact solutions for nondiffracting beams. I. The scalar theory. J Opt Soc Am A. 1987;4:651–4. 10.1364/JOSAA.4.000651Search in Google Scholar

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

[18] 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. 10.2528/PIERL11090505Search in Google Scholar

[19] Atkinson KE . An introduction to numerical analysis. 2nd Edition. New York: John Wiley & Sons; 1989. p. 257–8. Search in Google Scholar