Time Dependent Models of Magnetospheric Accretion onto Young Stars

TL;DR

This paper uses one-dimensional time-dependent simulations to show how smooth density variations in the inner disk can create shocks in the accretion column, explaining short-term variability in young stars.

Contribution

It demonstrates that smooth, day-scale density variations can induce shocks in magnetospheric accretion flows, linking simple models to observed stellar variability.

Findings

Shocks form in the accretion column due to density fluctuations.

Shocks cause rapid increases in accretion luminosity.

The model explains short-term variability in young stars.

Abstract

Accretion onto Classical T Tauri stars is thought to take place through the action of magnetospheric processes, with gas in the inner disk being channeled onto the star's surface by the stellar magnetic field lines. Young stars are known to accrete material in a time-variable manner and the source of this variability remains an open problem, particularly on the shortest (~ day) timescales. Using one-dimensional time-dependent numerical simulations that follow the field line geometry, we find that for plausibly realistic young stars, steady-state transonic accretion occurs naturally in the absence of any other source of variability. However, we show that if the density in the inner disk varies smoothly in time with ~ day long time-scales (e.g., due to turbulence) this complication can lead to the development of shocks in the accretion column. These shocks propagate along the accretion column and ultimately hit the star, leading to rapid, large amplitude changes in the accretion rate. We argue that when these shocks hit the star the observed time-dependence will be a rapid increase in accretion luminosity followed by a slower decline and could be an explanation for some of the short period variability observed in accreting young stars. Our one-dimensional approach bridges previous analytic work to more complicated, multi-dimensional simulations, and observations.

Equations

$$\begin{array}{r}\mathbf{B}=\frac{B_{\text{mbox}oct}}{2}{\xi }^{-5}[(5{\cos }^2\theta -3)\cos \theta \hat{\mathbf{r}}\\ +\frac{3}{4}(5{\cos }^2\theta -1)\sin \theta \hat{\mathbf{\theta }}]\\ +\frac{B_{\text{mbox}dip}}{2}{\xi }^{-3}(2\cos \theta \hat{\mathbf{r}}+\sin \theta \hat{\mathbf{\theta }})\end{array}$$

$$\Gamma \equiv \frac{B_{\text{mbox}oct}}{B_{\text{mbox}dip}}.$$

$$p=-\frac{1}{4}{\xi }^{-4}\Gamma (5{\cos }^2\theta -3)\cos \theta -{\xi }^{-2}\cos \theta ,$$

$$q=\frac{1}{4}{\xi }^{-3}\Gamma (5{\cos }^2\theta -1){\sin }^2\theta +{\xi }^{-1}{\sin }^2\theta .$$

$$\begin{array}{r}\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0\\ \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\frac{1}{\rho }\nabla P-\nabla \Psi -\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathbf{r})\\ \frac{\partial }{\partial t}(\frac{e}{\rho })+(\mathbf{u}\cdot \nabla )(\frac{e}{\rho })=-\frac{P}{\rho }\nabla \cdot \mathbf{u}\end{array}$$

$$P=(\gamma -1)e$$

$$a=\sqrt{\gamma \frac{k_{b}T}{m_H}}$$

$$\alpha \equiv \frac{\rho }{{\rho }_1}u\equiv \frac{u_p}{a_1}ϵ\equiv \frac{e}{a_1^2{\rho }_1}$$

$$\frac{\partial \alpha }{\partial t}$$

$$\frac{\partial u}{\partial t}$$

$$\displaystyle\frac{\partial}{\partial t}\Big{(}\frac{\epsilon}{\alpha}\Big{)}$$

$${\dot{M}}_{global}=f(4\pi R_{\ast }^2)\frac{\dot{M}}{A}$$

$${\rho }_{disk}(t)={\rho }_1[A(\sin (\frac{2\pi }{P}t)+1)+1].$$

$${\rho }_{disk}(t)={\rho }_1(1+A{e}^{-\frac{t^2}{2{\sigma }^2}})$$

$$V_{ff}=\bigg{(}\frac{2GM}{R_{*}}\bigg{)}^{1/2}\bigg{(}1-\frac{R_{*}}{R_{disk}}\bigg{)}^{1/2}$$

$$u^{n+1}_{i}=u^{n}_{i}-\frac{\Delta t}{h^{a}_{p,i}}\frac{2}{(\alpha^{n}_{i}+\alpha^{n}_{i+1})}\bigg{(}\frac{P^{n}_{i}-P^{n}_{i+1}}{p^{b}_{i}-p^{b}_{i+1}}\bigg{)}-\frac{\Delta t}{h^{a}_{p,i}}\bigg{(}\frac{b}{(\xi_{i}^{a})^{2}}\bigg{)}\bigg{(}{(\xi_{i}^{a}})^{-5}\frac{f^{a}_{i}\cos{\theta^{a}_{i}}}{H}\bigg{)}+\omega\Delta t\xi^{a}_{i}\sin{\theta^{a}_{i}}(\hat{x}\cdot\hat{p}^{a}_{i}).$$

$$\omega=\big{(}\frac{\Omega R_{*}}{a_{1}}\big{)}^{2}.$$

$$\hat{x}\cdot \hat{p}=\frac{{\xi }^{-5}\sin \theta \cos \theta (f+g)}{\sqrt{H}}$$

$$\begin{array}{r}\frac{\partial h_{x,i}^a}{\partial p}=\frac{h_{x,i}^b-h_{x,i-1}^b}{p_i^b-p_{i-1}^b}\\ \frac{\partial h_{x,i}^b}{\partial p}=\frac{h_{x,i}^a-h_{x,i-1}^a}{p_i^a-p_{i-1}^a}.\end{array}$$

$$h_p={\xi }^5[f^2{\cos }^2\theta +g^2{\sin }^2\theta {]}^{-1/2}$$

$$h_q=\frac{{\xi }^4}{\sin \theta }[f^2{\cos }^2\theta +g^2{\sin }^2\theta {]}^{-1/2}$$

$$h_{\varphi }=\xi \sin \theta$$

$$\begin{array}{r}f\equiv \Gamma (5{\cos }^2\theta -3)+2{\xi }^2\\ g\equiv \frac{3}{4}\Gamma (5{\cos }^2\theta -1)+{\xi }^2\end{array}$$

$$H\equiv {\left|\nabla p\right|}^2={\xi }^{-10}[f^2{\cos }^2(\theta )+g^2{\sin }^2(\theta )].$$

$$u^2=-W[-{\lambda }^2(q\xi {)}^{-6}(4-3q\xi )\exp (-2\epsilon -\frac{2b}{\xi }-\omega q{\xi }^3)]$$

$$\epsilon =\frac{1}{2}{\lambda }^2-\frac{3}{2}bq.$$

$$\ln \lambda -\frac{{\lambda }^2}{2}=3\ln q{\xi }_s-\frac{1}{2}\ln (4-3q{\xi }_s)-\frac{1}{2}+b(-\frac{3}{2}q+\frac{1}{{\xi }_s}+\frac{1}{2}{q}^4{\xi }_s^3)$$

$${\lambda }^2=-W[-\exp (-6\ln q{\xi }_s-\ln (4-3q{\xi }_s)-1+b(3q+\frac{2}{{\xi }_s}+q^4{\xi }_s))]$$

$$3{\xi }_s(\frac{8-5q{\xi }_s}{4-3q{\xi }_s})+3\omega q{\xi }_s^4=2b.$$

Full text

††slugcomment: Accepted for Publication in ApJ, March 13th, 2017

1. INTRODUCTION

Excess IR emission in classical T Tauri stars (CTTS) is attributed to optically thick circumstellar disks composed of dust and gas that form alongside the star from the collapsing natal molecular cloud (e.g. Mendoza V.,, 1966; Shu et al.,, 1987). Typical full (or primordial) circumstellar dust disks have inner radii at the dust sublimation radius, roughly 0.05-0.5 AU (Natta et al.,, 2001; Muzerolle et al.,, 2003), and extend out to tens or hundreds of AU (e.g. O’dell and Wen,, 1994; Isella et al.,, 2009; Andrews et al.,, 2009). The gaseous portion of the disk is thought to extend further in to where the stellar magnetic fields are strong enough to disrupt the disk (e.g., Koenigl,, 1991; Livio and Pringle,, 1992; Hartmann,, 1998; Muzerolle et al.,, 2003). This occurs roughly at the co-rotation radius, where the Keplerian rotational velocity is equal to the stellar rotation rate (typically 5-10 stellar radii). Although polarization measurements and modeling of T Tauri stars have shown the magnetic field is complicated near the photosphere (Valenti and Johns-Krull,, 2004; Donati et al.,, 2011; Gregory and Donati,, 2011), the large-scale magnetic structure of the star is dominated by its lowest order multipole components.

In addition to excess IR emission, CTTS exhibit UV luminosities above photospheric levels from shocks caused by material accreting onto the star. (e.g., Gullbring et al.,, 1998; Muzerolle et al.,, 1998, 2001, 2003; Ingleby et al.,, 2011, 2015). In the canonical magnetospheric accretion scenario, disk material travels from the magnetic truncation radius, along the magnetic field lines, and falls at super sonic free-fall speeds on stellar surface (Koenigl,, 1991; Calvet and Hartmann,, 1992; Shu et al.,, 1994; Hartmann et al.,, 1994). Since the material is supersonic, it shocks as it collides with the star, and heats the gas to temperatures of order $\sim 10^{5}-10^{6}$ K, which in turn leads to UV and soft X-ray emission (Calvet and Gullbring,, 1998). The accretion rate onto the star can be measured by characterizing and modeling this excess emission (Gullbring et al.,, 1998; Muzerolle et al.,, 1998; Ardila et al.,, 2013).

CTTS are known to be variable, hosting stochastic and quasi-periodic changes in luminosity that occur over a vast range of timescales. Measurements exhibiting variability have been taken across the electromagnetic spectrum, e.g., bursts of X-rays and increased UV luminosity as a result of rapid increases in the accretion rate (Ingleby et al.,, 2015), periodic decreases in optical flux explained as warps in the disk variably blocking the star and starspots (Bouvier et al.,, 1999; Alencar et al.,, 2010), as well as ‘seesaw’ variability in the infrared due to changes of the inner disk wall height (Muzerolle et al.,, 2009; Espaillat et al.,, 2011). Furthermore, equivalent widths and shapes of spectral lines have been observed to change (e.g., Johns and Basri,, 1995). Longer period changes, on the timescales of months and years, have been observed in CTTS and explained as variable star spot coverage and perhaps hints of stellar activity cycles analogous to the solar cycle (Cohen et al.,, 2004). All of the aforementioned processes can occur simultaneously, making it difficult to isolate the driving force for a given source of variability.

In this paper, we focus on the innermost region of the disk and variability due to changes in the accretion rate that occur on day to week timescales. This type of variability is faster than what is expected from viscous processes since the viscous timescale in the inner disk is on the order of several years (Hartmann,, 1998). Proposed sources of variability on these timescales include Rayleigh-Taylor instabilities that produce variable “tongues” of material that flow between magnetic field lines (Kulkarni and Romanova,, 2008) and magnetorotational instabilities that can cause changes in the amount of material available for accretion as a function of time (Romanova et al.,, 2012). A monitoring campaign of NGC 2264 observed multiple accretion signatures that appear to fit these models well (Stauffer et al.,, 2014), providing motivation for continued modeling of the inner regions of the disk.

The hydrodynamics of gas traveling along magnetic accretion columns has been studied analytically, assuming steady-state flows (e.g., Hartmann et al.,, 1994; Li,, 1996; Koldoba et al.,, 2002). Under the assumption that the accretion columns are well described by a polytropic equation of state ($P\propto\rho^{1+1/n}$), analytic work suggests a transonic steady state solution that approaches free-fall velocities may not be possible within the accretion columns if the effective polytropic index, $n$ is smaller than a critical value (Adams and Gregory,, 2012). These authors developed a coordinate system that follows magnetic field lines aligned with the stellar rotation axis under the assumption that the magnetic field is strong enough to dominate the flow. This allows for the transformation of the system to a one-dimensional (1D) problem with semi-analytic steady state solutions for the isothermal case where $n\rightarrow\infty$ The general form for the predicted constraint on the minimum value of $n$ that allows these solutions to exist is $n>\ell+3/2$, where $\ell$ is the highest order multipole component of the magnetic field near the star. This constraint only holds in the inner limit where $r\rightarrow 0$. The predicted absence of this type of solution under certain conditions may imply the possibility of spontaneous variability generated within accretion columns, even in the case of a constant accretion rate feeding the columns. However, this prediction cannot be tested analytically, and we must appeal to time-dependant numerical calculations.

Expanding upon this analytic work, we present a 1D hydrodynamic simulation of accretion onto a CTTS using a coordinate construction that follows the accretion columns from the gas disk to the star. We investigate two possible sources of short timescale accretion variability: inherent variability from the inability to reach free-fall speeds suggested by Adams and Gregory, (2012) and variability driven by turbulence/instabilities at the inner disk edge. In this work, we also study the hydrodynamics within the accretion column and derive possible observable quantities.

This paper is structured as follows. The construction of the magnetospheric accretion system and our numerical method used to simulate the system are presented in §2 and §3. We then present simulations testing for spontaneous variability and discuss their implications in §4. Explicit time dependence is then included, the results of which are shown in §5. We discuss implications of our simulation on the current understanding of variability within the inner disk and possibilities for observations in §6 and conclude in §7.

2. MODEL CONSTRUCTION

The magnetic fields of CTTS are strong near the star, in some cases reaching up to $6\mbox{kG}$ at the photosphere (Valenti and Johns-Krull,, 2004), and as such, we assume that the material in the accretion columns follow the stellar magnetic field lines, which are unperturbed by the accretion flow and co-rotate with the star (see Koldoba et al.,, 2002; Adams and Gregory,, 2012, for a discussion of the validity of this approximation). The magnetic truncation radius, and thus the inner edge of the disk, is assumed to be co-located with the co-rotation radius for simplicity. The large scale magnetic field of CTTS is dominated by the dipole component, but some CTTS have significant octupole components (Donati et al.,, 2011). In this work, we investigate two field geometries: a pure dipole field and a field formed from the sum of a dipole and an octupole component. Finally, we take the magnetic field to be constant in time and its axis to be aligned with the rotation axis.

This assumption does prevent the magnetic field lines from becoming twisted from differential rotation, which is an important part of angular momentum transfer within the inner disk (Ustyugova et al.,, 2006) and could perhaps be a source of variability itself. In regards to the transfer of angular momentum, the timescale in which the rotation period of the star changes ($\sim 10^{6}$ years) (Ustyugova et al.,, 2006) is many orders of magnitude longer than the roughly week long timescales this paper focuses on. Thus, we assume the stellar rotation rate is constant.

2.1. Coordinate System

Under the assumptions outliend above, the system can be reduced to a 1D problem using a set of coordinates that follows current-free, axis-symmetric magnetic field lines that co-rotate with the star. Because the magnetic field configuration is current free, and thus curl-free, it can be described as the gradient of scalar fields, or coordinates, denoted here as $(p,q)$. The coordinate $p$ describes the distance along a given field line, while $q$ describes which field line is being traced. These coordinates are orthogonal; $\nabla p$ points along the field line and $\nabla q$ points in the poloidal direction, perpendicular to $\nabla p$. The origin of the system is located at the center of the star. A complete formulation of this coordinate system and its associated scale factors can be seen in Adams and Gregory, (2012). For our numerical calculations, a grid was constructed along the magnetic field lines (the $p$ coordinate). A magnetic field containing dipole and octupole contributions is given by

[TABLE]

where $\xi$ the radius from the center of the star given in dimensionless stellar radius units, $\xi\equiv r/R_{*}$ and $\theta$ the colatitude. $B_{dip}$ and $B_{oct}$ are the strengths of the dipole and octupole terms in the multipole expansion of the magnetic field at the stellar surface. The parameter $\Gamma$ is defined as the ratio of the coefficients of the composite multipole field,

[TABLE]

The coordinates $q$ and $p$ can be written in terms of the standard spherical polar coordinates. The coordinate $p$, which traces along individual magnetic field lines, can be written as

[TABLE]

while the orthogonal coordinate $q$ is given by

[TABLE]

Two configurations of magnetic fields were used in our simulation, $\Gamma=0$ and $\Gamma=10$. Examples of the field lines for the different geometries are shown in Figure 1. The addition of the octupole components has a substantial impact on the regions of the accretion column at small radii. Near the surface of the star, the cross-sectional area of the accretion column is greatly reduced in the dipole plus octupole case compared to the pure dipole configuration, which will lead to increased levels of compression of the flow. Near the surface of the disk, the effect of adding octupole terms is minimal due to the steep radial fall off of the octupole component, causing the truncation radius to be set almost purely by the strength of the dipole field alone (Gregory et al.,, 2016).

2.2. Hydrodynamic Formulation

As the magnetic field is assumed to be strong enough to force material to only flow along field lines, the magnetic term in the force equation can be eliminated. The angular velocity vector is defined such that $\mathbf{\Omega}=\Omega\hat{z}$, and is aligned with the dipole and octupole moments. Transforming into a non-inertial frame introduces centrifugal force terms. Because the fluid is forced to flow solely along the field lines with no $\hat{\phi}$ components in this frame, the Coriolis term can be eliminated in the 1D force equation. After these simplifications, the continuity, force, and energy equations take the form

[TABLE]

where $\rho$ is the density, $e$ is the internal energy, $\mathbf{u}$ is the fluid velocity $\Psi$ is the gravitational potential and $P$ is pressure. Pressure, written in terms of an effective adiabatic index $\gamma=(n+1)/n$ and the ideal gas law, takes the form

[TABLE]

The sound speed $a$, derived from the adiabatic equation of state and assuming an ideal hydrogen gas, is written as

[TABLE]

where $k_{b}$ is the Boltzmann constant, and $m_{H}$ is the mass of a hydrogen atom, and $T$ is the gas temperature.

The governing equations can be written in dimensionless 1D forms along the co-rotating field line. Reference quantities were used to scale $\rho$, $e$, and the velocity along the field line, $u_{p}$, written as

[TABLE]

where $a_{1}$ and $\rho_{1}$ are the sound speed and density measured at the inner edge of the disk. Distances have been written in terms of the stellar radius ($\xi$) and time has been scaled using the stellar crossing time for a velocity $a_{1}$. This allows us to write the governing equations as

[TABLE]

where $\omega$ is a dimensionless parameter that measures stellar rotation, $\omega=\Omega R_{*}/a_{1}$. The gravitational potential is represented using the dimensionless quantity $b$, which measures the depth of the gravitational potential well relative to the sound speed, written as $b=GM_{*}/R_{*}a_{1}^{2}$. More discussion of this dimensionless formulation of this system is shown in (Adams and Gregory,, 2012). The curvilinear coordinate scale factor for the $p$ coordinate is denoted $h_{p}$, while $f$ and $H$ are ancillary functions of $\xi$ and $\theta$. The functional forms for these terms, along with an analytic expression for $(\hat{x}\cdot\hat{p})$ are shown in the appendix.

3. NUMERICAL METHOD

The hydrodynamic problem is solved using the method of finite differencing under a zeus-style framework using a time-explicit, multi-step solution procedure (see Stone and Norman, 1992). This style of simulation uses a staggered mesh framework in which the density, internal energy and gravitational potential are zone centered while the velocity is face centered. The grid is constructed with ‘ghost’ cells at both ends of the simulation domain, which are used to set boundary conditions. The velocity, density, and internal energy are updated using operator splitting in two sub-steps: the ‘source’ update and the ‘transport’ update.

Our numerical method is based on Owen and Adams, (2016), who simulated magnetically controlled isothermal outflows from hot Jupiters. In this work, the code has been modified to include the energy equation in the system of hydrodynamic equations. In addition, we have modified the geometry of the problem from a dipole plus vertical field to either a pure dipole or a dipole plus octupole magnetic field configuration and added the effects of rotation. More discussion of the changes can be found in the Appendix.

3.1. Sub-Steps

The first step is the source update, in which the velocity and the internal energy are updated by solving finite-difference versions of the governing differential fluid equations. The velocity is updated during the source update due to gravitational forces and pressure gradients along the accretion column. Because the system is in a non-inertial rotating reference frame, the velocity is also updated due to the centrifugal terms that appear in the force equation. The full difference equations are presented in the Appendix. In order to capture shocks an artificial viscosity is introduced using a full stress tensor formalism. The full stress tensor is need due to our complicated geometry. The full details of this artificial viscosity tensor are shown in Appendix A of Owen and Adams, (2016). Finally, the internal energy of the system is updated due to compressional heating in a time-centered manner (Stone and Norman,, 1992).

The second step is the transport update in which finite-difference approximations to the integral advection equations are used to update the velocity, density and internal energy. Reconstruction at the cell boundaries is performed in a second order manner with a van Leer limiter (see van Leer,, 1977).

3.2. Boundary Conditions

At the stellar boundary we apply standard outflow boundary conditions (see Stone and Norman,, 1992). At the disk boundary we used two different approaches. For the steady simulations described in §4 the density and internal energy were held fixed, while for the simulations with explicit time variability described in §5 the density was varied according to a prescribed time-dependent function and the internal energy was varied such that the ghost cells maintained constant entropy.

3.3. Simulation Parameters

The simulations were performed in a dimensionless manner, however, we present the results here in dimensional form in order to compare to real systems. Since the flows are scale-free, the specific values of the density and temperature do not matter inside the simulation.

The material at the foot-point in the disk is assumed to have a density of $3\times 10^{11}~{}\mbox{cm}^{-3}$, this is taken from the analytic accretion solutions of Adams and Gregory, (2012) for a standard accretion rate of $1\times 10^{-8}$ M*☉* yr*-1* and temperature of $10,000~{}\mbox{K}$ due to heating from UV radiation. Because the cross-sectional area of the accretion column is a function of chosen multipole geometry, the ‘standardized $\dot{M}$’ values reported in this work are mass fluxes at the star scaled to the equivalent mass flux through an area of one $\mbox{cm}^{2}$ at the surface of the disk. Mass accretion rates are scaled this way to facilitate comparisons between the two simulated geometries which would otherwise be drawing material from unequal sized magnetic footprints. Converting to the more observationally interesting quantities of the mass flux at the surface of the star can be done by dividing the reported scaled mass flux by the dimensionless cross section of the accretion column, $A$ (see Figure 1). For the two geometries discussed in this paper, $\Gamma=0$ and $\Gamma=10$, the ratio of the cross-sectional area of the column at the stellar surface to the disk is $5.14\times 10^{-4}$ and $4.18\times 10^{-5}$ respectively. The global stellar accretion rate in terms of a hot spot filling factor, $f$, and this standardized mass flux can be written as

[TABLE]

For reference, a young accreting star with $f=0.03$, $R=1.5R_{\odot}$ with a mass flux of $\dot{M}=1.0\times 10^{-7}\mbox{g\,cm}^{-2}$ under the pure dipole construction will have an accretion rate of $1.3\times 10^{-9}M_{\odot}\mbox{yr}^{-1}$ or $1.5\times 10^{-8}M_{\odot}\mbox{yr}^{-1}$ under the dipole plus octupole configuration. Both accretion rates are plausible for T Tauri stars (Ingleby et al.,, 2013).

We can determine the degree to which the flow is magnetically controlled by considering the magnetic plasma $\beta$ – the ratio of the gas pressure to the magnetic pressure. For a CTTS with a surface magnetic field strength $B_{\rm dip}$ = 3kG (e.g., Johns-Krull,, 2007) and an inner disk edge of 10 $R_{\ast}$, the parameter $\beta\ll 1$ several scale heights up in the disc where the accretion flow is launched. This finding validates our assumption that the magnetic field is unperturbed by the flow. In the midplane of disc, however, the density is larger and the magnetic plasma $\beta$ can approach unity. As result, MHD effects could be important in the midplane and hence in setting the location of the truncation radius.

The simulation domain is separated into to 1024 cells (resolution tests are discussed below). We fix the ghost cell used to set boundary conditions at the disk to a height of $0.05R_{*}$ above the mid-plane of the disk. Because gravitational potential and velocity gradients are larger near the star, a non-linear grid spacing is adopted to ensure the accuracy of the simulation. The grid spacing is scaled using a power law such that the resolution near the star is higher than near the disk.

The simulation was initialized with a very small negative velocity. Density throughout the simulation is initially set to the density at the disk. The internal energy of the simulation, $e$, is set such that the entropy is constant throughout the simulation. Transient behavior from initialization dissipates over the length of several stellar sound crossing times and the choice of initial flow profile is not important. Time steps were chosen within the simulation by enforcing a maximum Courant number of 0.5. Additionally, the time-step was not allowed to increase by more than 30% between subsequent steps.

4. RESULTS: SEARCHING FOR SPONTANEOUS VARIABILITY

As introduced above, previous analytic work suggested that a constraint can be placed on the index $n$ for steady-state solutions that pass smoothly through the sonic point and reach nearly free-fall speeds near the star (Adams and Gregory,, 2012). If $n$ does not satisfy this constraint, the system may be unable to reach this type of solution. This effect could generate spontaneous variability within the accretion column, and could therby explain some of the short-term variability observed in TTS. This variability would be predicted to occur when the polytropic index $n$ is sufficiently small (corresponding to a stiff equation of state), i.e., $n<\ell+3/2$, where $\ell$ is the highest order multipole with significant contributions to the magnetic field and in the inner limit where $r\rightarrow 0$. For the case where the stellar magnetic field is a pure dipole ($\ell=1$, $\Gamma=0$), the critical value for $n$ is $2.5$.

To test this prediction, a set of simulations were performed above and below the critical value of the polytropic index with other parameters chosen to be representative of those for CTTS. The simulation set spanned $b=[200,350,500]$ and $n$ in the range 2.0–4.7 with a step size of 0.3 (although, lower $n$’s and $b$’s were explored – see §4.1). Transient effects from initialization dissipated over the course of a few sound crossing times ($\sim\,$days). After these effects subsided, the simulations converged to a steady state and the flow was able to smoothly pass through the sonic point and reach free fall speeds for all combinations of $b$ and $n$ in this set. The condition of steady-state was confirmed by differencing the output velocity information between time steps, yielding accelerations in the regime of numerical rounding error. In order to be sure this is not a numerical issue, additional tests of a steady state flow are described in §4.1. Snapshots of an example simulation with values of $n=2.0$ with $b=500$ converging to a steady state solution can be seen in Figure 2. The absolute value of the residuals from differencing the density at a chosen point and the density at the same point in the final time step in the simulation is also shown in Figure 2. Velocity and density approach the final value with oscillatory behavior with a linearly decaying amplitude in semi-log space as time progresses, indicating the simulation is indeed converging towards a steady state solution.

A second set of simulations was used to simulate systems with octupole + dipole stellar magnetic field geometry using the same distribution of values for $b$ and $n$. The addition of octupole components increases the amount of compression that occurs near the star, further heating the gas and raising pressure. Spontaneous variability for this field line geometry is predicted to occur when $n<4.5$. The relative contribution of the octupole moment compared to the dipole moment for the magnetic field, $\Gamma$, was fixed to a value of 10 for this set of simulations. These simulations were analyzed using the same method as the pure dipole case, and returned similar results; the simulations converged to a steady state after the initial transient behavior subsided. Again, the velocity was able to reach free fall speeds for each combination of $b$ and $n$.

4.1. Subsonic Simulations

A subset of more extreme simulations with low masses and stiff equations of state were unable to reach a steady state solution because the flow failed to remain supersonic after the initial transient phase subsided. The depth of the stellar potential well ($b$) determines where the sonic point is located along the column, and if it is small enough, the sonic point will lie within the star. In this case, slow ($|u|<0.1v_{ff}$) time-variable solutions dominate the flow (analogous to the known time-variable “breeze” solutions Velli, 1994; Del Zanna et al., 1998; Owen and Adams, 2016). These solutions occur when $b$ is small, the compression of the flow is high (from octupole magnetic field components), and the flow is heated significantly ($n$ is small) and the sonic point criterion cannot be satisfied in the region between the inner edge of the disc and the stellar surface.

The original criterion for non-free-fall flow has the form $n<\ell+3/2$ and thus does not depend on the depth $b$ of the gravitational potential well. On the other the hand, the above results show departures from free-fall flow with sufficiently small values of $b$. We can understand this behavior as follows. The incoming flow fails to reach free-fall speeds when the pressure term does not decrease fast enough relative to the gravitational term as the radius decreases. Under conditions where $n<\ell+3/2$, the pressure term will always dominate in the limit $\xi\to 0$. As we have shown, however, the small dynamic range of the flow (the disc edge is only an order of magnitude larger than the stellar radius) does not generally allow for the flow to reach the regime (small $\xi$) where the pressure dominates, so the flow retains nearly free-fall speeds. As shown in equations (24 – 26) of Adams and Gregory, (2012), the depth $b$ of the gravitational potential well appears as a coefficient in the terms under comparison. As a result, when $b$ is large, the gravitational term dominates over a larger dynamic range (in $\xi$) and free-fall speeds can be reached even when $n<\ell+3/2$. For small values of $b$, the gravitational term is less effective at enforcing free-fall speeds and departures can arise (as discussed above).

Nonetheless, these time-variable subsonic solutions are unlikely to present themselves in the region of parameter space realized by CTTS. Furthermore, observed signatures associated with shocks and high velocity flows at the surface of the star also indicate these solutions are not occurring. Although we reject these subsonic solutions as realistic models for CTTS, they may be relevant for systems with lower masses, e.g., planets and brown dwarfs, should they possess strong octupole (or higher order) magnetic field components. We caution this conjecture needs to be investigated in more detail before any predictions can be made.

4.2. Tests

Numerous tests were performed to check that the simulations were producing accurate results to ensure that the lack of spontaneous variability in the transonic solutions was not a spurious result from a numerical error. In any steady-state flow the Bernoulli potential is constant; the Bernoulli potential of the flow was calculated and was found to converge to a nearly constant value after initial transient effects had dissipated. Figure 3 shows the difference between the potential at the disk and the Bernoulli potential along the field line for the simulation shown in Figure 2. Although there are deviations from a constant value along the field line, they are small and the magnitude of the deviations is in part dependant on the resolution of the simulation. Under a resolution of 1024 cells, the maximum deviations are on the order of $0.1\%$. The largest deviations from a constant potential generally occur near the star where the velocity of the flow is largest and the divergence of the magnetic field lines is steepest. A second simulation is also shown in Figure 3 in which the number of cells has been doubled to 2048, yielding maximum deviations from potential at the disk on the order of $0.02\%$.

The construction of our numerical scheme does not actively ensure that entropy is conserved to machine precision, as we explicitly work with internal energy. This provides another method to test the simulation, as the level of entropy in each cell should converge to a constant value once transient behavior has subsided under an adiabatic equation of state. The entropy of each cell was calculated and then normalized to the value of entropy at the inner edge of the disk. Entropy remained roughly constant throughout the system with the largest deviations on the order of $1.5\%$. The bottom panel of Figure 3 shows the normalized entropy throughout the simulation shown in Figure 2. The same simulation run at the higher resolution is also shown in that panel with significantly smaller deviations, indicating that these deviations from constant entropy will continue to decrease with increasing resolution.

Finally, an analytic solution exists for the isothermal case. This provides a way to test the portions of the simulation not associated with the energy fluid equation. The results for a simulation using dipole geometry run under an isothermal equation of state is compared against the analytic solution in the Appendix. As shown in Figure A1, the simulation closely agrees with the analytic result.

4.3. Lack of Variability

All of the transonic simulations run with a constant gas disk density converged to steady solutions with no indication of spontaneous variability, or lack of near free-fall velocities close to the stellar surface. The tests of the non-isothermal simulations in addition to the comparison of an isothermal flow against the analytic result suggest that the simulation is behaving physically. Instead, the assumptions that went into the prediction of a constraint on the possible values for $n$ by Adams and Gregory, (2012) may be insufficient for realistic systems. In particular, those authors note that the assumption of working in the limit where $\xi\rightarrow 0$ may not hold where the dynamic range of the system is not sufficiently large. For example in CTTs the radial range is only about a factor of 10 from the stellar surface to the disc’s edge. The constraint placed on the polytropic index is based on how the thermal pressure compares to the ram pressure in the innermost regions of the system. If the flow is to reach free-fall speeds, the ram pressure must dominate over the thermal pressure, which is the case for all of the supersonic simulations. In the cases where $n$ has a value lower than the predicted constraint, the dynamic range of $\xi$ does not appear to be large enough for the density term to dominate since the minimum value of $\xi$ is necessarily fixed at 1 (the stellar surface).

5. RESULTS: DISK DRIVEN VARIABILITY

Following the inability to explain short-term accretion variability through spontaneous variability in the accretion column in the transonic solutions, we ran sets of simulations using time variable boundary conditions. 3D magnetohydrodynamic simulations have shown that density in the inner disk can be highly variable due to turbulence (Romanova et al.,, 2012). As an approximation to these naturally occurring density perturbations, the density at the inner edge of the disk was modulated manually. Two driving functions were used, sinusoidal oscillations and Gaussian pulses, both on top of a baseline density. These driving functions were chosen for convenience, allowing us to study a well posed problem, but are unlikely to represent the true density variations in real systems, which are likely to be stochastic. The internal energy, and thus the pressure, at the boundary is set such that the cells at the boundary are varied adiabatically or isothermally in the isothermal models. The simulation is initialized in the same manner as before and is allowed to evolve for 15 days before the driving functions are initiated to avoid transients that occur upon initialization. It is shown in the bottom panel of Figure 2 that relative departures from the steady state solution are on the order of $10^{-2}$ after this amount of time.

In these models, we choose the dimension-less parameters to correspond to a typical CTTS. These models were constructed with the following stellar properties: $M_{*}=0.5\,M_{\astrosun}$, $R_{disk}=10\,R_{*}$, and $R_{*}=1.5\,R_{\astrosun}$. This gives us a value of $b\sim 750$. Again, the temperature at the disk is assumed to be 10000K and at a density of $3\times 10^{11}\mbox{cm}^{-3}$. The parameter space formed by $\Gamma$, $n$, the amplitude of the perturbation $A$, and the timescale of the perturbation was explored to examine the behavior of the accretion column under a variety of conditions.

5.1. Isothermal Sinusoidal Driving

In the first set of simulations with driven accretion, the accretion column is assumed to be isothermal and the density at the inner edge of the disk was forced to oscillate sinusoidally with a fixed period $P$, and amplitude $A$, written as

[TABLE]

Simulations were run with P = [0.5, 1, 2, 4, 8] days, $\Gamma=[0,10]$ and $A=[1,2,3]$.

The behavior of the flow fell into three regimes depending on the driving period of the simulation. Simulations of isothermal flows with driving periods shorter than about a day (roughly the stellar sound crossing time) showed smaller amplitudes than the driving function near the stellar surface. This is because the flow can not react fast enough to the changes in density in the outer boundary, i.e. the flow cannot adjust on a time-scale shorter than the sound crossing time. The over-densities in the accretion column are mostly time averaged out by the time that they collide with the star. The mass accretion rate of a simulation displaying this behavior with $P=0.5$, $A=3$, and $\Gamma=0$. is shown by the orange line in Figure 4.

The second regime of behavior was observed in simulations with periods roughly one day or longer. As the density is modulated at the inner edge of the disk, weak shocks form and then propagate along the field line. Figure 4 shows two simulations in this regime in green and purple. These shocks form because the velocity of the flow is slow in the outer regions of the accretion column, allowing material to build up during periods of high density in the disk, and then fall at once in the form of a higher density shock onto the star. The period of the accretion rate oscillations at the star matches the driving period of the boundary condition for all sets of simulation parameters. Finally, if the driving period becomes much longer than the flow time-scale (tens of days), the flow evolves steadily along the steady-state solutions.

In these driven accretion simulations, density is the quantity manually controlled at the outer boundary, which allows the mass loss rate to change as a function of model parameters. The averaged mass accretion rate at the star for all of the isothermal simulations is shown in Figure 5. Simulations with larger amplitudes show higher mass accretion rates, which is expected because the mean density in the outer regions of the accretion column is higher for larger values of $A$. Systems driven with shorter periods in most cases have lower accretion rates. This is in part due to sharp peaks in the density causing pressure gradients both inward and outward, that prevent matter from flowing onto the columns, and in some cases can even lead to flow away from the star. As the driving period increases, the pressure gradient become shallower. Simulations run with magnetic fields comprised of dipole and octupole components show lower accretion rates than simulations run under the pure dipole field, in agreement with the analytic steady state solutions presented by Adams and Gregory, (2012). In this construction of the magnetic field lines with $R_{disk}=10R_{*}$, the net centripetal and gravitational force along the field line is away from the star in the outer regions of the accretion column for both the pure dipole and octupole + dipole configurations. This contributes to slow velocities in these regions.

5.2. Non-isothermal Sinusoidal Driving

Another set of simulations with sinusoidal driving density perturbations at the inner edge of the disk were run with values of n = [3.5, 4.5, 5.5] using the same set of periods and amplitudes as the isothermal case. These simulations exhibit much of the same behavior as the isothermal simulations. Similar to the isothermal case, weak shocks form in the outer regions of the accretion column in simulations with longer periods ($P>1$ day). Figure 6 shows a snapshot of a simulation with a shock traveling along the accretion column. The shock speed for the inset plot was found using the Rankine–Hugoniot jump condition derived from the continuity equation using densities and velocities measured on either side of the shock.

The time averaged accretion rate for the non-isothermal simulations is shown in Figure 7. Simulations with both short $0.5$d and longer $4.0$d oscillations are plotted to show both regimes of shocking and non-shocking behavior. For all simulations, the average accretion rate decreases as the accretion column becomes closer to isothermal. Changing $A$, $\Gamma$ and $P$ reveals the same patterns as the isothermal models shown in Figure 5.

5.3. Gaussian Density Pulses

To isolate the effects that a single accretion burst would have on conditions near the stellar surface and test how long it takes the accretion column to relax from a perturbation, simulations were run with boundary conditions that manually vary the density at the disk-magnetosphere footprint as a Gaussian in time. This causes a density pulse to travel along the magnetic field line until it collides with the star, which for appropriate choices of the time-scale is a shock. The density at the inner edge of the disk during the pulse is described as

[TABLE]

where $A$ is the amplitude of the driving Gaussian and $\sigma$ represents the timescale of the pulse.

Figure 8 shows the properties of the flow measured near the surface of the star as a function of time for several values of $n$, including the mass accretion rate, temperature, density, and velocity. The Gaussian density pulses shown have $\sigma=1$ and $A=3$. As discussed in section 5.1, the amount of material available for accretion onto the star is a function of the parameters describing the system, even if the time averaged density in the boundary remains the same.

The temperature, density, and velocity were measured near the star during the passage of the density pulse for a set of simulations with $A=[1,2,3]$, $\Gamma=[0,10]$, and $n=[3.5,4.5,5.5]$, the results of which are shown in Figure 9. The initial values of $T$, number density and $u$ for the steady state solutions, described in Section 4, are also shown in gray. As expected, the polytropic index is important in determining the temperature of the gas. As the polytropic index increases, the temperature of the gas decreases for all simulations and approaches the temperature of the upper layers of the gas disk (10000K). The density during and before the pulse is decreased in the simulations with higher polytropic indexes, but the effect is less dramatic than the effect on the temperature within the column. The velocity near the star gradually increases with increasing values of $n$, but does not quite reach the free-fall velocity.

The geometry of the magnetic field lines also plays a critical role in the environment near the star. The field lines converge faster near the star for the octupole + dipole geometry, causing a higher level of compression at small radii compared to the pure dipole configuration. For the most isothermal simulations shown here ($n=5.5$) while the pulse is passing through, the temperature difference between the pure dipole and the octupole + dipole configurations is on the order of 15000K. During steady state, the difference is roughly 8000K. Density increases by roughly an order of magnitude while the pulse is passing through. The changes in velocity are quite small between the two field configurations compared to the changes in the density and temperature as the thermal energy is relatively insignificant compared to the kinetic energy of the flow. The largest change in the velocity between the two configurations for any of the simulations during the passage of the pulse is only $\approx 1.2\%$ of the free-fall velocity.

The amount of time it takes the column to relax to its initial state provides a measure of the timescale between discernible accretion events. This timescale was found to be a function of the polytropic index and the geometry of the magnetic field lines. Reminiscent of the bottom panel in Figure 2, the absolute value of the residuals between density in a cell near the stellar surface and the same cell at the final time step after the Gaussian pulse has passed can be enclosed within a power-law envelope. The envelope was found for each set of parameters within the sample by taking a finite temporal derivative of the density at a cell near the star and fitting the residuals with a power law at the times of the zero crossings.

The slope of the fit, $\beta$, and the half-life of the residuals is shown for all of the simulations in this sample in Figure 10. The half-life for simulations under the octupole plus dipole magnetic field configuration is $\approx 10\%$ longer than those under a pure dipole configuration. The polytropic index also plays a role in how quickly the accretion column relaxes. Simulations with smaller values of $n$ are more rapidly damped than more isothermal simulations. For comparison, two fully isothermal simulations with $A=3$ were run for $\Gamma=0$ and $\Gamma=10$, and had $t_{1/2}=2.3,3.1$ days respectively. The amplitude of the initial Gaussian pulse plays a limited role in the relaxation timescale, which suggests a single slope/half-life may be sufficient for describing the behavior an accretion column undergoing events with similar timescales while the parameters governing the system remain constant.

6. DISCUSSION

We have discovered a new hydrodynamic mechanism that can lead to short-period, rapid changes in the accretion rate onto the star that may explain some of the short period accretion variability observed in CTTS. We find that smooth time-dependent variability in the inner disk that feeds the magnetospheric accretion column is amplified along the column. Under certain conditions these variations can shock in the accretion column, leading to rapid, large amplitude changes in the accretion rate onto the star. This is a promising mechanism to explain some of the variable behaviour observed in CTTS, as the inner disk that feeds the accretion column is likely to be variable. For example, density variations due to turbulence at the inertial scale have eddy turn-over times $\sim\Omega^{-1}$ (e.g., Fromang and Papaloizou,, 2006), which at the disk-magnetosphere boundary gives a time-scale of order a day with amplitudes of order unity. In further studies, we plan to model the exact observation signatures these shocks will produce; however, below we discuss our results in the context of recent observations and speculate about the connection between observations and our models.

A monitoring campaign of CTTS in NGC 2264 revealed variability in the light curves of most of the stars with significant UV excesses is dominated by short duration accretion bursts (Stauffer et al.,, 2014). For all but one of these moderately accreting stars (which was systematically different than the others), the median rise and fall times were 0.4 and 0.5 days respectively, measured between the beginning or end of the burst and the center value. Qualitatively, this is comparable to the rise/fall behavior shown by some of the driven accretion scenarios presented here, e.g., the isothermal model driven using the sinusoidal boundary condition with a 1.0 day period shown in Figure 4. We also find shorter rise times in our simulations because of the shocking behavior presented in §5.1. Although further investigation is required, we can speculate that the small asymmetry in the observed rise and fall times is because of the weak shocking behavior followed by a slower relaxation period. It is important to note that the boundary conditions used in these simulations are continuous, and the build up of material occurs within the active region of the simulation. Most of the accretion burst dominated light curves presented by Stauffer et al., (2014) did not show periodic structure and instead bursted stochastically, providing inspiration for future simulations with density fluctuations at the outer boundary distributed randomly in time with varying amplitudes. In the case of the variability at the disk being turbulence driven, one naturally expects a range of time-scales rather than a simple periodic signal.

In these simulations, the velocity approaches but never reaches the free-fall velocity measured from the inner edge of the disk

[TABLE]

because of the pressure support along the column. This is relevant for popular methods of calculating mass accretion rates that assume the material reaches $V_{ff}$ as it strikes the surface of the star (e.g., Calvet and Gullbring,, 1998; Gullbring et al.,, 1998). Because the kinetic energy flux ($\rho u^{3}$) has such steep scaling with $u$, the difference between $V_{ff}$ and the velocities shown in Figure 9 can result in calculated accretion rates that differ by up to $\sim 10\%$. Differences in velocities between simulations are all within $\sim 0.01V_{ff}$ of each other, so the effect on the calculated accretion rate will be a roughly constant multiplicative factor. As this is a nearly systematic effect for all models that assume $u\rightarrow V_{ff}$, this effect could easily be accounted for in future models of accretion shocks, although it does not drastically improve the accuracy of accretion rate measurements. Typical uncertainties in current measurements of accretion rates are generally much larger than 10$\%$. For example, uncertainties in $A_{V}$ alone can introduce a factor of 2 or larger (Ingleby et al.,, 2015).

As shown in Figure 9, the addition of octupole components to the magnetic field drastically increases the density in the accretion column near the stellar surface. The overall average accretion rate for this field configuration in our simulations is decreased slightly because of marginally slower velocities (see Figure 9 and Figure 7) and a smaller column cross section (Figure 1). The increased density in accretion columns with dipole plus octupole fields leads to higher kinetic energy fluxes. Higher kinetic energy fluxes lead to shock emissions that peak at shorter wavelengths (Calvet and Gullbring,, 1998). Combined with the smaller column cross section, this suggests that stars with significant octupole components are more likely to show evidence of hotter accretion spots with more limited surface coverage. Because the fractional changes between the initial and final velocity are small while pulses are passing through, changes in the energy flux at the surface of the star closely trace the density profile as a function of time. The thermal component of the energy flux within the column is essentially negligible for all the simulations shown here as the flow is highly super-sonic.

Line emission associated with magnetospheric accretion columns has been observed (e.g., Johns and Basri,, 1995). The simulation in its current form does not explicitly include radiative energy losses, which could play a role in determining the dynamics occurring within the accretion column, rather they are packaged up into in the polytropic index. As a first guess as to its importance, we can compare the optically thin cooling timescale to the flow-timescale. We note since this uses optically thin cooling it will provide a maximal estimate of the cooling rate, if any of the lines were to become optically thick, then the cooling rate would be reduced. Figure 11 shows the flow timescale, the time it takes for material to flow over a distance where the density changes significantly compared the cooling timescale derived from optically thin, solar metallicity, equilibrium cooling curves (Schure et al.,, 2009). The cooling is faster near the surface of the star where the magnetic field lines converge and material within the columns has been heated. Radiative cooling would likely drastically change the temperature structure seen in Figure 8. Figure 11 reveals cooling may also be significant near the magnetic footprint where the column connects with the disk where the flow velocities are slower. Cooling in the outer region of the accretion column would cause the flow to be more isothermal, affecting the formation process of the small shocks predicted in some of the driven accretion scenarios presented here. That said, discontinuities are present within the isothermal simulations presented here.

Because the cooling timescale is proportional to $\rho^{-1}$ for optically thin cooling and density in the simulation is normalized, the cooling timescale could be shifted up and down by changing the value for number density at the disk. The cooling curves shown here were generated assuming a density at the disk of $3\times 10^{11}\mbox{cm}^{-3}$, yielding a surface mass flux of $1.26\times 10^{-7}\mbox{g cm}^{-2}$. The cooling timescale curve spans such a large amount of timescales that even if the density at the disk was decreased by a factor of 10, radiative cooling would still be important within the inner and outer regions of the simulation. Because the mass accretion rate, and thus the density at the disk, can be highly variable, a range of values of the polytropic index are possible even for a fixed location along the accretion column. If the accretion rate is very low, e.g., $\dot{M}\approx 10^{-11}M_{\astrosun}/\mbox{yr}$, a more adiabatic equation of state is appropriate for the magnetic footprint at the disk ($n\rightarrow 1.5$ for a monatomic ideal gas). For a larger accretion rate, the gas is best approximated by an isothermal equation of state (n $\rightarrow\infty$). Using the simple $\rho^{-1}$ scaling for cooling timescales, estimates of when these regimes occur can be found using Figure 11 by scaling the cooling timescale curves inversely with density. The shape of the cooling curve is also itself a function of the polytropic index, so Figure 11 only provides a rough guide. Although the inclusion of radiative cooling is a logical next step in the development of this simulation, it would be a significant expansion on the current framework and thus will be left for future work.

Finally, another source of heating or cooling that could arise is from MHD waves that are explicitly ignored in our treatment. Under an assumed magnetic field of 3kG the Alfvén crossing time is significantly shorter than the flow crossing time ($\sim 20$ times shorter). This suggests that Alfvén (and fast) waves111In the low $\beta$ flows we consider here the slow wave is essentially the same as sound waves. are likely too fast to generate the sort of variability that creates these shocks; however, they may be able to transport energy and hence lead changes in the thermodynamics of the flow, which in principle could generate variability. Full MHD simulations are required to investigate the importance of MHD waves on the dynamics and thermodynamics.

The goal of this set of simulations was to expand 1D analytic steady state solutions into dynamic solutions allowing us to gain physical insight into their results. The simplicity of the 1D approach of these simulations also allows us to easily computed observable quantities. However, the limitation of 1D comes at the cost of the loss of information associated with 2D and 3D simulations. This includes the naturally occurring instabilities at the disk boundary (e.g., Romanova et al.,, 2012), which we have artificially replicated with driving functions. We also lose information about anything that occurs in the toroidal direction. That said, the use of a 1D system allows for a higher degree of control over the system and allowed for the connections between theory and simulation. One possible future advancement that could be introduced within this 1D construction is the misalignment of the magnetic field and the stellar rotation axis.

7. CONCLUSIONS

We have performed 1D hydrodynamic simulations of accretion onto young stars under multipole magnetic field configurations aligned with the rotation axis, using both dipole and dipole plus octupole configurations. Static and time dependent inflow boundary conditions were used at the inner edge of the disk to investigate the possibility of spontaneous accretion variability and explore the effects of a time-variable accretion rate in the disk. Physical parameters that govern the flow were varied within the range expected in real accreting CTTS. The major conclusions of this work are summarized below.

We found no spontaneous variability is generated within the column and the flow approaches steady state free-fall solutions for transonic solutions with values of $n$ lower than the constraint suggested by Adams and Gregory, (2012). The analytic constraint on $n$ holds only in the limit of small radius $\xi\to 0$. In contrast, we found it is likely that realistic young stellar systems lack the dynamic range in radius necessary for the constraint to come into effect. 2. 2.

By introducing time variability in the density at the inner edge of the disk, weak shocks can form in the outer regions of the accretion column and propagate towards the star resulting in rapid, high-amplitude changes in the accretion rate in systems with variability periods longer than roughly a day (of order the sound crossing time). Variability driven with periods less than a day do not appear to shock, and the oscillations in density are averaged out, as the flow cannot respond fast enough to the changes. With variability timescales much longer than the flow timescale, the flow profiles closely follow the steady solutions but slowly vary between them and no shocks form. 3. 3.

The rapid increases in accretion rate seen in our simulations are qualitatively comparable to short duration accretion burst dominated variability observed in monitoring campaigns (Stauffer et al.,, 2014). Specifically, the shocks in the accretion columns naturally result in fast accretion rate rise times as the shock hits the star, and slower fall times as the flow relaxes to steady-state on the order of the sound crossing time (Figure 4). 4. 4.

Simulations were run featuring single time dependent Gaussian density pulses to reveal how the accretion column reacts during both quiescent and high accretion states. Decreasing the polytropic index leads to higher temperatures and higher mass accretion rates. Magnetic fields with dipole plus octupole components lead to much higher densities near the surface of the star with only slightly lower velocities compared to pure dipolar fields. 5. 5.

We find the outer and inner regions of the accretion column are likely to be isothermal while the central portion may not be. Optically thin cooling timescales of material in the column indicate the temperature profile is more complex than a single adiabatic or isothermal equation of state can provide (Figure 11). Additional work to incorporate radiative cooling into this simulation is one possible next step to producing a more physically accurate representation of this system.

We have found short term variability is not spontaneously generated within the accretion column, at least for our 1D examination of the system under realistic parameters for CTTS. Instead, variability on $\sim$day timescales is likely driven by instabilities and turbulence near the magnetic footprints of the flow, which can be modulated by shocks forming in the accretion columns. The limitations of this study lie in the assumptions of a polytropic treatment of pressure and the 1D treatment of the accretion column. Despite these limitations, these simulations provide a readily interpretable link between analytic solutions, more complex simulations and observations. In the future we expect an approach of the kind developed in this work can be used to directly compute time-dependent observables, such as the UV excess.

8. ACKNOWLEDGEMENTS

We thank the anonymous referee who’s comments helped improve the manuscript. CER and CCE acknowledge that this material is based upon work supported by the National Science Foundation under Grant No. AST-1455042. JEO is supported by NASA through Hubble Fellowship grant HST-HF2-51346.001-A, awarded by the Space Telescope Science Institute (which is operated for NASA under contract NAS5-26555). FCA acknowledges support from the NASA Exoplanets Research Program NNX16AB47G, and from the University of Michigan.

Appendix A NUMERICAL DETAILS

A.1. Finite-Difference Equations

In this section we describe the finite difference scheme used to solve the system of hydrodynamic equations and discuss the associated scale factors and ancillary functions. The majority of the difference equations and their implementation were detailed in Owen and Adams, (2016); here we detail the changes that have been made to include the magnetospheric accretion geometry and the energy equation. Variables are subscribed with the location on the grid (i.e. $i-1$). Time dependant variables ($u$, $\alpha$, $\epsilon$, and $P$) are superscribed with the time step (i.e $n+1$). Time independent variables are superscribed with $a$ or $b$, depending on if the variable is face centered or zone centered respectively (see also Stone and Norman,, 1992). Recall in this staggered scheme, $\alpha$ and $\epsilon$ (and thus $P$) are zone centered while $u$ is face centered (e.g Stone and Norman,, 1992; Owen and Adams,, 2016).

During the source step, the velocity is updated due to forces from pressure, gravity, rotation and viscosity. The velocity is first updated using the following finite-differenced version of the momentum equation:

[TABLE]

where $\omega$ is a dimensionless parameter that measures rotation,

[TABLE]

The component of $x$ along the field line, $\hat{x}\cdot\hat{p}$, can be written as a function of $\xi$, $\theta$, and the ancillary functions $f$, $g$ & $H$ as

[TABLE]

Expressions for the ancillary functions and scale factors are shown below. The velocity and internal energy are then updated using an artificial viscosity tensor. More discussion on this tensor and the velocity update itself are included in Appendix A of Owen and Adams, (2016). The source step is finished with the addition of the compressional heating update for the internal energy, performed using the explicit finite difference expression shown in (Stone and Norman,, 1992).

Derivatives of the scale factors for both the face centered and zone centered grids were computed by finite differencing using the following expressions

[TABLE]

Owen and Adams, (2016) showed at sufficient resolution, the derivatives of the scale factors found by finite differencing are accurate enough for these geometric constructions for use in these simulations. This is useful since analytically obtaining these derivatives can be cumbersome.

The functional forms of the scale factors for $p$, $q$ and $\phi$ can be written as

[TABLE]

[TABLE]

[TABLE]

where $f$ and $g$ are ancillary functions defined as

[TABLE]

.

These ancillary functions recur often when using these coordinates and can be used in both the pure dipole case ($\Gamma=0$) and the octupole plus dipole case. These functions can be combined to write $\mathinner{\!\left\lvert\nabla p\right\rvert}^{2}$, defined as $H$:

[TABLE]

Finally, the transport step, which solves for changes in the time dependant variables due to fluid advection, follows the finite difference equations described in Stone and Norman, (1992).

A.2. Analytic vs. Simulation Comparison

In order to be sure our numerical method is working as expected, we can compare our code to an analytic solution that exists for isothermal accretion. A relatively simple solution for the velocity structure of the accretion column exists for the case of isothermal flow along dipole magnetic field geometry can be obtained using the procedure of (Cranmer,, 2004; Adams and Gregory,, 2012). Following Adams and Gregory, (2012), an analytic solution for the velocity $u$ in the isothermal case can be written as

[TABLE]

where $W$ is the Lambert $W$ special function (see also Cranmer,, 2004), $\xi$ is the dimensionless radius, $b$ is the dimensionless depth of the potential well and $\varepsilon$ is an integration constant written as

[TABLE]

We note the Lambert $W$ function has two real branches, the ‘[math]’ branch corresponds to sub-sonic flow and the ‘$-1$’ branch to super-sonic flow. This solution for the velocity assumes the location where the field line intercepts the disk is the co-rotation radius such that $\omega=bq^{3}$. The integration constant $\lambda$ can be written in the transcendental form

[TABLE]

where $\xi_{s}$ is the sonic point radii. Using the Lambert W function, the solution for $\lambda$ takes the form

[TABLE]

where the sonic point $\xi_{s}$ is calculated using a numeric solution to the matching condition for the sonic point. The matching conditions is written as

[TABLE]

An isothermal dipole simulation with $b=250$, a disk radius of $\xi_{\mbox{d}}=10$ and 1024 cells was run using a constant density outer boundary condition and outflow inner boundary conditions. The output from this simulation was then compared against the analytic solution shown above. For this value of $b$, the sonic point is at a radius of 8.71 stellar radii. A comparison of the analytic solution and the simulation result can be seen in Figure A1.