Time-dependent Models of Magnetospheric Accretion onto Young Stars

Under the assumptions outlined previously, 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; ∇p points along the field line and ∇q points in the poloidal direction, perpendicular to ∇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 & 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

$$\begin{eqnarray}{\boldsymbol{B}} & = & \displaystyle \frac{{B}_{{\rm{oct}}}}{2}{\xi }^{-5}[(5{\cos }^{2}\theta -3)\cos \theta \,\hat{{\boldsymbol{r}}}\\ & & +\displaystyle \frac{3}{4}(5{\cos }^{2}\theta -1)\sin \theta \,\hat{{\boldsymbol{\theta }}}]\\ & & +\displaystyle \frac{{B}_{{\rm{dip}}}}{2}{\xi }^{-3}(2\cos \theta \,\hat{{\boldsymbol{r}}}+\sin \theta \,\hat{{\boldsymbol{\theta }}}),\end{eqnarray} \tag{ 1 }$$

where ξ is the radius from the center of the star given in dimensionless stellar radius units, ξ ≡ r/R_* and θ 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 Γ is defined as the ratio of the coefficients of the composite multipole field,

$$\begin{eqnarray}&&{\rm{\Gamma }}\equiv \displaystyle \frac{{B}_{{\rm{oct}}}}{{B}_{{\rm{dip}}}}.\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&p=-\displaystyle \frac{1}{4}{\xi }^{-4}{\rm{\Gamma }}(5{\cos }^{2}\theta -3)\cos \theta -{\xi }^{-2}\cos \theta ,\end{eqnarray} \tag{ 3 }$$

while the orthogonal coordinate q is given by

$$\begin{eqnarray}&&q=\displaystyle \frac{1}{4}{\xi }^{-3}{\rm{\Gamma }}(5{\cos }^{2}\theta -1){\sin }^{2}\theta +{\xi }^{-1}{\sin }^{2}\theta .\end{eqnarray} \tag{ 4 }$$

Two configurations of magnetic fields were used in our simulation, Γ = 0 and Γ = 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 with 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).

Zoom In Zoom Out Reset image size

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 ${\boldsymbol{\Omega }}={\rm{\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

$$\begin{eqnarray}\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{u}}) & = & 0\\ \displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{u}} & = & -\,\displaystyle \frac{1}{\rho }{\rm{\nabla }}P-{\rm{\nabla }}{\rm{\Psi }}-{\boldsymbol{\Omega }}\times ({\boldsymbol{\Omega }}\times {\boldsymbol{r}})\\ \displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{e}{\rho }\right)+({\boldsymbol{u}}\cdot {\rm{\nabla }})\left(\displaystyle \frac{e}{\rho }\right) & = & -\,\displaystyle \frac{P}{\rho }{\rm{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 5 }$$

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

$$\begin{eqnarray}&&P=(\gamma -1)e.\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&a=\sqrt{\gamma \displaystyle \frac{{k}_{{\rm{b}}}T}{{m}_{{\rm{H}}}}},\end{eqnarray} \tag{ 7 }$$

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 ρ, e, and the velocity along the field line, u_p, and can be written as

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{\rho }{{\rho }_{1}}\quad u\equiv \displaystyle \frac{{u}_{p}}{{a}_{1}}\quad \epsilon \equiv \displaystyle \frac{e}{{a}_{1}^{2}{\rho }_{1}},\end{eqnarray} \tag{ 8 }$$

where a₁ and ρ₁ are the sound speed and density measured at the inner edge of the disk. Distances have been written in terms of the stellar radius (ξ), and time has been scaled using the stellar crossing time for a velocity a₁. This allows us to write the governing equations as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \alpha }{\partial t}=-\displaystyle \frac{\alpha u}{{h}_{p}}\left[\displaystyle \frac{1}{\alpha }\displaystyle \frac{\partial \alpha }{\partial p}+\displaystyle \frac{1}{u}\displaystyle \frac{\partial u}{\partial p}+\displaystyle \frac{1}{{h}_{q}{h}_{\phi }}\displaystyle \frac{\partial }{\partial p}({h}_{q}{h}_{\phi })\right]\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial t} & = & -\,\displaystyle \frac{1}{{h}_{p}}\left[u\displaystyle \frac{\partial u}{\partial p}-\displaystyle \frac{1}{\alpha }\displaystyle \frac{\partial P}{\partial p}-\displaystyle \frac{b}{{\xi }^{7}}\displaystyle \frac{\cos \theta f}{H}\right]\\ & & +\omega \sin \theta (\hat{x}\cdot \hat{p})\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{\epsilon }{\alpha }\right) & = & -\,\displaystyle \frac{u}{{h}_{p}}\left[\displaystyle \frac{\partial }{\partial p}\left(\displaystyle \frac{\epsilon }{\alpha }\right)-\displaystyle \frac{(\gamma -1)}{u}\displaystyle \frac{\epsilon }{\alpha }\displaystyle \frac{\partial u}{\partial p}\right]\\ & & -\displaystyle \frac{u}{{h}_{p}}\left[\displaystyle \frac{(\gamma -1)}{{h}_{q}{h}_{\phi }}\displaystyle \frac{\epsilon }{\alpha }\displaystyle \frac{\partial }{\partial p}({h}_{q}{h}_{\phi })\right],\end{eqnarray} \tag{ 11 }$$

where ω is a dimensionless parameter that measures stellar rotation, $\omega ={\rm{\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 the dimensionless formulation of this system is shown in Adams & Gregory (2012). The curvilinear coordinate scale factor for the p coordinate is denoted h_p, while f and H are ancillary functions of ξ and θ. The functional forms for these terms, along with an analytic expression for $(\hat{x}\cdot \hat{p})$, are shown in the Appendix.

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 needed, due to our complicated geometry. The full details of this artificial viscosity tensor are shown in Appendix A of Owen & Adams (2016). Finally, the internal energy of the system is updated due to compressional heating in a time-centered manner (Stone & 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).

At the stellar boundary, we apply standard outflow boundary conditions (see Stone & Norman 1992). At the disk boundary, we used two different approaches. For the steady simulations described in Section 4, the density and internal energy were held fixed, while for the simulations with explicit time variability described in Section 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.

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 × 10¹¹ cm⁻³; this is taken from the analytic accretion solutions of Adams & Gregory (2012) for a standard accretion rate of 1 × 10⁻⁸ M_⊙ yr⁻¹ and temperature of 10,000 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 1 cm² at the surface of the disk. Mass accretion rates are scaled this way to facilitate comparisons between the two simulated geometries that 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, Γ = 0 and Γ = 10, the ratio of the cross-sectional area of the column at the stellar surface to the disk is 5.14 × 10⁻⁴ and 4.18 × 10⁻⁵, respectively. The global stellar accretion rate in terms of a hot spot filling factor, f, and this standardized mass flux can be written as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{global}}=f(4\pi {R}_{* }^{2})\displaystyle \frac{\dot{M}}{A}.\end{eqnarray} \tag{ 12 }$$

For reference, a young accreting star with f = 0.03, R_∗ = 1.5 R_⊙ with a mass flux of $\dot{M}=1.0\times {10}^{-7}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ under the pure dipole construction, will have an accretion rate of 1.3 × 10⁻⁹ M_⊙ yr⁻¹ or 1.5 × 10⁻⁸ M_⊙ yr⁻¹ 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 β—the ratio of the gas pressure to the magnetic pressure. For a CTTS with a surface magnetic field strength B_dip = 3 kG (e.g., Johns-Krull 2007) and an inner disk edge of 10R_*, the parameter β ≪ 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 β can approach unity. As a 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 (with resolution tests discussed later). We fix the ghost cell used to set boundary conditions at the disk to a height of 0.05R_* above the midplane 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.

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| \lt 0.1{v}_{{ff}}$) time-variable solutions dominate the flow (analogous to the known time-variable "breeze" solutions; Velli 1994; Del Zanna et al. 1998; Owen & Adams 2016). These solutions occur when b is small, the compression of the flow is high (from octupole magnetic field components), 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\lt {\ell }+3/2$ and thus does not depend on the depth b of the gravitational potential well. On the other the hand, the previous 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\lt {\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 ξ) where the pressure dominates, so the flow retains nearly free-fall speeds. As shown in Equations (24)–(26) of Adams & 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 ξ), and free-fall speeds can be reached even when $n\lt {\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 previously).

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.

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%.

Zoom In Zoom Out Reset image size

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 are compared against the analytic solution in the Appendix. As shown in Figure 12, the simulation closely agrees with the analytic result.

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 & Gregory (2012) may be insufficient for realistic systems. In particular, those authors note that the assumption of working in the limit where $\xi \to 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 with 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 ξ does not appear to be large enough for the density term to dominate, since the minimum value of ξ is necessarily fixed at 1 (the stellar surface).

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

$$\begin{eqnarray}&&{\rho }_{{\rm{disk}}}(t)={\rho }_{1}\left[\displaystyle \frac{A}{2}\left(\sin \left(\displaystyle \frac{2\pi }{P}t\right)+1\right)+1\right].\end{eqnarray} \tag{ 13 }$$

Simulations were run with P = [0.5, 1, 2, 4, 8] days, Γ = [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 cannot react fast enough to the changes in density in the outer boundary (i.e., the flow cannot adjust on a timescale shorter than the sound crossing time). The overdensities 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 Γ = 0 is shown by the orange line in Figure 4.

Zoom In Zoom Out Reset image size

The second regime of behavior was observed in simulations with periods of roughly 1 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 timescale (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 composed 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 & 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.

Zoom In Zoom Out Reset image size

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 jump condition derived from the continuity equation, using densities and velocities measured on either side of the shock.

Zoom In Zoom Out Reset image size

The time averaged accretion rate for the non-isothermal simulations is shown in Figure 7. Simulations with both short 0.5-day and longer 4.0-day 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, Γ, and P reveals the same patterns as the isothermal models shown in Figure 5.

Zoom In Zoom Out Reset image size

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 timescale is a shock. The density at the inner edge of the disk during the pulse is described as

$$\begin{eqnarray}&&{\rho }_{\mathrm{disk}}(t)={\rho }_{1}(1+{{Ae}}^{-\tfrac{{t}^{2}}{2{\sigma }^{2}}}),\end{eqnarray} \tag{ 14 }$$

where A is the amplitude of the driving Gaussian and σ 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 σ = 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.

Zoom In Zoom Out Reset image size

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], Γ = [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 (10,000 K). 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.

Zoom In Zoom Out Reset image size

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 with 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 15,000 K. During steady state, the difference is roughly 8000 K. 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 with the changes in the density and temperature, as the thermal energy is relatively insignificant compared with 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 ≈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, β, 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 ≈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 Γ = 0 and Γ = 10, and had t_1/2 = 2.3 and 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 of an accretion column undergoing events with similar timescales, while the parameters governing the system remain constant.

Zoom In Zoom Out Reset image size

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 & 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, α, , 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 & Norman 1992). Recall in this staggered scheme, α and (and thus P) are zone centered, while u is face centered (e.g., Stone & Norman 1992; Owen & 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:

$$\begin{eqnarray}{u}_{i}^{n+1} & = & {u}_{i}^{n}-\displaystyle \frac{{\rm{\Delta }}t}{{h}_{p,i}^{a}}\displaystyle \frac{2}{({\alpha }_{i}^{n}+{\alpha }_{i+1}^{n})}\left(\displaystyle \frac{{P}_{i}^{n}-{P}_{i+1}^{n}}{{p}_{i}^{b}-{p}_{i+1}^{b}}\right)\\ & & -\displaystyle \frac{{\rm{\Delta }}t}{{h}_{p,i}^{a}}\left(\displaystyle \frac{b}{{({\xi }_{i}^{a})}^{2}}\right)\left({({\xi }_{i}^{a})}^{-5}\displaystyle \frac{{f}_{i}^{a}\cos {\theta }_{i}^{a}}{H}\right)\\ & & +\omega {\rm{\Delta }}t{\xi }_{i}^{a}\sin {\theta }_{i}^{a}(\hat{x}\cdot {\hat{p}}_{i}^{a}),\end{eqnarray} \tag{ 16 }$$

where ω is a dimensionless parameter that measures rotation,

$$\begin{eqnarray}&&\omega ={\left(\displaystyle \frac{{\rm{\Omega }}{R}_{* }}{{a}_{1}}\right)}^{2}.\end{eqnarray} \tag{ 17 }$$

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

$$\begin{eqnarray}&&\hat{x}\cdot \hat{p}=\displaystyle \frac{{\xi }^{-5}\sin \theta \cos \theta (f+g)}{\sqrt{H}}.\end{eqnarray} \tag{ 18 }$$

Expressions for the ancillary functions and scale factors are shown as follows. 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 & 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 & 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:

$$\begin{eqnarray}\displaystyle \frac{\partial {h}_{x,i}^{a}}{\partial p} & = & \displaystyle \frac{{h}_{x,i}^{b}-{h}_{x,i-1}^{b}}{{p}_{i}^{b}-{p}_{i-1}^{b}}\\ \displaystyle \frac{\partial {h}_{x,i}^{b}}{\partial p} & = & \displaystyle \frac{{h}_{x,i}^{a}-{h}_{x,i-1}^{a}}{{p}_{i}^{a}-{p}_{i-1}^{a}}.\end{eqnarray} \tag{ 19 }$$

Owen & 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 ϕ can be written as

$$\begin{eqnarray}&&{h}_{p}={\xi }^{5}{[{f}^{2}{\cos }^{2}\theta +{g}^{2}{\sin }^{2}\theta ]}^{-1/2}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{h}_{q}=\displaystyle \frac{{\xi }^{4}}{\sin \theta }{[{f}^{2}{\cos }^{2}\theta +{g}^{2}{\sin }^{2}\theta ]}^{-1/2}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{h}_{\phi }=\xi \sin \theta ,\end{eqnarray} \tag{ 22 }$$

where f and g are ancillary functions defined as

$$\begin{eqnarray}f & \equiv & {\rm{\Gamma }}(5{\cos }^{2}\theta -3)+2{\xi }^{2}\\ g & \equiv & \displaystyle \frac{3}{4}{\rm{\Gamma }}(5{\cos }^{2}\theta -1)+{\xi }^{2}.\end{eqnarray} \tag{ 23 }$$

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

$$\begin{eqnarray}&&H\equiv {| {\rm{\nabla }}p| }^{2}={\xi }^{-10}[{f}^{2}{\cos }^{2}(\theta )+{g}^{2}{\sin }^{2}(\theta )].\end{eqnarray} \tag{ 24 }$$

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 & Norman (1992).

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 that exists for the case of isothermal flow along dipole magnetic field geometry can be obtained using the procedure of Cranmer (2004) and Adams & Gregory (2012). Following Adams & Gregory (2012), an analytic solution for the velocity u in the isothermal case can be written as

$$\begin{eqnarray}&&{u}^{2}=-W\left[-{\lambda }^{2}{(q\xi )}^{-6}(4-3q\xi )\exp \left(-2\varepsilon -\displaystyle \frac{2b}{\xi }-\omega q{\xi }^{3}\right)\right],\end{eqnarray} \tag{ 25 }$$

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

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{1}{2}{\lambda }^{2}-\displaystyle \frac{3}{2}{bq}.\end{eqnarray} \tag{ 26 }$$

We note the Lambert W function has two real branches, the "0" branch corresponds to subsonic flow and the "−1" branch to supersonic flow. This solution for the velocity assumes the location where the field line intercepts the disk is the co-rotation radius such that ω = bq³. The integration constant λ can be written in the transcendental form

$$\begin{eqnarray}\mathrm{ln}\lambda -\displaystyle \frac{{\lambda }^{2}}{2} & = & 3\mathrm{ln}q{\xi }_{s}-\displaystyle \frac{1}{2}\mathrm{ln}(4-3q{\xi }_{s})-\displaystyle \frac{1}{2}\\ & & +b\left(-\displaystyle \frac{3}{2}q+\displaystyle \frac{1}{{\xi }_{s}}+\displaystyle \frac{1}{2}{q}^{4}{\xi }_{s}^{3}\right),\end{eqnarray} \tag{ 27 }$$

where ξ_s is the sonic point radii. Using the Lambert W function, the solution for λ takes the form

$$\begin{eqnarray}{\lambda }^{2} & = & -\,W\phantom{\rule{}{3.8ex}}[-\exp (-6\mathrm{ln}q{\xi }_{s}-\mathrm{ln}(4-3q{\xi }_{s})-1\\ & & +\left.\left.b\left(3q+\displaystyle \frac{2}{{\xi }_{s}}+{q}^{4}{\xi }_{s}\right)\right)\right],\end{eqnarray} \tag{ 28 }$$

where the sonic point ξ_s is calculated using a numeric solution to the matching condition for the sonic point. The matching condition is written as

$$\begin{eqnarray}&&3{\xi }_{s}\left(\displaystyle \frac{8-5q{\xi }_{s}}{4-3q{\xi }_{s}}\right)+3\omega q{\xi }_{s}^{4}=2b.\end{eqnarray} \tag{ 29 }$$

An isothermal dipole simulation with b = 250, a disk radius of ξ_d = 10, and 1024 cells was run using a constant density outer boundary condition and outflow inner boundary conditions. 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 12.

Zoom In Zoom Out Reset image size