Interstellar radiation as a Maxwell field: improved numerical scheme and application to the spectral energy density

The interstellar radiation field (ISRF) in a galaxy is an electromagnetic (EM) field in a very high vacuum, hence it should be a solution of the Maxwell equations. However, the existing models for the ISRF do not take into account the full EM field with its six components coupled through the Maxwell equations. Consider, for example, the model of Chi & Wolfendale [1]. It assumes an axisymmetric distribution of the volume emissivities $j_{i}(\lambda,\rho,z)$ of four stellar components $(i)\ (i=1,...,4)$: $j_{i}$ decreases exponentially with both the distance $\rho$ to the galactic axis and the altitude $z$ over the galactic central disk. The contribution of component $(i)$ to the energy density of the ISRF at some position $(\rho^{\prime},z^{\prime})$ and wavelength $\lambda$ is obtained by integrating $j_{i}(\lambda,\rho,z)g/l^{2}$ over the whole galactic volume. Here $l$ is the distance between the studied position and the running point in the galactic volume; $g$ describes the dust absorption and is obtained by integrating the visual extinction per unit path length over the linear path joining the studied position and the running point in the galactic volume. Other models, e.g. by Mathis, Mezger and Panagia [2], Gordon et al. [3], Robitaille [4], Popescu et al. [5], are based on similar principles: all of these models consider quantities such as the stellar emissivity and luminosity, and the dust opacity, and they evolve the light intensity emitted by the stars by taking into account (in addition to the distance) the radiative transfer, in particular by dust absorption/ reemission. Clearly, those models do not produce an EM field, hence even less one that would be a solution of the Maxwell equations.

In a recent work [6], we proposed a model applicable to the relevant ideal case of an axisymmetric galaxy, and that provides for the ISRF such an exact solution of the Maxwell equations — a feature which, as discussed above, and to the best of our knowledge, appears to be fully new. This is indeed needed to study the relevance of a possible candidate for dark matter that emerges [7] from an alternative, scalar theory of gravity. However, it is also of astrophysical interest independently of the latter, since, as we noted, the ISRF must be an exact Maxwell field and this condition is not fulfilled by the existing models. As a step in checking the model proposed in Ref. [6], its application to predict the variation of the spectral energy density (SED) in our Galaxy has been subjected to a first test [8]. To this purpose, the model has been adjusted by asking that the SED predicted for our local position in the Galaxy coincide with the SED determined from spatial missions by Henry, Anderson & Fastie [9], Arendt et al. [10], Finkbeiner et al. [11], and Porter & Strong [12]. It has been found in that most recent work [8] that the spatial variation of the SED thus obtained with our model does not differ too much in magnitude from that predicted by the recent radiation transfer model of Ref. [5], but that the SED predicted by our model: (i) is extremely high on the axis of the Galaxy — i.e., on the axis of the axial symmetry that is assumed for the model of the Galaxy; (ii) has rather marked oscillations as function of the wavelength; and (iii) seems to decrease more slowly when the altitude $z$ increases (or rather when $\left|z\right|$ increases), as compared with the radiation transfer model.

The aim of this paper is to present an improved numerical scheme to operate that “Maxwell model of the ISRF”, and to apply this improved scheme to check the findings (i)–(iii) above. Section 2 provides a summary of the model. Section 3 describes the improvement of the numerical scheme. In Sect. 4, we check whether the model really predicts extremely high values of the SED on the axis of the Galaxy. Section 5 studies the spatial variation of the SED and compares it with results of the literature. In Sect. 6, asymptotic expansions are used to interpret the findings of the foregoing section. The Conclusion section 7 is followed by Appendix A, which discusses the relation between the discrete and continuous descriptions of the SED.

This model has been presented in detail in Ref. [6]. An axisymmetric galaxy is modelled as a finite set of point-like “stars”, the azimuthal distribution of which is uniform. Those points ${\bf x}_{i}\ (i=1,...,i_{\mathrm{max}})$ are obtained by pseudo-random generation of their cylindrical coordinates $\rho,\phi,z$ with specific probability laws, ensuring that the distribution of $\rho$ and $z$ is approximately that valid for the star distribution in the galaxy considered, and that the set $\{{\bf x}_{i}\}$ is approximately invariant under azimuthal rotations of any angle $\phi$ [6]. In the present work, as in Refs. [6, 8], $16\times 16\times 36$ triplets $(\rho,z,\phi)$ were thus generated, so that $i_{\mathrm{max}}=9216$, and the distribution of $\rho$ and $z$ is approximately that valid for the star distribution in the Milky Way.

The ISRF is also assumed axisymmetric, and thus depends only on $\rho$ and $z$. Since we want to describe, not the field inside the sources and in their vicinity, but instead the smoothed-out field at the intragalactic scale, we search for a solution of the source-free Maxwell equations. In the axisymmetric case, any time-harmonic source-free Maxwell field is the sum of two Maxwell fields: (i) one deriving from a vector potential having just the axial component $A_{z}$ non-zero, with $A_{z}$ obeying the standard wave equation, and (ii) one deduced from a solution of the form (i) by EM duality [13]. We consider for simplicity a model ISRF that has a finite frequency spectrum $(\omega_{j})_{j=1,...,N_{\omega}}$, hence we may apply the foregoing result to each among its time-harmonic components $(j)$, and then just sum these components. Moreover, we envisage the ISRF as being indeed an EM radiation field, thus excluding from consideration the purely magnetic part of the interstellar EM field [14]. Hence the ISRF is made of “totally propagating” EM waves, i.e., ones without any “evanescent” component [6, 16]. Specifically, we assume that the two scalar potentials $A_{j\,z}$ and $A^{\prime}_{j\,z}$ that define the decomposition (i)-(ii) of each time-harmonic component $(j)$, mentioned above, are themselves totally propagating. In that case, both $A_{j\,z}$ and $A^{\prime}_{j\,z}$ have the explicit form [15, 16]:

with $\omega_{j}$ the angular frequency, $K_{j}:=\omega_{j}/c$, $J_{0}$ the first-kind Bessel function of order $0$, and where $S_{j}$ is some (generally complex) function of $k\in[-K_{j},+K_{j}]$. For a totally propagating, axisymmetric EM field, but otherwise fully general, the two potentials $A_{j\,z}$ and $A^{\prime}_{j\,z}$ may be different, i.e., may correspond with different “spectra” in Eq. (1), say $S_{j}$ and $S^{\prime}_{j}$ [13].

To determine these potentials, that is, to determine the spectrum functions $S_{j}$, we use a sum of potentials emitted by the “stars”. We assume that every “star”, each at some point ${\bf x}_{i}$, contributes to the global potential $A_{j\,z}$ of a given frequency $\omega_{j}$ ($j=1,...,N_{\omega}$) by a spherically symmetric scalar wave of the same frequency $\omega_{j}$, whose emission center is its spatial position ${\bf x}_{i}$ — in order that all the directions starting from the star be equivalent. Thus, consider time-harmonic spherically symmetric solutions of the wave equation that have a given angular frequency $\omega$. It is easy to check by direct calculation that they can be either an outgoing wave, an ingoing wave, or the sum of an ingoing wave and an outgoing one, and that, up to an amplitude factor, the following is the only outgoing wave:

Clearly, only that outgoing solution is relevant here, given that the point-like “stars” must be indeed sources of radiation.  ¹¹1 The outgoing wave (2) can also be characterized [17] as being the only time-harmonic spherically symmetric solution of the wave equation that satisfies the Sommerfeld radiation condition [18]. However, the Sommerfeld condition aims precisely at selecting a boundary condition in order to find only “physical”, i.e., outgoing solutions for the Helmholtz equation. The latter equation applies to general time-harmonic solutions of the wave equation. In the spherically symmetric case, the time-harmonic solutions are easy to find and the outgoing solutions are immediate to recognize. Thus, the contributions of the $i$-th star to the potentials $A_{j\,z}$ and $A^{\prime}_{j\,z}$ can differ only in amplitude, since both must be a multiple of

where $\ K_{j}:=\frac{\omega_{j}}{c},\qquad r_{i}:=\left|{\bf x}-{\bf x}_{i}\right|$. But there is no apparent physical reason to affect different amplitudes to the contribution of the $i$-th star to $A_{j\,z}$ and to $A^{\prime}_{j\,z}$, hence we assume both of them to be equal to $\psi_{{\bf x}_{i}\,\omega_{j}}$. To determine the global potentials $A_{j\,z}$ and $A^{\prime}_{j\,z}$ $(j=1,...,N_{\omega})$, that generate the axisymmetric model ISRF with a finite frequency spectrum $(\omega_{j})$, the sum of the spherical potentials (3) emanating from the point stars is fitted to the form (1). As noted in Ref. [6], this is not assuming that the ISRF is indeed the sum of the radiation fields emitted by the different stars (which is not correct, due to the radiation transfers) — because (a) the equalities (4), (20) or (21) below are not exact equalities but ones in the sense of the least squares, and (b) nothing is really assumed regarding the EM field of the “star” itself, in particular we actually do not need to assume that it has the form (i)-(ii) above (e.g. the one corresponding with two equal potentials $A_{i\,j\,z}=A^{\prime}_{i\,j\,z}=\psi_{{\bf x}_{i}\,\omega_{j}}$).

In the previous works [6, 8], this fitting was done for all frequencies at once. That is, the following least-squares problem was considered:

where the sign $\cong$ indicates that the equality is in the sense of the least squares (the arguments of the functions varying on some spatio-temporal grid $G$), and where the numbers $w_{j}>0$ are the weights affected to the different frequencies. In view of the axial symmetry, the spatial position ${\bf x}$ is restricted to the plane $\phi=0$, so ${\bf x}={\bf x}(\rho,z)$ and

Since the contributions of the $i$-th star to $A_{j\,z}$ and to $A^{\prime}_{j\,z}$ have both been assumed to be equal to $\psi_{{\bf x}_{i}\,\omega_{j}}$, there is no possibility to distinguish between $A_{j\,z}$ and $A^{\prime}_{j\,z}$, either — whence $\psi_{\omega_{j}\ S_{j}}=A_{j\,z}=A^{\prime}_{j\,z}$ on the r.h.s. of (4). The unknowns of the problem are the spectrum functions $S_{j},\quad j=1,...,N_{\omega}$. We determine $S_{j}$ by the (generally complex) values

where

with $\delta_{j}:=2K_{j}/N$, is a regular discretization of the interval $[-K_{j},+K_{j}]$ for $k$ in the integral (1). Calculating those integrals with the “Simpson $\frac{3}{8}$ composite rule”, (4) becomes the computable least-squares problem

with

The $S_{nj}$ ’s are the solved-for parameters in the least-squares problem (8). In Eq. (8), $N$ must be a multiple of $3$, and in Eq. (9) we have

Part (i) of the decomposition of the model ISRF then obtains as follows [6]:

with  $R_{n}=\sqrt{K_{0}^{2}-k_{n}^{2}}$  and

(Here $k_{n}$ and $a_{n}$ ($0\leq n\leq N$) are as $k_{nj}$ and $a_{nj}$ in Eqs. (7) and (10), replacing $K_{j}$ by $K_{0}=\frac{\omega_{0}}{c}$, with $\omega_{0}$ some (arbitrary) reference frequency.) Since we assume $A_{j\,z}=A^{\prime}_{j\,z}$ for the global potentials generating the model ISRF, part (ii) of its decomposition is deduced from the first part by the EM duality:

It follows from this and from (13) that the model ISRF, sum of these two parts, has the components (14)–(16), and that the other components are just

Equation (4) may be split into the different frequencies (marked by the index $j$), simply by removing the sum on $j$ from both sides of either equation. The same is true for Eq. (8). Naturally, also the weight $w_{j}$ may then be removed from the l.h.s., by entering the inverse $1/w_{j}$ into the unknown spectrum function $S_{j}$ on the r.h.s. Equation (8) thus becomes

At this point, one notes that both $\psi_{{\bf x}_{i}\,\omega_{j}}$ [Eq. (3)] and $f_{nj}$ [Eq. (9)] have the same dependence on time, $\exp(-\mathrm{i}\omega_{j}t)$, which we can hence remove also, to obtain a least-squares problem with merely the spatial variables $\rho$ and $z$:

where $G^{\prime}=\{(\rho_{m},z_{p}),\ 1\leq m\leq N_{\rho},\ 1\leq p\leq N_{z}\}$ is the spatial grid, and

The separation, into the different frequencies, of the fitting of the sum of the potentials emitted by the “stars”, is consistent with the linearity of the wave equation and the Maxwell equations. Moreover, the elimination of the time variable from the fitting represents an appreciable gain in computing time. We recall that, for the EM field in a galaxy, the arguments of the Bessel function $J_{0}$ and the angular exponential, e.g. in Eq. (22), have the huge magnitude $\left|{\bf x}\right|/\lambda\sim 10^{25}$, which enforces us to use a precision better than quadruple precision in the computer programs, thus leading to slow calculations [6]. Note that the “separate fitting”, i.e. the least-squares problem (21), is not exactly equivalent to the “grouped fitting”, i.e. the least-squares problem (8) (this will be confirmed by the numerical results below): the two are slightly different ways of adjusting the global potentials (1). ²²2 If Eq. (20), or equivalently Eq. (21), were an exact equality instead of being an equality in the sense of the least squares, then of course it would imply Eq. (8) (with $w_{j}\equiv 1$) as an exact equality. However, equations (13)–(19) apply with the separate fitting as well — although the relevant values $S_{nj}$ are different. The separate fitting is more appropriate, because solutions corresponding with different frequencies behave independently in the Maxwell equations, and each frequency can be treated with more precision by considering it alone. Indeed a very important point is that, by switching to the separate fitting, we improve the situation regarding the “overfitting”, i.e., we decrease the ratio $R$ of the number of parameters $N_{\mathrm{para}}$ to the number of data $N_{\mathrm{data}}$: now, for each value of the frequency index $j$, we have to solve the least-squares problem (21), with $N_{\mathrm{para}}=N+1$ unknown parameters and $N_{\mathrm{data}}=N_{\rho}\times N_{z}$ data (the “data” are the values of the l.h.s. of (21) on the spatial grid $G^{\prime}$). Whereas, with the formerly used grouped fitting, we had to solve just one least-squares problem (8) with $N_{\mathrm{para}}=(N+1)\times N_{\omega}$ unknown parameters and $N_{\mathrm{data}}=N_{t}\times N_{\rho}\times N_{z}$ data.

On the other hand, through the processes of radiative transfer, there are indeed transfers of radiation intensity from some frequency domains to other ones, e.g. the interaction with dust leads to a transfer from higher to lower frequencies (see e.g. Fig. 3 in Ref. [3]). But these processes are not directly taken into account by the present model: not any more with the grouped fitting than with the separate fitting. They are indirectly taken into account through the adjustment of the energy density [8], which we briefly recall now.

The time-averaged volumic energy density of an EM field having a finite set of frequencies, $(\omega_{j})_{j=1,...,N_{\omega}}$, is given by [8]

where the complex numbers $C^{(q)}_{j}({\bf x})\ (q=1,...,6)$ are the coefficients in the expansion, in time-harmonic functions, of each among the six components of the EM field:

and where $\alpha_{q}=\epsilon_{0}$ for an electric field component, whereas $\alpha_{q}=\epsilon_{0}c^{2}$ for a magnetic field component (here $\epsilon_{0}$ is the vacuum permittivity, with $\epsilon_{0}=1/(4\pi\times 9\times 10^{9})$ in SI units). For an axisymmetric EM field, it is enough to consider the plane $\phi=0$, thus ${\bf x}={\bf x}(\rho,z)$, and we have

Using in that case the decomposition (i)-(ii), the expressions of three among the $C^{(q)}_{j}$ coefficients follow directly from the expressions (14)–(16) of the corresponding components of the EM field [8]. Moreover, in the special subcase (18) considered here, the other components are given by (19), whence in the same way the three remaining $C^{(q)}_{j}$ coefficients.

Now note that, in the least-squares problem (21), that we use to determine the values $S_{nj}$ allowing to compute the EM field (14)–(16) and (19), no data relative to the intensity of the fields emitted by the point-like “stars” has been used until now. Hence, we may multiply the l.h.s. of (21) by some number $\xi_{j}>0$, thus obtaining now new values $S^{\prime}_{nj}=\xi_{j}S_{nj}\ (n=0,...,N)$ as the solution of (21). ³³3 We might even affect different weights $\xi_{ij}$ to the radiations of frequency $\omega_{j}$ emitted by the different stars $(i)$, in order to account for different luminosities. However, given that our aim is to determine the spectra $S_{j}$, each of which characterizes according to Eq. (1) the axisymmetric radiation of frequency $\omega_{j}$ in the galaxy, we feel that this would not likely change the results very significantly. Therefore, to adjust the model, we determine the numbers $\xi_{j}\ (j=1,...,N_{\omega})$ so that the values $u_{j}({\bf x}_{\mathrm{loc}})$ of the SED for our local position ${\bf x}_{\mathrm{loc}}$ in the Galaxy and for the frequencies $\omega_{j}$, as given by Eq. (23), coincide with the measured values, as determined from spatial missions. We take the measured local values $f_{\bf x_{\mathrm{loc}}}(\lambda_{j})$ as plotted in Ref. [12] (see Appendix A), and we take $\rho_{\mathrm{loc}}=8$ kpc and $z_{\mathrm{loc}}=0.02$ kpc, see e.g. Ref. [19]. The model thus adjusted then allows us to make predictions: in particular, predictions of the spatial variation of the SED in the Galaxy. Such predictions may then be compared with predictions of the mainstream models of the ISRF, which models are very different from the present model.

In the foregoing work [8], the same adjustment just described was used in the framework of the “grouped fitting” (i.e. the least-squares problem (8)). A surprising result was that found for the values of the maximum of the energy density $u_{j}({\bf x})$ in the Galaxy — thus, owing to the axial symmetry (25), for the values of

found for the different spatial grids investigated, all having $\rho$ varying regularly from $\rho_{0}=0$ to $\rho_{\mathrm{max}}\simeq 10$ kpc and $z$ varying regularly from $z_{0}=0$ or $z_{0}=-z_{\mathrm{max}}$ to $z_{\mathrm{max}}\leq 1\,$kpc. ⁴⁴4  Precisely: $\rho_{0}:=\rho_{m=1}$ ; $\rho_{\mathrm{max}}:=\rho_{m=N_{\rho}}$ with, in this subsection, $\rho_{\mathrm{max}}=10\,\mathrm{kpc}\times\frac{N_{\rho}-1}{N_{\rho}}$ ; $z_{0}:=z_{p=1}$ and $z_{\mathrm{max}}:=z_{p=N_{z}}$ with, in this paper, $z_{\mathrm{max}}=1\,$kpc and $z_{0}=-z_{\mathrm{max}}$. These maximum values, which are always found at $\rho=0$, thus on the axis of symmetry, are extremely high for lower wavelengths $\lambda_{j}$, with $u_{j\mathrm{max}}\simeq 10^{27}\,\mathrm{eV/cm}^{3}$. Moreover, the value of $u_{j}(\rho=0,z)$ depends little on $z$ in the domain investigated.

These surprisingly high values occur in a larger or smaller range of wavelengths, depending on the settings of the calculation. Therefore, the question arises whether these extremely high values are a physical effect or a numerical artefact. However, the dependence on the settings is governed by the “amount of overfitting”: less overfitting increases the range of the high values [8]. This makes it plausible that the high values might be a true prediction of the model. We will now try to check whether this is indeed the case.

This model’s prediction for the spatial variation of the SED in the Galaxy was investigated, using again the separate fitting and the adjustment of the local SED on the measured values (both being described in Sect. 3). It was shown by using two different types of representations.

First, we plotted the SED at four different points in the Galaxy, and we compared the results with those obtained by Popescu et al. [5], who used a radiation transfer model built by them. (Their model also assumes axisymmetry.) Figures 8–11 show this comparison, our model being used here with $N=192$ and $N_{\omega}=76$. (Other choices of parameters that we tried gave similar figures.) It can be seen that the predictions of the present model do not differ very significantly from those of the radiation transfer model of Ref. [5]. The main difference is that our calculations oscillate somewhat strongly with the wavelength. The comparison of Figs. 8–11 here with Figs. 2-5 in Ref. [8] shows that the difference between the results of the two models is significantly smaller now than it was in our previous work, in which the calculations were based on the grouped fitting [8]: the difference in $\mathrm{log}_{10}(u_{j})$ between the results of our model and Ref. [5] is here $\lesssim 1$, whereas it went beyond $3$ and even $4$ in the previous calculations. However, the new calculations oscillate with the wavelength also at higher wavelengths. Whereas, when the grouped fitting was used, there was virtually no oscillation for $\lambda\gtrsim 10\,\mu\mathrm{m}$ at the two positions at $\rho=8\,\mathrm{kpc}$. (There were oscillations in the whole range of $\lambda$ for the two positions at $\rho=1\,\mathrm{kpc}$.) In order to check if those calculations inside the spatial domain of small values of the SED could be “polluted” by the extremely high values on the galaxy’s axis, we investigated the effect of doing the fitting on a “shifted” grid with $\rho\geq 1\,\mathrm{kpc}$. This did not lead to less oscillations. The general reason for these oscillations may be simply that this model takes fully into account the wave nature of the EM field.

Second, we plotted the radial and vertical profiles of the radiation fields at three wavelengths close to the ones considered in Fig. 7 of Popescu et al. [5] (“K, B , UV”). Figures 12 and 13 show these profiles as they are calculated at points $(\rho,z)$ belonging to the “logarithmic” grid on which the fitting was done for this calculation (see the legend). Those profiles of the radiation fields obtained with the present model on the fitting grid are relatively similar to those that are predicted with the very different model of Ref. [5], both in the levels and in the rough shape of the profiles. The most important difference is seen on the vertical profiles of Fig. 13: according to the Maxwell model of the ISRF, the energy density level decreases more slowly when the altitude $z$ increases — or even, for the $\lambda=2.29\mu\mathrm{m}$ radiation at $\rho=0.059\,$kpc or the $\lambda=0.113\mu\mathrm{m}$ radiation at $\rho=7.5\,$kpc, the level of the SED does not decrease in the range considered for $z$. A similar lack of decrease is found on the radial profiles of Fig. 12, for the $\lambda=2.29\mu\mathrm{m}$ radiation, either at $z\simeq 0$ or at $z=1.25\,\mathrm{kpc}$. Using that same fitting done on a logarithmic grid, we also calculated and plotted the radial and vertical profiles of the same radiations, but this time for regularly spaced values of $\rho$ (or respectively $z$), and in a larger range for $\rho$ (or respectively $z$), Figs. 14 and 15. The radial profiles of Figs. 12 and 14 are consistent, although, in contrast with Fig. 12, Fig. 14 plots the SED at points $(\rho,z)$ which were not involved in the fitting, and which moreover involve an extrapolation to a larger range for $\rho$ as compared with the fitting. The vertical profiles of Fig. 15, which also correspond with points which were not involved in the fitting, and also involve an extrapolation to a larger range for $z$ as compared with the fitting, show an oscillating behaviour without any tendency to a decrease at large $z$.

To help understanding the behaviours just noted, in this section we study the asymptotic behaviour of the expressions of the components of the EM field and of the SED, as they are given by the Maxwell model of the radiation field. The expressions (14)–(16) that are implemented in the numerical model, are deduced from the exact integral expressions of the EM field for a given angular frequency $\omega$, after the summation over the frequencies, and after the discretization (7) is done. Hence, we begin with the exact integral expressions of the EM field for a given angular frequency. These expressions, which are valid for any totally propagating, axisymmetric, time-harmonic EM field, are (Eqs. (13)–(15) in Ref. [6]):

where $K:=\omega/c$ — the other components being obtained by the duality (18) from the components (41)–(43), with in the most general case an other spectrum function $S^{\prime}(k)$.

The dependence in $\rho$ of the components (41)–(43) is determined by that of the Bessel functions $J_{0}$ and $J_{1}$, and by the form of the integrals which involve them. At large $x$ we have the asymptotic expansion [20]

However, the argument of the Bessel functions in Eqs. (41)–(43) is $x=\rho\sqrt{K^{2}-k^{2}}$. Hence, as $\rho\rightarrow\infty$, $x$ does not tend towards $\infty$ uniformly, depending on the integration variable $k$: we even have $x\equiv 0$ independently of $\rho$, for $k=\pm K$. Therefore, it is not obvious to see if the integrals (41)–(43) do have an expansion at fixed $z$ as $\rho\rightarrow\infty$.

As to the behaviour at fixed $\rho$ and at large $z$: up to the real part, and for a fixed value of $\rho$, the components (41)–(43) are expressions of the form $e^{-\mathrm{i}\omega t}I(z)$, with

and where, specifically, $\,a=-K,\,b=+K$, and the phase function is simply $\ g(k)\equiv k$, which has no stationary point. (The regular function $k\mapsto f(k)$ depends on the component being considered, and also on $\rho$ as a parameter.) In that case, we have the asymptotic expansion [21]

So at large $z$ and for a fixed value of $\rho$, all components of any totally propagating, axisymmetric, time-harmonic EM field are order $\frac{1}{z}$ (unless the coefficient of $\frac{1}{z}$ in this expansion is zero, which is true only in particular cases — then the relevant component is higher order in $\frac{1}{z}$). This applies indeed to the part (i) of the decomposition (i)-(ii), that is given by Eqs. (41)–(43), but also to the part (ii), since it is obtained from (41)–(43) by applying the EM duality (18) (with, in the most general case, a different spectrum function $S^{\prime}(k)$). Hence the SED (23) is order $\frac{1}{z^{2}}$ at large $z$, for any fixed value of $\rho$ — when the $C^{(q)}_{j}({\bf x})$ coefficients correspond with the exact expressions (41)–(43). [The explicit expression of the coefficient of $\frac{1}{z^{2}}$, depending on $\rho$ , $K$, $S(K)$, $S(-K)$ (and, in the most general case, of the values $S^{\prime}(K)$, $S^{\prime}(-K)$ of the spectrum function $S^{\prime}$ corresponding to the part (ii) of the decomposition (i)-(ii)) might easily be obtained from (23), (41)–(43), and (46).] The foregoing result applies to a general spectrum function $S(k)$ (and $S^{\prime}(k)$). By summation on the frequency index $j$, it extends to an EM field having a finite set of frequencies.

Let us now investigate the asymptotic behaviour of the EM field and the SED, still in the totally propagating case with axial symmetry, but now after the summation over the frequencies and the discretization (7). After the discretization, each among the $C^{(q)}_{j}$ coefficients in the expansions (24) of the components of the EM field has the form [8]:

where $R^{\prime}\,_{n}{}^{(q)}>0\$ (except for $R^{\prime}\,_{0}{}^{(q)}$ and $R^{\prime}\,_{N}{}^{(q)}$, both of which turn out to be zero) are constant numbers, and where $G_{nj}(z)=\exp\left(\mathrm{i}\omega_{j}\,t\right)F_{nj}(t,z)$ is just the function $F_{nj}$ in Eq. (17) hereabove, deprived of its periodic time-dependence (and thus is a periodic function of $z$). Together with (44), Eq. (47) shows that, at a given value of $z$, we have $C^{(q)}_{j}=O(1/\sqrt{\rho})$ as $\rho\rightarrow\infty$. The SED for an EM field having a finite set of frequencies is given by Eq. (23). For any given frequency $(j)$, $u_{j}$ is a quadratic form of the $C^{(q)}_{j}$ coefficients, hence