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

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

Keywords: Disk galaxy, Maxwell equations, axial symmetry, exact solutions, numerical model, dark matter

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

Introduction

Apart from pure magnetic fields, which are thought to be produced by a galactic dynamo action [START_REF] Beck | Magnetic fields in the Milky Way and in galaxies[END_REF][START_REF] Chamandy | Galactic spiral patterns and dynamo action I: a new twist on magnetic arms[END_REF], 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. [START_REF] Porter | A new estimate of the galactic interstellar radiation field between 0.1µm and 1000µm[END_REF][START_REF] Maciel | The interstellar radiation field[END_REF]. 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 towards testing the following prediction [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF] 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. [START_REF] Arminjon | Interaction energy of a charged medium and its EM field in a curved spacetime[END_REF] summarizes Ref. [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_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. [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_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 [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF][START_REF] Arminjon | Interaction energy of a charged medium and its EM field in a curved spacetime[END_REF]. 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 [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF][START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF]. I.e., 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 [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF]. (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 able later to calculate the interaction energy predicted by the scalar theory [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF], 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 paper 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 Sect. [START_REF] Maciel | The interstellar radiation field[END_REF]. Finally, we present our conclusions.

2 Description of the model

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, Eq. ( 19), and the corresponding sum is given by Eq. (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 [START_REF] Maciel | The interstellar radiation field[END_REF], "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 the Introduction 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, scattering, etc.

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

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 , φ mpq , z mp )} (that is, i = i(m, p, q)), is got by random generation, as follows:

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

P (a < ρ m < b) = 1 h b a e -ρ h dρ, ( 1
)
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 [START_REF] Kent | Galactic structure from the Spacelab infrared telescope. II. Luminosity models of the Milky Way[END_REF][START_REF] Robin | The radial structure of the Galactic disc[END_REF][START_REF] Porcel | The radial scale length of the Milky Way[END_REF].

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

P (a < z mp < b) = 1 h z b a e -z hz dz, (2)
with h z = 0.2 kpc in the numerical computations. This also roughly corresponds to the Milky Way [START_REF] Schneider | Extragalactic astronomy and cosmology: an introduction[END_REF].

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

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

Explicit representation for axisymmetric EM radiation fields

In Ref. [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF], two classes of time-harmonic axisymmetric solutions of the Maxwell equations were introduced. Although the aim of that paper [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF] 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

A := Ψe z , or A z := Ψ, A ρ = A φ = 0. (3)
(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 [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF] and, in more detail, in Ref. [START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF] -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 ) :

B φ = - ∂A z ∂ρ , E φ = 0, (4)
E ρ = i c 2 ω ∂ 2 A z ∂ρ ∂z , B ρ = 0, (5)
E z = i c 2 ω ∂ 2 A z ∂z 2 + iωA z , B z = 0. (6)
In Ref. [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF], this class was defined only when the scalar wave Ψ has the following form:

ψ ω S (t, ρ, z) = e -iωt +K -K J 0 ρ √ K 2 -k 2 e ik z S(k) dk, (7)
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 1 We are indeed willing to describe the EM radiation field on a galactic spatial scale, thus excluding from consideration the (much smaller) stars of the galaxy, which are the primary source of that radiation field. Given any axisymmetric solution A z of the standard wave equation in some open domain for the spacetime variables t, ρ, z, Eqs. ( 4)-( 6) define a solution of the source-free Maxwell equations in that same subdomain [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF][START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF]. 2 In Ref. [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF], K was defined as K := 2ω/c instead.

3 This is the terminology used in Ref. [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]. We can give the following justification. Eq. ( 7) is an integral over k of "Bessel beams", each having the form [START_REF] Durnin | Exact solutions for nondiffracting beams. I. The scalar theory[END_REF]

D ω k k (t, x, y, z) := e i(kz-ωt) 2π 0 exp[ik (x cos φ + y sin φ)] dφ 2π (8)
≡ e i(kz-ωt) J 0 (k ρ)

(ρ := x 2 + y 2 ), (9)
in the form (7) [START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]. However, as noted in Ref. [START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF], 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]), Eq. ( 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) [START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF]. Therefore, given any frequency ω > 0 and any integrable function S(k), Eqs. ( 4)-( 6), applied with A z = ψ ω S given by Eq. ( 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 [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]:

E = cB, B = -E/c. (11)
In the work [START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_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 [START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF]. Let (A, E, B) be any time-harmonic axisymmetric solution

with k 2 + k 2 = K 2 := ω 2 /c 2 . (10)
Eq. ( 8) is the average over the angle φ of plane waves W φ having k x = k cos φ, k y = k sin φ, and k z = k, and all having the same angular frequency ω. However, in Eq. ( 8), k and k may be complex and are only subjected to the constraint [START_REF] Porcel | The radial scale length of the Milky Way[END_REF]. For any such pair (k, k ), Eq. ( 8) defines an axisymmetric solution of the scalar wave equation, and may be rewritten as [START_REF] Robin | The radial structure of the Galactic disc[END_REF]. In particular, in the case k real with |k| > K, we get k = ±i √ k 2 -K 2 , so that then Eq. ( 8) is an angular average of "evanescent" plane waves, i.e., the amplitude of which depends exponentially on x cos φ + y sin φ. Thus, if in Eq. ( 7) one replaces the integration bounds ±K by ±K , with K > K, then, depending on whether |k| ≤ K or |k| > K, D ω k k will be an average of either propagating plane waves (with a real wave vector) or "evanescent" ones. Whereas, leaving the bounds ±K means considering a totally propagating solution, i.e., as stated in Ref. [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF], one that does not contain "evanescent" waves.

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:

E = E 1 + E 2 , B = B 1 + B 2 . ( 12
)
Of course, as usual, it is implicit that, in Eqs. ( 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 [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF], we obtain (using the fact that dJ 0 /dx = -J 1 (x)):

B φ ω S = Re e -iωt +K -K √ K 2 -k 2 J 1 ρ √ K 2 -k 2 S(k) e ikz dk , (13)
E ρ ω S = Re -i c 2 ω e -iωt +K -K √ K 2 -k 2 J 1 ρ √ K 2 -k 2 ik S(k) e ikz dk , (14)
E z ω S = Re ie -iωt +K -K J 0 ρ √ K 2 -k 2 ω - c 2 ω k 2 S(k) e ikz dk , (15)
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 [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF]:

Ψ (ω j ) (S j ) (t, ρ, z) = Nω j=1 ψ ω j S j (t, ρ, z), (16)
where each ψ ω j S j is a time-harmonic solution having the form [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF]. 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 .

The case of spherical waves

If, in the axisymmetric time-harmonic solution [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF], one puts

S(k) ≡ c 2ω (-K < k < K := ω c ), (17)
then one gets a spherical time-harmonic solution of the scalar wave equation: This yields indeed [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF][START_REF] Gradshteyn | Table of Integrals, Series, and Products. 7th English Edition[END_REF]

ψ ω S≡ c 2ω (t, ρ, z) = e -iωt sinc (Kr), K := ω c , (18)
where sinc θ := sin θ θ , r := |x| = ρ 2 + z 2 . However, because sin Kr = (e iKr -e -iKr )/(2i), 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:

ψ ω (t, ρ, z) = e i(Kr-ωt) Kr , K := ω c . (19)
Obviously, this is a totally propagating solution. Nevertheless, unfortunately, it does not seem be amenable to the integral form [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF] with an analytical function S(k) -in contrast to Eq. [START_REF] Ala | A wavelet operator on the interval in solving Maxwell's equations[END_REF]. However, for a spherical source (as opposed to a sink), we must use Eq. ( 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:

Ψ (x i ) (ω j ) (S j ) (t, ρ, z) = imax i=1 Nω j=1 S j e i(K j r i -ω j t) K j r i , (20)
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 ij .

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.

(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 Eqs. ( 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) 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 [START_REF] Schneider | Extragalactic astronomy and cosmology: an introduction[END_REF].

Step 1 (especially) is delicate numerically, as we will see in the next section. Therefore, we shall focus in this paper 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

Precise object of the fitting

The fitting of the sum (20) is done to determine the "wave vector spectra" S j in Eq. ( 16):

Ψ (ω j ) (S j ) (t, ρ, z) = j ψ ω j S j (t, ρ, z), (21)
where [Eq. [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF] with ω = ω j and

K j = ω j /c] ψ ω j S j (t, ρ, z) = e -iω j t +K j -K j J 0 ρ K 2 j -k 2 e ik z S j (k) dk. ( 22
)
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 galactic distances wavelength

kpc µm 3 × 10 25 , (23)
which is the order of magnitude of the arguments of the Bessel function J 0 and the complex exponential in Eq. ( 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 jn -as proposed (for a very different problem) by Garay-Avendaño & Zamboni-Rached [START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]. Indeed, considering (to begin with) one axisymmetric time-harmonic solution [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF] of the scalar wave equation, this method leads to the following expansion {Eq. (8) in Ref.

[13]}:

ψ ω S (t, ρ, z) = 2Ke -iωt ∞ n=-∞ C n sinc(h ω n (ρ, z)), (24)
with K = ω/c and

h ω n (ρ, z) = K 2 ρ 2 + (zK + πn) 2 . ( 25
)
Introducing the wavelength λ = 2πc/ω = 2π/K, we have from (25):

(h ω n (ρ, z))) 2 = π 2 2ρ λ 2 + 2z λ + n 2 . ( 26
)
On the r.h.s. of Eq. ( 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 Eq. ( 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 [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF] of the scalar wave equation, and removing its dependence in t, a formal inverse Fourier transform gives us

S(k) = 1 2πJ 0 ρ √ K 2 -k 2 +∞ -∞ ψ(ρ, z)e -ikz dz. ( 27
)
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 [START_REF] Ala | A wavelet operator on the interval in solving Maxwell's equations[END_REF] corresponding to a spherical wave, with the frequency ω 0 being defined in Sect. 4 below: i) The spectrum S obtained by Eq. ( 27) depends on ρ, which it shouldn't. ii) For a given value of ρ, S(k)

varies with k ∈]-K, +K[ instead of being the constant S(k) ≡ c 2ω = 1 2K .
iii) For values of ρ in the investigated range (10 -3 kpc, ..., 10 2 kpc), S(k) is much smaller than that value 1 2K -by a factor of at least 10 5 to 10 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 [Eq. [START_REF] Jackson | Classical electrodynamics[END_REF] 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. [START_REF] Ala | A wavelet operator on the interval in solving Maxwell's equations[END_REF]). 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

S nj := S j (k nj ) (n = 0, ..., N, j = 1, ..., N ω ), (28)
where

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

+K j -K j f (k)dk = N n=0 a nj f (k nj ) + O 1 N 4 , (30)
where N must be a multiple of 3, and

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

a nj = 3 × (3/8) δ j otherwise. ( 33
)
Using the approximation (30) to calculate Ψ [Eqs. ( 21) and ( 22)], we get:

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

f nj (t, ρ, z) = ω j ω 0 a n J 0 ρ ω j ω 0 K 2 0 -k 2 n exp i ω j ω 0 k n z -ω j t , (35)
where the real numbers a n ≥ 0 and k n (0 ≤ n ≤ N ) are as a nj and k nj in Eqs. (31) and (29), replacing K j by K 0 = ω 0 c , so that

a nj = ω j ω 0 a n , k nj = ω 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", Eq. (20), by Eq. (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:

G = {(t l , ρ m , z p ), 1 ≤ l ≤ N t , 1 ≤ m ≤ N ρ , 1 ≤ p ≤ N z }. (36)
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 :

D J = imax i=1 Nω j=1 S j exp(i(K j r imp -ω j t l )) K j r imp , r imp := |x(ρ m , φ = 0, z p ) -x i | . (37)
(Recall that x(ρ, φ, z) is the spatial point with cylindrical coordinates (ρ, φ, z).) Similarly, we denote A J nj = f nj (t l , ρ m , z p ). The fitting of the sum (20) by Eq. (34) on the spatiotemporal grid G amounts to solving the linear system

N n=0 Nω j=1 A J nj S nj = D J (J = 1, ...J max ) (38)
in the sense of the least squares, hence getting the complex numbers S nj as the output. These calculations are implemented on a PC, using the Matlab language and software.

Quadruple precision is needed

The ratio in Eq. ( 23), thus a number of the order of 10 25 , gives the magnitude of the spatial argument K j r i in the complex exponential e iθ of Eq. ( 20) and the argument of the Bessel function J 0 in Eq. ( 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(10 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 Eq. ( 20) and the sign of J 0 ρ ω 2 j c 2 -k 2 in Eq. ( 22). Therefore, double precision (16 significant digits) is not enough: quadruple precision (32 significant digits) is neededand 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

Calculation of the EM field and its exact character

As with Eq. ( 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 Eqs. ( 13)- [START_REF] Jackson | Classical electrodynamics[END_REF]. And as we did with A z to obtain Eqs. (34) and (35), we use the approximation (30) to calculate the integrals in Eqs. ( 13)- [START_REF] Jackson | Classical electrodynamics[END_REF]. We thus get:

B φ (t, ρ, z) = N n=0 Nω j=1 R n J 1 ρ ω j ω 0 R n Re [F nj (t, z)] + O 1 N 4 , (39)
E ρ (t, ρ, z) = N n=0 Nω j=1 c 2 ω 0 k n R n J 1 ρ ω j ω 0 R n Re [F nj (t, z)] + O 1 N 4 , (40) E z (t, ρ, z) = N n=0 Nω j=1 c 2 ω 0 k 2 n -ω 0 J 0 ρ ω j ω 0 R n Im [F nj (t, z)] + O 1 N 4 , (41) with R n = K 2 0 -k 2 n and F nj (t, z) = ω j ω 0 2 a n exp i ω j ω 0 k n z -ω j t S nj . ( 42
)
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 Eqs. ( 21) and ( 22) with exact "wave-vector spectra" S j . Then Eqs. (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 [ [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]. Moreover, given e.g. the first four arguments k nj and the four corresponding function values S nj (n = 0, ..., 3), there exists one and only one polynomial P of degree ≤ 3, such that P (k nj ) = S nj (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 nj and the four corresponding function values S nj 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, Eqs. ( 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 Eqs. ( 21) and ( 22). Similarly, Eqs. ( 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, Eqs. ( 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 3rd-degree polynomials.

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 Eqs. ( 4)-( 6), with A z = ψ ω S≡ c 2ω the spherically-symmetric time-harmonic solution [START_REF] Ala | A wavelet operator on the interval in solving Maxwell's equations[END_REF] of the scalar wave equation. This yields (after taking the real part):

B φ (t, ρ, z) = ρ cos ωt r 2
sin Kr Kr -cos Kr , (43) 5 The continuity of S j results from the common value S nj at the common bound of any two successive subintervals. (The derivatives of S j are in general not continuous at the bounds of the subintervals, though.) Note that the integration formula (30) gives the same result whether it is applied to the whole interval, or successively to each of the subintervals, each covering three steps δ j , because the weights a nj = (3/8) δ j at the bounds of two successive subintervals add to give the weight a nj = 2 × (3/8) δ j at a multiple of three steps inside the whole interval.

E ρ (t, ρ, z) = c 2 ω ρz r 3 sin ωt - 3 cos Kr r + 3 Kr 2 -K sin Kr , (44)
E z (t, ρ, z) = c 2 ω sin ωt r 1 - 3z 2 r 2 cos Kr r + 3z 2 Kr 3 - 1 + K 2 z 2 Kr + Kr sin Kr r , (45
) with K := ω/c and r := ρ 2 + z 2 . The outputs of the (exact) Eqs. ( 43)-( 45) were compared to those of Eqs. ( 39)-( 41) applied with the relevant constant spectrum S(k) := c 2ω , thus in Eq. ( 42):

S n := S(k n ) = c 2ω (n = 0, ..., N ). (46)
(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 Eqs. ( 39)-(42).) Note that, with the exact spectrum S(k) := c 2ω , Eqs. ( 13)-( 15) provide just the same exact fields as do Eqs. ( 43)-(45). However, even with the exact spectrum values (46), Eqs. ( 39)-(41) provide only an approximation of those exact fields, due to the discrete integration (30).

Figs. 1 to 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 Eqs. ( 18) and ( 43)-(45), or "with the spectrum (46)", i.e., by Eqs. (34) and ( 39)-(41) applied with the spectrum values (46). The different scales were here: scale = 10 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/(10N 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/ν = 10 -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. 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 Eq. ( 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 Eqs. ( 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 Eqs. (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 4)-( 6) applied to the outgoing spherical wave [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]. That sum is similar to the sum Ψ (x i ) (ω j ) (S j ) given by Eq. ( 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 Eqs. ( 4)-( 6) with A z = ψ ω as given by Eq. ( 19):

B φω (t, ρ, z) = ρ r 2 sin ϕ + cos ϕ Kr , ϕ := Kr -ωt, (47)
E ρω (t, ρ, z) = c 2 ω ρz r 3 3 cos ϕ r + K - 3 Kr 2 sin ϕ , (48)
E zω (t, ρ, z) = c 2 ωr 3z 2 r 2 -1 cos ϕ r + - 3z 2 Kr 3 + 1 + K 2 z 2 Kr -Kr sin ϕ r , (49)
ρ i := (x -x i ) 2 + (y -y i ) 2 , ( 50
)
Z i := z -z i , (51)
r i := |x -x i | = (x -x i ) 2 + (y -y i ) 2 + (z -z i ) 2 = ρ 2 i + Z 2 i . (52)
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

e ρ i = e ρ (x; x i ) = ((x -x i )e x + (y -y i )e y ) /ρ i , (53)
e φ i = e φ (x; x i ) := e z ∧ e ρ (x;

x i ) = ((x -x i )e y -(y -y i )e x ) /ρ i , (54)
and e z . When each spherical source has the same finite frequency spectrum (ω j , S j ), as in Eq. ( 20), it is understood that the components B φ i , E ρ i and E z i involve the corresponding weighted sum, e.g.

B φ i = Nω j=1 S j B φ i ω j , (55)
where B φ i ω j is given by Eq. ( 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): B φ := B.e φ = i B φ i e φ i .e φ , (56)
E ρ := E.e ρ = i E ρ i e ρ i .e ρ , (57)
E z := E.e z = i E z i . (58)
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 Sect 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 Eqs. ( 56)-(57), and therefore also the B φ and E ρ components, are not defined if ρ = 0.

Results and discussion

A file of some 10 4 randomly generated "stars" (as described in Subsect. 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, Eq. (20), and in the analytical expression used to fit it, Eq. ( 21). They are regularly spaced and symmetric around a central frequency ω 0 , thus

ω j = ω 0 -∆ω + (j -1) ∆ω N inter , (j = 1, ..., N ω = 2N inter + 1), (59)
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 Eq. ( 20) have the form

S j ∝ ω j exp - 1 2σ 2 (ω j -ω 0 ) 2 (j = 1, ..., N ω ) (60)
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:

0 ≤ t < T 0 := λ 0 c , (61)
ρ 0 ≤ ρ < scale , ( 62
) -scale/20 < z < scale/20, (63)
with scale = 3.086 × 10 20 m 10 kpc and ρ 0 = scale/10 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 [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF][START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Schneider | Extragalactic astronomy and cosmology: an introduction[END_REF], or [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF][START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF][START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF]. The symbols ≤ and < in Eqs. (61)-(63) are meant to indicate that, e.g. for Eq. (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 Eqs. (62)-( 63). (This is just the same kind of variation as for the validation test of Sect. 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. Whereas, 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 Eqs. ( 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 Eqs. (56)-(58), or "from the model", i.e., from Eqs. (39)-(41), using in the latter case in Eq. ( 42) the spectrum values S nj obtained from fitting the sum (20) by Eq. (34), Eq. (38). Thus, in both cases, Figures 4 to 9 show, for the [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF][START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Schneider | Extragalactic astronomy and cosmology: an introduction[END_REF] 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 [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF][START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Schneider | Extragalactic astronomy and cosmology: an introduction[END_REF] 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 [START_REF] Arminjon | Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation[END_REF][START_REF] Arminjon | An explicit representation for the axisymmetric solutions of the free Maxwell equations[END_REF][START_REF] Garay-Avendaño | Exact analytic solutions of Maxwell's equations describing propagating nonparaxial electromagnetic beams[END_REF].

If the increase of the error with scale, found for the smaller scales investigated in the validation test of Subsect. 3.4 (up to 10 3 λ), would continue up to the scale relevant to a typical disc galaxy ( 3 × 10 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 [START_REF] Gradshteyn | Table of Integrals, Series, and Products. 7th English Edition[END_REF] was available, and the discretized spectrum values were taken from it, Eq. (46). In contrast, in this section, the discretized spectrum values S nj are now obtained from fitting the sum (20) by Eq. (34) over a spatiotemporal grid with galactic-like spatial dimensions, Eqs. (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:

i) 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 Eqs. (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 Eqs. (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 Eq. (34) on the spacetime domain (61)-(63), using e.g. the parameters described in the paragraph following Eqs. ( 61)-(63) -but using these spectrum values to calculate the fields on a shifted grid, with ρ starting at ρ 0 = 0 instead of ρ 0 = scale/106 .

ii) 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 Eqs. (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 [START_REF] Arminjon | On the equations of electrodynamics in a flat or a curved spacetime and a possible interaction energy[END_REF][START_REF] Zamboni-Rached | Structure of nondiffracting waves and some interesting applications[END_REF][START_REF] Schneider | Extragalactic astronomy and cosmology: an introduction[END_REF] 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 Eq. ( 3) from a spherically symmetric solution Ψ of the scalar wave equation. As we saw at the end of Subsect. 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 Eqs. ( 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 disc galaxy as the field calculated directly can be.

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 Eqs. ( 13)- [START_REF] Jackson | Classical electrodynamics[END_REF], 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 [START_REF] Porter | A new estimate of the galactic interstellar radiation field between 0.1µm and 1000µm[END_REF][START_REF] Maciel | The interstellar radiation field[END_REF].

Figure 1 :

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

Figure 2 :

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

Figure 3 :

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

Figure 8 :

 8 Figure 8: E z calculated directly or from the model; t = T 0 /5

This was advised to the author by Pavel Holoborodko, of Advanpix.

However, the same remark applies to the direct calculation, e.g. Eq. (56) for B φ , with B φ i given by Eq. (55), where B φ i ωj is given by Eq. (43) with ρ i , Z i , r i , ω j , K j in the place of ρ, z, r, ω, K.

Figure 4: B φ calculated directly or from the model; t = 2T 0 /5

Figure 5: B φ calculated directly or from the model; t = 4T 0 /5

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

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

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

Acknowledgements. I am grateful to Christian Boily and to Garrelt Mellema for useful remarks at the Computational Astrophysics Conference in Saint Petersburg, September 2019.