Kengo TOMIDA (Princeton / JSPS Fellow)
James M. Stone (Princeton)• Introduction
 Star-Disk Interaction
 15 years ago and now• Code Development
 MHD in Spherical Polar Coordinates• (Preliminary) Results• SummarySo farGoalThis Work• Toward understanding the main accretion phase• Small scale structure around a protostar is crucial• Modeling the accretion flow interacting with a protostar
Questions:
1. How does the gas in the disk accrete onto the protostar?
2. How does the stellar spin evolve?
a. MRI-driven turbulence in the presence of global magnetic fields(shearing-box simulations cannot be realistic enough)b. Interaction between the accretion flow and stellar fields

3. How does the mass ejection process (outflow) occur?
4. How can transitional disks form?
a. Accretion?
b. Wind / Outflow?
5. Observations of T Tauri stars, especially variability

Because I’m now in Princeton, I want
• something new, different and scientifically significant,• and to contribute to Athena Project.

Considering today’s situation in computing...
• There will be no drastic revolution in computing.
• Node performance will not be improved significantly.
⇒ I need something extremely well scalable.
Star-Disk Interaction with Athena• Athena did not have spherical polar coordinates.
• Only one active group, but scientifically important• No self-gravity, No radiation transfer → scalable
富田、輻射やめるってよ ※やめませんcountour=Bφstreamlines=Bp

weak Bvery
strong B～15 years ago: 2D MHD (Hirose et al. 1997, Miller & Stone 1997)
MRI was not well developed because of short simulation time.
Strong magnetic braking may dominate the evolution.

Romanova et al. 2012

3D cubed-sphere MHD simulations by Romanova et al. @ Cornell
Weak fields: Boundary-Layer accretion through the disk
Strong fields: Magneto-Spheric accretion through the funnel-flow
Key: Interaction between stellar fields and accretion driven by MRI.

• New CT-MHD solver in Spherical Coordinates for Athena• PLM and PPM (Note: Blondin & Lufkin 93 is wrong)• Cell-averaged magnetic fields• Limitation: VL integrator, uniform mesh-spacing• Adiabatic EOS with γ =5/3
rdisk =• Point-source gravity
rcorotation• Star: fixed + density ceiling• Other boundaries: outflow• Cold disk + Hot Corona• Stellar dipole B + uniform Bz• Interior: rigid-body rotation• Disk / outer envelope: Keplerian

X-point

WeakModerate

A bit strongμ (stellar dipole)

0.32

0.32

0.32

RF (Fmag / Fgrav @ r=1)

0.3

0.3

0.3

Bz (uniform field)

0.001

0.001

0.001

rX (“X-point”)

6.84

6.84

6.84

rd (disk = corot.)

33

3

RT (Tcorona / Tdisk)

100

100100Θ d (cs/vK)

0.175 (10°)

0.175 (10°)

0.175 (10°)

ρd101

0.3

β i (plasma β at rd)～1400～140～40

Based on Miller & Stone 1997 (Models Ia-c) and Romanova+ 2002
Color: Independent (Red) / Dependent (Black)
Background: Variable (Yellow) / Fixed (None)
Note: this system is scale-freeIn ideal MHD, numerical resistivity is dominant dissipation.
“Turbulent reconnection” (Lazarian & Vishniac 1999, Kowal et al. 2009)
Dissipation in turbulent media does not depend on the input resistivity
and resolution if turbulence is well resolved. (*this must be tested*)

In order to resolve Magneto Rotational Instability,2𝜋𝑣𝑎 2𝜋𝐻𝜆𝑀𝑅𝐼 ==Ω𝛽at least 𝜆𝑀𝑅𝐼 / Δ𝑧 > 10, 𝜆𝑀𝑅𝐼 / 𝑅Δ𝜑 > 20 (Hawley+ 2011)
𝐻/Δ𝑥 ≳ 30 (depends on β ) (Latter+ 2010, Hawley+ 2013 ※net zero)
Romanova+ 2012: 𝜆𝑀𝑅𝐼 /Δ𝑧～ 8 (thin disk) but 𝜆𝑀𝑅𝐼 /𝑅Δ𝜑 is far lower
In this work, 𝐻/Δ𝑥 ≳ 20 (low resolution), 40 (test) at r ～ 5
Nr x Nθ x Nφ =800 x 320 x 640 (with 800 cores)(Nr x Nθ x Nφ =1440 x 640 x 1280 (with 4000 cores) looks converged.)
Tiger cluster @ Princeton, Cray XC30-A @ NAOJ, BlueGene/Q @ KEK

Flat MPI Parallelization, 2nd-order, VL, HLLD, sphericalWeak scaling: 643 per core, Measured on XC30 at NAOJ
Efficient scalability > 95% with 12,000 coresweakmoderatestrong

color: density

Weak @ t= 100orb

Moderate @ t=70 orbr< 8 only (because of the huge file size > 8GB / step)
100 orbits at the surface ～ 20 orbits at r=3 (corot.), 4.5 orbits at r=8

Omegaweak @ t～ 100 orb
moderate @ t～70 orbit
strong @ t～12 orbit

Stellar magnetic fields are notdynamically significant in the
weak and moderate models.
Magnetic braking may work inthe strong case, but longer
simulation is required.

WeakModerate

Saturation levels are similar: β ~ 1-100
MRI develops from the center and propagates outward
λ𝑀𝑅𝐼 ∝ 𝛽−1/2

WeakModerate
Magnetic Stress due to turbulent B-field is dominant in the disk,
but magnetic braking by global magnetic fields is not negligible
in the models with stronger fields.

Miller &Stone 1997

Weak
Moderate
Toward the pole, magnetic pressure driven winds are launched.
Layered outflows from the disks are observed, but they are notsignificant (very sparse), and not related to the original X-point.
Different B-field, rotation or turbulence suppresses it?

High-resolution MHD simulations of Star-Disk interaction• Weak-field cases are OK; direct accretion through boundary-layer• Strong-field cases need more work; time-consuming, less stable(but magnetospheric regime with strong field is more interesting...)• Longer-term simulations are required• Outflows in old 2D simulations may be due to the initial conditions• Time-dependent analyses are needed
Accretion rate, Stellar Spin, radial profile, “Butterfly diagram”, etc.
• Non-uniform mesh spacing and mesh (de-)refinement may help.

• Do NOT trust any code just because it’s public or famous.

Every physical quantity is defined by volume average.
(Magnetic fields in CT scheme are surface average.)
𝑎 𝑑𝑉𝑉𝑖,𝑗,𝑘𝑎 𝑖,𝑗,𝑘 =𝑑𝑉
𝑉𝑖,𝑗,𝑘𝜕𝑡 𝑸 𝑑𝑉 = −𝑭 ∙ 𝒏𝑑𝑆 +𝜕𝑡 𝑩 𝑑𝑆 = −𝑺 𝑑𝑉

𝜀 𝑑𝑙

Discretize these equations appropriately, using numericalfluxes calculated by an adequate Riemann solver (ex: HLLD)Piecewise Linear Method: linear profile within a cell + limiter𝑎𝑖 𝑟 = 𝑎𝐿,𝑖 + 𝐴 𝑟 − 𝑟𝑖− = 𝑎𝑅,𝑖 − 𝐴 𝑟𝑖+ − 𝑟
2𝐶𝐿 𝐶𝑅e.g. 𝐴 = 𝐿 𝐶𝐿 , 𝐶𝑅 =(van Leer limiter)𝑎 𝑖 − 𝑎 𝑖−1
𝐶𝐿 =𝑟𝑖 −𝑟𝑖−1
𝑎 𝑖 − 𝑎 𝑖−1⇒ 𝐶𝐿 = 𝑟 − 𝑟𝑖𝑖−1

𝐶𝐿 +𝐶𝑅(in Cartesian coordinates)(in Curvilinear coordinates)

Simple enough, but some codes do not take care of this.
(Note:1. This form is general.

2. PPM is worse; Blondin & Lufkin 1993 is wrong!)

Common recipe to evaluate B-fields at cell center:surface-weighted average∆𝑟𝑖 (∆𝑆𝑖− 𝐵𝑖− + ∆𝑆𝑖+ 𝐵𝑖+ )𝐵𝑟 𝑖 =∆𝑉𝑖
where ∆𝑆𝑖− = 𝑟𝑖− 2 , ∆𝑉𝑖 =𝑟𝑖+ 3 −𝑟𝑖− 33

This scheme may look good as it conserves div B better,but is actually not monotonic: ∆𝑉𝑖 ≠ ∆𝑟𝑖 (∆𝑆𝑖− + ∆𝑆𝑖+ ).
Option 1 (not recommended): ∆𝑉𝑖 = ∆𝑟𝑖 ∆𝑆𝑖− + ∆𝑆𝑖+Option 2: take volume-average of a simple linear profile(𝑟𝑖+ − 𝑟 𝑖 )𝐵𝑖− + ( 𝑟 𝑖 − 𝑟𝑖− )𝐵𝑖+𝐵𝑟 𝑖 =∆𝑟𝑖