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

Physical and Artificial Resistivity

(in smoothed particle magnetohydrodynamics) James Wurster

1st Phantom Users Workshop Monash University, 20 February 2018

Ideal magnetohydrodynamics

2

dB dt = (B · r) v − B (r · v)

Ideal MHD

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

3

Ideal MHD

4

z [AU] x [AU]

• 5000

5000

• 5000

5000 x [AU]

• 5000

5000

Density (rendered) + Magnetic field lines Ideal MHD. Left: Initial conditions. Right: at ρmax = 10-9g cm-3

Ideal MHD: Artificial Resistivity

5

dB dt = (B · r) v − B (r · v) + r × ηart (r × B)

ηart ≈ 1 2αBvsigh

where

Ideal MHD: Artificial Resistivity

6

Ø Artificial resistivity (Tricco & Price, 2013) Ø Always applied if there is a gradient in the magnetic field (i.e. |∇B | > 0 )

dBi

a

dt =

• 1

Ωaρa X

b

mb h vi

abBj arj aWab (ha) Bi avj abrj aWab (ha)

i + dBi

a

dt

• art

dBi

a

dt

• art

= ρa 2 X

b

mbBi

ab

" αB

avsig,aˆ

rj

abrj aWab(ha)

Ωaρ2

a

+ αB

b vsig,bˆ

rj

abrj aWab(hb)

Ωbρ2

b

# vi

ab

= vi

a vi b

Bi

ab

= Bi

a Bi b

vsig,a = q c2

s,a + v2 A,a

αB

a

= min ✓ha |rBa| |Ba| , 1 ◆ |rBa| ⌘ v u u tX

i

X

j

• ∂Bi

a

∂xj

a

• 2

dBi

a

dt =

• 1

Ωaρa X

b

mb h vi

abBj arj aWab (ha) Bi avj abrj aWab (ha)

i + dBi

a

dt

• art

dBi

a

dt

• art

= ρa 2 X

b

mbαBvsig,abBi

ab

" ˆ rj

abrj aWab(ha)

Ωaρ2

a

+ ˆ rj

abrj aWab(hb)

Ωbρ2

b

# Bi

ab

= Bi

a Bi b

vsig,ab = |vab ⇥ ˆ rab| αB ⌘ 1

Ideal MHD: Artificial Resistivity

7

Ø Artificial resistivity (Price, et al, submitted) Ø Always applied for non-zero velocity Ø Less resistive that that from Tricco & Price (2013)

Ideal MHD: Artificial Resistivity

8 ηv× r t=0.5 t=1.0 ηa+ηb ηab 0.1 0.2 0.3 0.4 0.5 density

Ø Price et. al. (2017) artificial resistivity Ø Tricco & Price (2013) Ø Tricco & Price (2013) with alternate averaging Wurster, Bate, Price & Tricco (2017)

• X

b

" vsig,a = q c2

s,a + v2 A,a

αB

a

= min ✓ha |rBa| |Ba| , 1 ◆ v u

• vsig,ab

= |vab ⇥ ˆ rab| αB ⌘ 1

Ideal MHD: Artificial Resistivity

9 ηv× r t=0.5 t=1.0 ηa+ηb ηab 0.1 0.2 0.3 0.4 0.5 density

0.0007 0.0008 0.0009 0.001 0.0011 0.0012 0.0013 0.2 0.4 0.6 0.8 1 Emag Time ηv X r ηa + ηb ηab

Price et. al. (2017) Tricco & Price (2013) Tricco & Price (2013) w/alternate averaging

Wurster, Bate, Price & Tricco (2017)

Motivation: Non-ideal MHD

10

Orion Molecular Cloud HL Tau Ionisation fraction ~ 10-14 ~ 10-12

Non-ideal MHD

ØPartially ionised plasma ØNon-zero resistivity & conductivity ØIons, electrons & neutrals behaviour is environment-dependent

11

B ρ Ohmic Resistivity Hall Effect Ambipolar Diffusion (ion-electron drift) (ion-neutral drift)

Non-ideal MHD

Adapted from Wardle (2007)

log n log B

Ambipolar Diffusion (dissipative) Hall Effect (non-dissipative) Ohmic Resistivity (dissipative)

12

Non-ideal MHD: Hall effect

Image credit: Tsukamoto et al (2017); see also: Braiding & Wardle (2012a,b)

13

z [AU] x [AU]

• 5000

5000

• 5000

5000 x [AU]

• 5000

5000

Ideal vs non-ideal MHD

14

Density (rendered) + Magnetic field lines During first core phase. Left: ideal MHD. Right: non-ideal MHD

Ideal vs non-ideal MHD

15

Density (rendered) + Magnetic field lines During first core phase. Left: ideal MHD. Right: non-ideal MHD

z [AU] x [AU]

• 10

10

• 10

10 x [AU]

• 10

10

ρ ρ 5 10 15 20 25 5 10 15 20

• 20
• 15
• 10
• 5

log η (cm2 s-1) log nn(cm-3) log ρn (g cm-3) ηOR ηHE > 0 ηHE < 0 ηAD 16

Non-ideal MHD

NICIL: Wurster (2016) Marchand+ (2016) NICIL v1.2.3 is implemented in the current git version of Phantom

10-10 10-5 100 105 100 105 1010 1015 1020 1025

η (s) nH (cm-3)

First collapse Isothermal First core Adiabatic Second collapse Second core Adiabatic

ηAD ηΩ

• ηH

ηH

ηAD Duffin & Pudritz 2008 ηΩ Machida et al. 2007

Non-ideal MHD in Phantom: the NICIL library

17 17

Ø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

Implementation

18 18

ØContinuum equations ØSPMHD equations

dBi

a

dt = − 1 aρa

• b

mb

• vi

abBj a ∇j a Wab (ha)

• d
• d
• −Bi

avj ab∇j a Wab (ha)

• + dBi

a

dt

• non-ideal

DOR

a

= −ηOR Ja,

dBa dt

• non-ideal

= −ρa

• b

mb Da aρ2

a

× ∇aWab(ha) D

• ideal
• b
• a

+ Db bρ2

b

× ∇aWab(hb)

• ,

Wurster, Price & Ayliffe (2014)

DHE

a

= −ηHE Ja × ˆ Ba, DAD

a

= ηAD

• Ja × ˆ

Ba

• × ˆ

Ba.

dB dt = (B · r) v − B (r · v) + r × ηart (r × B) + r × ηOR (r × B) + r × ηHE h (r × B) × ˆ B i + r × ηAD nh (r × B) × ˆ B i × ˆ B

Implementation

19 19

Density Loop: do i = 1,N do j = 1,Nneigh 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 Ji × Bi and (Ji × Bi) × Bj do j = 1,Nneigh Calculate Jj × Bj and (Jj × Bj) × Bi Using j, calculate dB/dtnon-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 enddo

ØTimestepping: ØPhantom includes super-timestepping (Alexiades, Amiez & Gremaud 1996) ØRight: cpu-hours required for the 106 particle models with µ0=5 in Wurster, Price & Bate (2016) ØNon-ideal MHD is slightly slower for t < tff, and much slower for t > tff

102 103 104 105 106 107 108 0.2 0.4 0.6 0.8 1 1.2 cpu time [cpu-hours] Simluation time [tff] ideal MHD non-ideal MHD

Implementation

20 20

dtCourant = Cc h vsig dtnimhd = Cni h2 |η|

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 ∝h2, hence can be computationally expensive

21

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 Nicil’s git repository: https://bitbucket.org/jameswurster/nicil