The Biermann Catastrophe in Numerical MHD
Carlo Graziani, Petros Tzeferacos, Dongwook Lee, Donald Q. Lamb, Klaus Weide, Milad Fatenejad, Joshua Miller
SS UBMITTED TO THE A STROPHYSICAL J OURNAL
Preprint typeset using L A TEX style emulateapj v. 5/2/11
THE BIERMANN CATASTROPHE IN NUMERICAL MHD C ARLO G RAZIANI , P
ETROS T ZEFERACOS , D
ONGWOOK L EE , D ONALD
Q. L
AMB , K
LAUS W EIDE , M
ILAD F ATENEJAD , & J
OSHUA M ILLER
Flash Center for Computational Science, Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL, 60637
Submitted to the Astrophysical Journal
ABSTRACTThe Biermann Battery effect is frequently invoked in cosmic magnetogenesis and studied in High-EnergyDensity laboratory physics experiments. Generation of magnetic fields by the Biermann effect due to mis-aligneddensity and temperature gradients in smooth flow behind shocks is well known. We show that a Biermann-effectmagnetic field is also generated within shocks. Direct implementation of the Biermann effect in MHD codesdoes not capture this physical process, and worse, produces unphysical magnetic fields at shocks whose valuedoes not converge with resolution. We show that this convergence breakdown is due to naive discretization,which fails to account for the fact that discretized irrotational vector fields have spurious solenoidal componentsthat grow without bound near a discontinuity. We show that careful consideration of the kinetics of ion viscousshocks leads to a formulation of the Biermann effect that gives rise to a convergent algorithm. We note twonovel physical effects: a resistive magnetic precursor in which Biermann-generated field in the shock “leaks”resistively upstream; and a thermal magnetic precursor , in which field is generated by the Biermann effectahead of the shock front due to gradients created by the shock’s electron thermal conduction precursor. Botheffects appear to be potentially observable in experiments at laser facilities. We re-examine published studiesof magnetogenesis in galaxy cluster formation, and conclude that the simulations in question had inadequateresolution to reliably estimate the field generation rate. Corrected estimates suggest primordial field values inthe range B ∼ − G — − G by z = 3 . Subject headings: magnetohydrodynamics — plasmas — magnetic fields INTRODUCTIONDynamo theories of proto- and extra-galactic primordialmagnetic fields, which endeavor to explain how those fieldsachieved their current strength and structure, work by ampli-fying small initial seed fields by means of turbulent plasmamotions (Kronberg 1994; Kulsrud & Zweibel 2008). However,the induction equation of resistive magnetohydrodynamics(MHD), ∂ B ∂t = ∇ ××× (cid:26) u ××× B − ηc π ∇ ××× B (cid:27) , (1)always admits the solution B ( x , t ) = 0 . This simple observa-tion poses a problem for the generation of the required seedfields, as they cannot be created in ideal MHD starting from afield-free state.There have been several proposals for generating the re-quired seed fields from mechanisms such as primordial phasetransitions, or from processes occurring during inflation (seereviews in Widrow 2002; Widrow et al. 2012). The Biermannbattery effect (Biermann 1950) provides another popular res-olution of this problem (Kulsrud et al. 1997; Widrow 2002;Kulsrud & Zweibel 2008). The effect, which arises in conse-quence of the large difference in the electron and ion mass,is attributable to small-scale charge separation in the plasma.Pressure forces produce much larger accelerations of electronsthan of ions, and the relative acceleration of the two compo-nents results in charge separation that must be balanced by anelectric field E B ≡ − ( en e ) − ∇ P e , (2)where n e and P e are the electron number density and pressure,respectively. Since this field is not, in general, irrotational, it [email protected] can act as a source of magnetic field in the induction equation, ∂ B ∂t = ∇ ××× (cid:26) u ××× B − ηc π ∇ ××× B + cen e ∇ P e (cid:27) , (3)generating a non-zero B from an initially unmagnetized state.The Biermann battery effect has been successfully invokedin numerical simulations exploring the generation of seed fieldsin cosmological ionization fronts (Subramanian et al. 1994;Gnedin et al. 2000), protogalaxies (Davies & Widrow 2000),and Pop-III star formation (Xu et al. 2008). Magnetic fieldgeneration by the Biermann mechanism is also of significantinterest in direct-drive and indirect-drive inertial confinementfusion, where strong gradients behind the converging shockcan lead to dynamically important field strengths (Srinivasan& Tang 2013), and more generally in the field of High-EnergyDensity Physics, where the effects of field generation (Gregoriet al. 2012) and amplification (Meinecke et al. 2014) can beexamined in a laboratory setting at laser facilities, in exper-iments where large gradients are produced in strong plasmashocks (Fryxell et al. 2010; Tzeferacos et al. 2012, 2014).While the generation of magnetic fields by the Biermanneffect as a result of strong mis-aligned density and temperaturegradients in the smooth flow behind shocks is well known, weshow that there exists a previously unrecognized Biermanneffect due to the electron-ion charge separation that occurs within ion viscous shocks. This effect is a consequence ofthe kinetic theory of shock structure in plasmas, as we showbelow.A straightforward implementation of the Biermann effect infinite-volume Eulerian, purely Langrangian, and ALE codes,whether as a split source term or as a flux term, does not cap-ture this physical process, and worse, leads to non-convergentresults (Fatenejad et al. 2013). In symmetric situations suchas planar or spherical shocks, where no field should arise, a r X i v : . [ a s t r o - ph . GA ] J a n GRAZIANI ET AL.such implementations produce anomalous field generation thatgrows without bound with resolution (Fatenejad et al. 2013).This behavior is observed across a range of different MHDcodes (see the discussion of codes in Fatenejad et al. 2013).This is the Biermann catastrophe of numerical MHD. We showthat this failure is intimately related to the failure of such codesto correctly model the structure of the plasma shock.In a gasdynamic/MHD formulation, where the shocks aremodeled as zero-width discontinuities of the flow, the troublearises from the behavior of the Biermann flux, which is to sayfrom the electric field, Eq. (2), in the vicinity of a shock. Thegradient ∇ P e , which analytically speaking acquires a Dirac δ component at the shock, is ascribed a numerical magnitude thatgrows without bound at the shock with increasing resolution. Itis this divergence that is connected with failure of MHD codesto correctly predict shock-driven magnetic field generation insupposedly simple test cases, where the correct value of thegenerated field is zero. The investigation of this failure is oneof the central concerns of this article.That this failure was not previously recognized is a conse-quence of the history of the numerical study of the BiermannBattery effect as a source of cosmic seed fields, which has acurious feature with respect to shocks. On the one hand, theimportance of shocks to lifting barotropic constraints, thusmaking available the kind of non-aligned gradients required todrive the Biermann effect, has been widely recognized (Kul-srud et al. 1997; Davies & Widrow 2000). On the other hand,the direct effect of shocks – as opposed to that of their trail-ing downstream gradients – on this magnetogenesis has notreally been carefully examined. This is probably because thecomputational strategies that have been adopted both circum-vent difficulties at shocks and direct attention away from themagnetizing properties of shocks. For example, Kulsrud et al.(1997) and Gnedin et al. (2000) perform essentially hydrody-namic simulations, in which magnetic fields evolved by theinduction equation have no dynamical role. Davies & Widrow(2000) forgo the induction equation altogether, in favor of theequation of inviscous evolution of vorticity ω ∂ ω ∂t = ∇ ××× (cid:26) u ××× ω + 1 ρ ∇ P (cid:27) , (4)which, by comparison with Eq. (3), implies that inviscous,non-resistive plasmas satisfy the equation ∂ B ∂t = cM i e (1 + ¯ Z ) ∂ ω ∂t . (5)(Kulsrud & Zweibel 2008). On its face, this relation wouldappear to suggest that in an initially unmagnetized and irrota-tional plasma, the magnetic field is simply proportional to thevorticity, B = cM i e (1+ ¯ Z ) ω , where M i is the ion mass and ¯ Z isthe average ionization fraction.One difficulty with these approaches is that they entirelyneglect to treat the field generation within the shock itself , aswe noted above. There are indeed large, non-aligned gradi-ents at the simulated shocks, but these gradients are unphys-ical side-effects of the numerical strategies used to integratethe hydrodynamic equations, and consequently are resolution-dependent. It is therefore hopeless to expect convergence withresolution of the resulting magnetic fields.Furthermore, any magnetization generated by the shock it-self imprints itself on the magnetic field structure as the shockmoves through, leaving behind a substantial residue that is superposed on the smooth-flow Biermann-generated field. Itis essential to come to a correct understanding of the behaviorof the Biermann term at shocks, to have any confidence thatresults arrived at in this manner bear any resemblance to real-ity. In particular, the presumption that shock jump conditionon vorticity should bear a relation to the jump condition onmagnetization that preserves the proportionality between thetwo has been hypothesized but never demonstrated (Kulsrudet al. 1997; Davies & Widrow 2000), and seems in fact far-fetched: it would imply that magnetization is vorticity, that is,that the coincidence expressed by Eq. (5) in fact reflects a deepidentity, and that magnetic degrees of freedom of ideal MHDare somehow already contained in unmagnetized Eulerian hy-drodynamics, encoded in derivatives of the velocity field. Sucha claim is difficult to accept, and to even attempt to support itwould require an analysis of the modification brought by theBiermann Battery term to the Rankine-Hugoniot jump condi-tion on B , an analysis that we furnish for the first time in thispaper. We will return to a discussion of vorticity in § B ?3. Assuming the Biermann source term, Eq. (2) can beinterpreted in a sensible manner in the vicinity of aninviscid shock, are its predictions consistent with thepredictions of kinetic theory near a shock? Or doesthe Biermann term in MHD need to be flux-limited, inthe style of thermal and radiation source terms whosemisbehavior must be limited on kinetic theory groundswhen gradients grow too large?We address question (1) in §
2, abstracting the essential in-gredients of plasma shock theory required to construct a validMHD model of the Biermann effect due to shocks; whereas weHE BIERMANN CATASTROPHE 3address the remaining two questions in §
3, where we demon-strate that (2) in fact the Biermann source of Eq. (2) is math-ematically consistent and well-behaved near weak-solutiondiscontinuities, and (3) the Biermann source term in fact yieldsthe correct EMF across the shock, matching the prediction ofplasma shock theory (Zel’dovich, Y.B. and Raizer, Y.P. 2002;Amendt et al. 2009; Jaffrin & Probstein 1964).We clarify the origin of the Biermann catastrophe as a nu-merical effect attributable to the difficulty of discretizing thesource term of Eq. (2) in the vicinity of a shock. We showthat the numerical anomaly can be eliminated by leveragingthe continuity of the electron temperature T e across shocks – abenefit of the kinetic-theory connection. Reformulation of theBiermann source term in terms of T e allows the singularity tobe isolated, and the flux of magnetic field due to the Biermanneffect to be rewritten in a manifestly finite form suitable fortranslation into a convergent numerical algorithm.The Biermann effect is due to electron-ion charge separa-tion, and is sensitive to departure from thermal equilibriumbetween electrons and ions. Such a departure is precisely whatoccurs at shocks, so that a correct treatment of the effect atshocks necessarily requires that the disequilibrium be mod-eled. For this, a 2-temperature plasma model is mandatory.An interesting consequence of this observation is the fact thatthe Biermann effect is enhanced not only at shocks, but alsoat contact discontinuities, at ionization fronts, and at materialspecies boundaries, because the electron partial pressure isdiscontinuous at such surfaces, even though the total pressureis continuous there. This enhancement cannot be computed inan equilibrium treatment of the plasma.An additional point of this paper is to point out two noveland interesting effects associated with the Biermann effect inthe neighborhood of a shock. In § resistivemagnetic precursor is generated in resistive MHD, whereinmagnetic field generated by the Biermann effect in the shock“leaks” resistively from the shock into a region of the upstreamfluid whose physical extent is proportional to the resistivity.And in § thermal magnetic precursor isgenerated ahead of the shock through the Biermann effectby plasma motions generated by the shock’s electron thermalconduction precursor. Both effects are potentially observablein laboratory conditions at high-intensity laser facilities suchas Vulcan, Omega, and NIF. We discuss the observability ofthese effects at laser facilities. We show that appropriately-designed experiments at such facilities could currently observethe resistive magnetic precursor, providing a clean experimen-tal validation test of the Biermann effect in plasma shocks. Wealso argue that the smaller thermal magnetic precursor mightbecome observable in future experiments. REVIEW OF KINETIC THEORY OF PLASMA SHOCKSWe begin by reviewing some essential results from the ki-netic theory of shocks in plasmas. The basic theory of the fluidstructure of planar shocks in plasmas was set out in Shafranov(1957), while the electromagnetic structure of such shocks wasdiscussed in Jaffrin & Probstein (1964), and more recently inAmendt et al. (2009). An extremely lucid presentation of theseresults may be found in Chapter VII of Zel’dovich, Y.B. andRaizer, Y.P. (2002).There are three essential ingredients to be imported from thekinetic theory of plasma shocks in order to fashion a workingMHD model of the Biermann effect: the loss of thermal equi-librium between electrons and ions at shocks; the adiabaticbehavior of electrons, up to electron heat conduction; and the charge-separation-induced electric field across the shock front.We now review these in order, and conclude by presenting theMHD model that forms the basis for the rest of this paper.2.1.
Ion-Electron Disequilibrium
As discussed on p.36 of Drake (2006), a strong shock dis-turbs the thermal equilibrium between electrons and ions in aplasma. That equilibrium is maintained by electron-ion colli-sions, and operates over timescales τ ei that are long comparedto the shock-crossing time of a parcel of fluid entering theion viscous shock in consequence of the large ratio m i /m e (Spitzer 1962). As a result, it is essential to describe the fluidin terms of an additional degree of freedom – the electrontemperature, T e – with respect to the usual equilibrium MHDmodel.Fortunately, it is unnecessary to model the fluid using newinertial degrees of freedom to describe the electron fluid. Asingle inertial component for the fluid as a whole gives anadequate description of the fluid structure near the shock, and acompletely satisfactory description of the fluid in smooth flowregions, where electron-ion collisions restore local thermalequilibrium (Drake 2006; Zel’dovich, Y.B. and Raizer, Y.P.2002). This makes it easier to adapt existing MHD codes totreat the Biermann effect correctly, since it is much easier toadd a scalar degree of freedom for T e than it would be to dealwith two velocity fields, one each for the electron fluid and forthe ion fluid. 2.2. Nearly Adiabatic Electrons
The large ion-to-electron mass ratio, which we recall fromthe discussion preceding Eq. (2) is responsible for the chargeseparation that produces the Biermann effect, has the furtherconsequence that electron mobility is much higher than ionmobility, and that thermal conductivity due to electron-electroncollisions dominates the heat transport in the fluid. At the sametime, the forces between electrons and ions as the fluid crossesthe ion viscous shock front are effectively dissipation-free,as they are simply electrostatic fields generated by chargeseparation, and collisional dissipation processes are too slowto act during the shock-crossing timescale.These observations lead to the conclusion that while theions undergoing shock compression experience the usualthermodynamically-irreversible, entropy-generating processassociated with shocks, the electrons are compressed adia-batically by the electrostatic forces exerted upon them bythe ions, and consequently do not suffer entropy incrementsdue to shock compression. The electron entropy would be apassively-advected scalar, then, except for the dissipative effectof electron thermal conduction, and for the slow (relative totimescales relevant to the shock) effect of electron-ion colli-sional equilibration. Electron thermal conduction effectivelyrules out any sudden change in T e at the shock, since such achange would produce an enormous restoring heat flux to healthe discontinuity.The fluid structure that follows from these considerations isof a sudden discontinuous compression of the ions at the shock,accompanied by a smooth increase in T e , which is continuousthroughout the shock (in the Eulerian limit where the shockwidth tends to zero, T e acquires a discontinuous derivative inthe direction of the shock normal). The electron temperatureexhibits a thermal precursor that leads the shock. The size λ T of the precursor region may be calculated by balancingadvection against heat diffusion in the frame of the shock, and GRAZIANI ET AL.is about λ T ≈ κ e ρc v,e D , (6)where κ e is the electron thermal conductivity, c v,e is the elec-tron specific heat at constant volume per unit mass, and D isthe shock speed (Shafranov 1957; Zel’dovich, Y.B. and Raizer,Y.P. 2002). The effect of this precursor region on accretionshocks in galaxy clusters has been recently studied in Smithet al. (2013). 2.3. Electric Structure
The sudden change in the motion of ions entering the shocksheath, together with the more highly mobile motion of theelectrons in the vicinity of the shock, results in a zone ofcharge separation-driven electric field in the normal directionto the shock (Jaffrin & Probstein 1964; Amendt et al. 2009;Zel’dovich, Y.B. and Raizer, Y.P. 2002). By solving the Navier-Stokes equation across the shock, (treating the results as meansof kinetic-theory distributions, since a fluid picture is clearlyinvalid inside the viscous shock sheath), Jaffrin & Probstein(1964) calculated the full electric structure across the shock.We do not require the full structure here, since in the end wewish to adopt an Eulerian picture of the shock, in which thewidth of the shock is zero. For our purposes, the EMF acrossthe shock is all that is required. It is shown in Amendt et al.(2009) that the EMF E is given by E = k B T e ln (cid:18) n e,d n e,u (cid:19) = k B T e ln (cid:18) ρ d ρ u (cid:19) , (7)where the subscripts “ u ” and “ d ” denote “upstream” and“downstream”, respectively, and where we’ve assumed com-plete ionization to obtain the second line.This electric field is due to charge separation, and thus has acommon ancestry with the electric field E B in Eq. (2). How-ever, it is proper to the shock, rather than to the smooth-flowportion of the fluid. It is essential, for the purpose of physicalconsistency, to demonstrate that an MHD model implementingthe Biermann effect should reproduce the shock-crossing EMFgiven in Eq. (7). That this is possible is demonstrated in thenext section. 2.4. The MHD Model
We wrap up this section by exhibiting the Biermann-effect-laden MHD model implied by these kinetic-theory considera-tions. The dynamical system describing the model is given in conservation form as follows: ∂ρ∂t + ∇ · ( ρ u ) = 0 (8) ∂ρ u ∂t + ∇ · (cid:20) ρ uu − π BB + (cid:18) P + B π (cid:19)(cid:21) = 0 (9) ∂∂t (cid:20) ρ (cid:18) u + (cid:15) T (cid:19) + B π (cid:21) + ∇ · (cid:20) ρ u (cid:18) u + (cid:15) T + Pρ (cid:19) + c π (cid:16) − ( u /c ) ××× B + cη π ∇ ××× B + E B (cid:17) ××× B (cid:21) = 0 (10) ∂ρs e ∂t + ∇ · ( ρ u s e ) = − T e − ∇ · ( − κ e ∇ T e )+ ρc v,e T e τ ei ( T i − T e ) (11) ∂ B ∂t + ∇ · (cid:20) uB − Bu − c η π ∇ B (cid:21) + c ∇ ××× E B = 0 . (12)In these equations, P = P e + P i is the total material pressure, (cid:15) T is the total specific thermal energy per unit mass, and s e isthe specific electron entropy per unit mass.Eqs. (8-9) are the usual MHD equations of mass and momen-tum conservation. The energy conservation equation, Eq. (10)includes a term to account for resistive dissipation, and anotherto account for the Biermann effect.Eq. (11) expresses the conservation of electron entropy s e (up to heat conduction and electron-ion heat exchange), which,as we pointed out above, is required by the same approxima-tion m e /m i → that gives rise to the Biermann effect in thefirst place. This equation is introduced to describe the addi-tional degree of freedom required by the 2-temperature plasmatreatment.Finally, in the induction equation, Eq. (12), the Biermanneffect is expressed here as a curl, rather than in the conservationform – a divergence of a flux tensor – that it will assume later.These dynamical equations are supplemented here byperfect-gas equations-of-state for both the electrons and theions. We assume total ionization throughout what follows forthe sake of simplicity. THE BIERMANN EFFECT AT SHOCKS3.1.
Kinematics of the Biermann Effect at Shocks
By inspection of the Biermann source term in Eq. (3) wesee that it is proportional to ∇ n e ××× ∇ P e . It follows thatfield generation by the Biermann effect can only occur if thegradients of n e and P e are not aligned. This means that inshocks with planar, cylindrical, or spherical symmetry, thefield generation rate is zero. We must therefore treat non-symmetric shock surfaces in order to analyze non-trivial casesof field generation. Accordingly, in this subsection, we set upthe basic kinematics of the Biermann effect at general shocksurfaces.To describe the shock surface, we introduce a level function Ψ( x , t ) , and use it to define the shock surface Ψ( x , t ) = 0 .Note that the function Ψ is of no dynamical significance, butis rather simply a mathematical convenience for describing theshock surface. We will assume that the shock is moving inthe direction of the normal vector n ≡ ∇ Ψ / |∇ Ψ | , so that theregion Ψ( x , t ) > is upstream, whereas the region Ψ( x , t ) < is downstream. These definitions are illustrated in Figure 1.We denote the local shock speed along n by D . By con-sidering the motion of the level surface Ψ( x , t ) = 0 it is notHE BIERMANN CATASTROPHE 5 Figure 1.
Illustration of the kinematics of the Biermann effect at a shocksurface. difficult to show that D = − ∂ Ψ ∂t / |∇ Ψ | . (13)To describe a moving MHD discontinuity that coincideswith the surface Ψ( x , t ) = 0 , we will frequently express afield quantity X by the decomposition X = X u ( x , t )Θ(Ψ) + X d ( x , t )Θ( − Ψ) , where X u and X d are continuous functions,and where Θ(Ψ) is a Heaviside step function. After substitut-ing such expressions into field equations, we will find someterms proportional to the Dirac distribution δ (Ψ) , resultingfrom differentiation of the Heaviside functions. We will re-fer to the collected terms of this form as the “shock part” ofthe evolution equations. As will be seen, such terms embodyRankine-Hugoniot jumps, which can thus be efficiently ex-tracted for these non-symmetric shocks. This trick is alsouseful for deducing hydrodynamic flux terms, as will also beevident below.It is convenient to reformulate the Biermann source termusing the ideal equation of state to replace n e with T e , theelectron temperature. Assuming an ideal gas equation of state,the reformulated electric field is E B ≡ − ( k B /e ) T e ∇ ln P e , (14)so that the source term due to the Biermann effect in the induc-tion equation is ∂ B ∂t (cid:12)(cid:12)(cid:12)(cid:12) B = ( ck B /e ) ∇ T e ××× ∇ ln P e . (15)An important reason for this reformulation is that as we sawin §
2, in the presence of electron thermal conduction, T e iscontinuous at the shock , whereas n e is not (Shafranov 1957;Zel’dovich, Y.B. and Raizer, Y.P. 2002)). As we saw, thecontinuity of T e is a consequence of the high mobility ofelectrons relative to ions, and is therefore part and parcel ofthe same approximation that led to the Biermann source term(Eq. 2) in the first place. It is central to the developments thatfollow.The log electron pressure ln P e is discontinuous at the shock,and may be represented at time t by ln P e = ln P e,u Θ (Ψ) + ln P e,d Θ ( − Ψ) , (16)where P e,u ( x , t ) and P e,d ( x , t ) are continuous functions. The discontinuity of P e leads to a Dirac- δ singularity in ∇ P e . Using the distributional relation d Θ(Ψ) /d Ψ = δ (Ψ) ,we have ∇ ln P e = ln (cid:18) P e,u P e,d (cid:19) δ (Ψ) ∇ Ψ+Θ (Ψ) ∇ ln P e,u + Θ ( − Ψ) ∇ ln P e,d . (17)We may easily interpret the term in Eq. (17) that is proportionalto δ (Ψ) as the gradient proper to the shock, and the remainingterms as the gradient in smooth flow. Since we are interestedin the shock behavior, we define ∇ ln P e | Shock ≡ ln (cid:18) P e,u P e,d (cid:19) δ (Ψ) ∇ Ψ= ln (cid:18) ρ u ρ d (cid:19) δ (Ψ) ∇ Ψ , (18)where we’ve used the continuity of T e and the assumption ofcomplete ionization to replace the pressure ratio across theshock with the corresponding compression ratio ρ u /ρ d .The electron temperature T e is continuous across the shock,but its normal derivative is discontinuous. We may thereforerepresent T e by T e = T e,u Θ (Ψ) + T e,d Θ ( − Ψ) , (19)where T e,u ( x , t ) and T e,d ( x , t ) are continuous functions, andwhere continuity at the shock implies Ψ( x , t ) = 0 ⇒ T e,u ( x , t ) = T e,d ( x , t ) . (20)In virtue of the continuity of T e , the gradient ∇ T e is not burdened by a Dirac- δ , and we obtain ∇ T e = Θ (Ψ) ∇ T e,u + Θ ( − Ψ) ∇ T e,d . (21)Inserting Eqs. (18) and (21) into Eq. (15), we obtain ∂ B ∂t (cid:12)(cid:12)(cid:12)(cid:12) B,Shock = ( ck B /e ) ln (cid:18) P e,u P e,d (cid:19) δ (Ψ( x )) × (cid:20)
12 [ ∇ T e,u + ∇ T e,d ] ××× ∇ Ψ (cid:21) , (22)where we’ve made use of the distributional relation Θ(Ψ) δ (Ψ) = Θ(Ψ) dd Ψ Θ(Ψ)= 12 dd Ψ Θ(Ψ) = 12 dd Ψ Θ(Ψ)= 12 δ (Ψ) , (23)where to get the third line, we’ve used the fact that the Heavi-side function squares to itself, since its values are either 0 or 1.This derivation of Θ δ = δ/ is admittedly cavalier in its treat-ment of distributional quantities, and is presented for brevity.The result may also be obtained by a limiting procedure, inwhich the Heaviside function is represented by the limit of afamily of continuous functions – this is, after all, the “Eulerian”limit of the description of a shock as the viscosity goes to zero.Introducing the shock normal vector n , and using the factthat the tangential derivative of T e is continuous (so that n ××× GRAZIANI ET AL. ∇ T e,u = n ××× ∇ T e,d at the shock), we obtain ∂ B ∂t (cid:12)(cid:12)(cid:12)(cid:12) B,Shock = δ (Ψ) |∇ Ψ | × ( ck B /e ) ln (cid:18) P e,u P e,d (cid:19) ∇ T e ××× n . (24)Eq. (24) is amenable to a clarifying interpretation, which fol-lows from the recognition that if V is a region containing someportion Σ of the shock surface, and the differential element ofsurface area on Σ is dA , then (cid:82) V d x δ (Ψ( x )) |∇ Ψ | = (cid:82) Σ dA ,which is to say that δ (Ψ( x )) |∇ Ψ | d x is the differential areaelement on the shock surface. We may therefore interpret thequantity ( ck B /e ) ln (cid:16) P e,u P e,d (cid:17) ∇ T e ××× n in Eq. (24) as a field gen-eration rate per unit time per unit area of the shock surface.This interpretation already allows us to see that the Biermannsource term should give rise to finite, mathematically sensiblefield generation.To verify this, we may now obtain the field generation ratedue to the passage of the shock by inserting the above expres-sions into the induction equation, Eq. (3), neglecting resistivity.The equation may be written as ∂ B ∂t = ∇ · ( Bu − uB ) + ∂ B ∂t (cid:12)(cid:12)(cid:12)(cid:12) B . (25)In the neighborhood of the shock, we define B = B n n + B T u
Θ(Ψ) + B T d Θ( − Ψ) (26) u = u u Θ(Ψ) + u d Θ( − Ψ) (27) u u ≡ u T u + u nu n (28) u d ≡ u T d + u nd n , (29)where the “tangential” components satisfy B T u · n = 0 , B T d · n = 0 , u T u · n = 0 , u T d · n = 0 . In virtue of thesolenoidal condition ∇ · B = 0 , B n is continuous across theshock.We transform the advection term ∇ · ( Bu − uB ) in Eq. (25)by inserting Eqs. (26-27), grouping the Heaviside functions,and collecting the “shock” terms proportional to the Dirac- δ that arises from differentiating the Heaviside functions. UsingEq. (13) to replace ∂ Ψ /∂t , and using Eq. (24), we obtain δ (Ψ) |∇ Ψ | ( − D )( B T u − B T d ) = δ (Ψ) |∇ Ψ | × (cid:110) B n [ u T u − u T d ] (30) − [ u nu B T u − u nd B T d ]+ ck B e ln (cid:18) P e,u P e,d (cid:19) ∇ T e ××× n (cid:111) , (31)whence the Rankine-Hugoniot jump condition, generalized bythe Biermann flux is [( D − u n ) B T ] ud + B n [ u T ] ud + ck B e ∇ T e ××× n [ln P e ] ud = 0 . (32)Here, the notation [ . . . ] ud means the difference of the up-stream and downstream values at the shock location. The firsttwo terms in Eq. (32) comprise the usual jump condition forthe induction equation (see Chapter 7 of Gurnett & Bhattachar-jee 2005, for example). The final term is the contribution tothe jump condition from the Biermann effect, which is seen to produce a finite, well-defined discontinuity at the shock.We may obtain a useful result by considering the jump in B due to a shock advancing into a quiescent, unmagnetized fluid( B n = 0 , B T u = 0 , u u = 0 ). The downstream magnetic fieldis then given by B T d = ck B e ( D − u nd ) − ln (cid:18) P e,u P e,d (cid:19) ∇ T e ××× n = ck B eD ρ d ρ u ln (cid:18) ρ u ρ d (cid:19) ∇ T e ××× n , (33)where we’ve used the condition of mass continuity at the shock, ρ u D = ρ d ( D − u nd ) , to replace the term ( D − u nd ) , and wherewe’ve also used the continuity of T e to replace the pressureratio in the log with the density ratio.We conclude from the above development that the Bier-mann source term is mathematically well-defined even at weaksolution discontinuities, and yields definite finite predictionsto which a properly designed numerical algorithm should beexpected to converge.3.2. Physical Compatibility of the Biermann Source TermWith Plasma Shock Theory
In the Introduction, we raised the question of whether theBiermann source term, Eq. (2) behaves near a shock accordingto the predictions of kinetic theory, as summarized in §
2. Wenow explain the line of thinking behind this question, andanswer it definitively.The induction equation, Eq. (3) can be cast in conservativeform. In this form, the Biermann source term assumes theform of the divergence of a flux tensor, whose components arelinear in ∇ P e . It is clear that the construction of these fluxcomponents from the gradient of a discontinuous function isin some way associated with the numerical troubles that arisewith the Biermann effect near shocks.The key point here is that “trouble with a flux computedfrom a derivative” is actually a familiar situation from radia-tion diffusion, where the radiation flux, which is proportionalto the gradient of the radiation temperature, yields unphys-ical (superluminal) fluxes at regions where the temperaturechanges sharply. (see pp. 478-481 of Mihalas & Weibel-Mihalas 1999, for example). Another analogous situationarises with respect to electron thermal conduction, where thethermal flux from ∇ T e may yield unphysically-large transportat shocks, in consequence of the discontinuous change in T e (Mihalas & Weibel-Mihalas 1999, p. 302). In both of thesecases, a straightforward work-around is furnished by a flux lim-iter , in which a maximum flux deduced from a kinetic-theorypicture of the transport is used to smoothly cut off the flux inregions where gradients get large, while leaving the fluxes insmooth regions unaltered.This is the reason that the question of the validity of theBiermann flux at shocks is a natural one to ask. If it were thecase that the flux is just wrong when |∇ P e | gets too large, areasonable approach would be to use the limiting value forthe flux from kinetic theory (Eq. 7), and impose that flux inthe shock as a cutoff through some kind of flux limiter. Ifthe electric field due to charge separation given by Eq. (2)results in an EMF across a shock that exceeds this value, wecould cut it off at this value, and obtain approximate resultsanalogous to flux-limited approximation-based results fromdiffusion-limited radiation transport.Being thus prepared for bad news about the Biermann termfrom kinetic theory, we instead are met by a surprise uponHE BIERMANN CATASTROPHE 7examination of the MHD flux. The Biermann electric field at ashock is, according to Eqs. (14) and (18) E B = ( k B /e ) T e ln (cid:18) ρ d ( x ) ρ u ( x ) (cid:19) δ (Ψ( x )) ∇ Ψ . (34)Integrating this over a vanishingly short shock-crossing path l along the normal to the shock, we obtain the electromagneticwork done by an electron in the presence of the Biermannelectric field, E B = e (cid:90) l d x · E B = k B T e ln (cid:18) ρ d ρ u (cid:19) , (35)which is the same as the value inferred from kinetic theory,Eq.(7).We conclude from this argument that there is no analogy toflux-limited diffusion in the behavior of the Biermann sourceterm at shocks. The “plain vanilla” source gives the correctEMF across the shock, and needs no flux limiter to cut it offthere. It should be perfectly possible to construct a physicallyvalid algorithm to represent the Biermann MHD source termunmoderated by limiters, one that faithfully reproduces thepredictions of kinetic theory at a shock.3.3. Origin of the Numerical Biermann Catastrophe
As discussed in the Introduction, MHD simulations incor-porating the Biermann Battery effect have invariably pro-duced catastrophic non-convergent behavior as soon as dis-continuities develop. There are two (related) kinds of catas-trophes that have manifested themselves: simulations withspherical or planar symmetry, in which the charge-separationelectric field is irrotational and the magnetic field genera-tion rate should therefore be zero, develop a non-zero fieldat shocks with a field intensity that grows with increasingresolution (Fatenejad et al. 2013); and simulations of High-Energy-Density Physics (HEDP) laser experiments containingspatially-inhomogeneous shocks, similar to the (highly simpli-fied) simulations presented below in § E B of Eq. (15) is irrotational, in the sense that ∇ ××× E B = k B e ∇ T e ××× ∇ ln P e = 0 . (36)To what extent is this irrotational character preserved after T e and P e are separately discretized in a volume-based scheme?We consider the discretization of two scalar functions f ( x ) and g ( x ) that are assumed to have gradients that are ev-erywhere collinear, so that there exists some function α ( x ) such that ∇ g = α ∇ f . A pair of such functions satisfy ∇ g ××× ∇ f = 0 , so that the vector field g ∇ f is irrotational.The Taylor series of such a pair of functions are related toeach other by the collinearity constraint. We Taylor-expandboth functions to third order, average the resulting expansionover a grid of control volumes, and compute the “discretecurl” – that is, the circulation integral about a cell face – ofthe discretized “ g ∇ f ” vector field. We omit the details ofthis extremely lengthy calculation for the sake of brevity, andmerely present the conclusions that may be drawn therefrom.The result of the calculation is that the leading-order be-havior of the discrete curl of such a discretized vector field is O (∆ ) , where ∆ is the grid spacing. The coefficient multiply-ing this term is homogeneous of order 3 in the derivatives of f and g . What this means is that in regions of space wherethe two functions are smooth ( C or smoother), the coefficienthas a finite approximation in the discretization, and the dis-crete curl tends to zero as ∆ with increasingly fine resolution,yielding the correct curl of an irrotational vector field in thelimit.On the other hand, in the presence of a discontinuity, thecoefficient of the discrete curl has no finite approximation inthe discretization, but rather diverges as ∆ − when ∆ → ,reflecting the meaninglessness of the Taylor approximation inproximity to a discontinuity. As a consequence, the discretecurl as a whole diverges as ∆ − with vanishing ∆ in theneighborhood of a discontinuity.Setting f = ln P e , and g = T e , we recognize immediatelythe source of the discrete pathology described above. The dis-cretization of P e and T e produces a spurious solenoidal com-ponent of E B above and beyond any real, physically-correctsolenoidal component. This solenoidal anomaly is small andconverges to zero with increasing resolution in smooth flow.Near a discontinuity, however, the anomaly is not convergent,but rather diverges as ∆ − .This is the explanation of the numerical Biermann catastro-phe. It is, fundamentally, a completely predictable failure ofnaive, stencil-based approximations to the Biermann sourceterm (Eq. 2), which are not meaningful when Taylor-seriesapproximations to the fluid variables fail in the presence of adiscontinuity. The failure is irreducible, and doesn’t dependon whether the Biermann flux is added to other MHD fluxesor computed separately as a split term.3.4. The Cure for the Biermann Catastrophe
The above diagnosis of the origin of the numerical Biermanncatastrophe does not directly suggest a cure for the problem.That a cure should exist, however, is strongly suggested by thefact that the analytic treatment of the Biermann effect at shocksin § D by means of theRankine-Hugoniot relations, D × ( Φ L − Φ R ) = F ( Φ L ) − F ( Φ R ) . (37)Here, Φ is the state vector of conserved field quantities, Φ T ≡ [ ρ, ρ u T , ρ E , B T ] , that is, the mass, momentum, and energydensities and the magnetic field. F denotes the vector of fluxescorresponding to Φ . The subscripts “L” and “R” denote “Left”and “Right” states, respectively, and D is positive for a wavetraveling from Left to Right.(Toro 2009).We can now state a minimal requirement in order to assimi-late the Biermann flux to the other MHD fluxes used to updatethe fluid state: We need to establish the terms added by theBiermann effect to the flux vector F . Once these additivefluxes are determined, they can be added to the ideal MHDfluxes.The additional flux of B due to the Biermann effect is ob-tainable by solving for the shock part of the “Biermann-only”induction equation, Eq. (14). That is to say, we assume thedistributional forms of Eqs. (16), (19), and (26), plug them intoEq. (15), and keep only the terms proportional to the dirac- δ .The result may be obtained by inspection of Eq. (32), settingthe velocities u to zero: D ( B T u − B T d ) = − ck B e ln (cid:18) P e,u P e,d (cid:19) ∇ T e ××× n , (38)Comparing Eqs. (38) and (37), and using the continuity of ∇ T e ××× n , we obtain for the Biermann effect flux of B along n f B ( n ) ≡ − ck B e ln( P e ) ∇ T e ××× n . (39)This expression is manifestly well-defined at discontinuities.As in Eq. (32), from which Eq. (39) is derived, the only residueof the previously toxic gradient ∇ P e is contained in the benigndirection vector n , which is effectively −∇ P e / |∇ P e | at theshock. The effect of the cross product n ××× ∇ T e is to elim-inate the normal component of ∇ T e , while leaving only thetangential components. It is clear from the form of Eq. (39)that only the tangential components of B are subject to changeaccording to the Rankine-Hugoniot condition, as expected.It is convenient for the purposes of algorithmic implementa-tion to re-express Eq. (39) in tensor form, as an anti-symmetrictensor G ( B ) whose components G ( B ) ik express the flux in the coordinate direction k of the magnetic field component B i : G ( B ) ik ≡ − ck B e ln( P e ) (cid:88) j =1 (cid:15) ijk ∂T e ∂x j , (40)where (cid:15) ijk is the totally anti-symmetric Levi-Civita tensor. G ( B ) ik is, of course, the spatial part of the dual Maxwell tensordue to the Biermann effect. It is worth pointing out that thisexpression for the flux tensor differs from the “naive” flux G ( B,Naive ) ik ≡ + ck B e T e (cid:88) j =1 (cid:15) ijk ∂ ln( P e ) ∂x j , (41)obtainable directly from the Biermann electric field (Eq. 14)by a pure gradient, which has no effect on the induction equa-tion. It follows that Eq. (39) yields a general flux of B in thedirection n that can be used anywhere in the fluid, not only atdiscontinuities.The flux in Eq. (39) is not the only correction that mustbe applied to F . There is also a correction to be applied tothe energy part of the flux vector, in virtue of the Poyntingflux ( c/ π ) E B ××× B that arises in connection with the charge-separation electric field E B . This term is a bit puzzling atfirst: since E B is normal to the discontinuity, the Poyntingflux is tangential to the discontinuity, and it is not immediatelyobvious how such a flux is to be pressed into service in aRankine-Hugoniot relation.Nonetheless, the required flux may be inferred in a manneranalogous to the calculation by which we derived F B , above:We solve the “Biermann only” energy equation ∂ ( ρ E ) ∂t + ∇ · (cid:16) c π E B ××× B (cid:17) = 0 , (42)inserting the distributional forms of Eqs. (16), (19), (26), aswell as ( ρ E )( x , t ) = ( ρ E ) u Θ (Ψ) + ( ρ E ) d Θ ( − Ψ) . (43)After collecting terms proportional to the Dirac δ , we obtain D [ ρ E ] ud = ck B πe (cid:2) (( T e ∇ ln P e ) ××× n ) · B − ( ∇ ( T e ln P e ) ××× n ) · < B > + T e ln P e n · ∇××× < B > (cid:3) ud , (44)where in Eq. (44) we distinguish between B — the magneticfield to the left or right of the interface — and < B > ≡ ( B u + B d ) , the average of the upstream and downstreamfields at the interface. Comparing Eq. (44) and Eq. (37), weobtain for the energy flux along n due to the Biermann effect f ρ E ( n ) ≡ ck B πe (cid:2) (( T e ∇ ln P e ) ××× n ) · B − ( ∇ ( T e ln P e ) ××× n ) · < B > + T e ln P e n · ∇××× < B > (cid:3) . (45)Here we see a potential problem: the third term in Eq. (45)requires an average of the upstream and downstream currents ∇ ××× B . This is a difficulty if a true flux is to be extractedand pressed into service in a numerical scheme, since possiblediscontinuities in B makes such a term ambiguous. The reso-lution of the difficulty is that in general, there is always somefinite resistivity in real plasmas. As we will see below in § B , allowing us to replace < B > with B . We then obtain forthe flux of energy in the direction n f ρ E ( n ) = ck B πe ln P e [ ∇ ××× ( T e B )] · n . (46)To summarize, at a cell interface with a normal vector n , thehydrodynamic fluxes are to be adapted to the Biermann effectby adding to the implemented flux vector F an additional fluxvector F ( B ) ( n ) , given by F ( B ) = F ( B ) ρ F ( B ) ρ u F ( B ) ρ E F ( B ) B ≡ f ρ E ( n ) f B ( n ) . (47)We expect a discretization of the Biermann effect based onthese flux expressions to be convergent so long as the meshresolves the scales on which T e and B are continuous . Thescale at which T e is continuous is the scale length λ T of theelectron thermal conduction precursor region (Shafranov 1957;Jaffrin & Probstein 1964; Zel’dovich, Y.B. and Raizer, Y.P.2002), discussed in § λ T (and hence traverse its own thermal precursor). Itis therefore the scale to be resolved in order for the discretiza-tion of Eq. (39) to converge. At coarser resolutions, T e appearsas discontinuous at shocks as all the other fluid variables, andthe discretization discussed here will not yield converged re-sults for B . Only as the thermal precursor zone is resolved canconvergence be expected.In order for the discretization to produce convergent resultsfor energy, it is necessary that the simulation resolve the re-sistive length scale discussed in § TWO NOVEL PHYSICAL EFFECTS ARISING FROMSHOCK BIERMANNWe now exhibit two novel predictions of the theory of theBiermann effect at shocks. They are the existence of two mag-netic precursors to the shock wave, which lead the wave inregions whose extent depends on the upstream conditions: a resistive magnetic precursor , which arises due to magneto-generation in the shock “leaking” into the upstream fluid inconsequence of the presence of finite non-zero resistivity; anda thermal magnetic precursor , due to smooth-flow Biermann-effect magnetogeneration in the fluid motions set up by theelectron thermal precursor. We discuss these in turn.4.1.
The Resistive Magnetic Precursor
When the resistivity η in the induction equation, Eq. (3) isfinite, a new effect appears in the vicinity of a discontinuity: a resistive magnetic precursor travels ahead of the disconti-nuity, as the impulsively-generated field due to the Biermanneffect at the shock diffuses out. This effect is analogous tothe well-known electron thermal precursor that precedes aplasma shock, which was discussed in § ∂/∂t = 0 in the fieldequations. The induction equation then becomes ∇· ( Bu (cid:48) − u (cid:48) B )+ ∇· (cid:18) c η π ∇ B (cid:19) + ck B e ∇ T e ×××∇ ln P e = , (48)where u (cid:48) is the fluid velocity in the rest frame of the shock.We will keep only the shock part of the Biermann flux inEq. (48), that is, the right-hand side of Eq. (24). In effect,this amounts to assuming that the magnetic field generationfrom the Biermann effect is impulsive at the shock, so that thesmooth-flow Biermann effect is much smaller than the effectat the shock. This approximation is justified if the size λ T ofthe electron thermal precursor is much smaller than the sizeof λ B of the magnetic precursor, so that thermal gradients areunimportant except at the discontinuity. The plausibility ofthis condition will be verified in § Ψ( x , t ) = n · x , |∇ Ψ | = | n | = 1 for the discontinuity-tracing level function.We further adopt coordinates such that x is along n , so that n = e , and such that x is along ∇ T e ××× n . Away from thediscontinuity, we ignore all derivatives except for ∂/∂x .The presence of the resistive term changes the discontinuousstructure of the solution. This term is the divergence of theresistive flux of B . That flux behaves analogously to the Fick’slaw heat flux − κ ∇ T , in that it opposes gradients in B . Adiscontinuous B is disallowed, because it would result in aninfinite restoring resistive flux. B is therefore now continuous.By the Rankine-Hugoniot condition for transverse momentum(Gurnett & Bhattacharjee 2005, Chapter 7), it follows that u (cid:48) T is also continuous. Only u (cid:48) n and P e have solution discontinu-ities in this case. Assuming the upstream fluid is at rest in thelab frame, we therefore set u (cid:48) = u (cid:48) n e , with u (cid:48) n = ( − D ) Θ( x ) + ( u nd − D ) Θ( − x ) . (49)By the coordinate choice, assuming the far upstream fluid isunmagnetized, and in virtue of Eq. (48), we may set B = B e .Putting this all together, we obtain the following structureequation for the magnetic precursor: ddx (cid:20) c η π dBdx − u (cid:48) n B (cid:21) + δ ( x ) ck B e |∇ T e ××× n | [ln P e ] ud = 0 . (50)The meaning of this equation is that the precursor structureis determined by the balance of magnetic diffusion and mag-netic advection, in the presence of impulsive magnetic fieldgeneration due to the Biermann effect.0 GRAZIANI ET AL.In the upstream ( x > ) region, we may neglect gradientsin η (which depends on T e ). Eq. (50) becomes c η π d Bdx + D dBdx = 0 , (51)the solution of which, given lim x → + ∞ B ( x ) = 0 is B u ( x ) = B exp (cid:18) − πDc η x (cid:19) = B exp (cid:18) − x λ B (cid:19) , x > . (52)Here, B is an integration constant. We see that the precursorhas an exponential shape and a characteristic length λ B givenby λ B = c η πD . (53)Here, η = η ( T e, + ∞ ) is the value of the resistivity far upstreamof the discontinuity.We may obtain a relation for B from the jump condition atthe discontinuity, by integrating Eq. (50) across a vanishingly-small shock crossing path. The result is c η ( T e )4 π (cid:20) dBdx (cid:21) ud − [ u n ] ud B + ck B e |∇ T e ××× n | [ln P e ] ud = 0 , (54)where now η ( T e ) is evaluated at the shock. This conditiontogether with the upstream solution and boundary conditiondetermine B .The ideal MHD jump conditions may be recovered in thelimit of spatially-constant η with η → . In that case, thedownstream solution satisfying finite- B boundary conditionsat x → −∞ is B d ( x ) = B , and the first term in Eq. (54) isjust − DB . Eq. (33) follows immediately.4.2. The Electron Thermal Precursor
One further consequence of the presence of the thermalprecursor in T e described in § λ T , given byEq. (6). In general, the field intensity due to this effect can beexpected to be small compared to the field intensity due to theresistive magnetic precursor, since the gradients in the thermalmagnetic precursor are small compared to those available nearthe shock itself. In the constant-conductivity simulations thatwe discuss in § κ e ∼ T / e (Spitzer 1962).The actual structure of the thermal precursor is not the gentleexponential decay that one expects for a constant κ e , but ratherexhibits the steep gradient near its outer terminus shown inFigure 7.20 of Chapter VII of Zel’dovich, Y.B. and Raizer, Y.P.(2002). It is possible that this large gradient may be responsiblefor an enhancement of the Biermann effect at the terminusof the thermal conduction precursor zone. This interestingpossibility is beyond the scope of this paper. One consequence of the presence of a thermal magneticprecursor is that even in a non-resistive approximation, thenotion of an un-magnetized fluid upstream of the shock, onwhich, for example, Eq. (33) is based, is not strictly correct. Itis a valid approximation only when the Biermann generationrate in the thermal precursor is small compared to that of theshock. 4.3. Physical Length Scales
It is useful to establish some of the relevant length scalesunder conditions of interest. Below, we establish these scalesfor conditions prevailing in HEDP experiments. In §
5, we willestablish the scales characteristic of galaxy cluster formation.Using expression for η from Huba (2007) in the expressionfor the resistive magnetic precursor length λ B given in Eq. (53),we have the following expression for λ B in a fully-ionizedplasma: λ B = (4 / π ) − / c e m / e Z ln Λ D − ( k b T e ) − / = 16 . cm × Z × ln Λ10 × (cid:18) D cm s − (cid:19) − × (cid:18) k B T e eV (cid:19) − / , (55)where Λ is the usual Coulomb logarithm. Similarly, the ex-pression for the electron thermal precursor length scale λ T is(Zel’dovich, Y.B. and Raizer, Y.P. 2002, Chapter VII, § λ T = 25 κ e ( T e, ) ρc v,e D = 415 ξ ( k B T e ) / e m / e ZA ln Λ n i D = 5 . × − cm × ξ × Z − A − × (cid:18) ln Λ10 (cid:19) − × (cid:18) D cm s − (cid:19) − (cid:18) k B T e eV (cid:19) / (cid:16) n i cm − (cid:17) − , (56)where T e, is the electron temperature at the shock, and where ξ is a number in the range 1–2.The numbers in Eqs. (55) and (56) have been scaled toconditions that are routinely obtainable in large Laser facilitiessuch as Omega and NIF (see, for example Gregori et al. 2012).In these conditions, it is clear that the assumption λ T (cid:28) λ B ,required for the validity of the derivation of λ B , is easilysatisfied.It is also immediately apparent that the magnetic precursorlength scale is a macroscopic scale that can in principle be well-tailored to the physical dimensions of the interaction chambersof such facilities. A carefully-designed experiment, whichlaunches a shock into a relatively cold upstream plasma (notethe dependence of λ B on T e ) should in principle be able todetect the resistive precursor due to the Biermann effect.Note that as the value of T e upstream rises, λ B decreases as T − / e , while λ T increases as T / e . We therefore only needto increase the upstream temperature to about 7 eV for λ T to become comparable to λ B , and at warmer temperaturesthan this λ T becomes dominant. It is therefore plausible thata physical situation could be created in which the thermalHE BIERMANN CATASTROPHE 11magnetic precursor might be observable without interferencefrom the resistive magnetic precursor. THE BIERMANN EFFECT AND GALAXY CLUSTERFORMATIONIn the present section, we attempt to estimate the strengthof the seed magnetic fields that result from a correct treatmentof the Biermann effect at shocks in galaxy clusters at the timeof magnetogenesis ( z ∼ − ), and the physical length scalesof the resistive magnetic precursor and the electron thermalprecursor in this context.5.1. Proto-Galaxy Cluster Field Generation
As discussed in the introduction, the Biermann effect hasbeen invoked as the source of seed magnetic fields that can beamplified by the turbulent dynamo mechanism (Kulsrud et al.1997).Given the fact that in these studies, the gradients near shocksthat were used to calculate the Biermann effect field generationrates are artifacts of the hydrodynamic advance schemes, it ispossible – even likely – that the calculated field strengths areincorrect. We now attempt to estimate the typical strength ofmagnetic fields that are to be expected in early galaxy clusterformation, based on a corrected treatment of the Biermanneffect. We do not perform new simulations, but rather attemptto estimate the typical field strength based on published simu-lation data, in a preliminary effort to determine the reliabilityof field strengths inferred to date in studies of galaxy clusterformation.While this is not an easy task without actual simulation dataon-hand to analyze, the analysis work described in Miniatiet al. (2000), based partly on Λ CDM simulations described inCen & Ostriker (1994) provides enough information for us tomake a rough estimate of the typical field strength. The simu-lation in question has the properties H = 60 km s − Mpc − , Ω M = 0 . , Ω Λ = 0 . , Ω b = 0 . . We now give a rela-tively detailed description of our analysis of the information inMiniati et al. (2000), in order to make clear both the basis forthe estimated field strength and the considerable uncertaintythat attend these estimates.Our starting point is Eq. (33), which gives the post-shockfield strength due to the Biermann effect, assuming an un-magnetized upstream fluid. In order to use this equation toestimate field strengths, we need typical values of the shockspeed D , the compression ratio ρ d /ρ u , and the tangential gra-dient |∇ ⊥ ( k B T e ) | ≡ |∇ ( k B T e ) ××× n | that arise at primordialepochs.For the sake of exploiting the information available inMiniati et al. (2000), we will choose z = 3 as our fiducialprimordial epoch. Examining the upper-right panel of Figure 8in Miniati et al. (2000), we conclude that a not-atypical valueof |∇ ⊥ ( k B T e ) | is |∇ ⊥ ( k B T e ) | ≈ k B × K Mpc ≈ − erg cm − , (57)Neglecting momentarily the fact that the simulation certainlyfails to resolve the thermal precursor length scale λ T (Eq. 60below) the uncertainty in this estimate seems to be a factor of2 or so.Miniati et al. (2000) supplies shock speeds at redshift z = 0 (Figure 5 of Miniati et al. (2000)). These may be approximatelyshifted to z = 3 using Figure 10 of Miniati et al. (2000), whichshows that the kinetic energy processed by shocks increased by a factor of 15 between z = 3 and z = 0 . This suggeststhat temperatures increased by about that much (neglectingshock filling-factor differences between the two redshifts), andthat velocities increased by about a factor of 4. We had typicalshock temperatures of ∼ K at z = 3 in the data leadingto the estimate in Eq. (57). This corresponds to a temperature ∼ K at z = 0 . From the top-right panel of Figure 5 ofMiniati et al. (2000), this corresponds to a shock speed ofabout × cm s − at z = 0 , and hence D ∼ cm s − at z = 3 .Putting these numbers in Eq. (33), and assuming the strong-shock limit ρ d /ρ u = 4 appropriate to a γ = 5 / ideal gas, weobtain B = 2 . × − G × (cid:18) D cm s − (cid:19) − × (cid:18) |∇ ⊥ k B T e | − erg cm − (cid:19) . (58)This value should be regarded as a lower limit, since, asremarked above, the simulation does not resolve λ T . An upperbound on B can be obtained by replacing, in Eq. (57) the5 Mpc length scale estimated from Figure 8 of Miniati et al.(2000) with λ T from Eq. (60). This gives B (cid:46) − G;i.e., a value that is a factor of about higher. Obviously,this is a very conservative bound, as it assumes temperaturefluctuations of order unity over the diffusive scale λ T .It is clear from the very uncertain nature of the factors usedto construct this estimate that the field strength given in Eq. (58)is itself subject to considerable uncertainty. A few actual MHDsimulations with a correct implementation of the Biermanneffect seem required to establish how much field strength isin fact made available by the Biermann effect, for turbulentdynamo effects to amplify. Such simulations would be chal-lenging given the spatial scale λ T that needs to be resolved.5.1.1. Galaxy Cluster Formation Length Scales
We now estimate the physical length scales of the resistivemagnetic precursor and the electron thermal precursor for thecase of galaxy cluster formation using the physical parametersinferred in § z (cid:39) . As in § T e = 10 K as our refer-ence temperature and D = 10 cm s − as our reference shockspeed. The assumed cosmological parameter Ω b = 0 . corresponds to a proton number density n i ≈ − cm − at z = 0 , and hence to n i ≈ × − at z = 3 , which we takeas the density upstream of the shocks. With these parameters,the Coulomb logarithm is ln Λ ≈ . Setting Z = 1 , A = 1 ,we then have λ B = 8 . × − cm × ln Λ40 × (cid:18) D cm s − (cid:19) × (cid:18) T e K (cid:19) − / , (59)and λ T = 5 . kpc × ξ × ln Λ40 × (cid:18) D cm s − (cid:19) × (cid:18) T e K (cid:19) / × (cid:18) n i × − cm − (cid:19) − . (60)The much more tenuous and considerably hotter primordialplasma reverses the relative sizes of λ B and λ T comparedto the HEDP case considered above: λ B is now negligible,while λ T is an astrophysically-significant length. It is worth2 GRAZIANI ET AL.comparing λ T to λ e and λ i , the expected Spitzer mean freepaths of electrons and ions in the plasma, as a consistencycheck. This is the product of the mean collision time τ withthe thermal speed (3 k B T /m ) / . Using the expression for τ e from Huba (2007), we have λ e = 3 m / e ( k B T e ) / √ πn e ln Λ e × (cid:18) k B T e m e (cid:19) / = 3 √ √ π ( k B T e ) n e ln Λ e = 2 . × − kpc (cid:18) T e K (cid:19) × (cid:18) ln Λ40 (cid:19) − × (cid:18) n e × − cm − (cid:19) − . (61)The mean free path is independent of particle mass, and for Z = 1 it is the case that λ e = λ i . This length is the charac-teristic size of the viscous shock sheath, and is reassuringlyshorter than the semi-hydrodynamically scaled λ T .It is also worth considering whether the plasma is collisional,as we have implicitly assumed by using Spitzer-type transportcoefficients. The magnetic field strength estimated in § λ c of characteristic size λ c = (cid:18) k B T e m e (cid:19) / (cid:18) eBmc (cid:19) − = 41 . kpc (cid:18) T e K (cid:19) / × (cid:18) B × − G (cid:19) − . (62)Since λ c ∼ λ e , the flow is comfortably collisional for thechosen physical parameters. An increase in the estimate of B by two orders of magnitude or more would change thisconclusion; however, it should be noted that at least at the out-set of the process of cosmic magnetogenesis envisioned here,field strengths were probably small enough that the collisionalplasma approach is valid. Primordial field strengths that mayhave originated in early-universe phase transitions or in infla-tionary scenarios are poorly constrained by theory and obser-vation, but are believed to have plausible values B < − Gin all but a few scenarios (see Widrow 2002; Widrow et al.2012). It follows that initially, at least, the Spitzer-type trans-port coefficients that we have employed here are likely to beappropriate. This is in contrast to the apparent suppression ofelectron conductivity at later times (see, e.g. Markevitch et al.2003; Russell et al. 2012; Sanders et al. 2013; ZuHone et al.2013). NUMERICAL VERIFICATION6.1.
Implementation Using the FLASH Code
We have implemented the corrected Biermann Effect al-gorithm, described in § FLASHUsers Guide . We plan to include an implementation of thenew Biermann effect algorithm in a future release of the code.The 2-temperature MHD model that we use for our numer-ical tests is expressed in conservation form by the followingdynamical system: ∂ρ∂t + ∇ · ( ρ u ) = 0 (63) ∂ρ u ∂t + ∇ · (cid:20) ρ uu − π BB + (cid:18) P + B π (cid:19)(cid:21) = 0 (64) ∂∂t (cid:20) ρ (cid:18) u + (cid:15) T (cid:19) + B π (cid:21) + ∇ · (cid:20) ρ u (cid:18) u + (cid:15) T + Pρ (cid:19) +14 π (cid:18) − u ××× B + c η π ∇ ××× B (cid:19) ××× B + ck B πe ln P e [ ∇ ××× ( T e B )] (cid:21) = 0 (65) ∂ρs e ∂t + ∇ ( ρ u s e ) = − T e − ∇ · ( − κ e ∇ T e )+ ρc v,e T e τ ei ( T i − T e ) (66) ∂ B ∂t + ∇ · (cid:20) uB − Bu − c η π ∇ B + G ( B ) (cid:21) = 0 . (67)These equations differ from Eqs. (8)-(12) only in that theBiermann term in Eq. (65) replaces the one in Eq. (10) toreflect Eq. (46), while the Biermann term in Eq. (67) replacesthe one in Eq. (12) to reflect Eq. (40).FLASH advances the MHD equations using an Unsplit Stag-gered Mesh (USM) scheme that prevents the development ofmagnetic monopoles due to numerical noise in the inductionequation, and which is described in (Lee 2013).The thermal conduction and heat-exchange source termson the right of Eq. (66) are computed in a time-split manner,separately from the ideal MHD advance. In particular, thermalconduction advance is fully implicit, and works as describedin § FLASH Users Guide , while heat-exchangeoperates as described in § FLASH Users Guide .The remaining, advective part of Eq. (66) is implemented bytreating s e as a passively-advected mass scalar. This equationrequires an EOS implementation capable of using electronentropy as an input variable. For the current case of perfect-gas EOS and total ionization, the standard Sackur-Tetrodeequation for entropy is implemented in the EOS.The resistive terms in Eqs. (65) and (67) are not treatedimplicitly, but rather are added directly to the MHD fluxes, andare thus treated explicitly, placing diffusive stability limitationson the maximum timestep. The Biermann terms in Eqs. (65)and (67) are also added explicitly to the MHD fluxes. Theseflux modifications are performed in such a way as to preservethe solenoidal character of the magnetic field, by adding themto the ideal MHD Godunov fluxes before these are interpolatedto the edge-centered electric fields of the USM scheme, afterwhich each face-centered normal magnetic field component isupdated by the circulation integral of the electric field alongthe edges bounding the face (Lee 2013).We have not performed an analysis of the timestep limitation http://flash.uchicago.edu/site/flashcode/user support/ HE BIERMANN CATASTROPHE 13imposed by the Biermann effect, trusting instead that the smallmagnitude of the effect in the cases considered makes it un-likely that Biermann timestep constraints could dominate thosedue to the hyperbolic and diffusive terms. In this connection, itis noteworthy that the Biermann term, being quadratic in fluidderivative terms, makes no contribution to wave dispersionrelations obtained by linearization about uniform solutions,and therefore does not appear to affect ordinary wave speeds.Intuitively, one would therefore expect any timestep limita-tions due to the Biermann effect to be higher-order corrections,hardly affecting the numerical evolution of plasmas for whichthe Biermann term is small in magnitude. Nonetheless, afull analysis of the effect of the Biermann term on the hy-perbolic structure of the system of PDEs would be useful –possibly even necessary – for cases of intense field generation.Such an analysis would be somewhat complicated by the non-quasilinear structure of the Biermann term, necessitating anextension of the PDE system to quasilinear normal form.6.2.
The Simulations
The simulations presented below represent idealized situa-tions simplified for the sake of verifying the code, and for illus-trating the numerical and physical principles under discussion.In particular, the conductivity κ e , electron-ion equilibrationtimescale τ ei , and the resistivity η are all treated as constants.FLASH does have the capability to use Spitzer-type functionsof the thermal state for these parameters, but this would com-plicate the results unnecessarily without adding any real valueto these verification tests.In what follows, we assume fully-ionized Hydrogen – A = 1 , Z = 1 , adiabatic index γ = 5 / , c v,e = k B N A , where N A is Avogadro’s number and k B is Boltzmann’s constant. Sim-ulations are conducted in cylindrical coordinates, assumingazimuthal symmetry, and are therefore 2-dimensional. In ad-dition, we impose a reflection boundary on the R -axis (where R is perpendicular distance from axis of rotational symmetry,which is to say, R is the cylindrical radius) so that the domainrepresents a hemisphere of a solution with reflection symme-try about that axis. In every case, the domain is 4 cm radiuscylinder that extends 4 cm in the z direction from the R axis.The boundary conditions at R = 4 cm and at z = 4 cm areoutflow.All of the verification tests that we present are variationson the theme of a Sedov-esque explosion: an analytic self-similar Sedov solution is (1) modified to prevent temperaturesfrom diverging at the center by setting all variables constantinside some chosen radius; (2) smoothed with a Gaussiannear the shock, to ease instabilities that can otherwise appearnear the shock; (3) where required, distorted to an ellipsoidalprofile capable of producing real Biermann effect fields; and (4)scaled to velocities and pressures estimated to produce usablesimulation times given the resolutions and domain size in play.Values of shock ellipticity, κ e and η are then chosen to producereasonable thermal and magnetic precursor region sizes, andreasonable field generation rates, and values of τ ei are chosento be harmless. In every case, we run an initial model for atime without the Biermann effect, so as to allow transients inthe 2-temperature hydrodynamic variables to subside. Then werestart the calculation with Biermann field generation turnedon.Magnetic fields are always azimuthal in these verificationtests – since the gradients of T e and P e are always in the R − z plane the Biermann effect only generates non-zero field alongthe azimuthal direction. 6.2.1. Null Test: A Spherical Shock
It is not a simple matter to construct a non-trivial analyticverification solution of a plasma generating magnetic fieldthrough the Biermann effect. A trivial solution, on the otherhand, may be straightforwardly constructed by imposing sym-metry requirements that align the gradients ∇ P e and ∇ T e ,thus guaranteeing zero field generation. We use a sphericalSedov-like explosion as an example of such a test. This is infact the test that revealed the Biermann catastrophe in the firstplace (Fatenejad et al. 2013).We assume an initial shock radius of 1.15 cm, in an ambientmedium with ρ = 1 g cm − and P set to a negligible value.We choose an initial velocity scale inside the shock that leadsto a shock velocity that is on average about . × cm s − during the course of the simulation. Other physical parametersare κ e = 1 . × erg K − cm − s − and τ ei = 2 × − s. The initial solution is advanced for − s with nofield generation, then re-started and run for another − swith field generation turned on. We repeated this experimentwith the correct Biermann flux term of Eq. (40), with the naiveflux of Eq. (41), and with the time-split source term of Eq. (15),so as to compare the performance of the three algorithms.The results of this study are displayed in Figure 2. Thetop-left panel is a colormap view of the electron temperature T e . The dashed shock-crossing line across the bottom of thisfigure is the location of the data displayed in the top-rightpanel, which shows electron and ion temperatures. The shockis visible as the peak in T i . The thermal precursor is clearlyvisible, as is the continuous behavior of T e , which changesslope at the shock in accordance with standard plasma shocktheory (Zel’dovich, Y.B. and Raizer, Y.P. 2002, Chapter VII, § L -norm of the difference between the an-alytic and numerical solutions. The blue circles show theconvergence with resolution of the correct formulation. Thisconvergence behavior is in contrast to the behavior of the othertwo algorithms, which fail altogether to converge.6.2.2. Verification With An Ellipsoidal Shock
Next we exhibit simulations designed to produce non-zerovalues of the magnetic field strength. We start an ellipsoidalshock surface, obtained by distorting the spherical Sedov solu-tion. The initial configuration has a semi-major axis (alignedwith the R -axis) of 2 cm, and a semi-minor axis (aligned withthe z -axis) of 1.667 cm. The ambient medium is uniformwith ρ = 1 g cm − and P set to a negligible value. Wechoose a velocity scale that leads to shock velocities of theabout × cm s − . The higher shock speed is chosento produce somewhat intense magnetic field strengths. Tokeep the thermal precursor resolved, we increase the con-ductivity to κ e = 1 . × erg K − cm − s − . Wealso set τ ei = 10 − s. The initial solution is advanced for . × − s with no field generation, then re-started andadvanced to a simulation time t = 3 × − s with field gener-4 GRAZIANI ET AL. Figure 2.
Two-dimensional simulation of a spherical shock in a 2-temperature Hydrogen plasma. The images refer to the final timestep, at t = 2 × − s, and tothe highest resolution simulation. Top left: Electron temperature distribution. Top right: Electron and ion temperatures along the shock-crossing line shown at thebottom of the figure in the top left panel. Bottom left: B φ due to numerical noise. Bottom right: Total magnetic energy as a function of simulation resolution. ation turned on.The advance of the shock during the period when magneticfield is generated is illustrated by the two pressure colormapplots in the top panels of Figure 3. The lower-left panel of Fig-ure 3 shows the magnetic field intensity generated accordingto the correct flux of Eq. (40), ranging into the tens of Gauss.The solid red line in the figure shows the shock location, whilethe red arrows display the shock unit normal vector. The lower-right panel displays the difference between the correct and thenaive flux formulations, which is substantial at the outgoingshock surface.We do not have an analytical solution for the magnetic fielddistribution to compare to the output of these simulations.We do, however, have the relation between shock quantitiesexpressed by Eq. (32), which determines the jump conditionfor B . By locating points adjoining the shock, and computingthe local shock velocity at those points, we can then verify that Eq. (32) in fact obtains, to some accuracy limited by thenumerical approximation. This is in principle a non-trivialverification, since the code does not know about the jumpcondition Eq. (32) as such – it only knows the fluxes expressedby Eq. (40), which were inferred from the jump condition.Recovery of the jump condition thus constitutes a verificationthat the code correctly models the theory.To perform this verification, we first locate cells adjoiningthe shock by the method described in Appendix A, whichyields a 1-cell wide shock surface by fitting the mass, momen-tum, and energy Rankine-Hugoniot conditions to the neigh-boring data, using a speed-of-sound weighted inner-producton the state space, and treating the shock speed D as a freefit parameter. The shocked cells are the ones that fit the R-Hconditions with small fit residuals ( R < . ), together withsome other natural auxiliary consistency conditions describedin the Appendix. The fitted shock speed is then used in theHE BIERMANN CATASTROPHE 15 Figure 3.
Ellipsoidal shock. Top Left: Initial pressure distribution; Top Right: Final pressure distribution; Bottom Left: Final magnetic field distribution. The solidred line displays the location of the shock, while the red arrows display the direction of the shock normal; Bottom Right: Difference between final magnetic fielddistributions computed using the correct and incorrect Biermann flux. verification of Eq. (32).In the absence of normal component field B n , Eq. (32)becomes ( D − u d ) B d − ( D − u u ) B u + ck B e (cid:18) n z ∂T e ∂R − n R ∂T e ∂z (cid:19) [ln P e,d − ln P e,u ] ≡ a d − a u + b d − b u = 0 , (68)where we’ve defined “Advection” terms a u,d and “Biermann”terms b u,d by a d ≡ ( D − u d ) B d (69) a u ≡ ( D − u u ) B u (70) b d ≡ ck B e (cid:18) n z ∂T e ∂R − n R ∂T e ∂z (cid:19) ln P d (71) b u ≡ ck B e (cid:18) n z ∂T e ∂R − n R ∂T e ∂z (cid:19) ln P u . (72)The sum of terms in Eq. (68) is required to be zero. Ina discretized numerical PDE integration, this really meansthat the terms must cancel up to some truncation or roundingprecision, which is expressed relative to the largest of themagnitudes in play. We therefore define the “shock condition parameter” C by C ≡ a d − a u + b d − b u max ( | a d | , | a u | , | b d | , | b u | ) . (73)We calculate the value of C at cells along the shock front,at each of our four resolution levels, for both the correct fluxand the “naive” flux implementations of the Biermann effect.We display cumulative distributions of C in the top panels ofFigure 4. It is evident from these figures that the distribution of C is in fact centered near zero for both flux implementations.The width of the distributions behave very differently, however.In the case of the correct Biermann flux implementation, thedistributions get narrower with each refinement of resolution,providing some evidence of convergence with resolution tothe expected result. In the case of the naive flux implemen-tation, there is no such evidence of convergence. The sameobservations can be made of the plots in the middle panels ofFigure 4, which summarize the means and sample standarddeviations of the C -distributions, as a function of resolution.The convergence of the correct flux implementation, and theconvergence failure of the naive implementation, are clear inthese figures.The bottom panels of Figure 4 show the evolution of totalmagnetic energy in the domain as a function of time, for the6 GRAZIANI ET AL. Figure 4.
Left Panels: Correct Biermann flux term. Right Panels: “Naive” Biermann flux term. Top Panels: Cumulative distribution of the normalized magneticshock condition C , defined by Eq. (73), evaluated at points along the shock surface, for four different resolutions, illustrating convergence to the correct jumpcondition. Middle Panels: sample standard deviations of C , as a function of resolution. Bottom Panels: Total magnetic energy as a function of simulation time forthe four resolutions studied. HE BIERMANN CATASTROPHE 17
Figure 5.
Left Panel: Magnetic field distribution due to passage of shock, in the presence of finite resistivity. The simulation differs from the one shown inFigure 3 only by the presence of non-zero resistivity. The solid red line shows the location of the shock, while the vectors represent the unit normal to the shock.The magnetic precursor can be clearly seen ahead of the shock surface. The white dashed shock-crossing line illustrates the location of the data displayed in theright panel. Right Panel: Magnetic field strength magnitude along the white dashed shock-crossing line shown in the left panel. The position of the shock is shownby the vertical solid red line. The predicted exponential decay of the precursor is evident. different resolutions and for the two flux implementations. Thecorrect flux implementation appears to be converging (althoughnot in any strong sense converged) by this measure, even at latetime, whereas any evidence of convergence in magnetic fieldenergy is simply lacking for the naive flux implementation.6.2.3.
The Resistive Magnetic Precursor
In this set of simulations, we maintain the simulation pa-rameters described in § η to a finite positive value, η = 7 . × s − , cho-sen in conjunction with the other parameters to yield aneasily-discernible exponentially-decaying resistive precursordescribed by Eq. (52). We repeat the simulation strategy of § t = 4 . × − swith no field generation, then re-starting it and advancing to asimulation time t = 3 × − s with field generation (usingthe correct flux formulation) turned on. This advance timeis adequate for the establishment of the expected precursorstructure, as we demonstrate below.The left panel of Figure 5 shows a colormap of the distribu-tion of magnetic field strength across the domain. The shocklocation is shown by the solid red line, and the superposed vec-tors indicate the shock-normal direction. The magnetic field isevidently smoothed out by the resistivity, as can be seen by adirect comparison with the bottom-left panel of Figure 3. Theprecursor is already evident in this figure, as the substantialamount of magnetic field that has “leaked” ahead of the shock.The diagonal white dashed shock-crossing in the left panelof Figure 5 illustrates the line along which data was extractedto produce the right panel of Figure 5. This figure showsthe magnetic field values plotted along that line. The verticalsolid red line at X = 0 marks the location of the shock. Theexponential decay of the field is easily visible in this figure.At each shock location, the value of η can be combinedwith the local shock velocity to calculate λ Predicted, the ex-pected decay length of the precursor, according to Eq. (53).At the same time, the actual decay length λ Measured can
Figure 6.
Cumulative distribution of the measured magnetic precursor decaylength, normalized to the predicted length, at 849 points along the shocksurface. The mean is shown by the dotted line, while the standard deviation isillustrated by the errorbar. be measured directly by comparing the field strength at twosuitably-selected distances along the local shock normal. Wehave performed both calculations at each shock location, andcompared them. The results are displayed in Figure 6, wherewe plot the cumulative distribution of λ Measured /λ Predicted.The distribution appears to be consistent with a mean value of1, as expected, with some scatter. The scatter is not surprising,since the expression in Eq. (53) for λ Predicted is an approx-imation, assuming, as it does, propagation of magnetic field8 GRAZIANI ET AL.from a planar shock. Since the shock is, of course, not planar,the value of the Biermann generation rate changes apprecia-bly across the shock over distances not hugely different from λ Predicted itself. Under the circumstances, then, the level ofagreement between observation and prediction is satisfactory.6.2.4.
The Thermal Magnetic Precursor
In our final simulation, we explore the properties of thethermal magnetic precursor by returning to a non-resistivesimulation, similar to the ones discussed in § κ e is ten times larger: κ e =1 . × erg K − cm − s − . This broadens the thermalprecursor zone to about 0.2 cm, making it easier to discern.The results at the final time step are illustrated in Figure 7.The top-left panel shows the distribution of magnetic field inthe domain. The solid red line in the figure shows the locationof the shock, while the arrows show the direction of the shocknormal. In the top-right panel we see the electron and ion tem-peratures along the shock-crossing line shown in the previouspanel. The shock location coincides with the sharp drop in T i , while the thermal conduction precursor zone, which has anexponential structure for the constant conductivity model usedhere, may be clearly seen in the upstream behavior of T e .The magnetic field structure along the shock-crossing lineis shown in the bottom-left panel of the figure. The precursoris difficult to see in this plot, so the bottom-right panel shows | B | on a semi-log plot. The magnetic precursor structure isseen in this figure. It clearly has a much lower field intensitythan the fields generated by the Biermann effect at the shock.We emphasize again, however, that in part this is due to thehighly simplified constant-conductivity model adopted for thiswork. A true Spitzer-type conductivity, with a T / e depen-dence, would create much sharper gradients at the upstreamtermination of the precursor zone, potentially generating moreintense fields – relative to those generated at the shock – thanthose shown here. SUMMARY AND DISCUSSIONTo summarize our findings: Using the kinetic theory ofplasma shocks, we have given a description of the Biermanneffect within ion viscous shocks. Using this description, wehave shown that the convergence failures in the presence ofshocks of MHD codes implementing the Biermann effect isnot traceable either to a mathematical mis-specification of theBiermann term – which is well-defined despite appearances –or to a failure of the Biermann MHD source term to correctlyapproximate the expected physics. The failure is instead due tonaive discretization, a condition that we have shown is curableby exploiting the continuity of the electron temperature T e .We have described a convergent algorithm that incorporatesthe Biermann effect within shocks in numerical MHD, andimplemented it in the FLASH code. We have developed andused verification tests to demonstrate formal convergence usinga null solution from spherical shock, and to verify predictionsfor physical conditions near a shock using an ellipsoidal shock.Comparisons of the new algorithm, which provides a correctand accurate treatment of the Biermann effect within shocks,with the previous, naive one exhibit clearly desirable physicaland numerical properties present in the new algorithm thatwere previously wanting.We have noticed, described, and simulated two previouslyunrecognized physical effects: a resistive magnetic precursorthat leads the shock in the presence of non-zero resistivity, and that is analogous to the well-known thermal precursor causedby the presence of finite electron thermal conduction; and athermal magnetic precursor, produced by plasma motions inthe electron thermal conduction precursor. We have estimatedthe expected magnitude of both effects in conditions encoun-tered in laboratory experiments at laser facilities, and shownthat the characteristic length of the two effects can be mademacroscopic. In particular, the resistive magnetic precursor isphysically measurable in an experimental setup that measuresboth the position of the shock front and a time-series of the in-tensity of the magnetic field at some location initially ahead ofthe shock. Such an experiment would have to somehow ensurethat the upstream fluid remains unheated by shock-generatedradiation and by laser absorption, so as to keep the value of theresistivity high enough to sustain a macroscopically-scaled pre-cursor length. The thermal magnetic precursor may be morechallenging to observe, as it is intrinsically weaker than theresistive precursor (relative to the field intensity generated atthe shock). Further studies with Spitzer-type electron thermalconductivity are needed to determine whether sharp gradientsat the upstream termination of the electron thermal precursorcan give rise to sufficiently intense magnetic field strengths tobe experimentally measurable.We emphasize again that the requirement for convergencein B of the new Biermann effect algorithm is that the electronthermal precursor of the shock should be resolved, while therequirement for convergence in magnetic energy of the newalgorithm is that the resistive precursor should be resolved. Ofthese two requirements, the first is more essential, becauseneglect of the contribution of the Biermann effect to the fluxof magnetic energy is frequently a valid approximation. De-manding resolution of the electron thermal precursor is not atrivial requirement in realistic simulations. It certainly man-dates some kind of Adaptive Mesh Refinement (AMR) strategyto resolve the region near the shock. But irrespective of AMR-related efficiencies, merely resolving that length scale canconstrain the code to taking very short time steps due to theCourant-Friedrichs-Lewy (CFL) stability condition, unless animplicit advancement scheme is implemented for MHD. Thisissue is the subject of current study, but it lies beyond the scopeof the present paper.With the discussion of the modification due to the Biermanneffect of the Rankine-Hugoniot jump condition on B – Eq. (32)and text thereabout – it now seems worth returning to the ques-tion of vorticity ω , and the beguiling Equation (5) that suggestsits proportionality to B . Recall that much hangs on the jumpconditions – if the jump in ω maintains the same proportional-ity to the jump in B as do their respective evolution equations,then, subject to some cavils about boundary conditions, B issimply ω , at least until resistivity and/or viscosity manifestthemselves, and in effect, unmagnetized gasdynamics containswithin it the magnetic degrees of freedom of MHD, encodedin derivatives of the velocity field.It can be seen more clearly now why this implausible as-sertion is in fact false. In the first place, B and ω are ratherdifferent types of quantities from the point of view of kinetictheory, in that B is perfectly well-defined at kinetic scales –such as inside an ion viscous shock sheath – whereas ω can-not even be defined consistently in regions where the fluidpicture breaks down. The jump condition on B is a straightfor-ward consequence of the induction equation, and is essentiallykinematic in origin, as are the other Rankine-Hugoniot con-ditions. The jump condition on vorticity, on the other hand,cannot be inferred from the dynamical PDE equations alone,HE BIERMANN CATASTROPHE 19 Figure 7.
Ellipsoidal shock, with 10 × greater thermal conductivity. Top Left: Final magnetic field distribution. The solid red line shows the shock location, whilethe arrows show the direction of the shock normal. The white dashed shock-crossing line illustrates the location of the data displayed in the remaining panels; TopRight: Electron and ion temperatures along the white dashed shock-crossing line in the previous panel. The shock location and thermal precursor are evident;Bottom Left: B along the shock-crossing line; Bottom Right: Semilog plot of | B | along shock-crossing line. but depends in an essential manner on the equation of state(Kevlahan 1997). It would be possible, in principle, to infer a“flux” of ω from Eq. (4), and to construct a “jump condition”using that flux in the Rankine-Hugoniot conditions, Eq. (37).That jump condition would disagree with the one obtained in(Kevlahan 1997), as it is manifestly independent of the EOS.The EOS dependence of the jump condition on vorticity im-plies that even if there exists an EOS according to which thejump condition on ω preserves the required proportionalityto the Biermann-laden jump in B , that proportionality wouldcertainly not be preserved for any other EOS.In general, then, the passage of a shock certainly spoilsEq. (5). This should be no surprise, as the manipulations ofthe hydrodynamic equations of motion required to arrive atEq. (4) constitute precisely the sorts of prestidigitation withvector derivatives that break down at fluid discontinuities (seep. 216 of LeVeque 2002, for example). While Eq. (4) may be derived from the momentum equation on either side of theshock, the connection between B and ω is certain to be lost atthe shock itself, well before resistive and viscous effects canassert themselves. We conclude that the vorticity connectionis not a useful tool for studying magnetogeneration by theBiermann effect in the presence of shocks.In concluding, it is worth pointing out again that since shocksare not the only weak-solution discontinuities that arise inmulti-material MHD simulations, they are not the only loca-tions where potential convergence failures are corrected bythe new algorithm. In particular, contact discontinuities, mate-rial discontinuities, and ionization fronts all present potentialproblems for the naive Biermann algorithm, since they all rep-resent locations where P e changes discontinuously (despitethe continuity of total pressure P at such locations), and aretherefore all sites where the Biermann effect can be expected togenerate magnetic field. They all represent physical situations0 GRAZIANI ET AL.in which the new algorithm provides a correct treatment of theBiermann effect.We would like to thank Gianluca Gregori for discussions thathelped to focus this work on the relevant physical issues. Wewould also like to thank the anonymous referee for encouragingus to discuss more fully the astrophysical implications of thiswork; and Andrey Kravtsov, Brian O’Shea, and John ZuHonefor correspondence and conversations about the application of this work to the creation of magnetic fields by the Biermanneffect in galaxy clusters. This work was supported in part bythe National Science Foundation under grant AST-0909132;and in part by the U.S. DOE NNSA ASC through the ArgonneInstitute for Computing in Science under field work proposal57789. The software used in this work was developed in partby funding from the U.S. DOE NNSA-ASC and OS-OASCRto the Flash Center for Computational Science at the Universityof Chicago.APPENDIXSHOCK DETECTIONWe describe here our shock detection algorithm, which we used in the verification of the Biermann effect at shocks. Thisalgorithm has some benefits over shock-detection algorithms such as the ones described in Balsara & Spicer (1999); Miniati et al.(2000); Ryu et al. (2003), in that it furnishes an estimate of the shock speed directly from a single time slice of spatial data, weightsthe mass, momentum, and energy components of the Rankine-Hugoniot relations equally, yields a single-cell-wide shock interface,verifies that local characteristics converge on the shock surface, and does not impose arbitrary thresholds on physical quantitiessuch as compression ratios or Mach numbers.We start with a mesh of cells containing fluid values (for present purposes it is immaterial whether these are cell averages orpoint values). The algorithm seeks a set of “shocked” cells satisfying the following criteria: • The Shocked Surface Is One Cell Wide: Each shocked cell is the location of a maximum of |∇ P | in the direction of ∇ P ,where P is the total fluid pressure; • Fluid Quantities Change Correctly: The neighborhood of each shocked cell exhibits compression, acceleration, and heatingin the direction of ∇ P ; • Rankine-Hugoniot Conditions Obtain: The full Rankine-Hugoniot conditions on mass, momentum, and energy hold betweenthe fluid upstream and downstream of each shocked cell, where “upstream” means in the direction n ≡ −∇ P/ |∇ P | ; • Characteristics Converge On The Shock Surface: The characteristic convergence criterion holds at each shocked cell: theshock is supersonic upstream, and subsonic downstream.The result is an efficient and reliable algorithm that among other things yields an accurate estimate of the shock speed, whichis essential to our verification work on the Biermann effect. Note that we neglect the induction equation in the shock detectionalgorithm, and operate using only the material pressure – not the total pressure, which includes the magnetic pressure. For weaklymagnetized plasmas, such as the ones we consider in this paper, this creates no difficulty in identifying the correct shock surface.Some generalization would be required for significantly magnetized plasmas, particularly to the third and fourth conditions above.We now describe in somewhat greater detail the algorithm outlined above.
The Shocked Surface Is One Cell Wide
The first condition above is tantamount to insisting that the shocked cell be sited at a position where the gradient of P issteepest. In general, volume-based hydrodynamic advance schemes spread shocks over several cells. If we were not to imposethis condition, we would obtain a thickened shell of shocked cells, which would complicate the criterion (specified below) forclassifying neighboring cells as adjoining the shock upstream or downstream.The condition is very easy to enforce: in the immediate neighborhood of the cell being tested for “shockedness”, we check theadjoining cells closest to the direction of the normal vector n and its mirror image − n . This is done by forming the vector ofintegers ∆ i ≡ NINT [ n / max ( n , n , n )] (where the NINT function ascribes the nearest integer to each component of its realvector argument) and adding it to the cell discrete index vector i , to reach the cells at i ± ∆ i . If a stencil-based estimate of |∇ P | isgreater in the candidate cell than in the two so-chosen adjoining cells, then the condition is passed successfully. Fluid Quantities Change Correctly
An easy sanity check for “shockedness” of a cell is that cells downstream should unfailingly be (a) compressed, (b) accelerated(in the + n -direction), and (c) heated, with respect to cells upstream.We introduce a shock-width parameter h , such that the shock-adjoining cells “upstream” and “downstream” of the shocked cellare respectively located at ± h × ∆ × n relative to the cell under study, where ∆ is the length of a cell side. For FLASH with theHLL Riemann solver (Toro 2009) and 2nd-order reconstruction, an appropriate value of h is 1.5. We then simply verify that thefluid variables ρ , u , and e T (the total thermal energy) satisfy the conditions [ ρ ] du > , [ u · n ] du > , and [ e T ] du > . Rankine-Hugoniot Conditions Obtain
The Rankine-Hugoniot (R-H) conditions, which express flux conservation in the frame of the shock, have the form D × ( Φ d − Φ u ) = F ( Φ d ) − F ( Φ u ) . (A1)HE BIERMANN CATASTROPHE 21Here, Φ is the state vector of conserved field quantities, Φ T ≡ [ ρ, ρ u · n , ρ E ] , that is, the mass, normal momentum, and energydensities. F denotes the vector of fluxes corresponding to Φ . The subscripts “d” and “u” denote “downstream” and “upstream”states, respectively (Toro 2009).We fit R-H conditions to the upstream/downstream data defined using the shock-width parameter h above. In this fit, the shockspeed D is a free parameter to be adjusted to minimize a fit residual. Given a positive-definite inner product [Φ , Φ ] on the vectorstate space in Eq. (A1), we may define the normalized residual R as R ≡ [( D ∆ Φ − ∆ F ) , ( D ∆ Φ − ∆ F )][∆ F , ∆ F ] , (A2)where ∆ Φ ≡ Φ d − Φ u and ∆ F ≡ F ( Φ d ) − F ( Φ u ) . We may then easily minimize R with respect to D , obtaining D Min = [∆ Φ , ∆ F ][∆ Φ , ∆ Φ ] (A3) R Min = 1 − [∆ F , ∆ Φ ] [∆ F , ∆ F ] [∆ Φ , ∆ Φ ] . (A4)Given these definitions, we consider that the R-H conditions are satisfied if R Min is less than some small threshold value.It remains to define the inner product [ · , · ] used in these expressions. It is obvious that the naive unweighted sum-of-squaresof the components of vectors such as ∆ Φ is not an acceptable inner product, as it is dimensionally senseless. We require at aminimum a positive-definite metric that brings all vector components to the same physical dimensions, so we don’t wind up addingmass densities to energy densities, and so on. An additional desirable requirement is that all three components of Eq. (A1) shouldcontribute similar magnitudes to Eqs. (A3) and (A4), so that they are all weighted equally in the fit.A simple metric that accomplishes these requirements may be build using the local sound speed c s of the candidate shocked cell: [ Φ , Φ ] ≡ Φ T M Φ , (A5)where M ≡ (cid:34) c s −
00 0 c s − (cid:35) . (A6)This choice produces compatible dimensions and comparable magnitudes, because in and downstream of the shock we have c s ∼ ( k B T /m i ) / ∼ D , so that in the frame of the shock, the mass flux is ρD ∼ ρc s , the momentum flux has magnitude P ∼ ρc s , and the energy flux has magnitude ρD ∼ ρc s .We adopt this metric in computing Eqs. (A3) and (A4). We find that a threshold of R Min < . is adequate for identifyingcells satisfying the R-H conditions with few false-positives and false-negatives, in the simulations exhibited in this paper. Characteristics Converge On The Shock Surface
Our final criterion is easily stated and checked: shocks are supersonic upstream, and subsonic downstream, so that characteristicsfrom the family implicated in the shock converge on the shock (see Courant & Friedrichs 1948, pp. 141-146). This is an essentialstability criterion for the solution. Using the shock speed D Min calculated while fitting the R-H conditions to the data, thiscondition is simply u d · n + c s,d > D Min > u u · n + c s,u , (A7)where c s,d and c s,u are the downstream and upstream sound speeds, respectively, as determined with respect to the shock width h defined above. REFERENCESAmendt, P., Milovich, J., Wilks, S., Li, C., Petrasso, R., & S´eguin, F. 2009, Plasma Physics and Controlled Fusion, 51, 124048Balsara, D. S., & Spicer, D. S. 1999, Journal of Computational Physics, 149, 270Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65Bryan, G. L. et al. 2014, ApJS, 211, 19Cen, R., & Ostriker, J. P. 1994, ApJ, 429, 4Courant, R., & Friedrichs, K. 1948, Supersonic flow and shock waves, Pure and applied mathematics (Interscience Publishers)Davies, G., & Widrow, L. M. 2000, The Astrophysical Journal, 540, 755Drake, R. 2006, High-Energy-Density Physics: Fundamentals, Inertial Fusion, and Experimental Astrophysics, High pressure shock compression of condensedmatter (Springer)Dubey, A., Antypas, K., Ganapathy, M. K., Reid, L. B., Riley, K., Sheeler, D., Siegel, A., & Weide, K. 2009, Parallel Computing, 35, 512Fatenejad, M. et al. 2013, High Energy Density Physics, 9, 172Fryxell, B. et al. 2010, High Energy Density Physics, 6, 162Fryxell, B. et al. 2000, ApJS, 131, 273Gnedin, N. Y., Ferrara, A., & Zweibel, E. G. 2000, The Astrophysical Journal, 539, 505Gregori, G. et al. 2012, Nature, 481, 480Gurnett, D., & Bhattacharjee, A. 2005, Introduction to Plasma Physics: With Space and Laboratory Applications (Cambridge University Press)Huba, J. 2007, 2007 Plasma Formulary - Handbook of Data and Tables for Plasma Physics and Engineering (Wexford College Press)Jaffrin, M., & Probstein, R. 1964, Physics of Fluids, 7, 1658Kevlahan, N.-R. 1997, Journal of Fluid Mechanics, 341, 371Kronberg, P. P. 1994, Reports on Progress in Physics, 57, 3252 GRAZIANI ET AL.