Physical and Artificial Resistivity (in smoothed particle - PowerPoint PPT Presentation

Physical and Artificial Resistivity (in smoothed particle magnetohydrodynamics) James Wurster 1 st Phantom Users Workshop Monash University, 20 February 2018

Ideal magnetohydrodynamics d B = ( B · r ) v − B ( r · v ) d t 2

Ideal MHD Ø Fully ionised plasma Ø Zero resistivity & infinite conductivity Ø Ions & electrons are tied to the magnetic field 3

Ideal MHD 5000 z [AU] 0 -5000 -5000 0 5000 -5000 0 5000 x [AU] x [AU] Density (rendered) + Magnetic field lines 4 Ideal MHD. Left: Initial conditions. Right: at ρ max = 10 -9 g cm -3

Ideal MHD: Artificial Resistivity d B = ( B · r ) v − B ( r · v ) d t + r × η art ( r × B ) where η art ≈ 1 2 α B v sig h 5

Ideal MHD: Artificial Resistivity Ø Artificial resistivity (Tricco & Price, 2013) d B i + d B i � 1 h i X a v j a v i ab B j a r j a W ab ( h a ) � B i ab r j a � = a W ab ( h a ) � m b � d t Ω a ρ a d t � art b " # r j r j α B ab r j + α B ab r j d B i � a v sig ,a ˆ a W ab ( h a ) b v sig ,b ˆ a W ab ( h b ) ρ a X a � m b B i = � ab Ω b ρ 2 2 Ω a ρ 2 d t � a art b b v i v i a � v i = ab b B i B i a � B i = ab b q s ,a + v 2 = c 2 v sig ,a A ,a ✓ h a | r B a | ◆ α B = min , 1 a | B a | v 2 � � u ∂ B i tX X � a � u | r B a | ⌘ � � ∂ x j � � a i j 6 Ø Always applied if there is a gradient in the magnetic field (i.e. | ∇ B | > 0 )

Ideal MHD: Artificial Resistivity Ø Artificial resistivity (Price, et al, submitted) d B i + d B i � 1 h i X a v j a v i ab B j a r j a W ab ( h a ) � B i ab r j a � = a W ab ( h a ) � m b � d t Ω a ρ a d t � art b " # r j r j ab r j ab r j d B i � ˆ a W ab ( h a ) + ˆ a W ab ( h b ) ρ a X a � m b α B v sig ,ab B i = � ab Ω b ρ 2 2 Ω a ρ 2 d t � a b art b B i B i a � B i = ab b = | v ab ⇥ ˆ r ab | v sig ,ab α B 1 ⌘ Ø Always applied for non-zero velocity Ø Less resistive that that from Tricco & Price (2013) 7

Ideal MHD: Artificial Resistivity � density 0.5 Ø Price et. al. (2017) artificial resistivity t=0.5 t=1.0 � | v ab ⇥ ˆ = r ab | v sig ,ab α B 1 ⌘ 0.4 " � η v× r X � � Ø Tricco & Price (2013) � b q 0.3 s ,a + v 2 = c 2 v sig ,a A ,a ✓ h a | r B a | ◆ α B = min , 1 a | B a | η a + η b 0.2 v � � u Ø Tricco & Price (2013) with alternate averaging 0.1 8 Wurster, Bate, Price & Tricco (2017) η ab

Ideal MHD: Artificial Resistivity density 0.5 t=0.5 t=1.0 0.0013 Price et. al. (2017) η v X r Tricco & Price (2013) η a + η b η ab Tricco & Price (2013) w/alternate averaging 0.0012 0.0011 0.4 E mag 0.001 η v× r 0.0009 0.0008 0.3 0.0007 0 0.2 0.4 0.6 0.8 1 Time η a + η b 0.2 0.1 9 Wurster, Bate, Price & Tricco (2017) η ab

Motivation: Non-ideal MHD Orion Molecular Cloud HL Tau 10 Ionisation fraction ~ 10 -14 ~ 10 -12

Non-ideal MHD Ø Partially ionised plasma B ρ Ø Non-zero resistivity & conductivity Ø Ions, electrons & neutrals behaviour is environment-dependent Ohmic Resistivity Hall Effect Ambipolar Diffusion 11 (ion-electron drift) (ion-neutral drift)

Non-ideal MHD Ambipolar Diffusion (dissipative) Hall Effect log B (non-dissipative) Ohmic Resistivity (dissipative) log n 12 Adapted from Wardle (2007)

Non-ideal MHD: Hall effect 13 Image credit: Tsukamoto et al (2017); see also: Braiding & Wardle (2012a,b)

Ideal vs non-ideal MHD 5000 z [AU] 0 -5000 -5000 0 5000 -5000 0 5000 x [AU] x [AU] Density (rendered) + Magnetic field lines 14 During first core phase. Left: ideal MHD. Right: non-ideal MHD

Ideal vs non-ideal MHD 10 z [AU] 0 -10 -10 0 10 -10 0 10 x [AU] x [AU] Density (rendered) + Magnetic field lines 15 During first core phase. Left: ideal MHD. Right: non-ideal MHD

Non-ideal MHD log ρ n (g cm -3 ) ρ ρ -20 -15 -10 -5 0 η OR 25 10 5 First collapse First core Second Second core η HE > 0 Isothermal Adiabatic collapse Adiabatic η HE < 0 η AD η AD Duffin & Pudritz 2008 20 10 0 - η H η AD η H log η (cm 2 s -1 ) η (s) 15 10 -5 η Ω Machida et al. 2007 η Ω 10 10 -10 5 10 0 10 5 10 10 10 15 10 20 10 25 n H (cm -3 ) 0 5 10 15 20 log n n (cm -3 ) NICIL: Wurster (2016) Marchand+ (2016) 16 NICIL v1.2.3 is implemented in the current git version of Phantom

Non-ideal MHD in Phantom: the NICIL library Ø Phantom includes the NICIL code (Wurster 2016) Ø Publically available at https://bitbucket.org/jameswurster/nicil Ø When compiling, set NONIDEALMHD=yes Ø Realistic defaults are set; these will self-consistently calculate the non-ideal coefficients Ø Fully parameterisable Ø Primary parameters are included in Phantom ’s .in file Ø All parameters are included at the top of nicil.F90 Ø Important parameters that can be modified Ø Included non-ideal MHD terms (default = ohmic + Hall + ambipolar) Ø Ionisation source (default = cosmic rays + thermal) Ø Cosmic ray ionisation rate (default = 10 -17 s -1 ) Ø Elements that can be thermally ionised (cannot be modified through .in file) Ø Grain properties (default = fixed size of 0.1µm; alternate is MRN, but is slow) Ø Important values are summarised in the dump files and the .ev file Ø Can optionally preselect non-ideal MHD coefficients (preferably for tests only) Ø All coefficients and required variables are calculated at runtime 17 17

Implementation Ø Continuum equations d B = ( B · r ) v − B ( r · v ) d t + r × η art ( r × B ) + r × η OR ( r × B ) h i ( r × B ) × ˆ + r × η HE B nh i o ( r × B ) × ˆ × ˆ + r × η AD B B � Ø SPMHD equations + d B i d B i � � � 1 � a v j a � a − B i ab ∇ j v i ab B j a ∇ j a W ab ( h a ) = − m b a W ab ( h a ) � � � d t d t � a ρ a � non - ideal b ideal a b � � � d d � D a � � � � D b d B a � � + × ∇ a W ab ( h b ) , = − ρ a m b × ∇ a W ab ( h a ) � � b ρ 2 � a ρ 2 d t � b non - ideal a b J a × ˆ × ˆ D AD � � � D B a . D HE = − η HE J a × ˆ = η AD B a D OR B a , = − η OR J a , a a a 18 18 Wurster, Price & Ayliffe (2014)

Implementation Density Loop: do i = 1, N do j = 1, N neigh Using j , calculate density of i Using j , calculate current density, J = r × B , of i enddo Using new density of i , calculate η nimhd enddo Force Loop: do i = 1, N Calculate J i × B i and ( J i × B i ) × B j do j = 1, N neigh Calculate J j × B j and ( J j × B j ) × B i Using j , calculate d B / d t non-ideal of i enddo Calculate non-ideal timesteps enddo Step Loop: do i = 1, N Updated magnetic field of i , using ideal, non-ideal and artificial terms 19 19 enddo

Implementation Ø Timestepping: h = d t Courant C c 10 8 v sig ideal MHD non-ideal MHD 10 7 h 2 = d t nimhd C ni cpu time [cpu-hours] 10 6 | η | 10 5 10 4 Ø Phantom includes super-timestepping 10 3 (Alexiades, Amiez & Gremaud 1996) 10 2 Ø Right: cpu-hours required for the 10 6 particle 0 0.2 0.4 0.6 0.8 1 1.2 Simluation time [t ff ] models with µ 0 =5 in Wurster, Price & Bate (2016) Ø Non-ideal MHD is slightly slower for t < t ff , and much slower for t > t ff 20 20

Conclusions Ø Artificial resistivity is required to stabilised magnetohydrodynamics equations Ø Ideal MHD is a poor approximation for modelling molecular clouds or protoplanetary discs Ø Non-ideal MHD requires an assumption of chemistry Ø The non-ideal MHD coefficients are not dependent on neighbours Ø The non-ideal MHD contribution to the magnetic field evolution is dependent on neighbours Ø Non-ideal MHD introduces a diffusion timestep ∝ h 2 , hence can be computationally expensive j.wurster@exeter.ac.uk http://www.astro.ex.ac.uk/people/wurster/ Presentation available at http://www.astro.ex.ac.uk/people/wurster/files/spmhd_resistivity.pdf 21 Nicil ’s git repository: https://bitbucket.org/jameswurster/nicil