Electron-Positron Pair Production in Plasma's: Implications for Magneto-sonic Waves

Abstract

In the high energy (high temperature) regime, pair production can occur due to scattering reactions, high energy photons which exchange energy for the production of an e e + e e + e^(-)e^(+)e^-e^+-pair (electron positron) or other reaction types. In this limit the presence of e e e^(-)e^- and e + e + e^(+)e^+ modifies the behaviour of the plasma. The typical energy limit for pair production (mass-energy limit) corresponds to a temperature of T 10 10 T 10 10 T∼10^(10)T\sim10^{10}~K, while numerically solving for the number densities of electron and positrons as a function of temperature in this high energy limit suggests that pair production becomes important at a much lower temperature T 3 × 10 8 T 3 × 10 8 T∼3xx10^(8)T\sim 3\times 10^{8}~K. I will focus on the effects of pair production on a plasma in the fluid limit ( L λ m f p r L L λ m f p r L L≫lambda_(mfp)≫r_(L)L\gg \lambda_{mfp} \gg r_L where L L LL is the cell size, λ m f p λ m f p lambda_(mfp)\lambda_{mfp} is the mean free path, and r L r L r_(L)r_L is the Larmor radius). As the fluid limit ignores the motion of individual particles, the effect of pair production mostly comes into play by modifying the composition of material. I use the Magneto-hydrodynamics (MHD) code FLASH, which uses a composition dependent equation of state (EOS) (helmholtz EOS), to model the effects of pair production in the fluid. I demonstrate the effect of pair production by perturbing a one dimensional fast magneto-sonic (FMS) wave in the pair production density-temperature limit in a domain with periodic boundary conditions. By changing the composition (electron fraction) by a factor 2, the pressure of the fluid varies by 5 % 5 % ∼5%\sim 5\%, the adiabatic constant γ γ gamma\gamma varies by 10 % 10 % ∼10%\sim 10\% in opposing directions resulting in a variation in the FMS wave speed of 0.5 % 0.5 % ∼0.5%\sim 0.5\%.

1. Introduction

Pair production is the creation of an electron-positron pair from energy. This reaction can take several different forms such as pair production from a gamma ray
(1) N + γ N + e + e + , (1) N + γ N + e + e + , {:(1)N+gamma rarr N+e^(-)+e^(+)",":}\begin{equation} N + \gamma \rightarrow N + e^- + e^+, \end{equation}
(where γ γ gamma\gamma is a photon, N N NN is some background nucleus necessary for momentum conservation, e e e^(-)e^- are electrons, and e + e + e^(+)e^+ are positrons) or by the collision of two photons via the Breit-Wheeler process (Breit & Wheeler 1934)
(2) γ + γ e + e + . (2) γ + γ e + e + . {:(2)gamma+gamma rarre^(-)+e^(+).:}\begin{equation} \gamma + \gamma \rightarrow e^- + e^+. \end{equation}
It is also possible to create electron-positron pairs through the scattering of relativistic particles for example the reaction in which a relativistic electron scatters off a proton (Thorne & Blandford 2017)
(3) e + p + e + p + + e + e + , (3) e + p + e + p + + e + e + , {:(3)e^(-)+p^(+)rarre^(-)+p^(+)+e^(-)+e^(+)",":}\begin{equation} e^- + p^+ \rightarrow e^- + p^+ + e^- + e^+, \end{equation}
(where p + p + p^(+)p^+ are protons) or through other interactions. The rate at which pair production occurs was first estimated by Landau & Lifshitz (1934) in which pair production was considered as a result of collisions of various kinds. They found that the cross section for producing an electron-positron pair scales as the cube of the logarithm of the energy of the colliding particles partially setting the stage for the discussion of the energy limit where pair production becomes important to come in the following sections.
Pair production as a high energy process becomes important for some astrophysical plasmas of very high temperature such as the plasma composing the early universe or in select astrophysical transient events that produce the requisite temperature and densities such as the black hole accretion disks that form in the centre of massive rotating stars during core collapse (collapsar disks) (Woosley 1993).
In the fluid limit, pair production's main effect is modifying the composition of the plasma, which in turn modifies how the pressure relates to the density in the equation of state.
The paper will be structured as follows I will discuss the pair production limit in section 2, I will discuss the effect of pair-production in the fluid limit in section 3, and a numerical experiment designed to demonstrate the effects of pair production in the fluid limit and its results in section 4. Finally, I will summarize the results in section 5.

2. Pair Production Limit

We will first discuss the pair production energy limit followed by a discussion of where pair production becomes important in a plasma.

2.1. Energy Limit

The production of electron-positron pairs can occur by several different production reactions Equations [1]-[3] as well as others. Typically when discussing the energy limit for electron-positron pair production we consider the reaction in Equation [1]. In order to have energy conservation in the reaction, the incoming photon must have at minimum the mass-energy of the two produced particles. Any excess of that becomes kinetic energy of the created particles. So we can estimate this energy requirement:
(4) E e = m e c 2 = ( 9.11 × 10 31 k g ) ( 3 × 10 8 m / s ) 2 = 8.199 × 10 14 J 5.12 × 10 5 e V = 0.511 M e V (4) E e = m e c 2 = ( 9.11 × 10 31 k g ) ( 3 × 10 8 m / s ) 2 = 8.199 × 10 14 J 5.12 × 10 5 e V = 0.511 M e V {:(4){:[E_(e)=m_(e)c^(2)],[=(9.11 xx10^(-31)kg)(3xx10^(8)m//s)^(2)],[=8.199 xx10^(-14)J≃5.12 xx10^(5)eV],[=0.511MeV]:}:}\begin{equation} \begin{split} E_e &= m_e c^2 \\ &= (9.11\times 10^{-31}~\rm{kg})(3\times 10^8~\rm{m/s})^2 \\ &= 8.199\times 10^{-14}~\rm{J} \simeq 5.12\times 10^5~\rm{eV} \\ &= 0.511~\rm{MeV} \end{split} \end{equation}
where E e E e E_(e)E_e is the mass energy of the electron, m e m e m_(e)m_e is the mass of the electron, and c c cc is the speed of light. Thus as we need a single photon to produce both an electron and a positron, the true energy limit is
(5) E p p = 2 E e 1.022 M e V (5) E p p = 2 E e 1.022 M e V {:(5)E_(pp)=2E_(e)≃1.022MeV:}\begin{equation} E_{pp} = 2E_e \simeq 1.022~MeV \end{equation}
where E p p E p p E_(pp)E_{pp} is the e e + e e + e^(-)e^(+)e^-e^+-pair production energy limit. This energy corresponds to a plasma temperature of T 10 10 T 10 10 T∼10^(10)T\sim 10^{10} K, however this does not tell the whole story as photons produced within the plasma will have a distribution of energies (Planck's Law, Figure 1), some with the energy required to produce e e + e e + e^(-)e^(+)e^-e^+-pairs and others will not, thus pairs may play an important role in plasma at lower temperatures than this limit. Figure 1 shows the spectral radiance which can be thought of as a proxy for the distribution of photon energies produced by a black-body of temperature T = 10 9 T = 10 9 T=10^(9)T=10^9 K, as a function of wavelength ( 1 / E 1 / E prop1//E\propto 1/E). This lies below the temperature limit, however the distribution of energies present in the photons suggests that a subset of photons will have sufficient energy to produce pairs.

2.2. n-T Limit

Following the prescription in Pathria & Beale (2011), we can consider the role of e e + e e + e^(-)e^(+)e^-e^+-pairs as the universe cooled initially from temperatures of T 1 × 10 11 T 1 × 10 11 T∼1xx10^(11)T\sim 1\times 10^{11} K to much cooler than the pair production temperature limit.
image Figure 1 Spectral radiance as a function of wavelength for a black-body of temperature T = 10 9 T = 10 9 T=10^(9)T=10^9 K (black). The pair production energy limit is shown in red.
As the plasma about a second after the Big bang cools and approaches the cross over temperature for e e + e e + e^(-)e^(+)e^-e^+-pair production the relativistic energy is given by (Pathria & Beale 2011):
(6) ε k = ( c k ) 2 + ( m e c 2 ) 2 (6) ε k = ( c k ) 2 + ( m e c 2 ) 2 {:(6)epsi_( vec(k))=sqrt((ℏck)^(2)+(m_(e)c^(2))^(2)):}\begin{equation} \varepsilon_{\vec{k}} = \sqrt{(\hbar c k)^2 + (m_e c^2)^2} \end{equation}
(where ε k ε k epsi_( vec(k))\varepsilon_{\vec{k}} is the relativistic energy, \hbar is Planck's constant, and k k kk is the wavenumber) we can use this to calculate the density of states a e ( ε ) a e ( ε ) a_(e)(epsi)a_e(\varepsilon) that satisfy this energy (dispersion relation) by integrating over all of k k kk-space (assuming spherical symmetry):
(7) a e ( ε ) = 1 ( 2 π ) 3 d 3 k δ ( ε ε k ) = 8 π ( 2 π ) 3 0 d k k 2 δ ( ε ε k ) = ε ε 2 ( m e c 2 ) 2 π 2 ( c ) 3 (7) a e ( ε ) = 1 ( 2 π ) 3 d 3 k δ ( ε ε k ) = 8 π ( 2 π ) 3 0 d k k 2 δ ( ε ε k ) = ε ε 2 ( m e c 2 ) 2 π 2 ( c ) 3 {:(7){:[a_(e)(epsi)=(1)/((2pi)^(3))intd^(3)k delta(epsi-epsi _(k))],[=(8pi)/((2pi)^(3))int_(0)^(oo)dkk^(2)delta(epsi-epsi _(k))],[=(epsisqrt(epsi^(2)-(m_(e)c^(2))^(2)))/(pi^(2)(ℏc)^(3))]:}:}\begin{equation} \begin{split} a_e(\varepsilon) &= \frac{1}{(2\pi)^3} \int d^3k \delta(\varepsilon - \varepsilon_k) \\ &= \frac{8\pi}{(2\pi)^3} \int_0^\infty dk k^2 \delta(\varepsilon - \varepsilon_k) \\ &= \frac{\varepsilon\sqrt{\varepsilon^2-(m_ec^2)^2}}{\pi^2(\hbar c)^3} \end{split} \end{equation}
If we assume that e e + e e + e^(-)e^(+)e^-e^+-pairs are being produced in statistical equilibrium we have the reaction:
(8) N + γ N + e + e + (8) N + γ N + e + e + {:(8)N+gamma harr N+e^(-)+e^(+):}\begin{equation} N + \gamma \leftrightarrow N + e^- + e^+ \end{equation}
using the assumption of statistical equilibrium we have:
(9) μ + μ + = μ γ = 0 (9) μ + μ + = μ γ = 0 {:(9)mu_(-)+mu_(+)=mu _(gamma)=0:}\begin{equation} \mu_- + \mu_+ = \mu_\gamma = 0 \end{equation}
where μ x μ x mu _(x)\mu_x are the chemical potentials of each individual species, and thus
(10) μ + = μ . (10) μ + = μ . {:(10)mu_(+)=-mu_(-).:}\begin{equation} \mu_+ = - \mu_-. \end{equation}
As we would like to discuss the number densities of electrons and positrons as we want to analyze the pair production limit, it is instructive to compute the ratio of number densities n ± / n γ n ± / n γ n_(+-)//n_( gamma)n_\pm/n_\gamma Pathria & Beale (2011). To do this we have the integrals:
(11) n ± n γ = β m e c 2 a e ( ε ) / n γ e ( ε μ ± ) / ( k B T ) d x β (11) n ± n γ = β m e c 2 a e ( ε ) / n γ e ( ε μ ± ) / ( k B T ) d x β {:(11)(n_(+-))/(n_( gamma))=int_(betam_(e)c^(2))^(oo)(a_(e)(epsi)//n_( gamma))/(e^((epsi-mu_(+-))//(k_(B)T)))(dx)/(beta):}\begin{equation} \frac{n_\pm}{n_\gamma} = \int_{\beta m_ec^2}^\infty \frac{a_e(\varepsilon)/n_\gamma}{e^{(\varepsilon-\mu_\pm)/(k_BT)}} \frac{dx}{\beta} \end{equation}
where k B k B k_(B)k_B is the boltzmann constant, and β = 1 / ( k B T ) β = 1 / ( k B T ) beta=1//(k_(B)T)\beta = 1/(k_BT). If we define the number density of photons to be the number density of photons in a blackbody enclosure:
(12) n γ ( T ) = 2 ζ ( 3 ) π 2 ( k T c ) 3 (12) n γ ( T ) = 2 ζ ( 3 ) π 2 k T c 3 {:(12)n_(gamma)(T)=(2zeta(3))/(pi^(2))((kT)/(ℏc))^(3):}\begin{equation} n_{\gamma}(T)=\frac{2 \zeta(3)}{\pi^{2}}\left(\frac{k T}{\hbar c}\right)^{3} \end{equation}
where ζ ( 3 ) = 1.20206 ζ ( 3 ) = 1.20206 zeta(3)=1.20206\zeta(3) = 1.20206 and we use μ + = μ μ + = μ mu_(+)=-mu_(-)\mu_+ = -\mu_-, then we have the integrals:
(13) n n γ = 1 2 ζ ( 3 ) β m e c 2 x x 2 ( β m e c 2 ) 2 e x e β μ + 1 d x n + n γ = 1 2 ζ ( 3 ) β m e c 2 x x 2 ( β m e c 2 ) 2 e x e β μ + 1 d x (13) n n γ = 1 2 ζ ( 3 ) β m e c 2 x x 2 β m e c 2 2 e x e β μ + 1 d x n + n γ = 1 2 ζ ( 3 ) β m e c 2 x x 2 β m e c 2 2 e x e β μ + 1 d x {:(13){:[(n_(-))/(n_(gamma))=(1)/(2zeta(3))int_(betam_(e)c^(2))^(oo)(xsqrt(x^(2)-(betam_(e)c^(2))^(2)))/(e^(x)e^(-betamu_(-))+1)dx],[(n_(+))/(n_(gamma))=(1)/(2zeta(3))int_(betam_(e)c^(2))^(oo)(xsqrt(x^(2)-(betam_(e)c^(2))^(2)))/(e^(x)e^(betamu_(-))+1)dx]:}:}\begin{equation} \begin{split} \frac{n_{-}}{n_{\gamma}}&=\frac{1}{2 \zeta(3)} \int_{\beta m_{e} c^{2}}^{\infty} \frac{x \sqrt{x^{2}-\left(\beta m_{e} c^{2}\right)^{2}}}{e^{x} e^{-\beta \mu_{-}}+1} d x \\ \frac{n_{+}}{n_{\gamma}}&=\frac{1}{2 \zeta(3)} \int_{\beta m_{e} c^{2}}^{\infty} \frac{x \sqrt{x^{2}-\left(\beta m_{e} c^{2}\right)^{2}}}{e^{x} e^{\beta \mu_{-}}+1} d x \end{split} \end{equation}
As these integrals are both non-trivial they will need to be numerically solved. Also each depends on μ μ mu_(-)\mu_- which is an unknown here, so we will first solve for the electron chemical potential. In order to solve for this we will use the assumption of charge neutrality which says that the number density of electrons minus positrons is equal to the number density of protons.
(14) n n + = n p = ( 1 q ) n B (14) n n + = n p = ( 1 q ) n B {:(14)n_(-)-n_(+)=n_(p)=(1-q)n_(B):}\begin{equation} n_- - n_+ = n_p = (1-q)n_B \end{equation}
where n B n B n_(B)n_B is the number density of baryons, and q q qq is the ratio of neutrons to baryons and is given by:
(15) q = n n n B = 1 e β Δ ε + 1 . (15) q = n n n B = 1 e β Δ ε + 1 . {:(15)q=(n_(n))/(n_(B))=(1)/(e^(beta Delta epsi)+1).:}\begin{equation} q = \frac{n_n}{n_B} = \frac{1}{e^{\beta \Delta \varepsilon}+1}. \end{equation}
where Δ ε Δ ε Delta epsi\Delta \varepsilon is mass energy difference between neutron and proton. Thus we can compute the number density ratio :
(16) ( n n + ) n γ = = sinh ( β μ ) 2 ζ ( 3 ) β m e c 2 x x 2 ( β m e c 2 ) 2 cosh ( x ) + cosh ( β μ ) d x = ( 1 q ) η (16) n n + n γ = = sinh β μ 2 ζ ( 3 ) β m e c 2 x x 2 β m e c 2 2 cosh ( x ) + cosh β μ d x = ( 1 q ) η {:(16){:[((n_(-)-n_(+)))/(n_(gamma))=],[=(sinh(betamu_(-)))/(2zeta(3))int_(betam_(e)c^(2))^(oo)(xsqrt(x^(2)-(betam_(e)c^(2))^(2)))/(cosh(x)+cosh(betamu_(-)))dx],[=(1-q)eta]:}:}\begin{equation} \begin{split} &\frac{\left(n_{-}-n_{+}\right)}{n_{\gamma}}= \\ &=\frac{\sinh \left(\beta \mu_{-}\right)}{2 \zeta(3)} \int_{\beta m_{e} c^{2}}^{\infty} \frac{x \sqrt{x^{2}-\left(\beta m_{e} c^{2}\right)^{2}}}{\cosh (x)+\cosh \left(\beta \mu_{-}\right)} d x \\ &=(1-q) \eta \end{split} \end{equation}
Where η η eta\eta is the ratio of baryons to photons in the universe and was estimated by WMAP (Bennett et al. 2003) to be η 6 × 10 10 η 6 × 10 10 eta≃6xx10^(-10)\eta \simeq 6\times 10^{-10}. Given this we can make the substitution z = x / ( β m e c 2 ) z = x / ( β m e c 2 ) z=x//(betam_(e)c^(2))z = x/(\beta m_ec^2), and factor out cosh ( β μ ) cosh ( β μ ) cosh(betamu_(-))\cosh(\beta \mu_-) from the denominator such that we have the form:
(17) ( n n + ) n γ = = tanh ( β μ ) 2 ζ ( 3 ) 1 z z 2 1 cosh ( β z ) cosh ( β μ ) + 1 d x = ( 1 q ) η (17) n n + n γ = = tanh β μ 2 ζ ( 3 ) 1 z z 2 1 cosh ( β z ) cosh β μ + 1 d x = ( 1 q ) η {:(17){:[((n_(-)-n_(+)))/(n_(gamma))=],[=(tanh(betamu_(-)))/(2zeta(3))int_(1)^(oo)(zsqrt(z^(2)-1))/((cosh(beta z))/(cosh(betamu_(-)))+1)dx],[=(1-q)eta]:}:}\begin{equation} \begin{split} &\frac{\left(n_{-}-n_{+}\right)}{n_{\gamma}}= \\ &=\frac{\tanh \left(\beta \mu_{-}\right)}{2 \zeta(3)} \int_{1}^{\infty} \frac{z \sqrt{z^{2}-1}}{\frac{\cosh (\beta z)}{\cosh \left(\beta \mu_{-}\right)}+1} d x \\ &=(1-q) \eta \end{split} \end{equation}
Equation [16] for μ μ mu_(-)\mu_-. This is a tricky numerical integral and care must be taken to avoid overflowing cosh ( X ) cosh ( X ) cosh(X)\cosh(X) and exp ( X ) exp ( X ) exp(X)\exp(X) numerically. To do this I implemented limiters on cosh ( β z ) / cosh ( β μ ) cosh ( β z ) / cosh ( β μ ) cosh(beta z)//cosh(betamu_(-))\cosh(\beta z)/\cosh(\beta\mu_-) such that if either β z 1 β z 1 beta z≫1\beta z \gg 1 or β μ 1 β μ 1 betamu_(-)≫1\beta \mu_-\gg 1, then cosh ( X ) 0.5 e X cosh ( X ) 0.5 e X cosh(X)rarr0.5*e^(X)\cosh(X) \rightarrow 0.5\cdot e^{X} or if both β z 1 β z 1 beta z≫1\beta z \gg 1 and β μ 1 β μ 1 betamu_(-)≫1\beta \mu_-\gg 1 then cosh ( X ) / cosh ( Y ) = e X Y cosh ( X ) / cosh ( Y ) = e X Y cosh(X)//cosh(Y)=e^(X-Y)\cosh(X)/\cosh(Y) = e^{X-Y}. Also, to avoid other numerical difficulties, the integral is performed from one to some large number beyond convergence of the integrand with zero.
As a result we have a numerical solution for μ ( T ) μ ( T ) mu_(-)(T)\mu_-(T). Finally, we use this to calculate the electron and positron number density ratios (with photon number density) by numerically solving the integrals in equation [13] with similar limiters to the above integral. Finally we multiply each result by the number density of photons (Equation [12]) to get the number densities of electron and positrons as a function of temperature (Figure 2). In the figure we can see how the number density of electrons and positrons vary with temperature. On the right, the red vertical line represents the temperature equivalent ( T = E p p / k B T = E p p / k B T=E_(pp)//k_(B)T = E_{pp}/k_B) of the pair production energy limit. Pair production however is not a process that is switched on and of at that energy, that is simply the energy that a single photon needs to have to produce a pair. Within the plasma there will be photons with a variety of energies, thus even at lower temperatures than that limit we will see significant pair production. In fact we see that the electron and positron number densities remain very similar as pair production remains important above a temperature of T 3 × 10 8 T 3 × 10 8 T∼3xx10^(8)T\sim 3\times 10^{8} K, where the electron number density levels off, and the positron number density dips. This result maintains the key features (leveling off of n e n e n_(e)n_e at 3 × 10 8 3 × 10 8 ∼3xx10^(8)\sim 3\times 10^{8} K and monotonic simultaneous increase of the number densities of electrons and positrons at higher temperatures) of the result found in Pathria & Beale 2011, however the spike in number densities at 6 × 10 8 6 × 10 8 ∼6xx10^(8)\sim 6\times 10^8 K, and 8 × 10 9 8 × 10 9 ∼8xx10^(9)\sim 8\times 10^{9} K I attribute to numerical errors as in order to perform this numerical integral several limiters were necessary, as described previously, and the solver has an increase in error at these points.
image Figure 2 Number density of electrons (solid) and positrons (dashed) as a function of temperature. The red line represents the photon energy limit in terms of temperature for e e + e e + e^(-)e^(+)e^-e^+-pair production. The number densities begin to diverge at 3 × 10 8 3 × 10 8 ∼3xx10^(8)\sim 3\times 10^{8} K where pair production becomes unimportant so the positron number density decreases significantly while the electron number density levels off.

3. Fluid Limit

While electron-positron pair production has the potential to cause interesting consequences in plasmas on small length scales , my research predominantly makes use of the Magneto-hydrodynamics code FLASH. FLASH models plasmas in the fluid limit ( L λ m f p r L L λ m f p r L L≫lambda_(mfp)≫r_(L)L\gg \lambda_{mfp} \gg r_L where L L LL is the cell size, λ m f p λ m f p lambda_(mfp)\lambda_{mfp} is the mean free path, and r L r L r_(L)r_L is the Larmor radius), so we consider motion on length scales much longer than the mean free path which in turn is much larger than the Larmor radius. This means we ignore small scale effects such as the motions of individual particles within the fluid. Instead fluid is discretized into individual "cells" within some domain we consider. The fluid in these cells are evolved using the equations of compressible ideal magneto-hydrodynamics (written non-dimensionally assuming heat conduction is negligible, the fluid is non-viscous, and assuming we have no gravitational potential) (Fryxell et al. 2000, Dubey et al. 2009):
(18) ρ t + ( ρ v ) = 0 ρ v t + ( ρ v v B B ) + p = 0 ρ E t + ( v ( ρ E + p ) B ( v B ) ) = 0 B t + ( v B B v ) = 0 (18) ρ t + ( ρ v ) = 0 ρ v t + ( ρ v v B B ) + p = 0 ρ E t + ( v ( ρ E + p ) B ( v B ) ) = 0 B t + ( v B B v ) = 0 {:(18){:[(del rho)/(del t)+grad*(rho vec(v))=0],[(del rho( vec(v)))/(del t)+grad*(rho vec(v) vec(v)- vec(B) vec(B))+gradp_(**)=0],[(del rho E)/(del t)+grad*( vec(v)(rho E+p_(**))- vec(B)( vec(v)* vec(B)))=0],[(del( vec(B)))/(del t)+grad*( vec(v) vec(B)- vec(B) vec(v))=0]:}:}\begin{equation} \begin{split} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) &=0 \\ \frac{\partial \rho \vec{v}}{\partial t} + \nabla \cdot (\rho \vec{v} \vec{v} - \vec{B} \vec{B}) + \nabla p_* &= 0 \\ \frac{\partial \rho E}{\partial t} + \nabla \cdot(\vec{v}(\rho E + p_*) - \vec{B}(\vec{v}\cdot \vec{B})) &= 0 \\ \frac{\partial \vec{B}}{\partial t} + \nabla \cdot(\vec{v}\vec{B} - \vec{B}\vec{v}) &= 0 \\ \end{split} \end{equation}
where
(19) p = p + B 2 2 E = 1 2 v 2 + ϵ + 1 2 B 2 ρ (19) p = p + B 2 2 E = 1 2 v 2 + ϵ + 1 2 B 2 ρ {:(19){:[p_(**)=p+(B^(2))/(2)],[E=(1)/(2)v^(2)+epsilon+(1)/(2)(B^(2))/(rho)]:}:}\begin{equation} \begin{split} p_* &= p + \frac{B^2}{2} \\ E &= \frac{1}{2}v^2 + \epsilon + \frac{1}{2} \frac{B^2}{\rho} \\ \end{split} \end{equation}
Here ρ ρ rho\rho is density, t t tt is time, v v vec(v)\vec{v} is the fluid velocity, B B BB is the magnetic field, and p p p_(**)p_* is the total pressure (fluid thermal pressure p p pp plus magnetic pressure B 2 / 2 B 2 / 2 B^(2)//2B^2/2). This system of equations is closed by an Equation of state:
(20) P = P ( ϵ , ρ ) (20) P = P ( ϵ , ρ ) {:(20)P=P(epsilon","rho):}\begin{equation} P = P(\epsilon,\rho) \end{equation}
So how would pair production affect a plasma in the fluid limit? Well, if we choose a composition dependent equation of state:
(21) P = P ( ϵ , ρ ) P ( ϵ , ρ , Y e ) (21) P = P ( ϵ , ρ ) P ( ϵ , ρ , Y e ) {:(21)P=P(epsilon","rho)rarr P(epsilon","rho","Y_(e)):}\begin{equation} P = P(\epsilon,\rho) \rightarrow P(\epsilon,\rho,Y_e) \end{equation}
such as the helmholtz equation of state that depends on the composition in terms of an electron fraction, the creation of pairs will modify the composition of the fluid thus changing the electron fraction, and modifying the regime of the Equation of state we use to determine the pressure from the internal energy and density.

3.1. Helmholtz Equation of State

The Helmholtz equation of state (Timmes & Swesty 2000 includes contributions of radiation, electron pressure and positron pressure as well as nuclei and coulomb corrections. This makes it useful in regimes where electrons and positrons are relativistic or degenerate, or if radiation may contribute significantly to the pressure. The total pressure and internal energy are computed as:
(22) P t o t = P r a d + P i o n + P e l e + P p o s + P c o u l ϵ t o t = ϵ r a d + ϵ i o n + ϵ e l e + ϵ p o s + ϵ c o u l (22) P t o t = P r a d + P i o n + P e l e + P p o s + P c o u l ϵ t o t = ϵ r a d + ϵ i o n + ϵ e l e + ϵ p o s + ϵ c o u l {:(22){:[P_(tot)=P_(rad)+P_(ion)+P_(ele)+P_(pos)+P_(coul)],[epsilon_(tot)=epsilon_(rad)+epsilon_(ion)+epsilon_(ele)+epsilon_(pos)+epsilon_(coul)]:}:}\begin{equation} \begin{split} P_{tot} &= P_{rad} + P_{ion} + P_{ele} + P_{pos} + P_{coul} \\ \epsilon_{tot} &= \epsilon_{rad} + \epsilon_{ion} + \epsilon_{ele} + \epsilon_{pos} + \epsilon_{coul} \\ \end{split} \end{equation}
Where the "rad" component assumes a blackbody in local thermodynamic equilibrium:
(23) P r a d = a T 4 3 ϵ r a d = 3 P r a d ρ (23) P r a d = a T 4 3 ϵ r a d = 3 P r a d ρ {:(23){:[P_(rad)=(aT^(4))/(3)],[epsilon_(rad)=(3P_(rad))/(rho)]:}:}\begin{equation} \begin{split} P_{rad} &= \frac{aT^4}{3} \\ \epsilon_{rad} &= \frac{3P_{rad}}{\rho} \\ \end{split} \end{equation}
the "ion" component is is assumed to act like an ideal gas with γ = 5 / 3 γ = 5 / 3 gamma=5//3\gamma=5/3, and the electrons and ions are treated as non-interacting fermi-gases, which in the fluid limit is a reasonable assumption.
The number densities of electrons and positrons are derived by integrating the Fermi-Dirac (FD) distribution across a range of energies, and as they are used in the Helmholtz equation of state, they are written as:
(24) n = 8 π 2 h 3 m e 3 c 3 α 3 / 2 [ F 1 / 2 ( η , α ) + F 3 / 2 ( η , α ) ] n + = 8 π 2 h 3 m e 3 c 3 α 3 / 2 [ F 1 / 2 ( η 2 / α , α ) + α F 3 / 2 ( η 2 / α , α ) ] (24) n = 8 π 2 h 3 m e 3 c 3 α 3 / 2 F 1 / 2 ( η , α ) + F 3 / 2 ( η , α ) n + = 8 π 2 h 3 m e 3 c 3 α 3 / 2 F 1 / 2 ( η 2 / α , α ) + α F 3 / 2 ( η 2 / α , α ) {:(24){:[n_(-)=(8pisqrt2)/(h^(3))m_(e)^(3)c^(3)alpha^(3//2)[F_(1//2)(eta,alpha)+F_(3//2)(eta,alpha)]],[n_(+)=(8pisqrt2)/(h^(3))m_(e)^(3)c^(3)alpha^(3//2)[F_(1//2)(-eta-2//alpha,alpha):}],[{:+alphaF_(3//2)(-eta-2//alpha,alpha)]]:}:}\begin{equation} \begin{split} n_- &= \frac{8\pi \sqrt{2}}{h^3} m_e^3 c^3 \alpha^{3/2}\left[\mathcal{F}_{1/2}(\eta,\alpha) + \mathcal{F}_{3/2}(\eta,\alpha)\right] \\ n_+ &= \frac{8\pi \sqrt{2}}{h^3} m_e^3 c^3 \alpha^{3/2}\left[\mathcal{F}_{1/2}(-\eta-2/\alpha,\alpha)\right. \\ &\left.+ \alpha \mathcal{F}_{3/2}(-\eta-2/\alpha,\alpha)\right] \\ \end{split} \end{equation}
where η = μ / ( k B T ) η = μ / ( k B T ) eta=mu//(k_(B)T)\eta = \mu/(k_BT), and α = 1 β m e c 2 α = 1 β m e c 2 alpha=(1)/(betam_(e)c^(2))\alpha = \frac{1}{\beta m_e c^2}. Here F k ( η , α ) F k ( η , α ) F_(k)(eta,alpha)\mathcal{F}_k(\eta,\alpha) are the Fermi-Dirac integrals given by:
(25) F k ( η , α ) = 0 x k ( 1 + 0.5 α x ) 1 / 2 d x exp ( x η ) + 1 (25) F k ( η , α ) = 0 x k ( 1 + 0.5 α x ) 1 / 2 d x exp ( x η ) + 1 {:(25)F_(k)(eta","alpha)=int_(0)^(oo)(x^(k)(1+0.5 alpha x)^(1//2)dx)/(exp(x-eta)+1):}\begin{equation} \mathcal{F}_k(\eta,\alpha) = \int_0^\infty \frac{x^k (1+0.5\alpha x)^{1/2}dx}{\exp(x-\eta)+1} \end{equation}
charge neutrality requires that n , m a t t e r = n n + n , m a t t e r = n n + n_(-,matter)=n_(-)-n_(+)n_{-,matter} = n_--n_+ solving this equation numerically using Equation [24] for η η eta\eta the re-scaled chemical potential, then the Fermi-Dirac integrals are evaluated, the number densities can be solved for using Equation [24], and the partial pressures using:
(26) P = 2 K 3 m e c 2 α 5 / 2 [ F 3 / 2 ( η , α ) + α 2 F 5 / 2 ( η , α ) ] P + = 2 K 3 m e c 2 α 5 / 2 [ F 3 / 2 ( η 2 / α , α ) + α 2 F 5 / 2 ( η 2 / α , α ) ] (26) P = 2 K 3 m e c 2 α 5 / 2 F 3 / 2 ( η , α ) + α 2 F 5 / 2 ( η , α ) P + = 2 K 3 m e c 2 α 5 / 2 F 3 / 2 ( η 2 / α , α ) + α 2 F 5 / 2 ( η 2 / α , α ) {:(26){:[P_(-)=(2K)/(3)m_(e)c^(2)alpha^(5//2)[F_(3//2)(eta,alpha)+(alpha)/(2)F_(5//2)(eta,alpha)]],[P_(+)=(2K)/(3)m_(e)c^(2)alpha^(5//2)[F_(3//2)(-eta-2//alpha,alpha):}],[{:+(alpha)/(2)F_(5//2)(-eta-2//alpha,alpha)]]:}:}\begin{equation} \begin{split} P_- &= \frac{2K}{3} m_ec^2 \alpha^{5/2}\left[ \mathcal{F}_{3/2}(\eta,\alpha) +\frac{\alpha}{2}\mathcal{F}_{5/2}(\eta,\alpha)\right] \\ P_+ &= \frac{2K}{3} m_ec^2 \alpha^{5/2}\left[ \mathcal{F}_{3/2}(-\eta-2/\alpha,\alpha)\right. \\ &\left.+\frac{\alpha}{2}\mathcal{F}_{5/2}(-\eta-2/\alpha,\alpha)\right] \\ \end{split} \end{equation}
So in summary, the addition of electrons and positrons to a plasma in the fluid modifies the composition, and thus modifies the electron fraction of the plasma. Changing the electron fraction, the number density of electrons and positrons will change which modifies the chemical potential. Finally with the modified chemical potential, the partial pressures from the positrons and electrons will change (Equation [26]) and this can modify macroscopic quantities such as the wave speeds.

4. Test Case

In order to test the effect of pair production on a plasma in the fluid limit we need to find a macroscopic effect (length scales much larger than the plasma gyro-radius) that changes due to the inclusion or exclusion of pairs in our plasma. In order to illustrate this I have chosen to model the Fast-Magnetosonic Wave in the MHD fluid code FLASH. As we derived in Assignment 2, the phase speed of the fast magnetosonic wave is given by:
(27) v p 2 = 1 2 ( v A 2 + c s 2 ) + 1 2 [ ( v A 2 c s 2 ) 2 + 4 v A 2 c s 2 sin 2 θ ] 1 / 2 (27) v p 2 = 1 2 ( v A 2 + c s 2 ) + 1 2 ( v A 2 c s 2 ) 2 + 4 v A 2 c s 2 sin 2 θ 1 / 2 {:(27)v_(p)^(2)=(1)/(2)(v_(A)^(2)+c_(s)^(2))+(1)/(2)[(v_(A)^(2)-c_(s)^(2))^(2)+4v_(A)^(2)c_(s)^(2)sin^(2)theta]^(1//2):}\begin{equation} v_p^2 = \frac{1}{2}(v_A^2 + c_s^2) + \frac{1}{2}\left[(v_A^2 - c_s^2)^2 + 4v_A^2c_s^2\sin^2\theta\right]^{1/2} \end{equation}
where c s c s c_(s)c_s is the sound speed, v A v A v_(A)v_A is the alfven speed and θ θ theta\theta is the angle between the wave vector and the background magnetic field. For a constant composition plasma, γ = constant γ = constant gamma="constant"\gamma = \text{constant} we have:
(28) c s = γ P ρ v A = B 0 2 4 π ρ 0 (28) c s = γ P ρ v A = B 0 2 4 π ρ 0 {:(28){:[c_(s)=sqrt(gamma(P)/( rho))],[v_(A)=sqrt((B_(0)^(2))/(4pirho_(0)))]:}:}\begin{equation} \begin{split} c_s &= \sqrt{\gamma \frac{P}{\rho}} \\ v_A &= \sqrt{\frac{B_0^2}{4\pi \rho_0}} \\ \end{split} \end{equation}
Thus by changing the composition of the plasma (in pair-production regime or out of it), the value of γ γ gamma\gamma will change and thus the sound speed will change. This is then reflected in the Fast magneto-sonic wave speed. The test case discussed in section 4.1 is based on the test case implemented in athena++, however that is limited to a gamma-law equation of state and thus cannot account for the effect of composition. In the fluid code FLASH we can use a composition dependent equation of state (helmholtz) for MHD simulations, thus I will implement and modify the test case from athena++ into FLASH.

4.1. Fast-Magnetosonic Wave Simulation

In order to model a fast magnetosonic wave in the MHD Fluid code FLASH, I implement the test case for a 1-3 dimensional Linear wave in MHD in the fluid code athena++ (Stone et al. 2020). Following the description in Stone et al. (2008), the equations of magneto-hydrodynamics can be linearized to form a vector equation of the form:
(29) t U = A x U (29) t U = A x U {:(29)del _(t) vec(U)= vec(A)del _(x) vec(U):}\begin{equation} \partial_t \vec{U} = \vec{A} \partial_x \vec{U} \end{equation}
Where U U vec(U)\vec{U} is the vector of conserved variables:
(30) U = ( ρ ρ v x ρ v y ρ v z E b y b z ) (30) U = ρ ρ v x ρ v y ρ v z E b y b z {:(30) vec(U)=([rho],[rhov_(x)],[rhov_(y)],[rhov_(z)],[E],[b_(y)],[b_(z)]):}\begin{equation} \vec{U} = \begin{pmatrix} \rho \\ \rho v_x \\ \rho v_y \\ \rho v_z \\ E \\ b_y \\ b_z \\ \end{pmatrix} \end{equation}
where b = B / 4 π b = B / 4 π vec(b)= vec(B)//sqrt(4pi)\vec{b} =\vec{B}/\sqrt{4\pi} and the matrix A A vec(A)\vec{A} is given by:
(31) A = ( 0 1 0 0 0 0 0 v x 2 + γ v 2 / 2 X ( γ 3 ) v x γ v y γ v z γ b y Y v x v y v y v x 0 0 b x 0 v x v z v z 0 v x 0 0 b x A 51 A 52 A 53 A 54 γ v x A 56 A 57 ( b x v y b y v x ) / ρ b y / ρ b x / ρ 0 0 v x 0 ( b x v z b z v x ) / ρ b z / ρ 0 b x / ρ 0 0 v x ) (31) A = 0 1 0 0 0 0 0 v x 2 + γ v 2 / 2 X ( γ 3 ) v x γ v y γ v z γ b y Y v x v y v y v x 0 0 b x 0 v x v z v z 0 v x 0 0 b x A 51 A 52 A 53 A 54 γ v x A 56 A 57 ( b x v y b y v x ) / ρ b y / ρ b x / ρ 0 0 v x 0 ( b x v z b z v x ) / ρ b z / ρ 0 b x / ρ 0 0 v x {:(31){:[ vec(A)=],[([0,1,0,0,0,0,0],[-v_(x)^(2)+gamma^(')v^(2)//2-X^('),-(gamma-3)v_(x),-gamma^(')v_(y),-gamma^(')v_(z),gamma^('),-b_(y)Y^(')],[-v_(x)v_(y),v_(y),v_(x),0,0,-b_(x),0],[-v_(x)v_(z),v_(z),0,v_(x),0,0,-b_(x)],[A_(51),A_(52),A_(53),A_(54),gammav_(x),A_(56),A_(57)],[(b_(x)v_(y)-b_(y)v_(x))//rho,b_(y)//rho,-b_(x)//rho,0,0,v_(x),0],[(b_(x)v_(z)-b_(z)v_(x))//rho,b_(z)//rho,0,-b_(x)//rho,0,0,v_(x)])]:}:}\begin{equation} \begin{split} \vec{A} = \\ \tiny{\begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ -v_x^2 + \gamma'v^2/2-X' & -(\gamma-3)v_x & -\gamma'v_y & -\gamma'v_z & \gamma' & -b_y Y' \\ -v_xv_y & v_y & v_x & 0 & 0 & -b_x & 0 \\ -v_x v_z & v_z & 0 & v_x & 0 & 0 &-b_x \\ A_{51} & A_{52} & A_{53} & A_{54} & \gamma v_x & A_{56} & A_{57} \\ (b_xv_y-b_yv_x)/\rho & b_y/\rho & -b_x/\rho & 0 & 0 & v_x & 0 \\ (b_xv_z-b_zv_x)/\rho & b_z/\rho & 0 & -b_x/\rho & 0 & 0 & v_x \\ \end{pmatrix} } \end{split} \end{equation}
where v 2 = v v v 2 = v v v^(2)= vec(v)* vec(v)v^2 = \vec{v}\cdot \vec{v},
(32) A 51 = v x H + γ v x v 2 / 2 + b x ( b x v x + b y v y + b z v z ) / ρ v x X A 52 = γ v x 2 + H b x 2 / ρ A 53 = γ v x v y b x b y / ρ A 54 = γ v x v z b x b z / ρ A 56 = ( b x v y + b y v x Y ) A 57 = ( b x v z + b z v x Y ) X = [ ( b y , L b y , R ) 2 + ( b z , L b z , R ) 2 ] / ( 2 ( ρ L + ρ R ) ) Y = ρ L + ρ R 2 ρ (32) A 51 = v x H + γ v x v 2 / 2 + b x b x v x + b y v y + b z v z / ρ v x X A 52 = γ v x 2 + H b x 2 / ρ A 53 = γ v x v y b x b y / ρ A 54 = γ v x v z b x b z / ρ A 56 = b x v y + b y v x Y A 57 = b x v z + b z v x Y X = b y , L b y , R 2 + b z , L b z , R 2 / 2 ρ L + ρ R Y = ρ L + ρ R 2 ρ {:(32){:[{:A_(51)=-v_(x)H+gamma^(')v_(x)v^(2)//2+b_(x)(b_(x)v_(x)+b_(y)v_(y)+b_(z)v_(z))//rho-v_(x)X^('):}],[A_(52)=-gamma^(')v_(x)^(2)+H-b_(x)^(2)//rho],[A_(53)=-gamma^(')v_(x)v_(y)-b_(x)b_(y)//rho],[A_(54)=-gamma^(')v_(x)v_(z)-b_(x)b_(z)//rho],[A_(56)=-(b_(x)v_(y)+b_(y)v_(x)Y^('))],[A_(57)=-(b_(x)v_(z)+b_(z)v_(x)Y^('))],[X=[(b_(y,L)-b_(y,R))^(2)+(b_(z,L)-b_(z,R))^(2)]//(2(sqrt(rho_(L))+sqrt(rho_(R))))],[Y=(rho_(L)+rho_(R))/(2rho)]:}:}\begin{equation} \begin{gathered} \small{A_{51}=-v_{x} H+\gamma^{\prime} v_{x} v^{2} / 2+b_{x}\left(b_{x} v_{x}+b_{y} v_{y}+b_{z} v_{z}\right) / \rho-v_{x} X^{\prime}} \\ A_{52}=-\gamma^{\prime} v_{x}^{2}+H-b_{x}^{2} / \rho \\ A_{53}=-\gamma^{\prime} v_{x} v_{y}-b_{x} b_{y} / \rho \\ A_{54}=-\gamma^{\prime} v_{x} v_{z}-b_{x} b_{z} / \rho \\ A_{56}=-\left(b_{x} v_{y}+b_{y} v_{x} Y^{\prime}\right) \\ A_{57}=-\left(b_{x} v_{z}+b_{z} v_{x} Y^{\prime}\right) \\ X=\left[\left(b_{y, L}-b_{y, R}\right)^{2}+\left(b_{z, L}-b_{z, R}\right)^{2}\right] /\left(2\left(\sqrt{\rho_{L}}+\sqrt{\rho_{R}}\right)\right) \\ Y=\frac{\rho_{L}+\rho_{R}}{2 \rho} \end{gathered} \end{equation}
plus we have:
(33) γ = ( γ 1 ) X = ( γ 2 ) X Y = ( γ 2 ) Y H = ( E + P + b 2 / 2 ) / ρ (33) γ = ( γ 1 ) X = ( γ 2 ) X Y = ( γ 2 ) Y H = E + P + b 2 / 2 / ρ {:(33){:[gamma^(')=(gamma-1)],[X^(')=(gamma-2)X],[Y^(')=(gamma-2)Y],[H=(E+P+b^(2)//2)//rho]:}:}\begin{equation} \begin{split} \gamma^{\prime}&=(\gamma-1) \\ X^{\prime}&=(\gamma-2) X \\ Y^{\prime}&=(\gamma-2) Y \\ H&=\left(E+P+b^{2} / 2\right) / \rho \\ \end{split} \end{equation}
This matrix equation can be solved to yield eigenvectors and eigen values of the matrix equation. From these solutions, we can choose the appropriate eigenvector corresponding to the Fast magneto-sonic wave in our case. We can use this result to define the necessary perturbations to the conserved variables necessary to launch a Fast magneto-sonic wave in our plasma. If we define a perturbation to the conserved variables of the form (Gardiner & Stone 2005):
(34) U = U 0 + A R sin ( k x ) (34) U = U 0 + A R sin ( k x ) {:(34) vec(U)= vec(U)_(0)+A vec(R)sin(kx):}\begin{equation} \vec{U} = \vec{U}_0 + A\vec{R}\sin(kx) \end{equation}
where A A AA is the perturbation amplitude, R R RR is the right eigenvector corresponding to the wave mode we are selecting, and we define k k kk to be the wavenumber such that there is one wavelength across the domain we define. In Athena++ they use a toy model with round numbers and a gamma law equation of state.

4.2. Implementation in FLASH

To implement this in FLASH I first implemented the gamma-law test case from athena++ and compared the results between the two codes. One note is that athena++ implements conserved variables whereas FLASH uses primitive variables, thus a conversion is needed.
Next I chose two reference points to compare Fast Magneto-sonic waves. We want to compare the effect of e e + e e + e^(-)e^(+)e^-e^+-pair production on a high temperature and density gas, so for reference (as is relevant for my topic) I will use the approximate temperature and density of a collapsar disk (MacFadyen 2003):
(35) T 10 10 K ρ 10 9 g c m 3 (35) T 10 10 K ρ 10 9 g c m 3 {:(35){:[T∼10^(10)K],[rho∼10^(9)gcm^(-3)]:}:}\begin{equation} \begin{split} T&\sim 10^{10}~\rm{K} \\ \rho &\sim 10^{9}~\rm{g~cm}^{-3} \\ \end{split} \end{equation}
As the electron fraction of our plasma is defined:
(36) Y e = number~of~electrons number~of~baryons (36) Y e = number~of~electrons number~of~baryons {:(36)Y_(e)=("number~of~electrons")/("number~of~baryons"):}\begin{equation} Y_e = \frac{\text{number~of~electrons}}{\text{number~of~baryons}} \end{equation}
when pair production becomes important, the number of electrons in our plasma will increase, corresponding to an increase in Y e Y e Y_(e)Y_e, so for our two test cases we want to compare two states, one with a lower electron fraction corresponding to little pair production, and another with a higher electron fraction corresponding to more pair production. For demonstration purposes we choose two values of electron fraction within the limits of the Helmholtz equation of state ( Y e = 0.25 Y e = 0.25 Y_(e)=0.25Y_e = 0.25 for little pair production and Y e = 0.5 Y e = 0.5 Y_(e)=0.5Y_e = 0.5 for larger pair production).
To set the remainder of our initial conditions, I will compare to the initial conditions proposed for the Gamma-law equation of state test case where (Stone et al. 2020):
(37) ρ 0 = 1 P 0 = 1 / γ γ = 5 / 3 v x , 0 = v y , 0 = v z , 0 = 0 b x , 0 = 1 b y , 0 = 2 b z , 0 = 1 2 E 0 = P 0 ( γ 1 ) ρ 0 c s = 1 v A , x = 1 v F M S , x = 2 (37) ρ 0 = 1 P 0 = 1 / γ γ = 5 / 3 v x , 0 = v y , 0 = v z , 0 = 0 b x , 0 = 1 b y , 0 = 2 b z , 0 = 1 2 E 0 = P 0 ( γ 1 ) ρ 0 c s = 1 v A , x = 1 v F M S , x = 2 {:(37){:[rho_(0)=1],[P_(0)=1//gamma],[gamma=5//3],[v_(x,0)=v_(y,0)=v_(z,0)=0],[b_(x,0)=1],[b_(y,0)=sqrt2],[b_(z,0)=(1)/(2)],[E_(0)=(P_(0))/((gamma-1)rho_(0))],[c_(s)=1],[v_(A,x)=1],[v_(FMS,x)=2]:}:}\begin{equation} \begin{split} \rho_0 &= 1 \\ P_0 &= 1/\gamma \\ \gamma &= 5/3 \\ v_{x,0} = v_{y,0} = v_{z,0} &= 0 \\ b_{x,0} &= 1\\ b_{y,0} &= \sqrt{2} \\ b_{z,0} &= \frac{1}{2} \\ E_0 &= \frac{P_0}{(\gamma-1)\rho_0} \\ c_s &= 1 \\ v_{A,x} &= 1 \\ v_{FMS,x} &= 2 \\ \end{split} \end{equation}
The values are chosen such that the speed of the fast magneto-sonic wave is well separated from the sound speed. In our case we can attempt to replicate this while using a background density motivated by approximate collapsar disk values ρ 0 = 1 × 10 9 g / c m 3 ρ 0 = 1 × 10 9 g / c m 3 rho_(0)=1xx10^(9)g//cm^(3)\rho_0 = 1\times 10^{9}~g/cm^3, we use the initial temperature of T 0 = 10 10 T 0 = 10 10 T_(0)=10^(10)T_0 = 10^{10} K. Using this temperature, density and composition the EOS is used to calculate the initial pressure (Table [1]) and γ γ gamma\gamma values. We scale each magnetic field by a constant K = P 0 / ( b 0 b 0 ) K = P 0 / ( b 0 b 0 ) K=P_(0)//( vec(b_(0))* vec(b_(0)))K = P_0/(\vec{b_0}\cdot \vec{b_0}) where b 0 b 0 vec(b_(0))\vec{b_0} is the magnetic field vector from the gamma-law case and P 0 P 0 P_(0)P_0 is the background pressure from the helmholtz case, such that the alfv'en speed in the x-direction is comparable to the sound speed, and the relative magnitudes remain the same:
(38) b x , 0 = K G b y , 0 = 2 K G b z , 0 = K 2 G (38) b x , 0 = K G b y , 0 = 2 K G b z , 0 = K 2 G {:(38)b_(x,0)=K"G"b_(y,0)=sqrt2K"G"b_(z,0)=(K)/(2)"G":}\begin{equation} b_{x,0} = K~\text{G} b_{y,0} = \sqrt{2} K~\text{G} b_{z,0} = \frac{K}{2}~\text{G} \end{equation}
where K = 4 13 P 0 K = 4 13 P 0 K=(4)/(13)P_(0)K = \frac{4}{13}P_0. We once again set the background speeds to be zero v x , 0 = v y , 0 = v z , 0 = 0 v x , 0 = v y , 0 = v z , 0 = 0 v_(x,0)=v_(y,0)=v_(z,0)=0v_{x,0} = v_{y,0} = v_{z,0} = 0. Summarized in Table 1.
Case ρ 0 ρ 0 rho_(0)\rho_0 (g/cm 3 3 ^(3)^{3}) T 0 T 0 T_(0)T_0 (K) P 0 P 0 P_(0)P_0 (dyn/cm 2 2 ^(2)^{2}) γ γ gamma\gamma v i , 0 v i , 0 v_(i,0)v_{i,0} (cm/s) b x , 0 b x , 0 b_(x,0)b_{x,0} (G) b y , 0 b y , 0 b_(y,0)b_{y,0} (G) b z , 0 b z , 0 b_(z,0)b_{z,0} (G) Y e Y e Y_(e)Y_e v F M S v F M S v_(FMS)v_{FMS} (cm/s)
pair 10 9 10 9 10^(9)10^{9} 10 10 10 10 10^(10)10^{10} 8.51 × 10 26 8.51 × 10 26 8.51 xx10^(26)8.51\times 10^{26} 1.39 0 0 00 K K KK 2 K 2 K sqrt2K\sqrt{2}K K / 2 K / 2 K//2K/2 0.5 0.5 0.50.5 1.26829 × 10 9 1.26829 × 10 9 1.26829 xx10^(9)1.26829\times 10^{9}
no_pair 10 9 10 9 10^(9)10^{9} 10 10 10 10 10^(10)10^{10} 8.19 × 10 26 8.19 × 10 26 8.19 xx10^(26)8.19\times 10^{26} 1.50 0 0 00 K K KK 2 K 2 K sqrt2K\sqrt{2}K K / 2 K / 2 K//2K/2 0.25 0.25 0.250.25 1.26163 × 10 9 1.26163 × 10 9 1.26163 xx10^(9)1.26163\times 10^{9}
$ represents all directions ( i = x , y , z i = x , y , z i=x,y,zi=x,y,z), γ γ gamma\gamma is obtained from the Helmholtz EOS within the code based on the composition Y e Y e Y_(e)Y_e, and K = 4 13 P 0 K = 4 13 P 0 K=(4)/(13)P_(0)K = \frac{4}{13}P_0. }
Table 1: Background physical quantities values. $v_{i,0
We numerically solve the vector equation (Equation [29]) for the right eigenvector (used for scaling amplitudes of conserved quantity perturbations) in each of our cases as is done in Stone et al. (2020) and find:
(39) R p p = ( 0.935913 g cm 3 1.18768 × 10 9 g cm 2 s 1 8.8529 × 10 7 g cm 2 s 1 3.12997 × 10 7 g cm 2 s 1 2.17141 × 10 18 erg g cm 3 12511.4 G 4423.43 G ) R n o p p = ( 0.935913 g cm 3 1.16479 × 10 9 g cm 2 s 1 8.68223 × 10 7 g cm 2 s 1 3.06963 × 10 7 g cm 2 s 1 2.08849 × 10 18 erg g cm 3 12270.2 G 4338.16 G ) (39) R p p = 0.935913 g cm 3 1.18768 × 10 9 g cm 2 s 1 8.8529 × 10 7 g cm 2 s 1 3.12997 × 10 7 g cm 2 s 1 2.17141 × 10 18 erg  g cm 3 12511.4 G 4423.43 G R n o p p = 0.935913 g cm 3 1.16479 × 10 9 g cm 2 s 1 8.68223 × 10 7 g cm 2 s 1 3.06963 × 10 7 g cm 2 s 1 2.08849 × 10 18 erg  g cm 3 12270.2 G 4338.16 G {:(39){:[R_(pp)=([0.935913"g cm"^(-3)],[-1.18768 xx10^(9)"g cm"^(-2)s^(-1)],[8.8529 xx10^(7)"g cm"^(-2)s^(-1)],[3.12997 xx10^(7)"g cm"^(-2)s^(-1)],[2.17141 xx10^(18)"erg ""g cm"^(-3)],[12511.4G],[4423.43G])],[R_(nopp)=([0.935913"g cm"^(-3)],[-1.16479 xx10^(9)"g cm"^(-2)s^(-1)],[8.68223 xx10^(7)"g cm"^(-2)s^(-1)],[3.06963 xx10^(7)"g cm"^(-2)s^(-1)],[2.08849 xx10^(18)"erg ""g cm"^(-3)],[12270.2G],[4338.16G])]:}:}\begin{equation} \begin{split} R_{pp} &= \begin{pmatrix} 0.935913~\text{g cm}^{-3}\\ -1.18768\times 10^{9}~\text{g cm}^{-2} \rm{s}^{-1} \\ 8.8529\times 10^{7}~\text{g cm}^{-2} \rm{s}^{-1} \\ 3.12997\times 10^{7}~\text{g cm}^{-2} \rm{s}^{-1} \\ 2.17141\times 10^{18}~\text{erg }\text{g cm}^{-3} \\ 12511.4~\rm{G}\\ 4423.43~\rm{G} \end{pmatrix} \\ R_{no~pp} &= \begin{pmatrix} 0.935913 ~\text{g cm}^{-3}\\ -1.16479\times 10^{9} ~\text{g cm}^{-2} \rm{s}^{-1}\\ 8.68223\times 10^{7} ~\text{g cm}^{-2} \rm{s}^{-1}\\ 3.06963\times 10^{7} ~\text{g cm}^{-2} \rm{s}^{-1}\\ 2.08849\times 10^{18} ~\text{erg }\text{g cm}^{-3}\\ 12270.2 ~\rm{G}\\ 4338.16~\rm{G} \end{pmatrix} \end{split} \end{equation}
Using these values and a perturbation amplitude of A = ρ 0 / 10 6 A = ρ 0 / 10 6 A=rho_(0)//10^(6)A = \rho_0/10^6, I perturb a Fast Magneto-sonic wave in the two cases presented in Table 1 according to Equation [34]. The simulation is initialized in 1 Dimension with a domain size of one and periodic boundary conditions. This allows the left-going wave to propagate through the boundary and passing through the domain multiple times. The wave is evolved allowing it to pass through the domain 10 times. This is repeated for the two compositions. I measure the phase velocity by finding the x x xx-position of the density maximum at t = 0 t = 0 t=0t=0 s, then measuring the time the density within that cell returns to the maximal value. The phase speed is then given by:
(40) v F M S = x-domain~size time~for~wave~to~pass (40) v F M S = x-domain~size time~for~wave~to~pass {:(40)v_(FMS)=("x-domain~size")/("time~for~wave~to~pass"):}\begin{equation} v_{FMS} = \frac{\text{x-domain~size}}{\text{time~for~wave~to~pass}} \end{equation}
image Figure 3 Fast magneto-sonic wave shown in terms of density as a function of x position. Black represents the example where pair production is more important while the low electron fraction case (purple) represents where pair production is less important. From top to bottom we see the evolution of the waves in time, being initialized exactly the same until the low electron fraction case eventually leads the high electron fraction case (left-going wave) owing to the difference in phase speeds between the two scenarios.
As you can see in Figure 3 the fast magneto-sonic wave speed varies by Δ v F M S 7 × 10 6 Δ v F M S 7 × 10 6 Deltav_(FMS)∼7xx10^(6)\Delta v_{FMS} \sim 7\times 10^{6}~cm/s between two wave that are initialized given the same density and temperature while varying the electron fraction by a factor two. So despite the fact that we are in the fluid limit, and are ignoring the motion of individual particles, pair production modifies the composition of the material which in turn has implications on magneto-sonic waves through the sound speed component.
Keeping the density and temperature constant, while modifying the composition (electron fraction) by a factor two, I find that the pressure varies by 5 % 5 % ∼5%\sim 5\%, the adiabatic constant γ γ gamma\gamma varies by 10 % 10 % ∼10%\sim10\%, partially counteracting, and the FMS phase speed varies by 0.5 % 0.5 % ∼0.5%\sim 0.5\%.

5. Conclusion

I presented a discussion of the importance of e e + e e + e^(-)e^(+)-e^-e^+-pair production in plasmas, particularly focusing on how pair production is accounted for in the MHD code FLASH. To begin I discussed the limit in which pair production becomes important. The often quoted limit for pair production is an energy limit where a photon with energy in excess of the mass energy of the an electron positron pair, then this photon can produce a pair with any excess energy going into the kinetic energy of the produced particles. Using this as an energy limit and comparing the energy to k B T k B T k_(B)Tk_BT corresponds to a temperature of T 10 10 T 10 10 T∼10^(10)T\sim 10^{10}~K.
However in a plasma, pair production will not act as a switch being on or off, for example if a pair is produced by high energy photons within the plasma, photons can be produced by many different mechanisms, for example black-body radiation. With a roughly uniform temperature plasma, black-body photons will have a distribution of energies, some with the requisite energy to produce the pair and others without. So it is more useful to discuss in what regime does pair production become important. By numerically integrating the Fermi-Dirac distribution in combination with the density of states, and using charge neutrality I was able to show how electron and positron number densities vary with temperature in this extremely high energy regime. Electron and positron number densities converge at a temperature of T 3 × 10 8 T 3 × 10 8 T∼3xx10^(8)T\sim 3\times 10^{8}~K, suggesting that pair production becomes important at temperatures almost 2 orders of magnitude less than the energy limit.
As spatial scales of the fluid limit neglect the importance of individual particles how does pair production effect an MHD fluid in the fluid limit? Pair-production predominantly is important in terms of its ability to change the composition of the fluid. Increasing pair production increases the number of electrons relative to baryons which in other words increases the electron fraction of the fluid. FLASH has the option of using the Helmholtz equation of state, a composition dependent EOS that returns the pressure given the internal energy, density and electron fraction. By changing the composition of the fluid, how the pressure and density relate is modified as well. The Helmholtz EOS tracks the number density of electrons and positrons by integrating the Fermi-Dirac distribution, each contributes a partial pressure changing the total pressure of the fluid.
I then presented a test case in which I could demonstrate directly that the inclusion of pair production within the fluid has physical effects on macroscopic parameters visible in the fluid limit. This test case was perturbing a Fast magneto-sonic wave in one dimension within a plasma at densities and temperatures consistent with pair production (Table 1) and varying the composition (electron fraction) of the plasma to demonstrate the change in wave speed as a result of the composition change. I find that in the high density and temperature regime by changing the electron fraction by a factor of two, the pressure of the fluid varies by 5 % 5 % ∼5%\sim 5\%, the adiabatic constant γ γ gamma\gamma varies by 10 % 10 % ∼10%\sim 10\%, partially counteracting one another resulting in a fast magneto-sonic wave speed that varies by 0.5 % 0.5 % ∼0.5%\sim 0.5\%. The accumulation of phase shift in the density profile is visible in Figure 3. Unfortunately by only setting the density and temperature (varying electron fraction), this test results in a change in both P P PP and γ γ gamma\gamma, so the change in the FMS speed is less apparent then if we could fix one of these quantities. This however would not ensure both cases fell within the pair production regime by re-scaling the pressure.
In summary pair production is relevant physical effect in high energy astrophysical plasmas, in the fluid limit its main contribution comes in the way of modifying the composition of the plasma. We can demonstrate this effect in the fluid limit through fast magneto-sonic waves whose phase speed depends on the relation between pressure and density.

6. References

Bennett, C. L., Halpern, M., Hinshaw, G., et al. 2003, The Astrophysical Journal Supplement Series, 148, 1, doi: 10.1086/377253
Breit, G., & Wheeler, J. A. 1934, Physical Review, 46, 1087, doi: 10.1103/PhysRev.46.1087
Dubey, A., Antypas, K., Ganapathy, M. K., et al. 2009, J. Par. Comp., 35, 512 , doi: DOI:10.1016/j.parco.2009.08.001
Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273, doi: 10.1086/317361
Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509, doi: https://doi.org/10.1016/j.jcp.2004.11.016
Landau, L., & Lifshitz, E. 1934, in Perspectives in Theoretical Physics, ed. L. PITAEVSKI (Amsterdam: Pergamon), 27–38, doi: https://doi.org/10.1016/B978-0-08-036364-6.50006-5
MacFadyen, A. I. 2003, in From Twilight to Highlight: The Physics of Supernovae, ed. W. Hillebrandt & B. Leibundgut, 97, doi: 10.1007/10828549 14
Pathria, R., & Beale, P. D. 2011, in Statistical Mechanics (Third Edition), third edition edn., ed. R. Pathria & P. D. Beale (Boston: Academic Press), 275–297, doi: https://doi.org/10.1016/B978-0-12-382188-1.00009-8
Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137, doi: 10.1086/588755
Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, The Astrophysical Journal Supplement Series, 249, 4, doi: 10.3847/1538-4365/ab929b
Thorne, K. S., & Blandford, R. D. 2017, Modern Classical Physics: Optics, Fluids, Plasmas, Elasticity, Relativity, and Statistical Physics (Princeton University Press)
Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
Woosley, S. E. 1993, ApJ, 405, 273, doi: 10.1086/172359

Recommended for you

Kaichao You
A New Paradigm for Exploiting Pre-trained Model Hubs
A New Paradigm for Exploiting Pre-trained Model Hubs
It is often overlooked that the number of models in pre-trained model hubs is scaling up very fast. As a result, pre-trained model hubs are popular but under-exploited. Here a new paradigm is advocated to sufficiently exploit pre-trained model hubs.
8 points
0 issues