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

Mayeul Arminjon

► To cite this version:

Mayeul Arminjon. An analytical model for the Maxwell radiation field in an axially symmetric galaxy. Open Physics, 2021, 19, pp. 77-90. 10.1515/phys-2021-0008. hal-02887969v2

HAL Id: hal-02887969

https://hal.science/hal-02887969v2

Submitted on 19 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An analytical model for the Maxwell radiation field in an axially symmetric galaxy

Mayeul Arminjon

Univ. Grenoble Alpes, CNRS, Grenoble INP, 3SR, F-38000 Grenoble, France

E-mail: Mayeul.Arminjon@3sr-grenoble.fr

Short title: *Maxwell radiation field in an axisymmetric galaxy*

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

**Keywords:** Disk galaxy; Maxwell equations; axial symmetry; exact solutions; numerical model; 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 towards 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]. 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 [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 able later
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 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. 4. Finally, we
present our conclusions.

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, 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 [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 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.

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

► An exponential distribution is assumed for $\rho > 0$, i.e., $n_\rho$ values $\rho_m$ ($m = 1, ..., n_\rho$) are got by quasi-random generation with a probability:

$$P(a < \rho_m < b) = \frac{1}{h} \int_a^b e^{-\frac{\rho}{h}} d\rho, \qquad (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 [8, 9, 10].

► Also, an exponential distribution is assumed for $z > 0$: for any $m = 1, ..., n_\rho$, 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) = \frac{1}{h_z} \int_a^b e^{-\frac{z}{h_z}} dz, \qquad (2)$$

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

► 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_\rho$ and $p = 1, ..., n_z$, we draw $n_\phi$ values $\phi_{mpq}$, with a uniform distribution between 0 and $2\pi$ (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 $\Psi(t, \rho, z) = e^{-i\omega t} \tilde{\Psi}(\rho, z)$ of the scalar wave equation, and one associates with
it a vector potential A by

$$A := \Psi e_z, \quad \text{or} \quad A_z := \Psi, \quad A_\rho = A_\phi = 0. \qquad (3)$$

(We shall denote by ($e_\rho, e_\phi, e_z$) the standard, point-dependent, direct orthonormal basis associated with the cylindrical coordinates $\rho, \phi, 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<fnref n="1" />):

$$B_{\phi} = -\frac{\partial A_z}{\partial \rho}, \quad E_{\phi} = 0, \qquad (4)$$

$$E_{\rho} = i \frac{c^2}{\omega} \frac{\partial^2 A_z}{\partial \rho \partial z}, \quad B_{\rho} = 0, \qquad (5)$$

$$E_z = i \frac{c^2}{\omega} \frac{\partial^2 A_z}{\partial z^2} + i\omega A_z, \quad B_z = 0. \qquad (6)$$

In Ref. [13], this class was defined only when the scalar wave $\Psi$ has the following form:

$$\psi_{\omega} S(t, \rho, z) = e^{-i\omega t} \int_{-K}^{+K} J_0 \left( \rho\sqrt{K^2 - k^2} \right) e^{ikz} S(k) dk, \qquad (7)$$

with $\omega$ the angular frequency, $K := \omega/c$,<fnref n="2" /> and $J_0$ the first-kind Bessel function of order 0. ($c$ is the velocity of light.) Any *totally propagating*,<fnref n="3" /> time-harmonic, axisymmetric solution of the scalar wave equation can be set

$$D_{\omega k k'} (t, x, y, z) := e^{i(kz - \omega t)} \int_{0}^{2\pi} \exp[ik'(x \cos \phi + y \sin \phi)] \frac{d\phi}{2\pi} \qquad (8)$$

$$\equiv e^{i(kz - \omega t)} J_0(k'\rho) \quad (\rho := \sqrt{x^2 + y^2}), \qquad (9)$$

<footnote n="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, \rho, z$, Eqs. (4)–(6) define a solution of the source-free Maxwell equations in that same subdomain [13, 14].</footnote>

<footnote n="2">² In Ref. [13], $K$ was defined as $K := 2\omega/c$ instead.</footnote>

<footnote n="3">³ This is the terminology used in Ref. [13]. We can give the following justification. Eq. (7) is an integral over $k$ of “Bessel beams”, each having the form [16]</footnote>
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 \in 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, \rho, z$, thus in the domain $\Omega$: $(t, z \in \mathbb{R}, \rho \in \mathbb{R}_+)$ (i.e. $\rho \ge 0$) [12, 14]. Therefore, given any frequency $\omega > 0$ and any integrable function $S(k)$, Eqs. (4)–(6), applied with $A_z = \psi_\omega S$ given by Eq. (7), define an axisymmetric solution of the source-free Maxwell equations in the whole spatiotemporal domain, $\Omega$.

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

$$E' = cB, \quad B' = -E/c. \qquad (11)$$

In the work [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 with

$$k^2 + k'^2 = K^2 := \omega^2/c^2. \qquad (10)$$

Eq. (8) is the average over the angle $\phi$ of plane waves $W_\phi$ having $k_x = k' \cos \phi$, $k_y = k' \sin \phi$, and $k_z = k$, and all having the same angular frequency $\omega$. However, in Eq. (8), $k$ and $k'$ may be complex and are only subjected to the constraint (10). For any such pair $(k, k')$, Eq. (8) defines an axisymmetric solution of the scalar wave equation, and may be rewritten as (9). In particular, in the case $k$ real with $|k| > K$, we get $k' = \pm i\sqrt{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 \phi + y \sin \phi$. Thus, if in Eq. (7) one replaces the integration bounds $\pm K$ by $\pm K'$, with $K' > K$, then, depending on whether $|k| \le K$ or $|k| > K$, $D_{\omega k k'}$ will be an average of either propagating plane waves (with a real wave vector) or “evanescent” ones. Whereas, leaving the bounds $\pm K$ means considering a totally propagating solution, i.e., as stated in Ref. [13], one that does not contain “evanescent” waves.
of the source-free Maxwell equations (whether totally propagating or not). There exist a unique solution ($\mathbf{E}_1, \mathbf{B}_1$) of the first class (4)-(6) and a unique solution ($\mathbf{E}'_2, \mathbf{B}'_2$) of the second class, both with the same frequency as has ($\mathbf{A}, \mathbf{E}, \mathbf{B}$), and whose sum gives just that solution:

$$ \mathbf{E} = \mathbf{E}_1 + \mathbf{E}'_2, \quad \mathbf{B} = \mathbf{B}_1 + \mathbf{B}'_2. \qquad (12) $$

Of course, as usual, it is implicit that, in Eqs. (4)-(6), $B_\phi$, $E_\rho$ 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 $dJ_0/dx = -J_1(x)):

$$ B_{\phi \omega S} = \mathcal{R}e \left[ e^{-i\omega t} \int_{-K}^{+K} \sqrt{K^2 - k^2} J_1 \left( \rho \sqrt{K^2 - k^2} \right) S(k) e^{ikz} dk \right], \quad (13) $$

$$ E_{\rho \omega S} = \mathcal{R}e \left[ -i \frac{c^2}{\omega} e^{-i\omega t} \int_{-K}^{+K} \sqrt{K^2 - k^2} J_1 \left( \rho \sqrt{K^2 - k^2} \right) ik S(k) e^{ikz} dk \right], \quad (14) $$

$$ E_{z \omega S} = \mathcal{R}e \left[ i e^{-i\omega t} \int_{-K}^{+K} J_0 \left( \rho \sqrt{K^2 - k^2} \right) \left( \omega - \frac{c^2}{\omega} k^2 \right) S(k) e^{ikz} dk \right], \quad (15) $$

where $K := \omega/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 ($\omega_j$) ($j = 1, ..., N_\omega$) is got by summing solutions of the form (7):

$$ \Psi_{(\omega_j) (S_j)} (t, \rho, z) = \sum_{j=1}^{N_\omega} \psi_{\omega_j S_j} (t, \rho, z), \quad (16) $$

where each $\psi_{\omega_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

$$S(k) = \frac{c}{2\omega} \quad (-K < k < K := \frac{\omega}{c}), \qquad (17)$$

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

$$\psi_{\omega} S \equiv \frac{c}{2\omega} (t, \rho, z) = e^{-i\omega t} \mathrm{sinc}(Kr), \quad K := \frac{\omega}{c}, \qquad (18)$$

where $\mathrm{sinc} \theta := \frac{\sin \theta}{\theta}$, $r := |\mathbf{x}| = \sqrt{\rho^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 $\omega$, and of course it is well known:

$$\psi'_{\omega} (t, \rho, z) = \frac{e^{i(Kr - \omega t)}}{Kr}, \quad K := \frac{\omega}{c}. \qquad (19)$$

Obviously, this is a totally propagating solution. Nevertheless, unfortunately, it does not seem amenable to the integral form (7) with an analytical function $S(k)$ — in contrast to Eq. (18). 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 $\mathbf{x}_i$, all sources having the same amplitude and the same frequency spectrum, (19) becomes:

$$\Psi_{(\mathbf{x}_i)} (\omega_j) (S'_j) (t, \rho, z) = \sum_{i=1}^{i_{\max}} \sum_{j=1}^{N_{\omega}} S'_j \frac{e^{i(K_j r_i - \omega_j t)}}{K_j r_i}, \qquad (20)$$

where $r_i := |\mathbf{x} - \mathbf{x}_i|$, $K_j := \omega_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 $\omega_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}$.
## 2.5 The model

* Step 1. Determine a relevant axisymmetric solution $\Psi$ 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_\rho, E_z, B_\phi)$ with $B_\rho = B_z = E_\phi = 0$, by summing its time-harmonic contributions given by Eqs. (4)–(6), with successively $A_z := \psi_{\omega_j s_j} (j = 1, \dots, N_\omega)$, where the functions $\psi_{\omega_j s_j} (j = 1, \dots, N_\omega)$ are the result of step 1.

* (2.2) Similarly, calculate the associated EM field of the second class, $(B'_\rho, B'_z, E'_\phi)$ with $E'_\rho = E'_z = B'_\phi = 0$, by using the EM duality (11).

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_\rho = B_z = E_\phi = 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 Eq. (16):

$$ \Psi_{(\omega_j) (S_j)} (t, \rho, z) = \sum_j \psi_{\omega_j S_j} (t, \rho, z), \quad (21) $$
where [Eq. (7) with $\omega = \omega_j$ and $K_j = \omega_j/c$]

$$ \psi_{\omega_j} S_j(t, \rho, z) = e^{-i\omega_j t} \int_{-K_j}^{+K_j} J_0 \left( \rho \sqrt{K_j^2 - k^2} \right) e^{ikz} S_j(k) \, dk. \quad (22) $$

Thus, we have one spectrum $S_j$ for each value of the index $j = 1, \dots, N_\omega$, the latter specifying the frequency $\omega_j$. To determine these spectra, several methods could be a priori envisaged. However, a difficulty comes from the huge ratio

$$ \frac{\text{galactic distances}}{\text{wavelength}} \approx \frac{\text{kpc}}{\mu\text{m}} \approx 3 \times 10^{25}, \quad (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 [13]. Indeed, considering (to begin with) one axisymmetric time-harmonic solution (7) of the scalar wave equation, this method leads to the following expansion {Eq. (8) in Ref. [13]}:

$$ \psi_{\omega} S(t, \rho, z) = 2Ke^{-i\omega t} \sum_{n=-\infty}^{\infty} C_n \text{sinc}(h_{\omega n}(\rho, z)), \quad (24) $$

with $K = \omega/c$ and

$$ h_{\omega n}(\rho, z) = \sqrt{K^2 \rho^2 + (zK + \pi n)^2}. \quad (25) $$

Introducing the wavelength $\lambda = 2\pi c/\omega = 2\pi/K$, we have from (25):

$$ (h_{\omega n}(\rho, z)))^2 = \pi^2 \left[ \left( \frac{2\rho}{\lambda} \right)^2 + \left( \frac{2z}{\lambda} + n \right)^2 \right]. \quad (26) $$

On the r.h.s. of Eq. (26), $\rho/\lambda$ and $z/\lambda$ have the huge magnitude (23). Therefore, the functions $h_{\omega n}$, hence also the functions $\text{sinc}(h_{\omega n})$ which form the basis in the expansion (24), are practically independent of $n$ for relevant values of the spatial variables $\rho$ 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 (7) of the scalar wave equation, and removing its dependence in $t$, a formal inverse Fourier transform gives us

$$S(k) = \frac{1}{2\pi J_0(\rho\sqrt{K^2 - k^2})} \int_{-\infty}^{+\infty} \psi(\rho, z)e^{-ikz}dz. \quad (27)$$

This should thus be independent of $\rho$ when $e^{-i\omega t}\psi(\rho, z)$ is indeed an axisymmetric time-harmonic solution, with frequency $\omega$, 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 $\omega_0$ being defined in Sect. 4 below: i) The spectrum $S$ obtained by Eq. (27) depends on $\rho$, which it shouldn’t. ii) For a given value of $\rho$, $S(k)$ varies with $k \in ]-K, +K[$ instead of being the constant $S(k) \equiv \frac{c}{2\omega} = \frac{1}{2K}$. iii) For values of $\rho$ in the investigated range ($10^{-3}$ kpc, ..., $10^2$ kpc), $S(k)$ is much smaller than that value $\frac{1}{2K}$ — by a factor of at least $10^5$ to $10^{10}$ at given $\rho$, and depending on $\rho$. 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. (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

$$S_{nj} := S_j(k_{nj}) \quad (n = 0, \dots, N, \quad j = 1, \dots, N_{\omega}), \quad (28)$$

where

$$k_{nj} = -K_j + n\delta_j \quad (n = 0, \dots, N) \quad (29)$$

is a regular discretization of the integration interval $[-K_j, +K_j]$ for $k$, corresponding with the frequency $\omega_j$, Eq. (22). (We remind that $K_j := \omega_j/c$; moreover, $\delta_j := 2K_j/N$ is the size of the discretization interval.) Using the so-called “Simpson $\frac{3}{8}$ composite rule”, integrals like the one in Eq. (22) are
approximated by discrete sums, as follows:

$$ \int_{-K_j}^{+K_j} f(k)dk = \sum_{n=0}^{N} a_{nj}f(k_{nj}) + O\left(\frac{1}{N^4}\right), \qquad (30) $$

where N must be a multiple of 3, and

$$ a_{nj} = (3/8) \delta_j \quad (n = 0 \text{ or } n = N), \qquad (31) $$

$$ a_{nj} = 2 \times (3/8) \delta_j \quad (\text{mod}(n, 3) = 0 \text{ and } n \neq 0 \text{ and } n \neq N), \qquad (32) $$

$$ a_{nj} = 3 \times (3/8) \delta_j \quad \text{otherwise.} \qquad (33) $$

Using the approximation (30) to calculate $\Psi$ [Eqs. (21) and (22)], we get:

$$ \Psi(t, \rho, z) = \sum_{n=0}^{N} \sum_{j=1}^{N_{\omega}} f_{nj}(t, \rho, z) S_{nj} + O\left(\frac{1}{N^4}\right), \qquad (34) $$

with

$$ f_{nj}(t, \rho, z) = \frac{\omega_j}{\omega_0} a_n J_0 \left( \rho \frac{\omega_j}{\omega_0} \sqrt{K_0^2 - k_n^2} \right) \exp \left[ i \left( \frac{\omega_j}{\omega_0} k_n z - \omega_j t \right) \right], \qquad (35) $$

where the real numbers $a_n \ge 0$ and $k_n$ ($0 \le n \le N$) are as $a_{nj}$ and $k_{nj}$ in Eqs. (31) and (29), replacing $K_j$ by $K_0 = \frac{\omega_0}{c}$, so that

$$ a_{nj} = \frac{\omega_j}{\omega_0} a_n, \quad k_{nj} = \frac{\omega_j}{\omega_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$, $\rho$ and $z$, that makes a regular three-dimensional grid of points of spacetime:

$$ G = \{(t_l, \rho_m, z_p), \ 1 \le l \le N_t, \ 1 \le m \le N_\rho, \ 1 \le p \le N_z\}. \qquad (36) $$

We group the indices $l, m, p$ as a single index $J = J(l, m, p)$ ($1 \le J \le J_{\text{max}} = N_t \times N_\rho \times N_z$). We denote the corresponding values of the sum (20) by $D_J$:

$$ D_J = \sum_{i=1}^{i_{\text{max}}} \sum_{j=1}^{N_\omega} S'_j \frac{\exp(i(K_j r_{\text{imp}} - \omega_j t_l))}{K_j r_{\text{imp}}} , \quad r_{\text{imp}} := |\mathbf{x}(\rho_m, \phi = 0, z_p) - \mathbf{x}_i|. \qquad (37) $$
(Recall that $\mathbf{x}(\rho, \phi, z)$ is the spatial point with cylindrical coordinates $(\rho, \phi, z)$.) Similarly, we denote $A_{J nj} = f_{nj}(t_l, \rho_m, z_p)$. The fitting of the sum (20) by Eq. (34) on the spatiotemporal grid $G$ amounts to solving the linear system

$$ \sum_{n=0}^{N} \sum_{j=1}^{N_{\omega}} A_{J nj} S_{nj} = D_{J} \quad (J = 1, \dots, J_{\max}) \qquad (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.

## 3.2 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\theta}$ 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\theta}$, oscillate around 0 with a pseudo-period which is of the order of unity (exactly $2\pi$, for $\cos\theta$ and $\sin\theta$). 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(\rho\sqrt{\frac{\omega_j^2}{c^2} - k^2})$ in Eq. (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”.<fnref n="4" />

<footnote n="4">⁴ This was advised to the author by Pavel Holoborodko, of Advanpix.</footnote>
### 3.3 Calculation of the EM field and its exact character

As with Eq. (21) and (22) for the $A_z$ potential: each of $B_\phi$, $E_\rho$, and $E_z$ is the sum of the time-harmonic components given by Eqs. (13)–(15). 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)–(15). We thus get:

$$B_{\phi}(t, \rho, z) = \sum_{n=0}^{N} \sum_{j=1}^{N_{\omega}} R_n J_1 \left( \rho \frac{\omega_j}{\omega_0} R_n \right) \mathcal{R}e[F_{nj}(t, z)] + O \left( \frac{1}{N^4} \right), \quad (39)$$

$$E_{\rho}(t, \rho, z) = \sum_{n=0}^{N} \sum_{j=1}^{N_{\omega}} \frac{c^2}{\omega_0} k_n R_n J_1 \left( \rho \frac{\omega_j}{\omega_0} R_n \right) \mathcal{R}e[F_{nj}(t, z)] + O \left( \frac{1}{N^4} \right), \quad (40)$$

$$E_z(t, \rho, z) = \sum_{n=0}^{N} \sum_{j=1}^{N_{\omega}} \left( \frac{c^2}{\omega_0} k_n^2 - \omega_0 \right) J_0 \left( \rho \frac{\omega_j}{\omega_0} R_n \right) \mathcal{I}m[F_{nj}(t, z)] + O \left( \frac{1}{N^4} \right), \quad (41)$$

with $R_n = \sqrt{K_0^2 - k_n^2}$ and

$$F_{nj}(t, z) = \left(\frac{\omega_j}{\omega_0}\right)^2 a_n \exp \left[ i \left( \frac{\omega_j}{\omega_0} k_n z - \omega_j t \right) \right] S_{nj}. \quad (42)$$

Suppose one starts from an exact solution, with a finite frequency spectrum ($\omega_j$), of the scalar wave equation: $A_z = \Psi_{(\omega_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 $\le 3$. This is because the remainder is proportional to $f^{(4)}(\xi)$ for some $\xi \in ]-K_j, +K_j[$ [19]. 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 $\le 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 [-Kj, +Kj] 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 knj and the four corresponding function values Snj 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 $\tilde{S}_j$ which continuously extends those N/3 polynomial functions to the whole interval [-Kj, +Kj].<fnref n="5" /> Therefore, Eqs. (34) and (35) actually define an exact solution $\tilde{\Psi}$ of the scalar wave equation, which corresponds with substituting the functions $\tilde{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 $\omega_j$, when in these contributions one considers, for the frequency $\omega = \omega_j$, the spectrum function $S = \tilde{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 $\tilde{S}_j$, which are piecewise 3rd-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 Eqs. (4)–(6), with $A_z = \psi_\omega S \equiv \frac{c}{2\omega}$ the spherically-symmetric time-harmonic solution (18) of the scalar wave equation. This yields (after taking the real part):

$$B_{\phi}(t, \rho, z) = \frac{\rho \cos \omega t}{r^2} \left( \frac{\sin Kr}{Kr} - \cos Kr \right), \quad (43)$$

<footnote n="5">⁵ The continuity of $\tilde{S}_j$ results from the common value $S_{nj}$ at the common bound of any two successive subintervals. (The derivatives of $\tilde{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 $\delta_j$, because the weights $a_{nj} = (3/8)\delta_j$ at the bounds of two successive subintervals add to give the weight $a_{nj} = 2 \times (3/8)\delta_j$ at a multiple of three steps inside the whole interval.</footnote>
$$E_{\rho}(t, \rho, z) = \frac{c^2}{\omega} \frac{\rho z}{r^3} \sin \omega t \left[ -\frac{3 \cos Kr}{r} + \left( \frac{3}{Kr^2} - K \right) \sin Kr \right], \quad (44)$$

$$E_z(t, \rho, z) = \frac{c^2}{\omega} \frac{\sin \omega t}{r} \left[ \left(1 - \frac{3z^2}{r^2}\right) \frac{\cos Kr}{r} + \left(\frac{3z^2}{Kr^3} - \frac{1+K^2z^2}{Kr} + Kr\right) \frac{\sin Kr}{r} \right], \quad (45)$$

with $K := \omega/c$ and $r := \sqrt{\rho^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) := \frac{c}{2\omega}$, thus in Eq. (42):

$$S_n := S(k_n) = \frac{c}{2\omega} \quad (n = 0, \dots, N). \qquad (46)$$

(We are considering the single angular frequency case: $N_\omega = 1$, hence the index $j$ is omitted since it takes only one value: $j = 1$; moreover, we set $\omega_1 := \omega_0 := \omega$ in Eqs. (39)–(42).) Note that, with the exact spectrum $S(k) := \frac{c}{2\omega}$, 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 $\Psi = A_z, B_\phi, E_\rho, 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\lambda$ ($n = 1, 2, 3$), with $\lambda := c/\nu$. For this test, the frequency was taken to be $\nu := \omega/(2\pi) = 100$ MHz, thus $\lambda = 3$ m. The average quadratic differences are in general evaluated on a regular three-dimensional spatio-temporal grid (36) for the variables $t, \rho, z$. Thus $t$ takes $N_t$ values between $t_0$ and $t_0 + (N_t - 1)\delta t$, $\rho$ takes $N_\rho$ values between $\rho_0$ and $\rho_0 + (N_\rho - 1)\delta \rho$, and $z$ takes $N_z$ values between $z_0$ and $z_0 + (N_z - 1)\delta z$, with $\delta t = T/N_t$, $\delta \rho = \text{scale}/N_\rho$, and $\delta z = \text{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/\nu = 10^{-8}$ s. So here the grid is two-dimensional in fact. Also, we took here $N_\rho = 14, N_z = 13$, and $\rho_0 = \delta \rho, z_0 = 5\delta z$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Average quadratic differences on a $(\rho, z)$ grid. Scale: 10$\lambda$</figcaption></figure>

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
<figure><img src="image_3.png" /><figcaption>Figure 2: Average quadratic differences on a $(\rho, z)$ grid. Scale: $10^2 \lambda$</figcaption></figure>

the source-free Maxwell equations by summing solutions obtained by Eqs. (4)–(6) applied to the outgoing spherical wave (19). That sum is similar to the sum $\Psi_{(x_i)}(\omega_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_{\phi i}$, $E_{\rho i}$ and $E_{z i}$. In the case of a single frequency $\omega$, one first applies Eqs. (4)–(6) with $A_z = \psi'_\omega$ as given by Eq. (19):

$$B_{\phi\omega}(t, \rho, z) = \frac{\rho}{r^2} \left( \sin \varphi + \frac{\cos \varphi}{Kr} \right), \quad \varphi := Kr - \omega t, \quad (47)$$

$$E_{\rho\omega}(t, \rho, z) = \frac{c^2 \rho z}{\omega r^3} \left[ 3 \frac{\cos \varphi}{r} + \left( K - \frac{3}{Kr^2} \right) \sin \varphi \right], \quad (48)$$

$$E_{z\omega}(t, \rho, z) = \frac{c^2}{\omega r} \left[ \left( \frac{3z^2}{r^2} - 1 \right) \frac{\cos \varphi}{r} + \left( -\frac{3z^2}{Kr^3} + \frac{1+K^2z^2}{Kr} - Kr \right) \frac{\sin \varphi}{r} \right], \quad (49)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Average quadratic differences on a $(\rho, z)$ grid. Scale: $10^3 \lambda$</figcaption></figure>

but with $\rho$, $z$ and $r$ being replaced by

$$ \rho'_{i} := \sqrt{(x - x_{i})^{2} + (y - y_{i})^{2}}, \qquad (50) $$

$$ Z_{i} := z - z_{i}, \qquad (51) $$

$$ r_{i} := |\mathbf{x} - \mathbf{x}_{i}| = \sqrt{(x - x_{i})^{2} + (y - y_{i})^{2} + (z - z_{i})^{2}} = \sqrt{\rho_{i}^{2} + Z_{i}^{2}}. \qquad (52) $$

The definition of $B_{\phi i}$, $E_{\rho i}$ and $E_{z i}$ is that the exact fields produced at $\mathbf{x}$ by the source at $\mathbf{x}_i$ are decomposed on the orthonormal direct basis made of

$$ \mathbf{e}'_{\rho i} = \mathbf{e}'_{\rho}(\mathbf{x}; \mathbf{x}_i) = ((x - x_i)\mathbf{e}_x + (y - y_i)\mathbf{e}_y) / \rho'_i, \qquad (53) $$

$$ \mathbf{e}'_{\phi i} = \mathbf{e}'_{\phi}(\mathbf{x}; \mathbf{x}_i) := \mathbf{e}_z \wedge \mathbf{e}'_{\rho}(\mathbf{x}; \mathbf{x}_i) = ((x - x_i)\mathbf{e}_y - (y - y_i)\mathbf{e}_x) / \rho'_i, \qquad (54) $$
and $e_z$. When each spherical source has the same finite frequency spectrum ($\omega_j, S'_j$), as in Eq. (20), it is understood that the components $B_{\phi i}, E_{\rho i}$ and $E_{zi}$ involve the corresponding weighted sum, e.g.

$$B_{\phi i} = \sum_{j=1}^{N_{\omega}} S'_{j} B_{\phi i \omega_j}, \qquad (55)$$

where $B_{\phi i \omega_j}$ is given by Eq. (47) with $\rho'_i, Z_i, r_i, \omega_j, K_j$ in the place of $\rho, z, r, \omega, K$. The total exact fields, sums of these different fields, are (reminding that $B_{\rho i} = B_{zi} = E_{\phi i} = 0$):

$$B_{\phi} := \mathbf{B} \cdot \mathbf{e}_{\phi} = \sum_{i} B_{\phi i} e'_{\phi i} \cdot \mathbf{e}_{\phi}, \qquad (56)$$

$$E_{\rho} := \mathbf{E} \cdot \mathbf{e}_{\rho} = \sum_{i} E_{\rho i} e'_{\rho i} \cdot \mathbf{e}_{\rho}, \qquad (57)$$

$$E_z := \mathbf{E} \cdot \mathbf{e}_z = \sum_i E_{zi}. \qquad (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_{\phi} = B_{\rho} = B_z = 0$. Note also that $e_{\rho}$ and $e_{\phi}$, hence also the scalar products in Eqs. (56)–(57), and therefore also the $B_{\phi}$ and $E_{\rho}$ components, are not defined if $\rho = 0$.

# 4 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 \times 16 \times 36$ triplets $(\rho, z, \phi)$.

The angular frequencies $\omega_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 $\omega_0$, thus

$$\omega_j = \omega_0 - \Delta\omega + (j - 1) \frac{\Delta\omega}{N_{\text{inter}}}, \quad (j = 1, \dots, N_{\omega} = 2N_{\text{inter}} + 1), \qquad (59)$$

where $\omega_0 = \frac{2\pi c}{\lambda_0}$, $\Delta\omega < \omega_0$. In the calculations, we took $N_{\text{inter}} = 5$ (hence $N_{\omega} = 11$), $\lambda_0 := 0.5 \, \mu\text{m}$, $\Delta\omega = \omega_0/2$. The weights $S'_j$ affected to the different
frequencies in Eq. (20) have the form

$$S'_{j} \propto \omega_{j} \exp \left[ -\frac{1}{2\sigma^{2}} (\omega_{j} - \omega_{0})^{2} \right] \quad (j = 1, \dots, N_{\omega}) \quad (60)$$

with $\sigma = \Delta\omega$, and are normed so that $\sum_j S'_j = 1$.

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

$$0 \le t < T_0 := \frac{\lambda_0}{c}, \quad (61)$$

$$\rho_0 \le \rho < \text{scale}, \quad (62)$$

$$-\text{scale}/20 < z < \text{scale}/20, \quad (63)$$

with $\text{scale} = 3.086 \times 10^{20} \text{ m} \approx 10 \text{ kpc}$ and $\rho_0 = \text{scale}/10^6$ in these calculations.

We have discretized that domain using grids (36) with $N_t = 4$, $N_\rho = 8$, $N_z = 7$ — in short (4, 8, 7) —, or (5, 12, 11), or (7, 14, 13). The symbols $\le$ 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 $\delta t = T_0/N_t$, so $t = (i_t - 1)\delta t$ ($i_t = 1, \dots, N_t$) — and similarly for $\rho$ 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 $\text{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 $\delta\rho/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_\phi, E_\rho, 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,
the field derives from exact solutions of the scalar wave equation by Eqs. (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 to 9 show, for the (5, 12, 11) spatiotemporal grid, the contour levels of the field components $B_\phi$, $E_\rho$, $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).

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\lambda$), would continue up to the scale relevant to a typical disc galaxy ($\simeq 3 \times 10^{25}\lambda$), then the relative quadratic errors between the field components calculated either directly or from the model (e.g. $|| \delta B_\phi || / || B_\phi ||$) 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, 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_\phi$, $E_\rho$ and $E_z$ also for $\rho = 0$ — which is not the case for the $B_\phi$ and $E_\rho$ components when one uses the direct calculation based on Eqs. (56)–(58), see after those equations. The model gives $B_\phi = E_\rho = 0$ for $\rho = 0$, because $J_1(0) = 0$. In contrast, for $E_z$, that model predicts very high values for $\rho = 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 fol-
lowing Eqs. (61)–(63) — but using these spectrum values to calculate the
fields on a shifted grid, with $\rho$ starting at $\rho_0 = 0$ instead of $\rho_0 = \text{scale}/10^6$.

ii) The maximum intensities (positive and negative), calculated either di-
rectly 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 = \lambda_0/c$.
However, to that large-scale spatial variation, is superimposed an oscillatory
variation at the very small scale of the main wavelength $\lambda_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. <fnref n="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 param-
eters. 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 Eq. (3) from a
spherically symmetric solution $\Psi$ 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 axisym-
metric solution of the source-free Maxwell equations. The set of point-like

<footnote n="6">⁶ However, the same remark applies to the direct calculation, e.g. Eq. (56) for $B_\phi$, with $B_{\phi i}$ given by Eq. (55), where $B_{\phi i \omega_j}$ is given by Eq. (43) with $\rho'_i, Z_i, r_i, \omega_j, K_j$ in the place of $\rho, z, r, \omega, K$.</footnote>
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 $\mathbf{x} = \mathbf{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.

# 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 Eqs. (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.** I am grateful to Christian Boily and to Garrelt Mellema for useful remarks at the Computational Astrophysics Conference in Saint Petersburg, September 2019.

References

[1] Beck R, Wielebinski R. Magnetic fields in the Milky Way and in galaxies. In: Planets, Stars and Stellar Systems. Oswalt TD, Gilmore G, editors, vol. 5, Dordrecht: Springer, 2013, pp. 641–723

[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–3589.

[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; 2005; Vol. 4, pp. 77-80.

[4] Maciel WJ. The interstellar radiation field. In: Astrophysics of the interstellar medium. Maciel WJ, editor. New York: Springer; 2013; Chapter 2, pp. 17-31.

[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–498.

[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; pp. 88-98.

[7] Arminjon M. Continuum dynamics and the electromagnetic field in the scalar ether theory of gravitation. Open Physics 2016;14:395–409.

[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–138.

[9] Robin AC, Crézé M, Mohan V. The radial structure of the Galactic disc. Astron. Astrophys. 1992;265:32–39.
[10] Porcel C, Garzón F, Jiménez-Vicente J. The radial scale length of the Milky Way. Astron. Astrophys. 1998;330:136-138.

[11] Schneider P. Extragalactic astronomy and cosmology: an introduction. Berlin: Springer; p. 55.

[12] Zamboni-Rached M, Recami E, Hernández-Figueroa H E. Structure of non-diffracting waves and some interesting applications. In: Localized waves; Hernández-Figueroa HE, Zamboni-Rached, Recami E, editors. Localized waves. Hoboken: John Wiley & Sons; 2008; pp. 43-77.

[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-4531.

[14] Arminjon M. An explicit representation for the axisymmetric solutions of the free Maxwell equations. Open Physics 2020;18:255-263.

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

[16] Durnin J. Exact solutions for nondiffracting beams. I. The scalar theory. J. Opt. Soc. Am. A 1987;4:651-654.

[17] Gradshteyn IS, Ryzhik IM. Table of Integrals, Series, and Products. 7th English Edition. Burlington (Mass., USA): Academic Press; 2007; §6.677, p. 722.

[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-140.

[19] Atkinson KE. An Introduction to Numerical Analysis. 2nd Edition. New York: John Wiley & Sons; 1989; pp. 257-258.
<figure><img src="image_5.png" /><figcaption>Figure 4: $B_{\phi}$ calculated directly or from the model; $t = 2T_{0}/5$</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: $B_{\phi}$ calculated directly or from the model; $t = 4T_{0}/5$</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: $E_ρ$ calculated directly or from the model; $t = 2T_0/5$</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: $E_ρ$ calculated directly or from the model; $t = 3T_0/5$</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: $E_z$ calculated directly or from the model; $t = T_0/5$</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: $E_z$ calculated directly or from the model; $t = 3T_0/5$</figcaption></figure>