Transport in Astrophysics: II. Diffusion with Advection in Expanding Nebulae

Lorenzo Zaninetti

Physics Department, Turin, Italy.

DOI: 10.4236/ijaa.2022.123012                

Abstract

The structure across an expanding shell in which drift and diffusion redistribute material requires careful consideration in order to correlate the surface brightness with the physical parameters. A new solution in terms of Fourier series is suggested for 1D diffusion in the presence of a drift velocity. The astrophysical parameters are chosen in agreement with the astronomical data for the planetary nebula A39. The new solution is then inserted into the existing theory for the astrophysical image which allows dealing with the intensity of radiation emitted in a spherical layer.

Keywords

Particle Diffusion, Random Walks

1. Introduction

In order to model the images of planetary nebulae (PN) and supernova remnants (SNR) we should parametrize the thickness of the expanding layer as well the involved physics. As an example, the number density of particles in the expanding emitting layer can be assumed constant or variable due to some diffusive process. A previous analysis of diffusion with drift has covered the stationary case [1] leaving the transient state as the subject of further research. A 1D diffusion with drift is governed by Fick’s second equation, which is a partial differential Equation (PDE) that, at the moment of writing, is a subject of research. We report two examples of the ongoing research on the above equation: an analytical and numerical way has been suggested by [2] and four different methods have been implemented by [3]. The intensity of the emitted radiation versus distance from the centre in astrophysical phenomena such PN and SNR presents a characteristic “U” shape which can be explained as the solution of the transport equation for the light in the case of an optically thin medium, starting with [4], or with a shell with constant emissivity [5].

This paper reviews the existing solutions for 1D diffusion with drift in Section 2, derives a new series solution for 1D diffusion with drift in Section 3, specifies in Section 4 the parameters to be adopted in an astrophysical environment such as that given by the planetary nebula (PN) A39 and applies in Section 5 the new results to the formation of an astrophysical image.

2. The Existing Solutions of Diffusion with Drift

We review 1D diffusion with drift in the framework of the mathematical theory of diffusion and asymmetric random walk with drift as developed in [1]. In order to avoid duplicates we adopt the following astrophysical units: pc for length and yr for time. In these units, the advection velocity v is expressed in pc/yr and the diffusion coefficient in pc²/yr.

2.1. 1D Diffusion with Drift, Stationary State

In one dimension and in the presence of a drift velocity, v, along the x-direction, a diffusion is governed by Fick’s second equation for the concentration, $C\left(x\right)$, see Equation (4.5) in [6],

$\frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x^2}-v\frac{\partial C}{\partial x}\mathrm{,}$ (1)

where v can take positive or negative values and D is the diffusion coefficient. The number density rises from 0 at x = a to a maximum value C_m at x = b and then falls again to 0 at r = c. The general solution to Equation (1) for a steady state is

$C\left(x\right)=A+B{\text{e}}^{\frac{vx}{D}}\mathrm{.}$ (2)

We now assume that v is negative: the solution is

$C\left(x\right)=A+B{\text{e}}^{-\frac{v}{D}x}\mathrm{.}$ (3)

The boundary conditions give

$C_{a\mathrm{,}b}\left(x\right)=\frac{C_m\left({\text{e}}^{-\frac{vx}{D}}-{\text{e}}^{-\frac{va}{D}}\right)}{-{\text{e}}^{-\frac{va}{D}}+{\text{e}}^{-\frac{vb}{D}}}a\le x\le b\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\text{downstreamside}\mathrm{,}$ (4)

and

$C_{b\mathrm{,}c}\left(x\right)=\frac{C_m\left({\text{e}}^{-\frac{vc}{D}}-{\text{e}}^{-\frac{vx}{D}}\right)}{-{\text{e}}^{-\frac{vb}{D}}+{\text{e}}^{-\frac{vc}{D}}}\mathrm{.}$ (5)

2.2. 1D Diffusion with Drift, Random Walk

Given a 1D segment of length side we can implement a random walk with

step-length $\lambda$ by introducing the numerical parameter $\text{NDIM}=\frac{\text{side}}{\lambda }$. We now report the adopted rules when the injection is in the middle of the grid:

1) The first of NPART particles is chosen.

2) The random walk of the particle starts in the middle of the grid. The probability of taking one step in the negative direction (downstream), is $p_1=\frac{1}{2}-\mu \times \frac{1}{2}$, and in the positive direction (upstream), $p_2=\frac{1}{2}+\mu \times \frac{1}{2}$, where $\mu$ is a parameter that characterizes the asymmetry ( $0\le \mu \le 1$ ).

3) When the particle reaches one of the two absorbing points, the motion starts again from (2) with a different diffusion pattern.

4) The number of visits is recorded on $\mathcal{M}$, a 1D grid.

5) The random walk terminates when all the NPART particles have been processed.

6) For the sake of normalization, the one-dimensional visitation or number density grid $\mathcal{M}$ is divided by NPART.

There is a systematic change of the average particle position along the x-direction:

$\langle dx\rangle =\mu \lambda \mathrm{,}$ (6)

for each time step. If the time step is $dt=\frac{\lambda }{v_{tr}}$ where $v_{tr}$ is the transport velocity, the asymmetry, $\mu$, that characterizes the random walk is

$\mu =\frac{v}{v_{tr}}.$ (7)

3. 1D Diffusion with Drift, Transient State

In 1D and in the presence of a drift velocity, v, along the radial direction the diffusion is governed by Fick’s second equation, see Equation (4.5) in [6],

$\frac{\partial }{\partial t}u\left(x,t\right)=D\left(\frac{{\partial }^2}{\partial x^2}u\left(x,t\right)\right)-v\left(\frac{\partial }{\partial x}u\left(x,t\right)\right).$ (8)

On assuming that the boundary conditions are $u\left(\mathrm{0,}t\right)=0$, $u\left(L\mathrm{,}t\right)=0$ and the profile of density at $t=0$ is $u\left(x\mathrm{,0}\right)=f\left(x\right)$ we have the following solution in terms of a Fourier series

$u\left(x,t\right)=N_0\frac{1}{L}\left(-2\left(\sinh \left(A\right)-\cosh \left(A\right)\right)\left(\cosh \left(B\right)+\sinh \left(B\right)\right)I_{fx}\left({\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\sin \left(\frac{\pi xn}{L}\right)}\right)\right),$ (9)

where

$A=\frac{\left(4D^2{\pi }^2{n}^2+L^2{v}^2\right)t}{4D{L}^2},$ (10)

$B=\frac{xv}{2D},$ (11)

and

$I_{fx}={{\displaystyle \int }}_{\text{}\text{}\text{}0}^{L}f\left(x\right)\left(\cosh \left(\frac{xv}{2D}\right)+\sinh \left(\frac{xv}{2D}\right)\right)\sin \left(\frac{\pi xn}{L}\right){\text{e}}^{-\frac{xv}{D}}\text{d}x,$ (12)

where N₀ is the number of particles at $x=0,t=0$. We now analyse three types of initial conditions. The first case to be analysed is a Gaussian density profile of the type

$u\left(x\mathrm{,0}\right)={\text{e}}^{-\frac{{\left(x-\frac{L}{2}\right)}^2}{b^2}}\mathrm{,}$ (13)

where b is an adjustable parameter. The integral (12) is

$\begin{array}{c}{I}_{\text{Gauss}}=-\frac{I}{4}b((\text{erf}\left(\frac{\left(2ID\pi n-Lv\right)b^2+2D{L}^2}{4bDL}\right)\\ +\text{erf}\left(\frac{\left(-2ID\pi n+Lv\right)b^2+2D{L}^2}{4bDL}\right)){\text{e}}^{\frac{\left(I\pi n{b}^{2}D-\frac{L{b}^{2}v}{2}+2D{L}^2\right)\left(ID\pi n-\frac{Lv}{2}\right)}{4D^2{L}^2}}\\ -{\text{e}}^{\frac{\left(ID\pi n+\frac{Lv}{2}\right)\left(I\pi n{b}^{2}D+\frac{L{b}^{2}v}{2}-2D{L}^2\right)}{4D^2{L}^2}}(\text{erf}\left(\frac{\left(2ID\pi n+Lv\right)b^2+2D{L}^2}{4bDL}\right)\\ +\text{erf}\left(\frac{\left(-2ID\pi n-Lv\right)b^2+2D{L}^2}{4bDL}\right)))\sqrt{\pi },\end{array}$ (14)

where $I=\sqrt{-1}$. In the case of a Gaussian profile of the density Figure 1 presents the solution for different values of the velocity.

Figure 1. The solution to Equation (9) in the case of a Gaussian profile of the density when $t=100\text{yr}$ (red line), $t=500\text{yr}$ (blue line) and $t=1000\text{yr}$ (green line). The other parameters are $n=100$, $v=3.521\times {10}^{-5}\frac{\text{pc}}{\text{yr}}$, $D=6\times {10}^{-7}\frac{{\text{pc}}^2}{\text{yr}}$, $L=0.228428\text{pc}$, $b=4\text{pc}$ and $N_0=1$.

The second case to be analysed is a parabolic density profile of the type

$u\left(x\mathrm{,0}\right)=\frac{4x\left(L-x\right)}{L^2}\mathrm{.}$ (15)

The integral (12) is

$\begin{array}{l}{I}_{\text{parabola}}=\frac{1}{{\left(4D^2{\pi }^2{n}^2+L^2{v}^2\right)}^3}\\ \times (-256L((2n\pi \left(-\frac{1}{2}{D}^{2}L{\pi }^2{n}^{2}v+D^3{\pi }^2{n}^2-\frac{1}{8}{L}^3{v}^3-\frac{3}{4}D{L}^2{v}^2\right)D\cos \left(\pi n\right)\\ +\sin \left(\pi n\right)\left(D^4{\pi }^4{n}^4+3D^{3}L{\pi }^2{n}^{2}v-\frac{1}{16}{L}^4{v}^4-\frac{1}{4}D{L}^3{v}^3\right))\cosh \left(\frac{Lv}{2D}\right)\\ +(-2n\pi \left(-\frac{1}{2}{D}^{2}L{\pi }^2{n}^{2}v+D^3{\pi }^2{n}^2-\frac{1}{8}{L}^3{v}^3-\frac{3}{4}D{L}^2{v}^2\right)D\cos \left(\pi n\right)\\ -\sin \left(\pi n\right)\left(D^4{\pi }^4{n}^4+3D^{3}L{\pi }^2{n}^{2}v-\frac{1}{16}{L}^4{v}^4-\frac{1}{4}D{L}^3{v}^3\right))\sinh \left(\frac{Lv}{2D}\right)\\ -2n\pi D\left(\frac{1}{2}{D}^{2}L{\pi }^2{n}^{2}v+D^3{\pi }^2{n}^2+\frac{1}{8}{L}^3{v}^3-\frac{3}{4}D{L}^2{v}^2\right))D^2).\end{array}$ (16)

In the case of a parabolic profile of the density Figure 2 presents the solution for different values of the drift velocity.

The third case to be analysed is a trigonometric density profile of the type

Figure 2. The solution to Equation (9) for the case of a parabolic profile of the density when $v={10}^{-6}\frac{\text{pc}}{\text{yr}}$ (red line), $v=2\times {10}^{-5}\frac{\text{pc}}{\text{yr}}$ (blue line) and $v=3.521\times {10}^{-5}\frac{\text{pc}}{\text{yr}}$ (green line). The other parameters are $n=100$, $t=1000\text{yr}$, $D=6\times {10}^{-7}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$, $L=0.228428\text{pc}$ and $N_0=1$.

$u\left(x\mathrm{,0}\right)=\cos \left(\frac{\left(\frac{L}{2}-x\right)\pi }{2L}\right)\mathrm{.}$ (17)

The integral (12) is

$\begin{array}{l}{I}_{\text{trigonometric}}\\ =\frac{1}{\left(4D^2{\pi }^2{n}^2-4D^2{\pi }^{2}n+D^2{\pi }^2+L^2{v}^2\right)\left(4D^2{\pi }^2{n}^2+4D^2{\pi }^{2}n+D^2{\pi }^2+L^2{v}^2\right)}\\ \times 8((\left({\pi }^3\left(-\frac{n^2}{2}+\frac{1}{8}\right)D^3-\frac{\left(n^2+\frac{1}{4}\right){\pi }^{2}vL{D}^2}{2}+\frac{L^{2}D\pi v^2}{8}-\frac{L^3{v}^3}{8}\right)\sin \left(\pi n\right)\\ \begin{array}{c}\stackrel{}{}\\ \\ \end{array}-\cos \left(\pi n\right)\pi \left({\pi }^2\left(n^2-\frac{1}{4}\right)D^2-\frac{DL\pi v}{2}+\frac{L^2{v}^2}{4}\right)nD)\cosh \left(\frac{Lv}{2D}\right)\end{array}$

$\begin{array}{l}+(\left({\pi }^3\left(\frac{n^2}{2}-\frac{1}{8}\right)D^3+\frac{\left(n^2+\frac{1}{4}\right){\pi }^{2}vL{D}^2}{2}-\frac{L^{2}D\pi v^2}{8}+\frac{L^3{v}^3}{8}\right)\sin \left(\pi n\right)\\ \begin{array}{c}\stackrel{}{}\\ \\ \end{array}+\cos \left(\pi n\right)\pi \left({\pi }^2\left(n^2-\frac{1}{4}\right)D^2-\frac{DL\pi v}{2}+\frac{L^2{v}^2}{4}\right)nD)\sinh \left(\frac{Lv}{2D}\right)\\ \begin{array}{c}\stackrel{}{}\\ \\ \end{array}+\left({\pi }^2\left(n^2-\frac{1}{4}\right)D^2+\frac{DL\pi v}{2}+\frac{L^2{v}^2}{4}\right)\pi nD)L\sqrt{2}D\mathrm{.}\end{array}$ (18)

The trigonometric solution as a function of x and t is presented as a surface plot in Figure 3.

4. Astrophysical Applications

We present the astronomical data on PN A39.

4.1. A39

The PN A39 is extremely round and therefore can be considered an example of spherical symmetry, see for example Figure 1 in [7] or Figure 4.

The radius of the shell in A39, $R_{shell}$, is

$R_{shell}=2.42\times {10}^{18}{\Theta }_{77}{D}_{21}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\text{cm}=\mathrm{0.78\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{pc}\mathrm{,}$ (19)

where ${\Theta }_{77}$ is the angular radius in units of 77" and $D_{21}$ the distance in units of 2.1 kpc, see [7].

Figure 3. The solution to Equation (9) as function of time and space for the case of a trigonometric profile of the density when $v=3.521\times {10}^{-5}\frac{\text{pc}}{\text{yr}}$. The other parameters are $n=100$, $t=1000\text{yr}$, $D=6\times {10}^{-7}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$, $L=0.228428\text{pc}$ and $N_0=1$.

Figure 4. Image of A39 through a blue-green filter that isolates the light emitted by oxygen atoms in the nebula at a wavelength of 500.7, the size of the image is ≈298" and the credit is given to WIYN/NOIRLab/NSF.

The expansion velocity has a range $\left[32\leftrightarrow \mathrm{37\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{\text{km}}{\text{s}}\right]$ according to [8] and the age of the free expansion is 23,000 yr, see [7]. In our astrophysical units we will use a $v=3.5215\times {10}^{-5}\frac{\text{pc}}{\text{yr}}$ which corresponds to the average value of the above range in velocity. The angular thickness of the shell is

$\delta r_{shell}=3.17\times {10}^{17}{\Theta }_{10}{D}_{21}\text{cm}=\mathrm{0.103\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{pc}\mathrm{,}$ (20)

where ${\Theta }_{10}$ is the thickness in units of 10.1" and the height above the galactic plane is 1.42 kpc, see [7]. For the sake of comparison with the observations, the length of some geometrical models will be expressed in arcsec. An occasional reader may ask whether A39 can be approximated by a sphere or not? One way to answer this question is to introduce the observational percentage of sphericity, ${\epsilon }_{\text{obs}}$, over the whole range of the polar angle $\theta$,

${\epsilon }_{\text{obs}}=100\left(1-\frac{{\displaystyle {\sum }_j{\left|r_{\text{obs}}-r_{\text{num}}\right|}_j}}{{\displaystyle {\sum }_j{r}_{\text{obs}}{}_{\mathrm{,}j}}}\right)\mathrm{,}$ (21)

where $r_{\text{num}}$ is the theoretical averaged radius of the nebula, $r_{\text{obs}}$ is the observed radius of the nebula, and the index j varies from 1 to the number of available observations. In the case of Figure 4, we have ${\epsilon }_{\text{obs}}=\mathrm{97\; <mi>\; \%\; </mi>}$ over 67 directions when the central star is assumed to be the source of the expansion.

4.2. Astrophysical Diffusion

Figure 5 presents $\mathcal{M}\left(x\right)$, the number of visits generated by the Monte Carlo simulation as well as the steady state mathematical solution with drift given by Formulas (4) and (5).

In the case of a transient diffusion, we present, in Figure 6, the parabolic solution that will be used in the theory of formation of the image.

5. Image

We review the radiative transfer equation and the formation of the image when a spherical object is emitting radiation in a thin layer with constant density and in the presence of a stationary diffusion with drift. We introduce two new models for the emission of a spherical PN: a purely geometrical model and a model based on a transient diffusion.

Figure 5. Number density in A39 of the 1D asymmetric random walk (full line), NDIM = 401, NPART = 200, side = 40 arcsec, $\lambda =0.1$ arcsec and $\mu =-0.013$. For astrophysical purposes $\mu$ is negative. The theoretical number density as given by Formulas (4) and (5) is presented when $u=1$, $C_m=1$, $a=60$ arcsec, $b=80$ arcsec, $c=100$ arcsec and $D=3.84$ (dotted line).

Figure 6. The solution to Equation (9) in the case of a parabolic profile of the density when $v=3.521\times {10}^{-5}\frac{\text{pc}}{\text{yr}}$ in shifted coordinates with zero at the centre of A39. The other parameters are $n=100$, $t=1000\text{yr}$, $D=6\times {10}^{-7}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$, $L=0.228428\text{pc}$ and $N_0=1$.

5.1. Transfer Equation

The transfer equation in the presence of emission only, see for example [9] or [10], is

$\frac{\text{d}{I}_{\nu }}{\text{d}s}=-k_{\nu }\zeta I_{\nu }+j_{\nu }\zeta \mathrm{,}$ (22)

where $I_{\nu }$ is the specific intensity, s is the line of sight, $j_{\nu }$ the emission coefficient, $k_{\nu }$ a mass absorption coefficient, $\zeta$ the mass density at position s and the index $\nu$ denotes the frequency of emission under consideration. The solution to Equation (22) is

$I_{\nu }\left({\tau }_{\nu }\right)=\frac{j_{\nu }}{k_{\nu }}\left(1-{\text{e}}^{-{\tau }_{\nu }\left(s\right)}\right)\mathrm{,}$ (23)

where ${\tau }_{\nu }$ is the optical depth at frequency $\nu$

$\text{d}{\tau }_{\nu }=k_{\nu }\zeta \text{d}s\mathrm{.}$ (24)

We now continue analysing the case of an optically thin layer in which ${\tau }_{\nu }$ is very small (or $k_{\nu }$ very small) and the density $\zeta$ is substituted with $C\left(s\right)$, our number density of particles. We now assume that the emissivity is proportional to the number density:

$j_{\nu }\zeta =KC\left(s\right)\mathrm{,}$ (25)

where K is a constant function. This can be the case for synchrotron radiation in the presence of an isotropic distribution of electrons with a power law distribution in energy, $N\left(E\right)$,

$N\left(E\right)dE=K_s{E}^{-{\gamma }_f}\mathrm{,}$ (26)

where $K_s$ is a constant. In this case the emissivity is

$j_{\nu }\rho \approx 0.933\times {10}^{-23}\alpha \left({\gamma }_f\right)K_s{H}_{\perp }^{\frac{{\gamma }_f+1}{2}}{\left(\frac{6.26\times {10}^{18}}{\nu }\right)}^{\frac{{\gamma }_f-1}{2}}\frac{\text{erg}}{\text{s}\cdot {\text{cm}}^3\cdot \text{Hz}\cdot {\text{rad}}^2}\mathrm{,}$ (27)

where $\nu$ is the frequency and $\alpha \left({\gamma }_f\right)$ is a slowly varying function of ${\gamma }_f$ which is of the order of unity and is given by

$\alpha \left({\gamma }_f\right)=2^{\left({\gamma }_f-3\right)/2}\frac{{\gamma }_f+7/3}{{\gamma }_f+1}\Gamma \left(\frac{3{\gamma }_f-1}{12}\right)\Gamma \left(\frac{3{\gamma }_f+7}{12}\right)\mathrm{,}$ (28)

for ${\gamma }_f\ge \frac{1}{2}$, see Formula (1.175) in [11]. Synchrotron emission is widely used to explain the radiation observed in SNR, see [12] - [17]. This non-thermal radiation continuum emission was also detected in a PN associated with a very long-period OH/IR variable star (V1018 Sco) [18], in the young PN IRAS 15103-5754 [19] and in IRAS 18062 + 2410 [20]. A discussion of thermal/non- thermal emission from PN’s can be found in [21] [22]. The intensity in the optically thin layer is therefore

$I_{\nu }\left(s\right)=K{\displaystyle {\int }_{s_0}^s}C\left(s\text{'}\right)\text{d}s\text{'}\text{optically}\text{thin}\text{layer}\text{linearcase}\mathrm{.}$ (29)

In Monte Carlo experiments, the number density is stored on the grid $\mathcal{M}$ and the intensity is

$I\left(i\mathrm{,}j\right)={\displaystyle \underset{k}{\sum }}\Delta s\times \mathcal{M}\left(i\mathrm{,}j\mathrm{,}k\right)\text{optically}\text{thin}\text{layer}\text{linearcase}\mathrm{.}$ (30)

where Δs is the spatial interval between the various values and the sum is performed over the interval of existence of the index k. The theoretical intensity is then obtained by integrating the intensity at a given frequency over the solid angle of the source.

5.2. Image Theory at Constant Density

The simplest model for the image is characterized by a radiative process with constant number density of emitting particles in a thin layer around the advancing sphere. We therefore assume that the number density C is constant in the spherical thin layer, and in particular rises from 0 at r = a to a maximum value C_m, remains constant up to r = b and then falls again to 0. The line of sight, when the observer is situated at the infinity of the x-axis, is the locus parallel to the x-axis which crosses the position y in the Cartesian x-y plane and terminates at the external circle of radius b, see Figure 7. The length of this locus is

$l_{0a}=2\times \left(\sqrt{b^2-y^2}-\sqrt{a^2-y^2}\right);0\le y<a$

$l_{ab}=2\times \sqrt{b^2-y^2};a\le y<b.$ (31)

The number density C_m is constant between two spheres of radii a and b and therefore the intensity of radiation is

$I_{0a}=C_m\times 2\times \left(\sqrt{b^2-y^2}-\sqrt{a^2-y^2}\right);0\le y<a$ (32)

Figure 7. The two circles (section of spheres) which include the region with constant density are represented by a full line. The observer is situated along the x direction, and three lines of sight, s₁, s₂ and s₃, are presented.

$I_{ab}=C_m\times 2\times \sqrt{b^2-y^2};a\le y<b.$ (33)

The ratio between the theoretical intensity at the maximum ( $y=a$ ) and at the minimum ( $y=0$ ) is given by

$\frac{I_{0a}=\left(y=a\right)}{I_{0a}=\left(y=0\right)}=\frac{\sqrt{b^2-a^2}}{b-a}.$ (34)

The quality of the fit is measured by ${\chi }^2$ :

${\chi }^2={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{{\left(T_i-O_i\right)}^2}{{\sigma }_i{}^2},$ (35)

where n is the number of points, $T_i$ is the theoretical value, $O_i$ is the observed value and ${\sigma }_i$ is the error with respect to the observed value.

A reduced merit function ${\chi }_{red}^2$ is evaluated by

${\chi }_{red}^2={\chi }^2/NF,$ (36)

where $NF\mathrm{<mo>\; =\; </mo>}n-k$ is the number of degrees of freedom, n is the number of elements and k is the number of parameters. A comparison of the observed data for A39 and the theoretical intensity of the simplest model are presented in Figure 8.

5.3. Image Theory in Stationary State with Drift

Figure 9 shows a spherical shell source of radius b between a spherical absorber of radius a and a spherical absorber of radius c. The number density rises from 0 at r = a to a maximum value C_m at r = b and then falls again to 0 at r = c. In order to continue, the 3D spherical diffusion is approximated by the 1D steady state diffusion with drift and the numbers density to be used are Formulas (4) and (5) once $r=\sqrt{x^2+y^2}$ is imposed. In this case, the operation of integration of the number density as given by Equation (29) can be performed only numerically, see Figure 10.

Figure 8. Cut of the mathematical intensity I (Formulas (32) and (33)), crossing the centre (full line) of A39 and real data (dotted line with some error bars). The parameters are $a=72.5^{″}$, $a=90.18^{″}$ and $C_m=1$. The number of data is 801 and for this model ${\chi }_{red}^2=10.84$ against ${\chi }_{red}^2=6.26$ for the rim model fully described in Jacoby et al. (2001).

Figure 9. The circular source inserted in the great box is represented by a dashed line, and the two absorbing boundaries by a full line. The observer is situated along the x direction, and three lines of sight are indicated.

Figure 10. Cut of the numerical intensity I in the case of steady state with drift which crosses the centre (full line) of A39 with parameters as in Figure 5 and real data (dotted line with some error bars). The number of data is 801 and for this model ${\chi }_{red}^2=20.96$ against ${\chi }_{red}^2=10.36$ of the rim model fully described in Jacoby et al. (2001).

5.4. Image Theory in a Geometrical Model

We assume that the number density rises from 0 at r = a to a maximum value C_m at r = b and then falls again to 0 at r = c following a parabolic behavior with vertex at r = b. The parabola that satisfies the above requirements is

$C\left(r;a,b,c\right)=-\frac{C_m\left(c-r\right)\left(a-r\right)}{\left(b-c\right)\left(a-b\right)},$ (37)

and we continue inserting $r=\sqrt{x^2+y^2}$. The integration which leads to the image, see Equation (29) and Figure 9, can be done and is

$\begin{array}{l}{I}_{0a}=C_m(3y^2\text{arcsinh}\left(\frac{\sqrt{c^2-y^2}}{y}\right)a+3y^2\text{arcsinh}\left(\frac{\sqrt{c^2-y^2}}{y}\right)c\\ -3y^2\text{arcsinh}\left(\frac{\sqrt{a^2-y^2}}{y}\right)a-3y^2\text{arcsinh}\left(\frac{\sqrt{a^2-y^2}}{y}\right)c\\ -3\sqrt{c^2-y^2}ca+\sqrt{c^2-y^2}{c}^2-4\sqrt{c^2-y^2}{y}^2-\sqrt{a^2-y^2}{a}^2\\ \begin{array}{c}\stackrel{}{}\\ \end{array}+3\sqrt{a^2-y^2}ac+4\sqrt{a^2-y^2}{y}^2);0\le y<a\end{array}$ (38)

$\begin{array}{l}{I}_{ac}=\frac{1}{3\left(b-c\right)\left(a-b\right)}\times C_m(3y^2\text{arcsinh}\left(\frac{\sqrt{c^2-y^2}}{y}\right)a\\ +3y^2\text{arcsinh}\left(\frac{\sqrt{c^2-y^2}}{y}\right)c-3y^2\text{arcsinh}\left(\frac{\sqrt{a^2-y^2}}{y}\right)a\\ -3y^2\text{arcsinh}\left(\frac{\sqrt{a^2-y^2}}{y}\right)c-3\sqrt{c^2-y^2}ca+\sqrt{c^2-y^2}{c}^2\\ -4\sqrt{c^2-y^2}{y}^2-\sqrt{a^2-y^2}{a}^2+3\sqrt{a^2-y^2}ac\\ \begin{array}{c}\stackrel{}{}\\ \end{array}+4\sqrt{a^2-y^2}{y}^2)\mathrm{<mo>\; ;\; </mo>}a\le y<c.\end{array}$ (39)

In the case of the geometrical model the theoretical intensity of A39 is presented as a cut, see Figure 11 and as a surface, see Figure 12.

5.5. Image Theory in a Transient State

Figure 9 can model the formation of the image in the case of a transient state developed in Section 3 once a is replaced by zero, b is replaced by L/2 and c is replaced by L. The solution considered is Equation (9) for a parabolic density profile, see Equation (15). In this case the operation of integration which leads to the image is done numerically and is presented as a cut, see Figure 13 and as an image, see Figure 14.

The effect of the insertion of a threshold intensity, $I_{tr}$, given by the observational techniques, will now be analysed. The threshold intensity can be parametrized by $I_{\max }$, the maximum value of intensity characterizing the ring: a typical image with a hole is visible in Figure 15 when $I_{tr}=I_{\max }/2$.

Figure 11. Cut of the numerical intensity I in the case of the geometrical model which crosses the centre (full line) of A39 ( $a=72^{″}$, $b=89^{″}$, $b=106^{″}$ ) and real data (dotted line with some error bars). The number of data is 801 and for this model ${\chi }_{red}^2=10.69$ against ${\chi }_{red}^2=5.65$ of the rim model fully described in Jacoby et al. (2001).

Figure 12. Surface of the theoretical intensity of A39 in the geometrical model, parameters as in Figure 11.

Figure 13. Cut of the numerical intensity I in the case of a transient state with drift and parabolic density profile which crosses the centre (full line) of A39 with parameters as in Figure 6 and real data (dotted line with some error bars). The number of data is 801 and for this model ${\chi }_{red}^2=7.92$ against ${\chi }_{red}^2=8.65$ of the rim model fully described in Jacoby et al. (2001).

Figure 14. Contour map of I particularized to simulate A39, parameters as in Figure 13.

Figure 15. The same as Figure 14 but with $I_{tr}=I_{\max }/2$.

6. Conclusions

PDE & Boundary Conditions

A new solution for the diffusion with drift has been derived as a Fourier series, see Equation (9). The above solution has been particularized for three types of density profiles at t = 0: Gaussian, parabolic and trigonometric.

Astrophysical Numbers

Once the time has been fixed in years and the velocity as 3.521 × 10⁻⁵ pc/yr we consequently choose the diffusion coefficient around 10⁻⁷ pc²/yr.

Image Theory

The image theory of a spherical emitter such as PN A39 has been modeled by an advancing layer with a constant density, a stationary advection with drift, an existing geometrical profile for the emitting layer and a time dependent diffusive model in presence of drift.

Plans for the Future

Formula (21) fixes the degree of sphericity of A39 at 97% but some asymmetries in the spatial dimensions and in the emissivity are present. An explanation of these anomalies requires introducing the concept of velocity as a function of the position angle, which was begun to be used in [1].

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Zaninetti, L. (2009) Scaling for the Intensity of Radiation in Spherical and Aspherical Planetary Nebulae. Monthly Notices of the Royal Astronomical Society, 395, 667-691. https://doi.org/10.1111/j.1365-2966.2009.14551.x
[2]	Mojtabi, A. and Deville, M.O. (2015) One-Dimensional Linear Advection— Diffusion Equation: Analytical and Finite Element Solutions. Computers & Fluids, 107, 189-195. https://doi.org/10.1016/j.compfluid.2014.11.006
[3]	Liu, D. and Wang, Y. (2015) High Order Numerical Solutions to Convection Diffusion Equations with Different Approaches Journal of Applied & Computational Mathematics, 4, 208-227.
[4]	Zaninetti, L. (2007) X-Ray Profiles in Symmetric and Asymmetric Supernova Remnants. Baltic Astronomy, 16, 251-285.
[5]	Opsenica, S. and Arbutina, B. (2011) Numerical Code for Fitting Radial Emission Profile of a Shell Supernova Remnant. Application. Serbian Astronomical Journal, 183, 75-85. https://doi.org/10.2298/SAJ1183077O
[6]	Berg, H.C. (1993) Random Walks in Biology. Princeton University Press, Princeton.
[7]	Jacoby, G.H., Ferland, G.J. and Korista, K.T. (2001) The Planetary Nebula A39: An Observational Benchmark for Numerical Modeling of Photoionized Plasmas. The Astrophysical Journal, 560, 272-286. https://doi.org/10.1086/322489
[8]	Hippelein, H. and Weinberger, R. (1990) The Expansion of Highly Evolved Planetary Nebulae. Astronomy & Astrophysics, 232, 129-134.
[9]	Rybicki, G. and Lightman, A. (1991) Radiative Processes in Astrophysics, Wiley-Interscience, New-York.
[10]	Hjellming, R.M. (1988) Radio Stars. In: Verschuur, G.L. and Kellermann, K.I., Eds., Galactic and Extragalactic Radio Astronomy, Springer-Verlag, New York, 381-438. https://doi.org/10.1007/978-1-4612-3936-9_9
[11]	Lang, K.R. (1999) Astrophysical Formulae. 3rd Edition, Springer-Verlag, New York.
[12]	Velázquez, P.F., Martinell, J.J., Raga, A.C. and Giacani, E.B. (2004) Effects of Thermal Conduction on the X-Ray and Synchrotron Emission from Supernova Remnants. The Astrophysical Journal, 601, 885-895. https://doi.org/10.1086/380755
[13]	Berezhko, E.G. and Volk, H.J. (2004) The Theory of Synchrotron Emission from Supernova Remnants. Astronomy & Astrophysics, 427, 525-536. (Preprint arXiv:astro-ph/0408121). https://doi.org/10.1051/0004-6361:20041111
[14]	Bamba, A., Yamazaki, R. and Hiraga, J.S. (2005) Chandra Observations of Galactic Supernova Remnant Vela Jr.: A New Sample of Thin Filaments Emitting Synchrotron X-Rays. The Astrophysical Journal, 632, 294-301. (Preprint arXiv:astro-ph/0506331). https://doi.org/10.1086/432711
[15]	Cassam-Chenai, G., Decourchelle, A., Ballet, J. and Ellison, D.C. (2005) Morphology of Synchrotron Emission in Young Supernova Remnants. Astronomy & Astrophysics, 443, 955-959. (Preprint arXiv:astro-ph/0507194). https://doi.org/10.1051/0004-6361:20052853
[16]	Bykov, A.M., Uvarov, Y.A. and Ellison, D.C. (2008) Dots, Clumps, and Filaments: The Intermittent Images of Synchrotron Emission in Random Magnetic Fields of Young Supernova Remnants. The Astrophysical Journal, 689, L133-L136. (Preprint 0811.2498). https://doi.org/10.1086/595868
[17]	Katsuda, S., Petre, R., Mori, K., Reynolds, S.P., Long, K.S., Winkler, P.F. and Tsunemi, H. (2010) Steady X-Ray Synchrotron Emission in the Northeastern Limb of SN 1006. The Astrophysical Journal, 723, 383-392. (Preprint 1009.0280). https://doi.org/10.1088/0004-637X/723/1/383
[18]	Cohen, M., Chapman, J.M., Deacon, R.M., Sault, R.J., Parker, Q.A. and Green, A.J. (2006) Radio Observations of the Planetary Nebula around the OH/IR Star OH354.88-0.54 (V1018 Sco). Monthly Notices of the Royal Astronomical Society, 369, 189-196. (Preprint arXiv:astro-ph/0603321). https://doi.org/10.1111/j.1365-2966.2006.10279.x
[19]	Suárez, O., Gómez, J.F., Bendjoya, P., Miranda, L.F., Guerrero, M.A., Uscanga, L., Green, J.A., Rizzo, J.R. and Ramos-Larios, G. (2015) Time-Variable Non-Thermal Emission in the Planetary Nebula IRAS 15103-5754. The Astrophysical Journal, 806, Article No. 105 (Preprint 1504.04277). https://doi.org/10.1088/0004-637X/806/1/105
[20]	Cerrigone, L., Umana, G., Trigilio, C., Leto, P., Buemi, C.S. and Ingallinera, A. (2017) Radio Variability and Non-Thermal Components in Stars Evolving towards Planetary Nebulae. Monthly Notices of the Royal Astronomical Society, 468, 3450-3460. (Preprint 1703. 06005). https://doi.org/10.1093/mnras/stx690
[21]	Hajduk, M., van Hoof, P.A.M., Sniadkowska, K., Krankowski, A., Błaszkiewicz, L., Dąbrowski, B., et al. (2018) Radio Observations of Planetary Nebulae: No Evidence for Strong Radial Density Gradients. Monthly Notices of the Royal Astronomical Society, 479, 5657-5677. (Preprint 1807.01055). https://doi.org/10.1093/mnras/sty1673
[22]	Hajduk, M., Van Hoof, P.A.M., Sniadkowska, K., Krankowski, A., Błaszkiewicz, L., Dabrowski, B., et al. (2019) Radio Continuum Spectra of Planetary Nebulae Galaxies, 7, Article No. 6. https://www.mdpi.com/2075-4434/7/1/6