Global 3D Simulations of Disc Accretion onto the classical T Tauri Star V2129 Oph
The magnetic field of the classical T Tauri star V2129 Oph can be modeled approximately by superposing slightly tilted dipole and octupole moments, with polar magnetic field strengths of 0.35 kG and 1.2 kG respectively (Donati et al. 2007, hereafter D07). Here we construct a numerical model of V2129 Oph incorporating this result and simulate accretion onto the star using a three-dimensional magnetohydrodynamic code. Simulations show that the disk is truncated by the dipole component and matter flows towards the star in two funnel streams. Closer to the star, the flow is redirected by the octupolar component, with some of the matter flowing towards the high-latitude poles, and the rest into the octupolar belts. The shape and position of the spots differ from those in a pure dipole case, where crescent-shaped spots are observed at the intermediate latitudes.
Simulations show that if the disk is truncated at the distance of which is comparable with the co-rotation radius, , then the high-latitude polar spot dominates, but the accretion rate obtained from the simulations (and from the accompanying theoretical calculations) is about an order of magnitude lower than the observed one. The accretion rate matches the observed one if the disk is disrupted much closer to the star, at . However, in that case the octupolar belt spots strongly dominate. In the intermediate case of , the polar spots are sufficiently bright, and the accretion rate is within the error bar of the observed accretion rate, and this model can explain the observations. However, an even better match has been obtained in experiments with a dipole field twice as strong compared with one suggested by D07.
The torque on the star from the disk-magnetosphere interaction is small, and the time-scale of spin evolution, years is longer than the years age of V2129 Oph. This means that V2129 Oph probably lost most of its angular momentum in the early stages of its evolution, possibly, during the stage when it was fully convective, and had a stronger magnetic field.
The external magnetic flux of the star is strongly influenced by the disk: the field lines connecting the disk and the star inflate and form magnetic towers above and below the disk. The potential (vacuum) approximation is still valid inside the Alfvén (magnetospheric) surface where the magnetic stress dominates over the matter stress.
keywords:accretion, accretion discs - magnetic fields - MHD - stars: magnetic fields.
Classical T Tauri stars (CTTSs) are young, low-mass stars and often show signs of a strong magnetic field (e.g., Basri, Marcy & Valenti 1992; Johns-Krull, Valenti & Koresko 1999; Johns-Krull 2007). The magnetic field plays a crucial role in disc accretion by disrupting the inner regions of the disc and channeling the disk matter onto the star (Ghosh & Lamb 1978, 1979; Königl 1991). Zeeman measurements of the surface magnetic fields of CTTSs show that the field is strongly non-dipole in many stars (Donati & Collier Cameron 1997; Donati et al. 1999; Jardine et al. 2002, 2006; Johns-Krull 2007). The magnetic field of the CTTS V2129 Oph has been recently measured with the ESPaDOnS and NARVAL spectropolarimeters (D07). Decomposition of the observed field into multipolar components has shown that the field of V2129 is predominantly octupolar ( kG), with a smaller dipole component ( kG) and a much smaller contribution from the higher multipole components. The large-scale magnetic field geometry has been calculated by extrapolating the surface magnetic field to larger distances by suggesting that the field is potential (vacuum approximation). The data presented in D07 correspond to the June 2005 observations. More recent measurements (July 2009) have shown that the field of V2129 Oph is larger: kG and kG (Donati et al. 2010, hereafter D10). Those authors suggest that the field of V2129 Oph varies in time due to the dynamo process. In this paper we assume that the field is fixed and corresponds to that given by D07. However, we do perform test experimental runs at larger values of the dipole field.
Analysis of the chromospheric spot distribution in V2129 Oph have shown that the polar, high-latitude spot dominates (D07, see also Jardine, Gregory & Donati 2008). D07 developed an approximate model where they investigated possible paths of matter flow in the potential field along closed field lines (see e.g., Gregory et al. 2006). They concluded that the disk should be disrupted by the dipole component far away from the star, at . In this case matter flows mainly along the field lines of the dipole component and hits the star at high latitudes which corresponds to observations. These estimates are useful first-order approximations to the processes in V2129 Oph. As a next step it is important to investigate matter accretion onto V2129 Oph in global MHD simulations, where the assumption that the field is potential can be dropped, and the interaction of the disk with the complex field of the star can be investigated.
In this paper, we developed a model star with parameters of V2129 Oph and with magnetic field distribution close to that observed in V2129 Oph (as described in D07). We approximated the field with superposition of the main (dipole and octupole) moments, and performed global 3D MHD simulations of matter flow around such a star. In the past, we were able to model accretion onto stars with quadrupole (Long, Romanova & Lovelace 2007, 2008), and octupole (Long, Romanova & Lamb, 2010a) magnetic field components. In this paper we investigate for the first time accretion onto a star with realistic magnetic field and parameters corresponding to CTTS V2129 Oph. In the accompanying paper, we investigate numerically accretion onto a model star with parameters close to BP Tau (Long et al., 2010b).
The goals of this paper are the following: (1) to investigate magnetospheric accretion onto this model star in global 3D MHD simulations, (2) to investigate numerically the magnetic field structure around the star and to compare it with the initial, potential field, (3) to compare the simulated and observed accretion spots, (4) to compare the results with results for a purely dipolar field, (5) to obtain the accretion rate from simulations and compare it with that obtained from observations, (6) to investigate the torque on the star.
Section 2 briefly describes the numerical model we use in our simulations. Section 3 compares the observed and modeled magnetic fields of V2129 Oph. Section 4 discusses the results of the simulations. Our conclusions are summarized in Section 5.
2 Numerical model and reference values
We solve the 3D MHD equations with a Godunov-type code in a reference frame rotating with the star using the “cubed sphere” grid. The model has been described earlier in a series of previous papers (e.g., Koldoba et al. 2002; Romanova et al. 2004a; Long, Romanova & Lamb 2010a). Hence, we describe it only briefly here.
Initial conditions. A rotating magnetic star is surrounded by an accretion disc and a corona. The disc is cold and dense, while the corona is hot and rarefied, and at the reference point (the inner edge of the disc in the disc plane), , and , where the subscripts ‘d’ and ‘c’ denote the disc and the corona. The disc and corona are initially in rotational hydrodynamic equilibrium (see Romanova et al. 2002 for details). An -type viscosity is incorporated into the code and helps regulate the accretion rate. The viscosity is nonzero only in the disc, that is, above a certain threshold density (which is in the dimensionless units discussed below). We use in this work.
Boundary conditions. At both the inner and outer boundaries, most of the variables are set to have free boundary conditions, . On the star (the inner boundary) the magnetic field is frozen into the surface of the star, that is, the normal component of the field, , is fixed, though all three magnetic field components may vary. At the outer boundary, matter is not permitted to flow back into the region. The free boundary condition on the hydrodynamic variables at the stellar surface means that accreting gas can cross the surface of the star without creating a disturbance in the star or the flow. This neglects the complex physics of the interaction of the accreting gas with the star, which is expected to produce a shock wave in the stellar atmosphere (e.g., Koldoba et al. 2008).
A “cubed sphere” grid is used in the simulations. The grid consists of concentric spheres, where each sphere represents an inflated cube. Each of the six sides of the inflated cube has an curvilinear Cartesian grid. The whole grid consists of cells. In our dipole+octupole model of V2129 Oph, the octupole component dominates and a high radial grid resolution is needed near the star. We achieve this by choosing the radial size of the grid cells to be 2.5 times smaller than the angular size in the region , while it is equal to the angular size in the outer region as in all our previous work. The typical grid used in simulations has grid cells.
Reference values. The simulations are performed in dimensionless variables . To obtain the physical dimensional values , the dimensionless values should be multiplied by the corresponding reference values as . These reference values include: mass , where is the mass of the star; distance , where is the radius of the star; velocity . The time scale is period of rotation at : .
To derive the reference values for the magnetic field, we take into account that the tilts of the dipole and octupole moments relative to the rotational axis are small, and hence we use the formulae for the axisymmetric case: and (see, e.g., Long, Romanova & Lamb 2010a). We suggest that the polar field: , . The reference magnetic moments for the dipole and octupole components are and respectively, where is the reference magnetic field. Hence, the dimensionless magnetic moments are: , . We take one of the above relationships; for example, the one for the dipole; to obtain
Hence, the reference magnetic field depends on the dimensionless parameter and the dipole magnetic field, . The reference density and the reference mass accretion rate also depend on these parameters:
We find the other reference values in a similar way: the torque ; energy flux, ; energy flux per unit area, ; temperature , where is the gas constant; and the effective blackbody temperature , where is the Stefan-Boltzmann constant. Tab. 1 shows the reference values which are fixed in all models. Tab. 2 shows the reference values for our main models which depend on the parameters and (see §4). In the paper, we show our results in dimensional units. Formulae 2, 3 and Tab. 2 help scale the results to other parameters.
In most of the models we change the dimensionless octupolar moment in same proportion as so as to keep the ratio fixed (models 1-3). To find the ratio between the dimensionless moments, we use the approximate formulae for aligned moments: and obtain the ratio (where ). In a few test cases, we increased the relative strength of the dipole field (models 4-5).
|Model||(G)||(G)||(g cm)||(yr)||(erg cms)||(g cms)|
3 Reconstructed magnetic field of V2129 Oph and the field used in 3D simulations
The magnetic field of V2129 Oph has been reconstructed from observations by D07. The authors decomposed the measured surface magnetic field into its spherical harmonic components. They then approximated the field as a superposition of low-order multipoles, adjusting the magnitudes of the multipole moments and their orientations to give the best match of the modeled field to the measured field.
In one of the best matches the surface magnetic field of V2129 Oph is approximated by a 0.35 kG dipole and a 1.2 kG octupole field, with corresponding magnetic moments tilted by and (towards phases 0.7 and 0.9) with respect to the rotational axis (D07).We choose this field distribution as the one for our dipole+octupole model. Higher-order multipoles and azimuthal component are also present, but they carry much less magnetic energy than the two main components.
Note that method of field reconstruction from observations has some uncertainties, and the error bars are different in different stars and for different components. The error bar is typically for the strongest field component (the octupole in the case of V2129 Oph) but larger for the weaker field components. In V2129 Oph, the dipole field is 3-4 times weaker than the octupole, and the error bar can be 3-4 times larger, that is . The error bars on the phases and tilts of the magnetic axes (about the rotational axis) also depend on the actual values: the phases are less accurate when the tilts are low. In the case of V2129 Oph, the error bar in the phase of the dipole axis is likely to be 0.1-0.2 rotation cycle, or larger, while the tilt of the dipole axis is accurate to no better than .
Fig. 1 shows the polar projection of the modeled (dipole plus octupole) field on the star’s surface. The style of the figures and the components of the field , and are similar to those reconstructed from observations by D07 (see their Figs. 12 and 13). Comparisons of our plots with the D07 plots show that the radial (the strongest) component of the magnetic field is very similar in both plots (see top left panel in Fig. 12 of D07). It shows a strong positive pole, slightly misplaced relative to the rotational axis, then a negative octupolar ring above the equator, and a positive octupolar ring below the equator. In the D07 plot, the octupolar rings are not smooth in the azimuthal direction, which shows the presence of the non-axisymmetric modes. However, this difference is not very important for the present modeling, because most of the magnetic energy is in the axisymmetric multipoles.
The surface distribution of the meridional component of the field, (Fig. 1, right panel) is quite similar to that in Fig. 13 of D07 (left bottom panel). In both cases we see an inner ring of positive polarity (red), and an outer ring of negative polarity (blue). The distribution of the azimuthal component, , shown in Fig. 1, is somewhat similar to that shown in Fig. 12 of D07 (left middle panel). However, our azimuthal field is weaker because we do not take into account the toroidal component of the field (see the right panel of Fig. 13 of D07 for the toroidal component). We suggest that the toroidal and the higher-multipole poloidal components may influence the spot shapes near the star and can be included in future modeling.
Fig. 2 shows the distribution of the modeled field on the surface of the star. The middle and right panels show strong northern (positive) and southern (negative) poles. The left panel shows that there are also northern (negative) and southern (positive) belts. One can see that the field is a typical octupole field tilted by (see also Long, Romanova & Lamb 2010a). The dipole component is present but is much smaller.
An important radius is the radius where the dipole and octupole components of the field are approximately equal. Below, we derive an approximate value of such a radius, assuming that the dipole and octupole moments are aligned with the rotational axis. The magnetic fields of the dipole and octupole vary with distance as (see Long, Romanova & Lamb 2010a, Eqn. 1): (i) in the equatorial plane: , , and (ii) in the polar direction: , . Equating and in the equatorial plane and evaluating and at the surface of the star (e.g., , we obtain the distance where the dipole and octupole components are equal:
Here, we took into account that kG and kG. That is, the octupolar field is expected to dominate at , whereas the dipole field is expected to dominate at larger distances.
Fig. 3 shows the distribution of the intrinsic magnetic field in our dipole+octupole model of V2129 Oph (at ). The left panel shows that the field near the star has an octupolar structure up to a radius of which is in agreement with the estimated above. The right panel shows that at larger distances the field has a more ordered, dipolar structure.
4 Modeling of accretion onto V2129 Oph
In simulating accretion onto CTTS V2129 Oph, we use its mass and radius derived from observations which are , (see D07 and references therein) and rotation period days (see Shevchenko & Herbst 1998). These values imply that the corotation radius . An estimate age of V2129 Oph is years (D07).
A number of simulation runs were performed. To change the size of the magnetosphere (which is equivalent to the changing of the dimensional accretion rate at the fixed magnetic field on the star), we varied the parameter from 0.5 to 5 and obtained a set of simulation runs where the disk stops at different distances from the star. We find from the simulations that at the inner disk radius is larger than the corotation radius and the disk moves away due to the “propeller” effect. At smaller parameters, the magnetosphere is large, , however, the accretion rate is much smaller than the observed one. We chose some intermediate case with as the main one. In that case the disk stops far away from the star (which is close to the corotation radius, ), and at the same time an accretion rate is not too low. We refer to this case as the main one. It corresponds to model 1 in Tab. 3. For this case we performed detailed analysis of different features of the flow, hot spots and the large-scale magnetic field evolution (see sections 4.1-4.3, 4.5, 4.6). In models 2 and 3 (see Tab. 3) we investigate additional cases, where the disk comes closer to the star and accretion rate is higher and is closer to the observed one. However, the low-latitude octupolar belt spots become stronger than the polar spot (see sec. 4.7). In models 4 and 5 we keep the same octupolar field, but increase the dipole magnetic field by a factor of 2 and 3 respectively. This helps to increase an accretion rate and at the same time octupolar belt spots are weak and the spot shape is closer to the observed one (see sec. 4.9). In model 6 we investigate accretion on to a star with a pure dipole field, and compare results with the main case (see sec. 4.4). Models 7-9 are theoretical models (see sec. 4.8).
4.1 Matter flow around the star
Here we discuss results of simulations in our main case (model 1 in Tab. 3). Initially, we place the inner edge of the disk at the radius . The small viscosity incorporated into the disk leads to inward matter flow. First, the disk moves inwards over rotations (Keplerian rotations at ). Later, the funnel streams form and the matter accretes onto the star quasi-stationary. We were able to reach time in this type of simulations. Fig. 4 shows the matter flow and selected field lines at . One can see that the disk is stopped by the dipole component of the field, and matter flows towards the star from two directions, which correspond to the directions of the dipole component tilt ( plane).
Fig. 5 shows the density distribution in three slices, , and . One can see that the disk is disrupted by the dipole component of the field far away from the star, at , and the funnel streams form at . The red line shows the surface where the matter and magnetic field stresses balance each other. The dominant stress is given by the -component of the stress tensor. The magnetospheric radius is then the radius where . Inside this radius the magnetic stress dominates (Ghosh & Lamb 1978, 1979; Long, Romanova & Lamb 2010a; see also Bessolaz et al. 2008 for different definitions of the inner disk radius). We obtain .
Closer to the star, the octupole component of the field becomes stronger, and each funnel stream splits into two parts, one of which flows towards the magnetic pole, and the other into an octupolar belt. This produces two types of accretion spots on the star. We discuss this in more detail in §4.3.
4.2 Magnetic field structure and the validity of the potential approximation
Next, we investigate the magnetic field evolution due to the disk-magnetosphere interaction. Fig. 6 (left two panels) shows the initial magnetic field distribution near the star (top panel) and in a larger region (bottom panel). These plots correspond to the potential approximation which is often used for extrapolating the surface magnetic field to larger distances (e.g., D07). As the simulations proceed, the plasma and magnetic field motions create currents and fields in the simulation region and the resulting magnetic field geometry strongly departs from the potential one. The dominant processes are connected with the differential rotation of the foot-points of the magnetospheric field lines which have different angular velocities at the star’s surface, in the disk, and in the corona. Differential rotation leads to conversion of the rotational energy into magnetic energy and to inflation of the field lines (e.g., Aly 1980; Lovelace et al. 1995; Romanova et al. 1998; Agapitou & Papaloizou 2000). The foot-points threading the star usually rotate faster than those threading the disc, and the magnetic field lines wrap around the rotational axis forming a magnetic tower (e.g. Lynden-Bell 1996). Formation of such a tower has been observed in different numerical simulations (e.g., Ustyugova et al. 2000; Kato, Hayashi & Matsumoto 2004; Romanova et al. 2004b, 2009; Ustyugova et al. 2006) and is a natural outcome of disk-magnetosphere interaction. The middle and right bottom panels of Fig. 6 show formation of such a tower. So, at large scales the magnetic field strongly departs from the potential one.
The situation is different near the star. Fig. 6 (top panels) shows that the field structure at is almost the same as at . This is the region , where the magnetosphere is not disturbed by the disk, and magnetic stress dominates over the matter stress, and hence the magnetic fields produced by plasma flow are insignificant compared with the main magnetic field of the star produced by currents flowing inside the star. The potential approximation is valid in this region.
We conclude that the potential approximation used in many works to recover the external magnetic field from the observed one (e.g. Jardine et al. 2002; D07) is sufficiently good in the vicinity of the star, where the magnetic stress dominates over the matter stress (inside the Alfvén surface, ), and where the magnetosphere is not disturbed by the accretion disc. However, the field strongly departs from the potential one at larger distances.
It is often suggested that in young stars a strong stellar wind stretches the stellar field lines in the radial direction, like in the solar corona, and this determines the geometry of the external magnetic field (e.g. Safier 1998; Gregory et al. 2006; Jardine, Gregory & Donati 2008). We should note that this solar-type scenario can dominate in the case of diskless, weak-line T Tauri stars. However, if a large-scale magnetic field truncates the accretion disk, then a different, magnetic tower geometry dominates. Such a tower can accelerate a small amount of coronal matter up to high velocities and can carry some angular momentum from the star to the corona. However, both processes are significant only in cases of rapidly rotating stars (e.g., Romanova et al. 2005; Ustyugova et al. 2006), and are less important in slowly rotating stars, like V2129 Oph (Romanova et al., 2009).
Gregory et al. (2008) noted that the dominance of the octupolar component near the star causes a smaller fraction of the stellar surface to be connected to open field lines compared with the case of a purely dipolar field. We observed that in V2129 Oph the amount of flux in open field lines depends on the position of the disk and on the strength of the dipole component which dominates at and which determines the amount of open flux. The field lines inside the inner disk are closed. All other field lines of the dipole component inflate and represent open magnetic flux. The octupolar component contributes towards open flux in polar region, but this flux is much smaller than the dipole open flux, as suggested by Gregory et al. (2008).
4.3 Accretion spots in V2129 Oph
D07 analyzed the brightness enhancements observed in V2129 Oph in the chromospheric lines, like in CaII IRT, and also in the HeI D3 line, which presumably forms in the shock wave near the star’s surface. The right panel of Fig. 9 of D07 shows the CaII IRT excess emission relative to the basal chromospheric emission, which can be interpreted as the chromospheric “hot spot”. Fig. 9 of D07 shows that the main part of the spot, the polar spot, is located at colatitudes with a center at . The position of this spot almost coincides with the position of the magnetic pole (see Fig. 12 top left panel of D07). In addition, non-axisymmetric spots of lower brightness are seen at higher colatitudes up to .
In our numerical model, matter flows towards the star in funnel streams. We suggest that all matter which approaches the star freely moves inward through the star, and hence we neglect the complex physics of the stream-star interaction. Instead, we calculate the matter and energy fluxes carried towards the star, and our accretion spot represents a slice in density/energy distribution across the funnel stream.
In reality, the stream interacts with the dense regions of the chromosphere and photosphere, and forms a shock wave where gas is heated, and the energy is radiated mainly behind the shock front (e.g., Lamzin 1995; Calvet & Gullbring 1998). The shock wave may oscillate due to cooling/heating instability behind the shock front (Koldoba et al., 2008), and the stream-star interaction may lead to formation of waves on the star’s surface and possible outflows (Cranmer 2008, 2009; see also Brickhouse et al. 2010). We suggest that the heating processes in the dense layers of star’s atmosphere may excite such lines as CaII IRT, HeI D3 and others, and hence we expect a correlation between the position of the chromospheric spots observed by D07, and accretion spots observed in our simulations.
Fig. 7 shows different views of the density spots. One can see that there are two main spots, the polar ones, which almost coincide with the octupolar magnetic poles. There are also two elongated spots at much lower latitudes where matter flows towards the magnetic octupolar belts.
Fig. 8 shows a polar view of the spots. The left panel shows the density distribution. The right panel shows the distribution of the energy flux which is calculated as: (here is the surface velocity of the star) and shows the location of the main energy deposition in the funnel stream. One can see that the polar density spot is located at colatitudes of and is centered at . The octupolar belt spots are azimuthally elongated and are located at . The energy spots have similar shapes and locations. We also calculated the distribution of the matter flux and temperature (calculated in the suggestion that all kinetic energy is converted into black body radiation). We observed, that their distribution is similar to that of the density and energy fluxes shown in Fig. 8. Any of these values can be used for analysis of accretion spots on the star. Note that the size of spots in our simulations depends on the chosen density/energy level, so that the densest/brightest parts of spots have a smaller size.
Comparison of the accretion spots obtained in our simulations with the CaII IRT chromospheric “spots” observed in V2129 Oph (Fig. 9 of D07) shows that the polar spot is present in both cases and has a similar location and shape. The simulated spot is somewhat more compact than the observed one (see however §3.6 where we discuss the dependence of the spot size on the density/energy levels). Also, the center of the observed spot is located at a slightly higher latitude than the simulated spot. In both cases, the polar spots are located near the magnetic poles seen in Fig. 1 of our paper (and Fig. 12 of D07), though the position of the simulated density spot is somewhat different in phase compared with the magnetic pole position.
Fig. 9 shows that in our simulations an accretion spot (shown in black contour lines) does not coincide with the maximum of the magnetic field (red region on the star). A spot is leading (about 0.15 in phase). Here we note that the spot is located between the dipole and octupolar poles. We suggest that the main accretion stream moves toward the dipole magnetic pole and hence the accretion spot is located close to this pole. However, an octupole field also influences the accretion stream, and as a result a spot is located between the dipole and octupolar magnetic poles. This may explain why we observe the phase difference of about 0.15 between the main (octupolar) magnetic pole and accretion spot. The uncertainty in the tilt and phase of the dipole field derived from observations is about 10%. If we assume that the phase difference between the dipole and octupolar poles is smaller, then we expect that the accretion spot will be closer to the octupolar magnetic pole (as observed by D07). There is also another possibility, which is that the funnel stream is dragged by the rapidly rotating disk and moves slightly ahead of the more slowly rotating magnetosphere, and this leads to the phase difference between the spot position and the magnetic pole (as shown in simulations by Romanova et al. 2003, 2004a).
4.4 Comparison with a pure dipole case (model 6)
We performed a special set of runs where the magnetic field of V2129 Oph is approximated with the dipole component only (zero octupolar field). The strength of the dipole field, 0.35 kG, and its tilt, , are the same as in the main dipole+octupole case. Figure 10 (top right) shows the matter flow and spots in the -plane in the case of a pure dipole field. One can see that matter flows in two narrow streams and the hot spots have the azimuthally elongated, crescent shapes, typical for stars with the dipole field (Romanova et al., 2004a; Kulkarni & Romanova, 2005). The spots are centered at a much lower latitude, , compared to spots in the dipole+octupole case.
The left panels show the cross-section and spots in our main dipole+octupole case for comparison. One can see that initially matter flows along the dipole field lines. However, closer to the star, the octupole component determines the flow pattern: the funnel stream becomes wider and splits into a polar stream and an equatorial belt stream. Note that the polar stream is redirected by the octupolar component to much higher latitudes compared with the stream in a pure dipole case.
We conclude that the high-latitude round accretion spots obtained in simulations of the dipole+octupole model are in much better qualitative agreement with the observed accretion spots on V2129 Oph (D07) than the lower-latitude crescent-shaped accretion spots obtained in the pure-dipole model.
In this work we have neglected the higher-order multipoles which can dominate near the stellar surface. We should note that they cannot change the position of the main accretion spot, but instead can only redistribute the density inside the main spot.
4.5 Area covered by hot spots
Here we calculate the area covered by the accretion spots. We choose a density , and calculate the area covered by the part of the spots with density higher than , and the fraction of the star covered by these spots: . We calculate this fraction for different values of . To separate the contribution from the polar spots, we calculate separately the fraction of the star covered by all spots (), and by polar spots only, . We also calculate the fraction of the star covered by spots in the pure dipole case, . Fig. 11 shows how the spot coverage varies with . One can see that the area covered by spots in the dipole case is larger than that in the dipole+octupole case if g cm (). The situation is reversed if the density level is smaller. Note that in earlier comparisons of accretion to stars with dipole and quadrupole fields, we observed a similar trend: at higher density levels, the spots in the pure dipole case are larger than those in the dipole plus quadrupole cases, while at lower density levels, the situation is reversed (Long, Romanova & Lamb, 2010a).
In another figure, Fig. 12, we choose several density levels and show the size of spots corresponding to these density levels. One can see that at the lowest density level, g cm (), the spots occupy , while at the highest density level, g cm () , they occupy . If only the polar spot is taken, then the coverage for the same density levels is smaller — and respectively. Fig. 11 shows that at even higher density levels, the spot coverage is very small. For example, at g cm, .
Thus, the spot size strongly depends on the density level. This density distribution shows the structure of the funnel streams: the density is higher in the middle of the stream, and decreases toward the periphery. The matter flux, kinetic energy flux and temperature have a similar spot-centered distribution and hence the spot coverage is different at different energy/temperature levels with the hottest spots occupying a smaller fraction of the star’s surface (see also Romanova et al. 2004a).
The size of the accretion spots reconstructed from observations, e.g., of the CaII line, also depends on brightness as seen in Figs. 9 and 11 of D07 and Fig. 6 of D10. The fraction of the star covered with each of visible spots is approximately 5% which is about 10% for both spots (observations take into account only the spots in the visible hemisphere, and do not take into account a probable antipodal spots; to compare simulations with observations, we multiply the observed fraction by a factor of two). One can see that the brightest parts of the observed spots cover a much smaller area of the star. The position and shapes of the observed and simulated spots are similar: both are round and are located at the co-latitude of .
Fig. 12 (middle panel) shows that at this spot coverage () the polar part of the spot strongly dominates, while the octupolar belt spots are weaker. We conclude that there is a good correspondence between the sizes and longitudinal positions of the main, brightest parts of the spots obtained in observations and simulations. However, there is a difference in shapes and phases of the lower-density/energy parts of the spots (see also §4.4 and §4.6).
It was predicted that the area covered by spots in the case of complex magnetic fields should be smaller compared with the spot area in the pure dipole case (e.g., Mohanty & Shu 2008; Gregory et al. 2008, see also Calvet & Gullbring 1998). In addition, Calvet & Gullbring (1998) derived from their model and observations of a number of CTTSs that their spot coverage is expected to be very small, less than a percent in most stars. However, our simulations show that the spot coverage is higher than that predicted by Calvet & Gullbring (1998). It is possible that these authors analyzed only the central, brightest parts of the spots. This problem requires further investigation.
|Observations of V2129 Oph||(G)||(G)||(Myr)||spots (at )|
|Model of V2129 Oph|
|model 1.||dip+oct 3D||1.5||0.33||350||1200||0.015||6.2||polar+small arcs|
|model 2.||dip+oct 3D||1.0||0.22||350||1200||0.029||4.3||polar+long arcs|
|model 3.||dip+oct 3D||0.5||0.11||350||1200||0.023||3.4||polar+rings|
|model 4.||dip+oct 3D||1.5||0.165||700||1200||0.029||5.9||polar|
|model 5.||dip+oct 3D||1.5||0.11||1050||1200||0.036||5.5||polar|
|model 6.||dip 3D||1.5||0.00||350||0||0.046||6.0||crescent|
|model 7.||dip (theory)||350||0||6.2||-|
|model 8.||dip (theory)||350||0||4.3||-|
|model 9.||dip (theory)||350||0||3.4||-|
4.6 Rotational modulation
D07 observed that the fluxes and amplitudes of the Zeeman signatures in the CaII IRT and HeI D3 emission lines exhibited strong rotational modulation over the rotational cycle of the star which is close to the 6.53 d period in V2129 Oph. The time-dependent longitudinal field derived from these emission lines also shows the rotational modulation (see Fig. 5 of D07). The emission in CaII IRT traces the regions of higher excitation of this line in the chromosphere, while the HeI D3 line is thought to form in the shock wave (Beristain, Edwards & Kwan, 1998). These features can be compared with the rotationally modulated accretion spots obtained in our simulations.
We calculate the rotationally modulated radiation from accretion spots towards the observer located at an angle relative to the rotational axis of the star. To calculate the radiation, we suggest that all the energy flux of the funnel stream is converted into isotropic black-body radiation in the shock wave (see details in Romanova et al. 2004a). We choose some late moment of time () when the position of the spots does not vary much and obtain the rotationally modulated radiation from such spots.
Fig. 13 compares the rotationally modulated longitudinal magnetic field derived from CaII IRT observations by D07 (blue line) with our rotationally modulated spot curves (red lines) which we calculated at the inclination angle (D07). One can see the gradual rise and sharper drop in both simulated and observed curves which may be a sign of the spots asymmetry. We calculated a similar rotationally modulated curve for V2129 Oph with the magnetic field approximated by a pure dipole field and found that the curve is of a sinusoidal shape (see the dashed black line in Fig. 13). We suggest that the non-sinusoidal shape may be connected with the presence of the low-latitude (octupolar belt) part of the spot in our simulations, and with the low-latitude components of the spots in the observations.
The observed and simulated curves have a phase difference of about 0.1-0.2 which probably reflects the phase difference between the simulated and observed spots.
4.7 Accretion rate onto V2129 Oph
The accretion rate in V2129 Oph has been estimated by D07 using the CaII IRT line to be M yr (using the approach of Mohanty, Jayawardhana & Basri 2005). Using fluxes from several emission lines and more recent empirical relations between line emission fluxes and accretion luminosities (Fang et al., 2009), a revised mass accretion rate was derived from the same spectra: , or: with the maximum and minimum values of and .
Below we compare the observed accretion rate with the rates obtained from simulations and derived from theoretical estimations. The dimensional accretion rate is
In an attempt to resolve this issue, we performed simulation runs at smaller values of the magnetospheric parameter, and (see models 2 and 3 in Tab. 3). In each case, we adjusted the dimensionless octupolar moment using the ratio: , thus fixing the ratio between the dipole and octupole moments. We observed from the simulations that for , the accretion rate is higher and it matches the minimum value derived from observations, (see model 2 in the Table). For , the accretion rate is sufficiently high to match the observed accretion rate (see model 5).
Fig. 14 (top panels) shows matter flow at different . One can see that the size of closed magnetosphere systematically decreases with . We determine (very approximately) the truncation radius as the radius at which the last closed (or, strongly disturbed) field line crosses the equatorial plane (shown as a dashed line on the plot). Then for , we obtain a set of radii, . Each of these radii has an error of about , because there are usually a few field lines which can be considered as slightly disturbed lines. Nevertheless, we see a clear trend of decreasing with .
Fig. 14 (bottom panels) shows that an area covered with the octupolar belt spots increases when decreases, because the disk is truncated closer to the star, where the octupolar component dominates. Note that in the case of , the polar spot dominates over octupolar ring spots if we take the brightest/densest of the stellar surface. However, in the case of , octupolar ring spots give the brightest part at any density/energy level. Such octupolar spots have not been observed by D07. So, the intermediate case (, model 2) can possibly explain the observations, though the truncation radius is somewhat smaller than the corotation radius, .
To further investigate this issue, we estimated the accretion rate using a theoretical approach. Prior simulations done in our group show that the formula for the Alfvén radius (Elsner & Lamb, 1977):
is approximately applicable for the truncation of the accretion disks, if (Long, Romanova & Lovelace 2005; see also Bessolaz et al. 2008 for different approaches for finding the inner disk radius). We can rewrite this formula in the convenient form:
Tab. 3 (models 7-9) shows theoretically derived values of for the truncation radii obtained in models 1-3. One can see that these accretion rates approximately correspond to those obtained in the 3D simulations.
Hence, both simulations and theory point to the fact that a star with an equatorial dipole magnetic field of 175 G can truncate the disk at the corotation radius only if the accretion rate is very small, smaller than that found from observations. Therefore, the disk should be truncated at smaller distances from the star to explain the accretion rate.
4.8 Test simulation runs at larger dipole field component (models 4 and 5)
To further resolve this issue, we suggest that the dipole component is larger than that reconstructed by D07 from observations. In fact, the July 2009 observations revealed a dipole component about three times stronger (D10). So, our test simulations may be relevant to some other moments of time, when the dipole field is stronger.
We performed two test simulation runs (models 4 and 5) where the octupolar component is fixed at kG, while the dipole component is twice and thrice as large: G () and G ().
Fig. 15 shows the results of these runs versus the main case. One can see that in all cases the disk stops at approximately same distance from the star (because in all cases). However, the dimensional accretion rate is proportional to and hence it is expected to have about 4 and 9 times higher accretion rate compared with the main case. Exact value also depends on the dimensionless parameter obtained from simulations. Tab. 3 (row 12) shows that the accretion rate in these two cases is sufficiently high. One can see that even the dipole field of G is sufficient in disrupting the disk at far distances from the star, and at the same time to give the observed accretion rate. Another advantage of simulations with a higher dipole component is that the main polar spot strongly dominates which gives a better match with observation of accretion spots in V2129 Oph.
|Model of V2129 Oph||(g cms)||(yr)|
|model 1. dip+oct 3D||1.5||0.33||0.006||rot. equilibrium|
|model 2. dip+oct 3D||1.0||0.22||+0.01||spin-up|
|model 3. dip+oct 3D||0.5||0.11||+0.015||spin-up|
|model 4. dip+oct 3D||1.5||0.165||+0.01||spin-up|
|model 5. dip+oct 3D||1.5||0.11||0.02||rot. equilibrium|
|model 6. dip 3D||1.5||0.00||0.003||rot. equilibrium|
4.9 Calculation of the torque and spin evolution
We also calculate the torque on the star from the disk-magnetosphere interaction. Matter coming from the disk in funnel streams can move faster or slower than the star and hence it creates either a positive or negative azimuthal field on the star, and correspondingly a positive (spin-up) or negative (spin-down) magnetic torque, . In addition, matter brings some positive angular momentum which spins the star up, . We calculated both the magnetic, , and matter, , torques, and found that the torque associated with the magnetic field is about 100 times larger than that associated with matter. This is typical in accretion through the funnel streams (see also Romanova et al. 2002, 2003). We calculated the torque for models 1-6 and the results are presented in Table 4. We observed from the simulations that the torque is small and it wanders around zero taking both positive and negative values in models 1, 5 and 6 (see Tab. 4). In models 2-4, the torque is slightly positive and the star spins up. We calculated the dimensional torque using corresponding reference values from Tab. 2.
Now, we can estimate the time-scale of spin evolution. The angular velocity of V2129 Oph is (at its period of days) , its angular momentum is , where . Taking and different values of for different models from Tab. 4, we obtain the characteristic time-scale as (see Tab. 4). One can see that the time-scale in all models is longer than the age of V2129 Oph which is estimated to be years (see D07). Hence, the torque observed in simulations is not sufficient to spin down the star to its present slow rotation.
If the star has lost it spin in the past due to disk-magnetosphere interaction, then we should suggest that in the past the accretion rate was higher (to ensure rapid spin-down), and the star’s magnetic field was stronger (to ensure disk truncation at a large distance from the star). In fact it is expected that in the past, V2129 Oph has been fully convective and its magnetic field has been amplified due to the dynamo mechanism. The “propeller” mechanism can be responsible for the efficient spinning down at this stage (e.g., Romanova et al. 2005; ?). The stellar winds can also carry angular momentum out at different stages of star’s evolution (e.g., Matt & Pudritz 2005, 2008; Cranmer 2008, 2009).
On the other hand, the July 2009 observations by D10 did show that the magnetic field of the star is stronger than in 2005 observations, and maybe the field of V2129 Oph is higher on an average than in the June 2005 observations. Note, however, that the torque is still somewhat small, even when the dipole field is three times stronger (see Tab. 4 for model 5).
Note that if a torque at the present epoch is small, and a star is not in the rotational equilibrium, then the truncation radius can be at any place (not at the corotation radius). In such a case one of the main condition which should be satisfied in comparisons of models with observations is the similarity of the spot shapes and positions.
5 Conclusions and Discussions
We performed global 3D simulations of disk accretion onto a star with parameters close those of T Tauri star V2129 Oph. The large-scale magnetic field of V2129 Oph is approximated by slightly tilted dipole and octupole magnetic moments with polar fields of 0.35 kG and 1.2 kG respectively. Test simulations for the case of a pure dipole field, smaller magnetospheres and a larger dipolar field component were performed for comparison. Below, we summarize our findings and compare our results with the D07 observations/model:
1. The field lines connecting the disk and the star inflate and form magnetic towers above and below the disk. The potential (vacuum) approximation is still valid inside the Alfvén (magnetospheric) radius, , where the magnetic stress dominates over the matter stress. At larger distances, the magnetic field distribution strongly departs from the potential one (see Fig. 6).
2. Simulations show that the disk is disrupted by the dipole component of the field and matter flows towards the star in two funnel streams. The flow in the stream is redirected by the octupole component with part of matter flowing towards the high-latitude pole, and another part into the low-latitude octupolar belt. The polar spots are located at high latitudes even if the disk is truncated close to the star (say, at ).
3. In the dipole case the spots are latitudinally-elongated, they have a typical crescent shapes and are located at lower co-latitudes (see also R04). The polar chromospheric spot observed in V2129 Oph (D07) is much closer in shape and position to the round, high-latitude polar spot obtained in the dipole+octupole model than the crescent-shaped spot obtained in the pure dipole model.
4. If the disk is disrupted at , like in our main case (model 1), then the accretion rate obtained from the simulations is about an order of magnitude smaller than that obtained from observations. If the disk is disrupted at (model 3), then the accretion rate has a good match, but the low-latitude octupolar belt spots strongly dominate, which is in contrast with observations. However, if the disk is disrupted at a somewhat larger distance from the star, (model 2), then the brightest part of the spot is located in the polar region, while the octupolar ring spots are much weaker. In this case the accretion rate is smaller than the observed accretion rate, but it is within the error bar given by D10 (see Tab. 3). Hence, model 2 can possibly explain observations of V2129 Oph. In that case, future, more sensitive observations should reveal the arcs of the octupolar belt spots.
5. Both simulations and theoretical estimations point to the fact that a weak dipole field of G cannot disrupt the disk at large distances from the star, unless the accretion rate is very low. In a test simulation with a dipole field twice as strong ( G, model 4), the accretion rate matches the observed accretion rate and the polar spots dominate. This model explains observations better than the other models.
6. We observed that the main (polar) accretion spot is shifted relative to the main (octupolar) magnetic pole (at the phase of 0.15) towards the dipole magnetic pole (see fig. 9). Such a shift may be connected with the fact that matter flows toward the dipole magnetic pole and is only partially redirected towards the octupolar pole. The dipole and octupole poles differ in phases, and this may lead to the observed phase difference. On the other hand, the funnel stream can be shifted forward/backward relative to the stellar surface and the phase of the accretion spot may differ from that of the magnetic pole (Romanova et al., 2003, 2004a).
7. The torque on the star is small and corresponds to a small accretion rate. The time-scale of spin evolution ( years) is longer than the estimated age of V2129 Oph of years. Hence, disk-magnetosphere interaction cannot be responsible for the star’s spin evolution at present epoch. Hence, we conclude that the star probably lost most of its angular momentum during the early stages of its evolution, for example, at the stage when a star has been fully convective and had a stronger, dynamo-generated magnetic field. The “propeller” mechanism could be responsible for the rapid spin-down (Romanova et al., 2005; Ustyugova et al., 2006). In addition, at all stages of evolution angular momentum flows into stellar winds (e.g. Matt & Pudritz 2005, 2008).
Among the different simulation runs presented here we prefer those runs where the disk is disrupted at a large distance from the star, within , in which the star is expected to be in the rotational equilibrium state (zero torque on average) (e.g., Long, Romanova & Lovelace 2005). Such a state is probable in slowly rotating CTTSs (Königl, 1991; ?; Romanova et al., 2002). However, the low values of torque obtained in our simulations raise the question of whether V2129 Oph is in the rotational equilibrium state, and hence the position of the inner disk is less important for the analysis.
A new set of observations of V2129 Oph done in July 2009 have shown that the octupole and dipole components are about 1.5 and 3 times larger compared with the June 2005 observations discussed in this paper (D10). If the field of V2129 Oph varies due to the dynamo mechanism, then (at a steady accretion rate) the disk will move closer to the star during periods when the field is weaker, and further away when it is stronger. Then, model 2 may explain the June 2005 observations, and models 4 or 5 (with dipole moments twice and thrice as large) the July 2009 observations. If the dipole component is large on an average, say, about (0.7-1)kG, then the accretion is sufficiently high to explain observations.
Jardine, Gregory & Donati (2008) investigated the properties of the funnel streams in V2129 Oph for a fixed magnetic field geometry. They chose a much higher accretion rate in their model M/yr, which is a typical accretion rate for mildly accreting CTTSs (see also D07, Eisner et al. 2005). The density they obtained in the funnel streams was about a hundred times higher compared with our model. At this accretion rate, however, the disk will be truncated very close to the star unless the dipole magnetic field is even higher than discussed above.
There is an interesting possibility that most of the disk matter may accrete directly through the equatorial plane due to magnetic interchange instabilities (Romanova, Kulkarni & Lovelace, 2008; Kulkarni & Romanova, 2008, 2009), while a smaller amount of matter climbs the magnetosphere and forms accretion spots. This model can possibly explain the difference between the low accretion rates obtained by D10 from spectral lines, and the much higher accretion rates derived from the disk measurements (e.g., Eisner et al. 2005).
We performed similar 3D MHD simulations of accretion onto another T Tauri star, BP Tau (Long et al., 2010b) where the dipole and octupole components of 1.2 kG and 1.6 kG strongly dominate over other multipoles (Donati et al., 2008). In that case the dipole component is much stronger than in V2129 Oph, and there is a better match between the simulated and observed accretion rates. Both V2129 Oph and BP Tau are examples of stars with complex fields where the dipole component determines the disk-magnetosphere interaction, while octupolar component shapes the spots.
In some T Tauri stars the magnetic field is more complex than in V2129 Oph, like in CTTSs CVCha and CRCha (Hussain et al., 2009). Modeling of such stars requires incorporation of even higher order multipoles, which can be done in future research.
Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and the NASA Center for Computational Sciences (NCCS) at Goddard Space Flight Center. The authors thank A.V. Koldoba and G.V. Ustyugova for the earlier development of the codes, S.A. Lamzin and R.V.E. Lovelace for discussions, and the referee for many useful questions and comments. The research was supported by NSF grant AST0709015, and funds of the Fortner Endowed Chair at Illinois. The research of MMR was supported by NASA grant NNX08AH25G and NSF grant AST-0807129.
- Aly, J.J. A&A, 86, 192
- Agapitou, V., Papaloizou, J. C. B., 2000, MNRAS, 317, 273
- Basri G., Marcy G.W., Valenti J.A., 1992, ApJ, 390, 622
- Beristain, G., Edwards, S., Kwan, J., 1998, ApJ, 499, 828
- Bessolaz N., Zanni C., Ferreira J., Keppens R., Bouvier J., 2008, A&A, 478, 155
- Brickhouse N.S., Cranmer S.R., Dupree A.K., Luna G.J.M., Wolk S., 2010, ApJ, 710, 1835
- Calvet N., Gullbring E., 1998, ApJ, 509, 802
- Cameron, A., Campbell, C, 1993, A&A, 274, 309
- Cranmer S. R., 2008, ApJ, 689, 316
- Cranmer S. R., 2009, ApJ, 706, 824
- Donati J.-F., Collier Cameron A. C., 1997, MNRAS, 291,1
- Donati J.-F., Collier Cameron A., Hussain G.A.J., Semel M., 1999, MNRAS, 302, 437
- Donati, J.-F., Jardine, M. M., Gregory, S. G., Petit, P., Bouvier, J., Dougados, C., Ménard, F., Cameron, A. C., Harries, T. J., Jeffers, S. V., Paletou, F. 2007, MNRAS, 380, 1297 (D07)
- Donati, J.-F., Jardine, M. M., Gregory, S. G., Petit, P., Paletou, F., Bouvier, J., Dougados, C., Mènard, F., Cameron, A. C., Harries, T. J., Hussain, G. A. J., Unruh, Y., Morin, J., Marsden, S. C., Manset, N., Auriére, M., Catala, C., Alecian, E. 2008, MNRAS, 386, 1234
- Donati J.-F., Bouvier J., Walter F.M., Gregory, S. G., Skelly, M.B., Hussain, G. A. J., Flaccomio, E., Argiroffi, C., Grankin, K.N., Jardine, M. M., Menard, F., Dougados, C., Romanova, M.M. & the MaPP collaboration,, 2010, MNRAS, in preparation (D10)
- Eisner J., Hillenbrandt L., White R., Akeson R., Sargent A., 2005, ApJ, 623, 952
- Elsner R. F., Lamb F. K., 1977, ApJ, 215, 897
- Fang, M., van Boekel, R., Wang, W., Carmona, A., Sicilia-Aguilar, A., Henning, Th. 2009, A&A, 504, 461
- Ghosh P., Lamb F. K., 1978, ApJ, 223, L83
- Ghosh P., Lamb F. K., 1979, ApJ, 232, 259
- Gregory S. G., Matt S. P., Donati J.-F., Jardine M., 2006, MNRAS, 371, 999
- Gregory S. G., Matt S. P., Donati J.-F., Jardine M., 2008, MNRAS, 389, 1839
- Gregory S. G., Jardine M., Gray, C.G., Donati J.-F., 2010, Rep. Prog. Phys., in press
- Hartmann, L., Hewett, R., Calvet N., 1994, ApJ, 426, 669
- Hussain G. A. J., Collier Cameron A., Jardine M. M., et al., 2009, MNRAS, 398, 189
- Jardine M. M., Wood K., Collier Cameron A., Donati J.-F., Mackay D.H., 2002, MNRAS, 336, 1364
- Jardine M. M., Collier Cameron A., Donati J.-F., Gregory S.G., Wood K., 2006, 367, 917
- Jardine M. M., Gregory S.G., Donati J. -F., 2008, MNRAS, 386, 688
- Johns-Krull C., Valenti J.A., Koresko C., 1999, ApJ, 516, 900
- Johns-Krull C. M., 2007, ApJ, 664, 975
- Kato Y., Hayashi M.R., Matsumoto R., 2004, ApJ, 600, 338
- Koldoba A. V., Romanova M. M., Ustyugova G. V., Lovelace R. V. E. 2002, ApJ, 576, L53
- Koldoba A. V., Ustyugova G. V., Romanova M. M., Lovelace R. V. E. 2008, MNRAS, 388, 357
- Königl, A. 1991, ApJ, 370, L39
- Kulkarni A. K., Romanova M. M. 2005, ApJ, 633, 349
- Kulkarni, A. K., Romanova, M. M. 2008, MNRAS, 386, 673
- Kulkarni, A. K., Romanova, M. M. 2009, MNRAS, 398, 701
- Lamzin S. A., 1995, A&A, 295, L20
- Long M., Romanova M. M., Lovelace R. V. E., 2005, ApJ, 634, 1214
- ———. 2007, MNRAS, 374, 436
- ———. 2008, MNRAS, 386, 1274
- Long M., Romanova M.M., Lamb F.K., 2010a, MNRAS, in press, (arxiv:0911:5455)
- Long, M., Romanova, M.M., Kulkarni, A.K., Donati, J.-F., 2010b, MNRAS, in press
- Lovelace, R.V.E., Romanova, M.M., Bisnovatyi-Kogan, G.S., 1995, MNRAS, 275, 244
- Lynden-Bell D. 1996, MNRAS, 279, 389
- Matt, S., Pudritz, R.E. 2005, ApJ, 632, L135
- Matt, S., Pudritz, R.E. 2008, ApJ, 678, 1109
- Mohanty, S., Jayawardhana, R. & Basri, G. 2005, ApJ, 626, 498
- Mohanty S., Shu F.H. 2008, ApJ, 687, 1323
- Romanova M. M., Kulkarni A. K., Lovelace R. V. E., 2008, ApJ, 673, L171
- Romanova M. M., Ustyugova G. V., Koldoba A. V., Chechetkin, V.M., Lovelace R. V. E., 1998, ApJ, 500, 703
- Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E., 2002, ApJ, 578, 420
- ———. 2004a, ApJ, 610, 920
- ———. 2004b, ApJ, 616, L151
- ———. 2005, ApJ, 635, L165
- ———. 2009, MNRAS, 399, 1802
- Romanova M. M., Ustyugova G. V., Koldoba A. V., Wick J. V., Lovelace R. V. E., 2003, ApJ, 595, 1009
- Safier, P. N., 1998, ApJ, 494, 336
- Shevchenko, V., Herbst, W. 1998, AJ, 116, 1419
- Ustyugova, G. V., Koldoba, A. V., Romanova, M. M., Lovelace, R. V. E., 2006, ApJ, 646, 304
- Ustyugova, Lovelace, R. V. E., Romanova, M.M., Li, H., Colgate, S.A., 2000, ApJ, 541, L21