text:
# New Tests of Magnetospheric Accretion in T Tauri Stars

## 1 Introduction

Classical T Tauri stars (CTTSs) appear to be roughly solar mass pre-main sequence stars still surrounded by disks of material which actively accrete onto the central star. It is within such disks that solar systems form. Understanding the processes through which young stars interact with and eventually disperse their disks is critical for understanding the rotational evolution of stars and the formation of planets. A key question is to understand how young stars can accrete large amounts of disk material with high specific angular momentum, yet maintain rotation rates that are observed to be relatively slow (e.g. Hartmann & Stauffer 1989, Edwards et al. 1994). Current models state that this process proceeds via magnetospheric accretion: stellar magnetic fields truncate the inner disk and channel accreting material onto the stellar surface, preferentially at high latitude. More generally in astrophysics, magnetized accretion onto a central source is a common theme in various contexts including accretion onto white dwarfs (e.g. AM Her stars – Wickramasinghe & Ferrario 2000), pulsars (pulsating X-ray sources – Ghosh & Lamb 1979), and possibly black holes at the center of active galactic nuclei (e.g. Koide, Shibata, & Kudoh 1999). Detailed tests of magnetospheric accretion theory are necessary for understanding CTTSs as well as these other classes of objects.

Support for magnetospheric accretion in CTTSs is significant. Current models can account for the relatively slow rotation of most CTTSs (Camenzind 1990; Königl 1991; Shu et al. 1994; Paatz & Camenzind 1996). Studies of the spectroscopic and photometric variability of CTTSs often suggest magnetically controlled accretion (e.g. Bertout, Basri, & Bouvier 1988; Bertout et al. 1996; Herbst, Bailer–Jones, & Mundt 2001; Alencar, Johns–Krull, & Basri 2001) with the magnetic axis inclined to the rotation axis in some cases (e.g. Kenyon et al. 1994, Johns–Krull & Basri 1995, Bouvier et al. 1999). Models of high resolution Balmer line profiles, computed in the context of magnetospheric accretion, reproduce many aspects of observed line profiles (Muzerolle, Calvet, & Hartmann 1998, Muzerolle et al. 2000). T Tauri stars are observed to be strong X-ray sources indicating the presence of strong magnetic fields on their surfaces (see Feigelson & Montmerle 1999 for a review), and CTTSs are observed to have strong surface magnetic fields (Guenther et al. 1999; Johns–Krull et al. 1999b; Johns–Krull et al. 2001), and strong magnetic fields have been observed in the formation region of the He I emission line at 5876 Å (Johns–Krull et al. 1999a; Johns–Krull & Valenti 2000), which is believed to be produced in a shock near the stellar surface as the disk material impacts the star.

Despite these successes, open issues remain. Current theoretical models assume a dipole geometry for the stellar magnetic field; however, Johns–Krull et al. (1999a) and Johns–Krull and Valenti (2001) show that the surface fields on BP Tau and TW Hya are not dipolar. On the other hand, it is expected that the dipole component of the field should dominate at distance from the star where the interaction with the disk is taking place, so this may not contradict current theory. Even in the case of the complex magnetic topology of the Sun, the dipole component appears to become dominant at $`2.5R_{}`$ or closer (e.g. Luhmann et al. 1998). For expected disk truncation radii of 3 – 6 $`R_{}`$ in CTTSs, this suggests the dipole component will govern the stellar interaction with the disk. On the other hand, Stassun et al. (1999) find no correlation between rotation period and IR signatures of circumstellar disks in a sample of 254 stars in Orion, leading them to question the validity of magnetically mediated disk locking which is intrinsic to magnetosphereic accretion theory. Bertout and Mennessier (1996) review points in favor of both the magnetospheric accretion model and the original axisymmetric, equatorial boundary layer model (based on the work of Lynden–Bell & Pringle 1974) which came before. While the magnetospheric model remains the favorite in the current literature, Bertout and Mennessier (1996) conclude that the existing data does not favor either model. Thus, the overall success of magnetospheric accretion as the best description of accretion onto CTTSs remains in some doubt.

The current equilibrium theories of magnetospheric accretion predict specific relationships between the stellar mass ($`M_{}`$), radius ($`R_{}`$), rotation period ($`P_{rot}`$), magnetic field ($`B_{}`$), and mass accretion rate ($`\dot{M}`$) (e.g. Johns–Krull et al. 1999b). Unfortunately, observation of all these quantities are only available for a handful of stars, with magnetic field measurements the most lacking. One could look for correlations among the other quantities to test the theories. Ghosh (1995) investigated the rotation and accretion properties of a number of T Tauri stars (TTSs) by assuming a dynamo relation for TTS magnetic fields in terms of Rossby number (convective turnover time divided by rotation period) since there were almost no magnetic field measurements for TTSs available at the time. Ghosh (1995) found that TTSs generally cluster below the limiting rotation rate as is to be expected; however, no clear correlations were seen in the data. Johns–Krull, Valenti, and Linsky (2001) have shown that TTS Rossby numbers generally place them in the region of saturated magnetic activity compared to main sequence and evolved stars, so it is not clear that parameterizing the stellar field in terms of the Rossby number is appropriate for these stars. Furthermore, TTSs do not show a large range of magnetic field properties as described below. While we are unaware of any other systematic studies along this line, Muzerolle, Calvet, and Hartmann (2001) note that current theory predicts a correlation between rotation period and mass accretion rate which has not been observed. Muzerolle et al. (2001) suggest that variations in the stellar magnetic field strength from star to star may account for the lack of correlation.

Johns–Krull et al. (2001) show that the mean photospheric magnetic fields of 6 TTSs range only between 2.1 and 2.7 kG, suggesting a surprising uniformity to the field strengths on these stars. The field on T Tau obtained by Guenther et al. (1999) also falls in this range. These observations are of the magnetic field averaged over the entire stellar surface, while for the current discussion we are most interested in the field in the actual accretion zones. Johns–Krull et al. (1999a) discovered a magnetic field in the He I 5876Å emission line on BP Tau from observations of circular polarization, arguing that the narrow component of this line likely forms in the postshock region as material accretes onto the stellar surface (see also Hartmann, Hewett, & Calvet 1994 and Lamzin 1995). The field detected is 2.5 kG, which is a lower limit to the true field due to the unknown angle between the field lines and the line of sight. Johns–Krull et al. (2001) also present similar spectropolarimetric observations of 4 CTTSs over 6 nights, finding fields in the He I line formation region that vary smoothly with time, probably due to stellar rotation which suggests the fields participating in the accretion flow are stable at least on a rotational time scale (see also Valenti, Johns–Krull, & Hatzes 2001). In all cases, the field determined from the He I polarization is equal to or lower than the mean photospheric field, as is expected if the fields are in fact the same but the He I measurement is reduced as a result of field line inclination to the line of sight.

As mentioned above, the theoretical relationships predicted by magnetospheric accretion theory involve many parameters, not just rotation period and mass accretion rate. Since all the exisiting observational data suggests relative uniformity of surface magnetic field strengths on CTTSs (at least in Taurus where the most detailed studies have taken place), we have examined several sources of data from the literature to look for the expected correlations assuming the field strength is in fact constant from star to star. This constancy of the magnetic field may represent some limit to the dynamo operating in these young stars. Certainly, the details of such a dynamo are far from well understood; however, Johns–Krull et al. (2001) show that within the context of other active stars, TTS are expected to show saturated magnetic activity based on their rotation and convective parameters, perhaps giving further justification to assuming a constant magnetic field. In §2 we review the predicted relationships between the various stellar and accretion parameters. In §3 we examine literature compilations of the requisite parameters for testing the predicted correlations. Section 4 presents the analysis of this data. In §5 we discuss our results, and our conclusions are given in §6.

## 2 Magnetospheric Accretion Predictions

Johns–Krull et al. (1999b) examined three analytic theories (Königl 1991, Cameron & Campbell 1993, Shu et al. 1994) of magnetospheric accretion, presenting equations for the stellar magnetic field strength as a function of $`M_{}`$, $`R_{}`$, $`P_{rot}`$, and $`\dot{M}`$. These studies are outgrowths of the pioneering work done by Ghosh and Lamb (1979) on accreting neutron stars. The work of Königl (1991) and Cameron and Campbell (1993) are very closely related to that of Ghosh and Lamb (1979) in that these papers only consider accretion onto the star. The Shu et al. (1994) work is different in that it is an integrated theory of accretion onto the central star plus a wind which is driven off the disk at the truncation point (the X-region). While coefficients representing such things as the efficiency of magnetic coupling to the disk vary from one study to the other, the predictions from Königl (1991) and Shu et al. (1994) depended in the same way on the above stellar and accretion parameters for an assumed dipolar magnetic field. Holding the stellar magnetic field constant, Königl (1991) and Shu et al. (1994) predict that

$$R_{}^3M_{}^{5/6}\dot{M}^{1/2}P_{rot}^{7/6}.$$
(1)
In the case of Cameron and Campbell (1993), the relationship is

$$R_{}^3M_{}^{2/3}\dot{M}^{23/40}P_{rot}^{29/24},$$
(2)
which is very similar in its dependencies to Königl (1991) and Shu et al. (1994). In fact, scatter plots created using equation (2) are virtually identical to those created using equation (1). Therefore, we only show plots based on equation (1); however, the results for equation (2) are cited in the tables and discussed in the text.

In a more detailed analysis of the magnetospheric accretion flow itself, Ostriker and Shu (1995) identified the importance of the notion of trapped flux in the Shu et al. (1994) theory. While still assuming a dipole stellar field, Ostriker and Shu noted that what their theory really predicted was the amount of magnetic flux trapped in the X-region. In their preferred model, this is 1.5 times the magnetic flux which would thread the disk exterior to the truncation point (the co-rotation radius in their model) in an undisturbed dipole field. For the system to remain in equilibrium, this is the flux that must be trapped in the X-region independent of the topology of the magnetic field on the stellar surface. Armed with this information, we can derive a new relationship among the observable parameters which is unique to the Ostriker and Shu (1995) treatment of the Shu et al. (1994) theory.

Re-writing equation (4.1) of Ostriker and Shu (1995) for the radius of the X-point, $`R_x`$, we have

$$R_x=\mathrm{\Phi }_{dx}^{4/7}\left(\frac{\mu _{}^4}{GM_{}\dot{M}_D^2}\right)^{1/7}$$
(3)
where $`\mathrm{\Phi }_{dx}^{4/7}`$ is a constant evaluated by Ostriker and Shu (1995) to be 0.923 in their preferred model, $`\mu _{}`$ is the dipole moment of the star, $`G`$ the gravitational constant, and $`\dot{M}_D`$ is the mass accretion rate through the disk. The general form of this equation is well known from many studies of accretion onto compact objects (see for example White & Stella 1988). In the Shu et al. (1994) theory, the disk mass accretion rate is not equal to the accretion rate onto the star (the quantity measured in the observations described below) since the theory also incorporates a magnetocentrifugally driven wind from the X-region; however, the two are related by a factor of order 2/3.

For a purely dipolar field, the magnetic flux threading the disk exterior to the point $`R_x`$ is

$$\varphi _x=\frac{2\pi \mu _{}}{R_x}.$$
(4)
In the preferred model of Ostriker and Shu (1995) the total flux trapped in the X-region is $`\varphi _t=1.5\varphi _x`$, of which a third participates in the accretion funnel flow of material onto the star (another third is in the dead zone, and the final third is part of the wind). Thus, the magnetic flux threading the X-region which directs disk material onto the star is $`0.5\varphi _x`$ independent of the exact topology of the field at the stellar surface, as long as the dipole dominates in the X-region. Equating this magnetic flux to that on the star in the accretion zones, we have

$$B_{}A_{}=0.5\varphi _x=\frac{\pi \mu _{}}{R_x}$$
(5)
where $`A_{}`$ is the area on the surface of the star on which the accreting disk material falls.

If we recall that $`R_x`$ in the Shu et al. (1994) theory is the co-rotation radius, where the Keplerian angular velocity in the disk is equal to the angular velocity of stellar rotation, we can solve equation (5) for $`\mu _{}`$, substitute the result into equation (3) and get

$`B_{}A_{}`$ $`=`$ $`2^{1/2}\pi ^{1/2}G^{1/2}\mathrm{\Phi }_{dx}M_{}^{1/2}\dot{M}_D^{1/2}P_{rot}^{1/2}`$ (6)
$`=`$ $`8.88\times 10^{24}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{1/2}\left({\displaystyle \frac{\dot{M}_D}{10^7M_{}\mathrm{yr}^1}}\right)^{1/2}\left({\displaystyle \frac{P_{rot}}{1\mathrm{day}}}\right)^{1/2}\mathrm{G}\mathrm{cm}^2`$
after some manipulation. In this paper, we assume $`B_{}`$ does not vary from star to star and $`A_{}=4\pi R_{}^2f_{acc}`$ where $`f_{acc}`$ is the filling factor of accretion zones, so equation (6) implies

$$R_{}^2f_{acc}M_{}^{1/2}\dot{M}^{1/2}P_{rot}^{1/2}$$
(7)
as a unique prediction of the Shu et al. (1994) theory as detailed by Ostriker and Shu (1995) and expanded here to incorporate a non-dipole field topology near the stellar surface. In equation (7), $`\dot{M}`$ is the mass accretion rate onto the star and $`\dot{M}=\alpha \dot{M}_D`$ where $`\alpha 2/3`$ since the wind carries away about one third of the disk accretion rate (Shu et al. 1994). Hartigan, Edwards, and Ghandour (1995) attempted to observationally determine the ratio of mass loss to mass accretion rate, finding a value of 0.01; however, more recent estimates of this ratio are $`0.1`$ (see review of Königl & Pudritz 2000) so the Shu et al. (1994) value seems quite reasonable given the current state of the observations.

In equation (7) above, the filling factor, $`f_{acc}`$, is the filling factor of the accretion spots on the stellar surface, which is not the same as the filling factor of the magnetic field. In the theories which use a strict dipole field geometry in which the magnetic axis is aligned with the rotation axis (the typical assumption), the magnetic filling factor is 1.0, but accretion onto the star occurs only in axisymmetric rings located at high latitude. Thus the filling factor of accretion zones is expected to be significantly less than 1.0, and its value is strongly influenced by the small scale topology of the magnetic field at the stellar surface. On the other hand, it is the bright accretion shock emission produced in these regions which is believed to give rise to the UV continuum excess and optical veiling emission seen in CTTSs (e.g. Calvet & Gullbring 1998). As a result, it is possible to observationally determine $`f_{acc}`$, and this has been done for reasonably large samples of stars by Valenti, Basri, and Johns (1993) and Calvet and Gullbring (1998). As discussed further below, accurate estimates of $`f_{acc}`$ are required to test the correlations predicted by equation (7), and these estimates can be difficult to make.

In this paper we focus on the above mentioned theories since they are largely analytic and permit the derivation of relatively simple expressions which show the dependence of stellar and accretion parameters on one another in the equilibrium situation. These expressions can be tested in a rather straightforward way using existing observational data. This is not to suggest that the above studies are the only such studies of magnetospheric accretion in CTTSs. There are a number of others, some of which derive their results from numerical studies of the governing equations (e.g. Paatz & Camenzind 1996), while some focus on the instabilities that can result between the star, its magnetosphere, and the disk which may launch winds and jets (e.g. Hayashi, Shibata, & Matsumoto 1996; Goodson, Winglee, & Böhm 1997; Miller & Stone 1997; Goodson, Böhm, & Winglee 1999; Ferreira, Pelletier, & Appl 2000). These studies are not directly testable in the fashion explored in this paper; however, important tests of many of these works may be possible with variability studies (see Goodson et al. 1999). Finally, we note that CTTSs are variable as has been mentioned in §1. This variability is not accounted for in the equilibrium theories examined in this paper, but is not necessarily evidence that these theories should be abandoned. The theories may be perfectly adequate to describe the average conditions of CTTSs. In this case, we expect to see the predicted correlations in the data, and the variability (as well as measurement error) should produce scatter about the mean relations. It is these correlations we are looking for here.

## 3 Data Used

The goal of this paper is test the proportionalities predicted in equations (1), (2), and (7) by using data from the literature to look for the predicted correlations. For consistency within each of the studies we employ, we use the stellar and accretion parameters from each individual study when testing for the predicted correlations instead of trying to come up with a single best estimate for each of the relevant parameters. The rotation periods are all taken from photometric monitoring studies which do not depend on these stellar parameters, so we apply the detected periods to all the accretion data sets.

Mass accretion rates onto large samples of CTTSs have been estimated by Valenti et al. (1993 – hereafter VBJ); Hartigan, Edwards, and Ghandour (1995 – hereafter HEG); Gullbring et al. (1998 – hereafter GHBC); Calvet and Gullbring (1998 – hereafter CG98); and White and Ghez (2001 – hereafter WG01). Of these papers, only VBJ and CG98 estimate the filling factor, $`f_{acc}`$, of accreting material on the stellar surface, which is required to test the relationship given in equation (7). CG98 do not separately estimate, $`T_{eff}`$, $`M_{}`$, and $`R_{}`$, instead using the values tabulated by GHBC which are reproduced in the tables here. For the 4 continuum stars (CW Tau, DG Tau, DL Tau, and DR Tau) from CG98 which are used in this study, CG98 assume $`R_{}=2R_{}`$ and $`M_{}=0.5M_{}`$, which we use here. In addition, CG98 do not tabulate $`\dot{M}`$ specifically, but the mass loss rate can be calculated from their equation (11) and the quantities tabulated in the paper. We note that in CG98, they do not perform a detailed fit to the data presented in GHBC, instead adopting the values from their grid of models which gives the best fit to the data as estimated by eye. However, this should be adequate for our purpose of looking for general correlations among different parameters. Rotation period data comes from the following papers in order of preference: Bouvier et al. 1995, Bouvier 1990, Bertout et al. 1996, and Bouvier et al. 1993. Tables 1 through 5 give the final compilation of stars and their relevant parameters for the 5 studies examined here. In a few cases (noted in the tables), required quantities such as the stellar radius are not tabulated in the original studies, but can be estimated from the tabulated stellar luminosity and effective temperature. For those studies which only tabulate spectral type, the spectral type – effective temperature calibration of Johnson (1966) is used in this paper to estimate $`T_{eff}`$.

Some of the stars used in this study are now known to be binary, while this was not known to be the case in earlier studies. In addition, pre-main sequence evolutionary tracks have been refined over time, which influences the derivation of stellar parameters. For many of the binaries, their separation is large enough (DK Tau, GK Tau, RW Aur) that they do not fall in the slit used by HEG, while WG01 explicitly take into account the binary nature of these stars. For the studies based on low to medium resolution spectroscopy (VBJ, GHBC, CG98), the slits used were wide enough to encompass both members of the binary systems. However, it is generally the case (DF Tau, DK Tau, GG Tau A, RW Aur) that the primary star dominates the optical and short wavelength light (Mathieu 1994; Ghez, White, & Simon 1997; Duchene et al. 1999; WG01) from which the accretion parameters are derived. As a result, the accretion rates and filling factors apply to the primary in the case of those stars with originally unknown companions. Inclusion of the WG01 results in this study will also serve to see if accurate accounting for binarity significantly changes the results. The use of different pre-main sequence tracks and differences in the extinction corrections used (as well as unknown binaries) between the different studies results in small (generally less than a factor of 2) changes in the stellar mass and radius from one study to the other, but these differences are small compared to the full range of parameters explored in the samples. We choose to use the samples as is since changes made in one parameter can also affect others (such as the mass accretion rate), so here we hope to maintain consistent sets of parameters from one study to the next as we look for general correlations among the parameters.

In addition to observational issues related to the influence of binaries, it is fair to consider the theoretical influence of binary interactions. This paper is concerned with the star-disk interaction, and the theories examined here have the disk accretion rate as a free parameter. Therefore, if CTTSs in binaries have different disk accretion rates compared to single CTTSs, this difference should naturally be accounted for. Theoretical effects studied in CTTS binary systems inlclude disk truncation and circumbinary replenishment (e.g. Artymowicz & Lubow 1994, 1996) and changes to disk accretion rates (e.g. Ostriker, Shu, & Adams 1992), but most authors consider the star-disk interaction to be the same as for single CTTSs (e.g. Armitage, Clarke, & Tout 1999). Observationally, WG01 find no difference in the accretion rate as a function of binary separation in the range 10 - 1000 AU, and Simon & Prato (1995) find no difference in disk lifetime between binary and single CTTSs. Thus, it is unlikely there is any real difference in the star-disk interaction for CTTSs in relatively wide binaries such as those considered here compared to single CTTSs. As a result, we include both single CTTSs and CTTSs in binaries in the current analysis, but we separate them by plot symbol for clarity. Stars identified as binary in WG01, Mathieu (1994), or Ghez, Neugebauer, and Matthews (1993) are labelled as binary; while stars indicated as single in WG01 are labelled so here.

As mentioned above, VBJ and CG98 are the only studies which determine the filling factor, a necessary quantity for testing equation (7). Generally speaking, the mass accretion rate is

$$\dot{M}=4\pi R_{}^2f_{acc}\rho v,$$
(8)
where $`\rho `$ is the density of the infalling material and $`v`$ is the velocity at which it is coming in. Therefore,

$$f_{acc}=\frac{\dot{M}}{4\pi R_{}^2\rho v}.$$
(9)
HEG, GHBC, and WG01 all give $`\dot{M}`$ and either $`R_{}`$ or the data necessary to derive it. One could either simply take single representative values for $`\rho `$ and $`v`$ or try to estimate them individually for each star. In either case though, the value of $`f_{acc}`$ determined depends explicitly on the mass accretion rate. Therefore, one would naturally expect a good correlation between the left and right side of equation (7) since $`\dot{M}`$ appears on both sides. We have estimated $`f_{acc}`$ from equation (9) and do indeed find good correlations for equation (7), but this cannot be viewed as support for our extended version of the Ostriker and Shu (1995) model, so we do not include these plots for the data of HEG, GHBC, or WG01.

In summary then, we use data from 5 observational studies to look for the correlations predicted in equations (1), (2), and (7). Of these, 3 studies do not independently determine the filling factor of accretion spots, $`f_{acc}`$, and therefore are only suitable to test the relationships of the simple (dipole) accretion theories embodied in equation (1) and (2). These three studies are HEG, GHBC, and WG01. The investigations of VBJ and CG98 do independently estimate $`f_{acc}`$ and are therefore suitable for testing the prediction of our modified Ostriker and Shu (1995) model, as well as the predictions of the simple models. Thus, VBJ and CG98 can be used to test equations (1), (2), and (7).

## 4 Analysis

### 4.1 Simple Accretion Theories

Figures 1a, 2, 3, 4a, and 5 show the results for the data from the 5 studies when looking for the correlation predicted by equation (1), which results from assuming a dipole magnetic field geometry with a constant field strength from star to star following Königl (1991), Shu et al. (1994 – and the unmodified Ostriker and Shu 1995 theory). In each figure, the quantity on the right hand side of the equation (1) is plotted versus the quantity on the left hand side of the equation. For simplicity in the plot units, we have expressed the stellar mass and radius in solar units for each star, and the mass accretion rate in units of $`10^7`$ M yr<sup>-1</sup>. As mentioned above, the scatter plots which result from equation (2) following Cameron and Campbell (1993) look almost identical to those shown in these figures and are therefore not shown separately. Filled circle plot symbols denote single CTTSs, and asterisks denote CTTSs in binary systems. In all the figures, the expected correlation is shown with a dashed line. For those plots which do show a significant correlation (see below), the best fit line is shown in a solid line. In general, the expected relationship for equation (1) is only weakly present, with the data of HEG showing the strongest correlation.

We quantify these results by computing the linear correlation coefficient (also known as Pearson’s $`r`$) and its associated false alarm probability, $`P_f`$ (Press et al. 1986). For each investigation, the associated correlation coefficient and false alarm probability are computed by using the logarithm of the quantities on each axis of the plots, thus giving the linear correlation of the scatter plots actually shown in these log-log plots. This is done to prevent the calculation of the correlation from being overly influenced by one or a few points with very large values on either axes. The resulting correlation coefficients and false alarm probabilities are given in Table 6. It is readily apparent from this table that when testing the correlations predicted by equations (1) or (2), only the HEG data produces a false alarm probability $`0.01`$, which is typically indicative of a real correlation. The data of HEG give $`P_f=0.009`$ in the case of equation (2), suggestive of a significant correlation. The best fit line (solid line in the plot) is determined from a standard least squares analysis, assuming each data point has the same uncertainty in the y-axis quantity. To estimate the uncertainty in the slope (predicted to be 1 by the theory), we then use the best fit line to estimate observed residuals for each point. The standard deviation of these residuals is then assumed to be the $`1\sigma `$ uncertainty of each point in the y-axis quantity, which in turn gives an uncertainty for the slope in the least squares fit. Given the number of potential uncertainties that go into estimating the error bars for the points in each figure, we feel this means of estimating the final uncertainty of the slope from the scatter plots is the least biassed. The slope and its uncertainty are given in Table 6 for that equation \[(1), (2), or (7) discussed below\] which produces the best correlation for each data set. For the case of HEG data and equation (2), the slope is $`0.88\pm 0.29`$ which is insignificantly different from the expected value of 1.0.

### 4.2 Modified Ostriker and Shu

#### 4.2.1 Initial Tests

Figures 1b and 4b show the results for the data from VBJ and CG98 respectively when looking for the correlation predicted by equation (7), which results from using the trapped flux concept of Ostriker and Shu (1995) while abandoning a pure dipole magnetic field geometry and assuming a constant stellar magnetic field strength. The plot symbols have the same meaning as given above, as do the solid and dashed lines. Both of these plots show much stronger correlations than present in the data when testing equations (1) and (2). This is also apparent in the correlation coefficients and false alam probabilities given Table 6. Also given in the table is the best fit slope and uncertainty calculated in the same way as described above. The expected slope is 1.0 and the VBJ data give this slope to within the uncertainties. The CG98 data give a slope substantially below this expected value. This low slope is largely driven by the two stars with the lowest values of $`R_{}^2f_{acc}`$; however, excluding them still gives a slope which is too low at the $`4\sigma `$ level. This gives some concern, but as mentioned before (and discussed further below), CG98 do not actually fit their models to the observations, indicating that their study may not be the best suited for strictly testing the predicted correlations. Nevertheless, it is certainly true that equation (7) produces far better correlations in the data than either equations (1) or (2).

#### 4.2.2 A Closer Look at VBJ and CG98

The data from both VBJ and CG98 appears to support the Ostriker and Shu (1995) theory as extended here to include non-dipole field topologies, while the data from all 5 observational studies at best only weakly supports the other magnetospheric accretion theories. As a result, it is important to look closely at VBJ and CG98 to see just how well they agree with equation (7). We also note from equation (8) that generally, $`\dot{M}`$ depends on $`f_{acc}R_{}^2`$, so there should be some correlation between these quantities. The observations primarily determine the accretion luminosity, $`L_{acc}`$, which is then used to determine the mass accretion rate either directly by applying an equation of the form

$$\dot{M}=\beta R_{}L_{acc}/GM_{}$$
(10)
where $`\beta `$ depends on whether accretion occurs through an equatorial boundary layer or through a magnetospeher; or a model is used in the case of CG98 to fit the accretion luminosity and the resulting parameters give the mass accretion rate. At the same time, the filling factor, $`f_{acc}`$, is usually determined from the equation

$$L_{acc}=F_{acc}4\pi R_{}^2f_{acc}$$
(11)
where $`F_{acc}`$ is the emission per unit area from the accretion zones and is determined by fitting the observations to a model, either a hydrogen slab model in the case of VBJ or a shock model for CG98. Thus, both the determination of $`\dot{M}`$ and $`f_{acc}`$ primarily depend on $`L_{acc}`$ and the value of $`R_{}`$. Therefore, it is important to see how well $`\dot{M}`$ is correlated with $`R_{}`$ (or $`R_{}^2`$, but when looking at the log of the quantities, it does not matter which power is considered when evaluating the correlation between $`\dot{M}`$ and $`R_{}`$) and $`f_{acc}`$ separately, as well as with the product $`R_{}^2f_{acc}`$ to see if the correlations seen in Figures 1b and 4b are driven by these potential correlations. Indeed, CG98 point out that $`f_{acc}`$ appears to be well correlated with $`\dot{M}`$ in their study, making this a real concern. On the other hand, we already know that $`\dot{M}`$ is not too well correlated with $`R_{}`$ since the top two panels in each figure do not show strong correlations.

In table 7 we give the correlation coefficients and resulting false alarm probabilities when comparing log($`\dot{M}^{1/2}`$) to log($`f_{acc}`$), log($`R_{}^2`$) and log($`f_{acc}R_{}^2`$) using the data from VBJ and CG98. These correlations are to be compared with the values in the last two columns of Table 6 which gives the correlation coefficient and false alarm probability when comparing all the parameters relevant to our extended Ostriker an Shu theory \[log($`M_{}^{1/2}\dot{M}^{1/2}P_{rot}^{1/2}`$) versus log($`f_{acc}R_{}^2`$)\]. In the case of VBJ, there is no correlation between the mass accretion rate and the filling factor, some between $`\dot{M}`$ and $`R_{}^2`$ and a real correlation between $`\dot{M}`$ and $`f_{acc}R_{}^2`$, but the full Ostriker and Shu (1995) comparison (Table 6) produces the best correlation. Thus, it appears that the inclusion of the rotation period and the mass do improve the results as predicted by the theory. The case is not quite the same for CG98, with the $`\dot{M}`$ showing the best correlation with $`f_{acc}R_{}^2`$, but the correlation based on equation (7) is nearly as high. It is certainly true that the CG98 shock model is the most physical of the two studies; however, CG98 do not actually fit the observed data, but pick by the eye the closest model from their grid to match the observations. It is also apparent from their Figure 6 that they do not generally match the observations (particularly the slope of the Paschen continuum) as well as VBJ, so it cannot be assumed that their values for various accretion parameters are more accurate.

We can look further at the results of VBJ. We substitute equation (10) into equation (7) to eliminate $`\dot{M}`$ and then use equation (11) to eliminate $`f_{acc}`$. Dropping constants and re-arranging, we then find

$$R_{}^{1/2}L_{acc}^{1/2}F_{acc}P_{rot}^{1/2}.$$
(12)
In the case of VBJ, $`F_{acc}`$ is determined by a model fit and its value is largely driven by the shape of the Paschen continuum and the size of the Balmer jump. Equation (12) above involves parameters which are independently determined to the greatest extent possible within the context of the VBJ study. The accretion luminosity is given by VBJ and we reproduce these in Table 1 for each star. Since we have access to the models used by VBJ, we compute $`F_{acc}`$ by integrating over wavelength the excess flux produced by the best fitting hydrogen slab model for each star and also give these values in Table 1. We note that VBJ do not estimate $`f_{acc}`$ using equation (11). Use of this equation with values of $`F_{acc}`$ from Table 1 to estimate $`f_{acc}`$ produce values which are systematically larger that those given by VBJ. The model of VBJ used the observations to determine the filling factor of the visible hemisphere ($`\pi R_{}^2`$) covered by the boundary layer (accretion zones). The systematic difference between the two estimates of $`f_{acc}`$ reflect the unknown forshortening of the emitting region if it is not located at disk center on the stellar surface. This depends on the exact topology of the accretion and the limb darkening that should be applied. The two estimates of $`f_{acc}`$ are extremely well correlated with each other and using values estimated from equation (11) produces a correlation in Figure 6b very similar to that shown. However, there are still worries that $`f_{acc}`$ and $`\dot{M}`$ may show an aritifical correlation, so testing the relationship given in equation (12) is a better test of our modified version of the Ostriker and Shu (1995) theory. Figure 6 shows the two sides of equation (12) plotted against one another, showing an excellent correlation. The resulting correlation coefficient is $`r=0.92`$ with an associated false alarm probability of $`P_f=1.5\times 10^7`$. A least squares analysis finds a best fit line (shown in solid) in the log-log plot with a slope of $`1.17\pm 0.12`$ which is not significantly different from the predicted slope of 1.0 (shown in the dashed line). We are unable to repeat this analysis for CG98 since we do not have access to their models. We can do a similar exercise for the data showing the highest correlation for the other magnetospheric accretion theories. Examining Table 6, the data from HEG may show a possible correlation when used to test equation (2). On the other hand, the HEG data shows a correlation between $`\dot{M}`$ and $`R_{}`$ with $`r=0.75`$ and $`P_f=3.1\times 10^4`$. Again, using equation (10) to eliminate $`\dot{M}`$ from equation (2) gives (after dropping constants)

$$R_{}^{97/40}M_{}^{11/120}L_{acc}^{23/40}P_{rot}^{29/24}.$$
(13)
HEG tabulate $`L_{acc}`$ for each star which we also give in Table 2. Using the HEG data to test equation (13) shows only a very weak correlation, with $`r=0.51`$ and $`P_f=0.031`$. The data are shown in Figure 6 and the slope is $`0.88\pm 0.36`$.

In summary, it appears both VBJ and CG98 are consistent with our extended version of the Ostriker and Shu (1995) theory, but the correlation given by CG98 may be the result of the generally expected correlation between $`\dot{M}`$ and $`f_{acc}R_{}^2`$. The data from VBJ actually show the best correlation when considering all the parameters required by the theory, and the data from this study shows an excellent correlation with the theory when one tries to remove the dependence of various parameters on each other. There is only weak support for the more general magnetospheric accretion theories of Königl (1991), Cameron and Campbell (1993), Shu et al. (1994), and Ostriker and Shu (1995) without the extensions described above for considering non-dipole field geometeries.

#### 4.2.3 Filling Factor Versus Stellar Radius

The above results suggest the Ostriker and Shu (1995) formulation can be shown to be consistent with the data if the strict dipole field geometry is dropped and it is assumed that the stellar magnetic field is basically constant from one CTTS to another. Ostriker and Shu (1995) state that one prediction of their work is that smaller CTTSs should have smaller filling factors of accretion spots on their surface. Thus we might expect $`f_{acc}`$ and $`R_{}`$ to be directly related to one another. On the other hand, the relationship given in equation (7) above clearly predicts that $`f_{acc}`$ and $`R_{}`$ should be inversely related to one another. This apparent contradiction deserves some comment.

This result can be understood in terms of our equation (6), which was derived assuming a dipole field geometry in the vicinity of the truncation point of the disk, but which maps back to an arbitrarily complex field topology at the stellar surface. This equation certainly holds for a pure dipole as used by Ostriker and Shu (1995). Remembering that $`A_{}=4\pi R_{}^2f_{acc}`$, we can re-write equation (6) as

$$B_{}R_{}^2f_{acc}=C_1M_{}^{1/2}\dot{M}^{1/2}P_{rot}^{1/2}$$
(14)
where $`C_1`$ is a constant, $`C_1=2^{5/2}\pi ^{1/2}G^{1/2}\mathrm{\Phi }_{dx}`$. If we now assume the field is a pure dipole as Ostriker and Shu (1995) do in their default model, then

$$B_{}=\mu _{}/R_{}^{3+ϵ}$$
(15)
where $`\mu _{}`$ is the dipole moment of the star and $`ϵ`$ is a small positive number ($`0<ϵ<1`$) which takes into account the fact that as the stellar radius decreases, the relevant field lines which thread the disk truncation region in the Ostriker and Shu (1995) model are closer to the pole where the field strength is stronger. Putting this into equation (14) above yields:

$$\mu _{}f_{acc}=R_{}^{1+ϵ}C_1M^{1/2}\dot{M}^{1/2}P_{rot}^{1/2}$$
(16)
after some re-arrangement. Holding other things constant, you expect a direct relationship between $`f_{acc}`$ and $`R_{}`$ as stated by Ostriker and Shu (1995). However, in our extension of the Ostriker and Shu (1995) theory to non-dipole field geometries at the stellar surface, equation (7) above clearly predicts an inverse relationship between $`f_{acc}`$ and $`R_{}`$ if the magnetic field strength at the stellar surface is held constant along with the stellar mass, rotation period, and the mass accretion rate. In a log-log plot of $`f_{acc}`$ versus $`R_{}`$, equation (7) predicts a linear relationship with a slope of -2.

Figure 6 shows the relationship between the accretion filling factor, $`f_{acc}`$, and the stellar radius, $`R_{}`$, for three data samples. The top panel uses the data given in Table 4 from CG98. This panel shows no correlation bewteen the two quantities. This may simply reflect the fact that the other relevant parameters ($`M_{}`$, $`\dot{M}`$, $`P_{rot}`$) must also be considered simultaneously; however, as mentioned above, the CG98 analysis does not strictly fit the observations, and in some case does not match the continuum slopes well. The middle panel of Figure 6 shows the VBJ data from Table 1. In this case an inverse correlation is observed: the correlation coefficient is $`r=0.72`$ with an associated false alarm probability of $`P_f=1.1\times 10^3`$. In a least squares analysis where we again use the distribution around the best fit line to estimate the uncertainty, the best fit slope is $`1.70\pm 0.41`$, which is not significantly different from the predicted value of $`2`$. The best fit is shown in a solid line, and the best fit with an assumed slope of $`2`$ is shown in the dashed line.

One can test this relation using significantly more data points by using the entire sample of VBJ. The required filling factors are taken from Table 5 of VBJ, and the stellar radii are calculated from the stellar luminosities given in Table 5 of VBJ and the spectral types given in Table 1 of VBJ. Again, to translate spectral type into effective temperature, we use the calibration given by Johnson (1966). This entire data set consists of 43 stars which are plotted in the bottom panel of Figure 6. This figure again shows an obvious correlation which is borne out statistically ($`r=0.69`$; $`P_f=9.3\times 10^7`$). The best fit line (solid) has a slope of $`1.65\pm 0.28`$, which again in not significantly different from the predicted value of $`2`$ (shown in the dashed line). We also see that in the entire data sample there is a large amount of scatter, particularly at the radius corresponding to the stars in the CG98 sample. Some of the scatter is likely due to the inherent variability of CTTS and some may also be due to the fact Figure 6 ignores the effects of other stellar parameters on the correlation. Figure 6 also illustrates the need for a large data sample when making comparisons of this sort. The fact that we do see the predicted correlation in the full data set of VBJ with a slope very close to the predicted value suggests that our extension of the Ostriker and Shu (1995) theory may be a realistic way to relate the magnetic field in the X-region to the stellar surface. On the other hand, we note that there is an observational bias that also predicts $`f_{acc}R_{}^2`$, which is shown in equation (11). For a larger star, a smaller filling factor of accretion shock emission is required to give the same area of emitting material, and hence the same luminosity, independent of exactly how the accreting material falls to the stellar surface. While the observed correlation is completely consistent with the predicition of our extended Ostriker and Shu (1995) model, it also suffers from a potential observational bias. Figure 6, based on equation (12), remains the strictest test of our modification to the Ostriker and Shu (1995) model, and it shows the best correlation found in this study.

## 5 Discussion

### 5.1 Results of This Study

The model of magnetically controlled accretion onto a compact object is encountered in many astrophysical contexts. In the case of CTTSs, such magnetospheric accretion appears to explain many of the properties of these stars; however, the general relationships between stellar and accretion parameters predicted by the theory has not been demonstrated observationally. In this paper, we use compilations of rotation periods and studies of accretion onto CTTSs from 5 sources (VBJ, HEC, GHBC, CG98, WG01) to test the general predictions of 4 analytic magnetospheric accretion models (Königl 1991, Cameron & Campbell 1993, Shu et al. 1994, and Ostriker & Shu 1995). All the theoretical studies assume a dipole configuration for the stellar magnetic field; however, Ostriker and Shu’s (1995) detailed analysis of the Shu et al. (1994) theory does not demand it. This treatment identifies the essential quantity of the trapped magnetic flux in the X-region, and this same amount of magnetic flux must thread the stellar surface in the accretion zones, independent of the exact magnetic topology on the stellar surface.

While magnetic field measurements for CTTSs are not numerous, those that are available show a surprising uniformity of the mean stellar field strength from star to star (Guenther et al. 1999, Johns–Krull et al. 1999b, Johns–Krull et al. 2001). Therefore, we have assumed a constant magnetic field from star to star and used data from the literature to look for the predicted correlations among stellar and accretion parameters which follow from the above mentioned theories. We suggest that a physical reason for the constancy of magnetic fields from one star to the other lies in a limit to the efficiency of the dynamo operating in these fully convective stars. Generally, the predicted correlations are absent or very weak in the data. This is also true of the strict Ostriker and Shu (1995) theory as applied to a pure dipole magnetic field geometry. On the other hand, Ostriker and Shu (1995) introduce the concept of the trapped flux, which allows one to abandon a strict dipole field geometry and still use the theory to make testable predictions (equations 6 and 7). In order to test these predictions, the filling factor, $`f_{acc}`$, of the accreting material on the stellar surface is required.

The studies of VBJ and CG98 are the only two considered here which present values for $`f_{acc}`$. For the stars in common in Tables 1 and 4, there is a large range in derived value of $`f_{acc}`$ between the two studies for some stars. However, it is well known that CTTSs are quite variable, showing both rotational modulation of accretion signatures (e.g. Vrba et al. 1993, Johns & Basri 1995, Bouvier et al. 1999) as well as longer timescale variations in their accretion properties. In fact, Ardila and Basri (2000) find variations in the filling factor on BP Tau with a full range from $`f_{acc}=0.000620.026`$, a factor (42) comparable to the extremes found between VBJ and CG98. This suggests that the filling factor on CTTSs can vary substantially, which may account for much of the difference in the filling factors found by VBJ and CG98. On the other hand, for the stars in Tables 1 and 4, the filling factors derived by VBJ are larger by a median factor of 3.5 than those found by CG98, which may reflect a small systematic difference in the two studies. For their default model, Ostriker and Shu (1995) predict a filling factor $`f_{acc}=0.06`$, comparable in value to many of the observed values. More reliable values of the filling factor would provide more stringent tests of the Ostriker and Shu (1995) model. Gómez de Castro and Lamzin (1999) show that analysis of ultraviolet (UV) emission lines in conjunction with shock models may result in tight constraints on accretion parameters. These results are promising, since the UV emission lines arise from the shock itself, as opposed to the photons being reprocessed through the photosphere as in the case of the optical excess studied by all five groups whose results are used in this paper. Unfortunately, the UV analysis has not yet proceeded to the point required to test the current magnetospheric accretion theories. The data in hand (Figures 1b and 4b) show much stronger correlations relative to the pure dipole accretion theories, with the VBJ data fully consistent with the predictions of equation (7). Uncertainties and variations in $`f_{acc}`$ (as well as the other relevant parameters including $`P_{rot}`$) likely contribute to much of the observed scatter in the plots.

A closer look at the VBJ and CG98 data show that the VBJ data does support the Ostriker and Shu (1995) model when as many of the interdependencies of the various parameters as possible are removed. The CG98 data suggests that its apparent support for the extended Ostriker and Shu (1995) model may be the result of the expected correlation between the mass accretion rate and the area of the star participating in the accretion flow. On the other hand, none of the other magnetospheric accretion theories are well supported by the existing data, making the our extension of the Ostriker and Shu (1995) description the most viable of the current analytic theories. We note though that this study by no means rules out the work of Königl (1991), Cameron and Campbell (1993), or Shu et al. (1994) which employ strictly dipole field geometries. It may be that these models are essentially correct, and it is the dipole component of the field that varies from star to star in the manner required by the theory. We do note, however, that for the case of BP Tau (Johns–Krull et al. 1999a) and TW Hya (Johns–Krull & Valenti 2001) that the measured dipole components on these stars are substantially below the theoretically predicted values.

### 5.2 Implications

The apparent success of the Ostriker and Shu (1995) trapped flux model once a dipole geometry is abandoned, and the lack of success for the other theories which explicitly assume a dipole field geometry, suggests the magnetic fields on CTTSs are not large scale dipoles. Johns–Krull et al. (1999b) tabulate dipole field predictions for Königl (1991), Cameron and Campbell (1993), and Shu et al. (1994). Measured field strengths (Basri, Mary, & Valenti 1992; Guenther et al. 1999, Johns–Krull et al. 1999b, Johns–Krull et al. 2001) on TTSs are in general agreement with these predictions; however, the dipole component of these fields appears to be smaller than the predictions by a factor of at least 10 in the case of Königl (1991) and Shu et al. (1994) based on the lack of detection of circular polarization in magnetically sensitive lines of CTTSs (Johns–Krull 1999a, Johns–Krull & Valenti 2001). The predicted dipole fields of Cameron and Campbell (1993) are weaker, but are still ruled out by the observations at the $`6\sigma `$ level, and the correlation analysis here does not generally support the Cameron and Campbell (1993) formulation. Independent of the theoretical preditions, the measured strength of the mean photospheric fields and the lack of observed circular polarization in the magnetically sensitive photospheric lines (Johns–Krull et al. 1999b, Johns–Krull & Valenti 2001) indicates that the surface field topology on at least two CTTSs (BP Tau and TW Hya) is not dominted by a dipole. On the other hand, the dipole component of the stellar field will fall off the slowest with distance from the stellar surface and is likely to dominate at the disk truncation radius. We again note that this generally happens in the case of the solar magnetic field at $`2.5R_{}`$ (Luhmann et al. 1998). This suggests that the material accreting onto the star generally does follow dipole-like magnetic field lines which explains the general success of line profile calculations done for dipole geometries (Hartmann, Hewett, & Calvet 1994; Muzerolle et al. 1998, 2001) as well as the circular polarization observed in the He I 5876 Å emission line of several CTTSs (Johns–Krull et al. 1999a, Johns–Krull & Valenti 2001).

It is interesting to briefly consider the study of Ferreira et al. (2000) which explores CTTS interactions with the surrounding disk for both dipole and higher order magnetic field components. This study is primarily concerned with the rotational spin down of young protostars into and through the TTS phase; however, the authors explicity ignore the torque associated with the surrounding accretion disk. As a result, the results may not be accurate. In all their numerical results, Ferreira et al. (2000) find a (sometimes steep) increase in the stellar field strength as the star ages from $`10^6`$ to $`10^7`$ years and beyond. The limited magnetic field observations available do not currently support this picture: Johns–Krull and Valenti (2001) find a mean magnetic field on TW Hya (age $`10^7`$ yr - Webb et al. 1999) which is well within the range found for CTTSs is Taurus. Additionally, Valenti and Johns–Krull (2001) report a mean field on a K3V Pleiades (age $`10^8`$ yr - Basri, Marcy, & Graham 1996) member which is less than those of TTS by about a factor of 2. More data is needed, particularly on young stars with an age of $`10^7`$ yr to see if any trends emerge.

We have just argued that the magnetic fields of CTTSs are not large scale dipoles and that the dipole component of CTTSs is significantly smaller than required by the magnetospheric accretion models of Königl (1991), Cameron & Campbell (1993), and Shu et al. (1994). Ostriker and Shu (1995) detail the theory of Shu et al. (1994) and show that the dipole assumtion is not necessary by introducing the idea of trapped flux. How does the field strength in the trapped flux model compare with the observed field strength? Solving equation (6) for the stellar magnetic field gives:

$$B_{}=14,500\left(\frac{M_{}}{M_{}}\right)^{1/2}\left(\frac{\dot{M}_D}{10^7M_{}\mathrm{yr}^1}\right)^{1/2}\left(\frac{P_{rot}}{1\mathrm{day}}\right)^{1/2}\left(\frac{R_{}}{R_{}}\right)^2\left(\frac{f_{acc}}{0.01}\right)^1\mathrm{G}$$
(17)
where we have set $`A_{}=4\pi R_{}^2f_{acc}`$. Since only the studies of VBJ and CG98 independently estimate the filling factor, $`f_{acc}`$, required in this equation, we only use these studies to estimate the stellar magnetic field strength, $`B_{}`$. Table 8 gives these estimates for the 3 stars observed by Johns–Krull et al. (2001) for which the calculation can be made. Also given in the table is the observed mean field strength from Johns–Krull et al. (2001). It is important to note that in the analysis of Johns–Krull et al. (1999b, 2001) evidence is found for a range of field strengths on the surface of CTTSs, and we are unable to specify the field strength only in the accretion zone. The He I polarization data of Johns–Krull and Valenti (2000) and Valenti et al. (2001) probably probes the field in the accretion zone at the stellar surface, but there is an unknown inclination of this field to the line of sight which reduces the observed field strength below the actual field value. With these caveats, the predicted and observed field strengths agree to within a factor of 2-3, with the predictions using the parameters of VBJ generally falling below those based on the parameters of CG98 in this sample of 3 stars. Except for the rotation period, these two studies find different values for all the relevant parameters in equation (11); however, the one which is most significant in producing the difference in the predicted field strengths is the filling factor, $`f_{acc}`$. CTTSs are know to vary due to rotational modulation. In addition, there are also likely to be real changes in the instantaneous accretion rate onto the star (e.g. Ardila & Basri 2000). Both processes could change the apparent filling factor of accreting material. As a result (and as mentioned above), it is not clear that the difference in the derived filling factor values between VBJ and CG98 is significant or simply the result of the variability of these stars.

Finally, we note that the relationships set forth in this paper are based on equilibrium calculations, while CTTSs do show substantial variation. As a result, we expect that the relationships should hold for the mean value of the relevant parameters (particularly the accretion rate), but probably do not describe the relationship between these parameters as they vary over relatively short timescales (compared to the evolutionary timescale of CTTSs and their disks). For example, Ardila and Basri (2000) study the accretion variability of the CTTSs BP Tau, finding that $`f_{acc}\dot{M}^2`$. At first, this appears to contradict equation (7). Rotational modulation may produce some of the variation observed by Ardila and Basri (2000). Aside from this, over the timescale of the variations in the instantaneous accretion rate studied by Ardila and Basri (2000), it is doubtful that the star reached a new equilibrium rotation rate which it would be expected to do. In other words, the torque balance on the star may also be varying, alternately trying to speed it up or slow it down, all the time varying around the equilibrium value. We would only expect equation (7) to hold for the mean values, as it seems to based on the correlation plots presented in this paper. Indeed, some of the scatter seen in these plots is likely to result from the intrinsic variability displayed by CTTSs.

## 6 Conclusions

The magnetospheric accretion model can explain many observed features of CTTSs. Equilibrium models of this process have relied primarily on dipole field geometries to explore the interaction of the star with the disk, which seems appropriate since the dipole component should dominate at large distance from the star. Observations suggest that the magnetic fields of CTTSs do not change much from star to star, though further observations are needed to confirm this. Under the assumption that the fields do not vary strongly, equilibrium theories predict specific correlations between stellar and accretion parameters which are generally not observed. This may mean that equilibrium theories are simply not adequate, or that a dipole field topology is too simple an assumption for relating observable parameters. On the other hand, we can use the notion of trapped flux introduced by Ostriker and Shu (1995) to extend this equilibrium theory to non-dipole field topologies at the stellar surface. The extension we make relies on assuming that the dipole component of the magnetic field does dominate at the X-region in the Ostriker and Shu (1995) model and that their analysis is accurate in terms of the accounting of how much magnetic flux is in the accretion flow, the dead zone, and the wind. This determines the magnetic flux in the accretion flow which we assume maps back to an equal amount of magnetic flux on the stellar surface. Again holding the magnetic field constant, new correlations are predicted bewteen stellar and accretion parameters. We argue that the data of VBJ are the best suited for looking for these correlations, and that they are apparent in the data. This suggests that equilibrium models are adequate to describe the general features of magnetospheric accretion in CTTSs once the dipole field assumption is dropped. Additional tests of our modified Ostriker and Shu (1995) model will confirm or refute this. Particularly valuable tests would be to determine all the relevant stellar and accretion parameters for stars with substantially different stellar properties than the Taurus stars used here. For example, detailed observations (including magnetic field measurements) of the more rapidly rotating CTTSs in Orion could provide critical tests of equation (6) and (7). Tests of these equations will also benefit from better estimates of accretion rates and filling factors which may result from the analysis of UV emission lines produced in the accretion shock (e.g. Gómez de Castro & Lamzin 1999). Until then, it appears that our modified Ostriker and Shu (1995) is a reasonable description of magnetospheric accretion onto CTTSs.

CMJ-K would like to acknowledge stimulating discussions with F. Shu and to thank him for pointing out the role of the trapped flux. We acknowledge the careful reading of the manuscript by S. Mohanty and the useful discussions which resulted. Communications with L. Prato and R. White also provided important input on several issues raised in the manuscript. We also wish to thank an anonymous referee for feedback provided on this paper. CMJ-K and ADG both acknowledge partial support from a NASA Origins grant NAG5-8098 to the Regents of the University of California.


format: text

text:
# Bistable Gradient Networks II: Storage Capacity and Behaviour Near Saturation

## I  Introduction

In this paper we consider the behaviour at high memory loading of a Hopfield-like attractor network of $`N`$ bistable elements, the bistable gradient network or BGNChinarov IWANN . This is a sequel to an earlier paper LowLoading  which considered the BGN at low loading. We compare the BGN to the deterministic Hopfield network (HN)Hopfield -Amit2 , examining the storage capacity and other key properties.

To begin, we review the BGN model and establish some notation. The BGN is described by the coupled differential equations

$$\frac{dx_i}{dt}=\frac{H}{x_i},$$
(1)
where $`x_i`$ are $`N`$ continuous-valued real variables (or components of an $`N`$-dimensional state vector $`𝐱`$) representing the outputs of the $`N`$ nodes of the network, and $`H`$ is the Hamiltonian

$$H=H_0+H_{int}=\underset{i=1}{\overset{N}{}}\left(\frac{x_i^4}{4}\frac{x_i^2}{2}\right)\frac{1}{2}\gamma \underset{i,j=1}{\overset{N}{}}w_{ij}x_ix_j.$$
(2)
The quantities $`w_{ij}`$ are the elements of a symmetric matrix of connection strengths and $`\gamma `$ is a control parameter determining the overall strength of the coupling among nodes. As in the Hopfield modelHopfield Little , the connections or synaptic weights are determined by the Hebb learning ruleHebb

$$w_{ij}=\frac{1}{N}\underset{\mu =1}{\overset{p}{}}\xi _i^\mu \xi _j^\mu \frac{p}{N}\delta _{ij}$$
(3)
where the N-dimensional vectors $`\xi ^\mu `$ represent a set of $`p`$ distinct *memory patterns* to be recognized by the network. We take these patterns to consist of binary elements $`\pm 1`$ only, and we assign them random values, thus introducing quenched disorder. The BGN’s key difference from the HN and from most of its continuous-valued relatives ContinuousHop  lies in the presence of the local quartic potential term $`H_0`$ in the Hamiltonian, which renders each node bistable.

The interaction term $`H_{int}`$ furnishes an input to each node given by $`h_i=\gamma _{\mu =1}^pw_{ij}x_j`$ so that the dynamical equation for each node is given by

$$\frac{dx_i}{dt}=x_ix_i^3+h_i.$$
(4)
The input may also be referred to as a “magnetic field” by analogy with Ising spin systems. In the absence of input each node has two stable fixed points at $`\pm 1`$, but a nonzero field shifts the two fixed points. If the critical magnitude $`h_c=\frac{2\sqrt{3}}{9}0.38`$ is exceeded then the field is strong enough to overcome the potential barrier in the quartic potential and there is then only one fixed point.

We define two sets of order parameters $`m_\mu `$ and $`b_\mu `$ by

$$m_\mu \xi ^\mu 𝐱=\frac{1}{N}\underset{i=0}{\overset{N}{}}\xi _i^\mu x_i$$
(5)
and

$$b_\mu \frac{1}{N}\underset{i=0}{\overset{N}{}}\xi _i^\mu \mathrm{s}ign(x_i).$$
(6)
$`m_\mu `$ are inner products or overlaps of the state vector with the memorized patterns, while $`b_\mu `$, the “bit overlaps,” encode information about sign agreements between the state vector and the stored patterns. For the purposes of this paper, we will for the most part be more interested in $`b_\mu `$ than in $`m_\mu `$, so where there is little risk of confusion we will drop the word “bit” and simply refer to $`b_\mu `$ as an “overlap.”

The degree of loading of the network’s memory can be parametrized by the ratio $`\alpha p/N`$. In the companion paper LowLoading  we examined the behaviour of the BGN in the low-loading limit $`\alpha 1`$. It was shown that the network can function as an associative memory and correct sign flip errors in a stored pattern as long as $`\gamma >\frac{1}{3}`$ . We found that the attractors of the BGN’s dynamics are readily classified into three categories which are separated from each other in energy. The lowest energy states are the *memory* or *retrieval states*, each of which corresponds to one of the memorized patterns. In addition to these there are higher-energy spurious attractors of two types. The *mixture* or *spin glass states* have partial overlaps with several patterns and thus lie close to the subspace spanned by the memory patterns. *Uncondensed states*, which have no counterpart in the HN, are states in which none of the fields acting on the nodes are strong enough to cause sign flips and the dynamics is therefore dominated by the local potential. They have energies per node close to $`0.25`$. The spin glass states are intermediate in energy between the memory states and the uncondensed states. In the range $`\frac{1}{3}\gamma 1`$, pattern recognition by the network was found to be highly selective; the input must be close to the stored pattern in order for the pattern to be fully restored. The uncondensed states are numerous (of order $`2^N`$) and fill most of the configuration space. The memorized patterns and their basins of attraction occupy isolated valleys among these states, but these valleys expand as $`\gamma `$ increases. For $`\gamma 1`$, on the other hand, the behaviour resembles that of the HN: there are no uncondensed states, and the memory states have large basins of attraction.

In this paper we turn to the case where $`p/N`$ is of order unity. We are interested in the maximum storage capacity, or the maximum number of patterns that can be successfully stored and retrieved, as well as changes in the network’s performance as this limit is approached. Earlier work with small networks IWANN  suggested that at least under certain conditions the BGN could store many more patterns than the HN while possessing few spurious attractors. It is known that the Hopfield storage limit of $`p0.14N`$ memory patterns can be exceeded if a more complicated learning algorithm is used Amit2 , but in the BGN case improved capacity is achieved with the familiar Hebb rule. Since the previous results Chinarov IWANN  were obtained with networks much too small to be of practical interest (e.g., $`N=5`$), we now examine larger networks, mainly through numerical simulations. (At the end of the paper we will return briefly to the small-network case.) We find that the high-loading behaviour, like that at low loading, depends strongly on $`\gamma .`$ For $`\gamma 1`$, a Hopfield-like first-order phase transition results in the destabilization of all memory states at a critical value of $`\alpha `$. For $`\gamma =2`$ this transition occurs at $`\alpha _c0.1`$, compared to $`\alpha _c0.14`$ for the HN. For $`\gamma =0.5`$, on the other hand, we find that it is possible for the stored patterns to remain stable even at loading factors of $`\alpha 0.3`$ and higher. Furthermore, there is no sudden blackout; instead, the performance degrades gradually as $`\alpha `$ increases. The price of this high capacity is that the retrieval of the patterns from corrupted versions may be imperfect.

The remainder of the paper is organized as follows: In section II we discuss in general terms the effects of crosstalk, or interference between different stored patterns. It is crosstalk which is responsible for limitations on storage capacity. We compare the effects of crosstalk in the BGN and the HN. This discussion provides a framework for interpreting our numerical results. In section III we examine numerically the stability of memory patterns as a function of $`\alpha `$. We find evidence of a first-order memory blackout phase transition in the BGN at high $`\gamma `$, but not at low $`\gamma `$. In section IV, we examine the effects of high loading on the basins of attraction for the memory states and on their retrieval from corrupted input. We see that increasing the loading markedly alters the energy landscape. In section V we comment briefly on the behaviour of smaller networks and on the relation between the previous small network results Chinarov IWANN  and those we have obtained for larger networks. We conclude with some discussion of the BGN results and possible directions for further study. We make some conjectures concerning the performance of the BGN in the presence of stochastic noise.

## II Pattern retrieval and crosstalk

In a given state $`𝐱`$, the input to the $`i`$-th node of the network can be expressed as

$$h_i=\gamma \underset{j=1}{\overset{N}{}}w_{ij}x_j=\gamma \underset{\mu =1}{\overset{p}{}}\xi _i^\mu m_\mu \gamma \frac{p}{N}x_i,$$
(7)
using the Hebb rule (3) and the definition (5) of $`m_\mu `$. In LowLoading , we showed that when $`\alpha 1`$ there are stable retrieval states in which $`b_\nu =1`$ for one particular $`\nu `$. For a state in which the $`\nu `$-th overlap is large, we can decompose the field as follows:

$`h_i`$ $`=\gamma \left(\xi _i^\nu m_\nu +{\displaystyle \underset{\mu \nu }{}}\xi _i^\mu m^\mu {\displaystyle \frac{p}{N}}x_i\right)`$ (8)
$`=\gamma \left(S_i+C_i{\displaystyle \frac{p}{N}}x_i\right).`$
We refer to the first term, $`S_i`$, as the “signal” term and the second term, $`C_i`$, as the “crosstalk” term. In the limit $`\frac{p}{N}0`$, the mutual overlaps between different patterns is small, so that $`m_\nu 1(\nu \mu )`$, and the last two terms in (8) can be neglected. The signal term is then dominant and there is a stable retrieval state given by $`𝐱=\sqrt{1+\gamma }\xi ^\nu `$ with $`m_\nu =\sqrt{1+\gamma }`$. We expect this solution to be approximately valid for small but nonzero values of $`\frac{p}{N}`$. For this case, the overlaps $`m_\mu (\mu \nu )`$ behave as Gaussian distributed random variables with zero mean and variance $`1/\sqrt{N}`$. Accordingly, the crosstalk term $`C_i\underset{\mu \nu }{}\xi _i^\mu m^\mu `$ in (8), being a sum of $`p`$ such quantities, is a random quantity with zero mean and variance $`\sqrt{p/N}`$. The third term, which arises from the subtraction of the diagonal elements, is of order $`p/N`$ and thus is generally smaller than the crosstalk term. We will neglect it for the moment.

In the absence of crosstalk, a retrieval state is not only linearly stable (i.e., stable against small perturbations) but, for $`\gamma >\frac{1}{3}`$, it is also stable against individual sign flips. The latter means that if a retrieval state is corrupted by changing the sign of one or a small number $`N`$ of nodes, then the dynamics will reverse the flipped sign and restore the pattern. This happens because, in the absence of crosstalk, each node experiences a field $`\gamma S_i=\gamma \xi _i^\nu m_\nu `$ which has the same sign as $`\xi _i^\nu `$ and for $`\gamma >\frac{1}{3}`$ that field is strong enough to overcome the potential barrier of the individual node. Now consider a given node (say, the $`i`$-th node) in the presence of crosstalk. The crosstalk field acting on the $`i`$-th node may be either aligned with or opposed to $`x_i`$. If it is aligned, then its effect on that node is to increase the magnitude of $`x_i`$, making it larger than $`\sqrt{1+\gamma }`$. If it is opposed to $`x_i`$, then its effect is to *decrease* the equilibrium magnitude of $`x_i`$. If the crosstalk term is large enough, then it may be sufficient to overcome the local potential barrier and reverse the sign of $`x_i`$, thus introducing a sign error into the pattern. This will occur only if

$$C_i+S_i<\frac{2\sqrt{3}}{9\gamma }.$$

By contrast, in the Hopfield model a sign error is introduced if

$$C_i+S_i<0.$$

Thus the relative strength of crosstalk required to introduce sign errors is greater for the BGN than for the HN, and the discrepancy is greatest for small values of $`\gamma `$. One might expect that this would make the BGN less vulnerable to crosstalk (and the memory states more stable) than the HN, especially at low $`\gamma `$, but this is not a foregone conclusion as the BGN’s dynamics include mechanisms which tend to amplify small initial overlapsLowLoading  and could conceivably also amplify crosstalk. Our numerical results confirm that the BGN is in fact less prone to crosstalk-induced errors at low values of $`\gamma `$, but not at high values.

If

$$\frac{2\sqrt{3}}{9\gamma }<C_i+S_i<\frac{2\sqrt{3}}{9\gamma }$$

the crosstalk will not be strong enough to reverse the sign of $`x_i`$ if initially $`x_i`$ is correctly aligned with $`\xi _i^\nu ,`$ but it will nonetheless destroy the stability against a sign flip. In other words, if $`x_i`$ is initially misaligned, it will not be corrected. We may say that the node is *“bistabilized”* but not destabilized. (Such an effect cannot occur in the HN where the nodes are not individually bistable.) Since the crosstalk is random, it will in general bistabilize some nodes and not others, with the result that the memory state will be stable against sign flips of certain nodes but not of others. Thus, even though a memory state may be linearly stable, the dynamics may only be able to correct some sign errors, not all. This contrasts with the low-loading case where crosstalk is negligible and there is a single threshold coupling strength, $`\gamma =\frac{1}{3},`$ above which any single sign error can be corrected. In general, crosstalk results in non-uniform behaviour among nodes, including different magnitudes of $`x_i`$ for different $`i`$.

## III Stability of the memory states and remanent overlap

To examine the stability of the memory states, we followed a procedure similar to that of reference Amit2 . Using an ensemble of realizations of the random patterns, we made a number of trials in which the network was initialized to the state $`𝐱=\xi ^\nu `$ for some pattern $`\nu `$. The initial bit overlap $`b_\nu `$ was thus equal to 1. We then allowed the state to evolve until an attractor was reached. In each case we measured the final energy as well as the bit overlap between the initial and final states, to which we refer as the *remanent* bit overlap $`b_r`$. We then constructed histograms showing the probabilities $`P(b_r)`$ for $`b_r`$ falling in intervals of width 0.05. Each of the histograms discussed in this section represents an ensemble of at least 200 initial conditions. If the memory state is entirely stable, as is the case at very low loading, then after the dynamics converges the final bit overlap will still be equal to $`1`$. For higher loading, however, the crosstalk fields may introduce sign errors so that $`b_r<1`$.

### III.1 HN and BGN with high $`\gamma `$

First consider the case of the HN. Figure 1 shows a series of $`P(b_r)`$ histograms for the HN at different values of $`\alpha =p/N`$. (The data here are our own, but the results are comparable to those given in Amit2 . We include them for comparison with BGN results.) At lower levels of loading, crosstalk introduces few errors, so $`b_r>0.95`$ in the large majority of cases. However, as $`\alpha `$ increases beyond the critical value $`\alpha _c0.148`$, the high-$`b_r`$ peak of the distribution begins to vanish and a second, broader peak begins to grow in the vicinity of $`b_r0.3`$. The states in the second peak are spin glass states. As shown in figure 2, this transition becomes sharper with increasing network size, and finite-size scaling analysis shows behaviour characteristic of a first order phase transition in the thermodynamic limit Amit2 . In the limit $`N\mathrm{}`$ the associative memory fails suddenly as the critical loading is exceeded— the remanent overlap drops abruptly from near 1 to 0.3. This nonzero value of the remanent overlap above the critical loading was noted in Amit2  and attributed to replica symmetry breaking, as the replica symmetric theory predicts that $`b_r`$ should drop to zero above the phase transition. This phenomenon is related to the non-zero remanent magnetization of a spin glass Mezard .

In the BGN with $`\gamma =2`$, a similar transition evidently occurs at $`\alpha _c0.1`$. As evidence, in figure 3 we show two series of histograms at increasing $`N`$, one below the suspected transition and one above. As in the HN case, the transition grows sharper with increasing network size. Below the critical loading, the high-$`b`$ peak remains robust as $`N`$ increases, while above the critical loading the high-$`b`$ peak shrinks with increasing network size and the low-$`b`$ peak grows. Two quantitative differences are that the critical loading $`\alpha _c`$ is lower in the BGN case, $`\alpha _c0.1`$, while the average remanent overlap above the critical loading is higher, near 0.45.

The first-order nature of the transition is confirmed by examining the energies of the final states. These energies and the overlaps are shown in a scatter plot in figure 4A. The spin glass states associated with the low-$`b_r`$ peak are clustered at energies below those of the retrieval states. The gap in energy between these two clusters corresponds to the latent heat of the phase transition.

### III.2 BGN at low $`\gamma `$

For $`\gamma =0.5`$, in contrast to $`\gamma =2`$, the BGN’s behaviour differs markedly from that of the HN. A series of $`P(b_r)`$ histograms for different values of $`p/N`$ is shown in figure 5. Two features are evident. First, the stored patterns remain stable with few errors even up to high levels of storage, $`\alpha =0.3`$. Second, there is no evidence of a discontinuous memory failure; rather, the retrieval quality as measured by $`b_r`$ appears to degrade continuously as $`\alpha `$ increases. No second peak appears in the histograms; instead, the high-$`b`$ peak first spreads and then begins to drift downward as errors accumulate.

At an intermediate value $`\gamma =1`$ the $`P(b_r)`$ histograms (figure 6) suggest a first-order transition, although the evidence is less pronounced than in the $`\gamma =2`$ case. A second peak appears above the transition, and the high-$`b_r`$ peak shows clear signs of shrinking as $`N`$ increases, but the tails of the two peaks overlap substantially . The greater overlap between the peaks comes about for two reasons. First, the high-$`b_r`$ peak above the critical loading is broader than in the $`\gamma =2`$ case. Second, the remanent magnetisation is much higher (i.e., the drop in $`b_r`$ at the critical point is much smaller.) The latent heat is also much smaller, as can be seen from the scatter plot of the energy (figure 4B). The critical loading, or storage capacity, is approximately 0.17, higher than for $`\gamma =2`$ and higher than for the HN.

## IV Attractors, basins and the energy landscape at high loading

In the previous section, we examined the trajectories of initial conditions corresponding to memorized patterns and determined whether these trajectories remain close to the pattern or move away from it. Such experiments, however, probe only one aspect of network performance. We are interested not only in the stability of memory states but also in the sizes of their basins of attraction. The function of associative memory depends on the ability of the dynamics to correct partially corrupted patterns. More generally, we are interested in the evolution of the energy landscape with increasing $`\alpha `$. To address these issues, we performed two additional sets of numerical experiments. First, we examined the attractors reached from a large number of *random* initial conditions to obtain a uniform sampling of the phase space and a broad picture of the energy landscape. Second, we examined the fate of initial conditions at specified Hamming distances from memory patterns. The latter set of experiments probes the landscape in the vicinity of the memory states. Results of similar experiments were given in ref. LowLoading  for the case of low loading.

### IV.1 Evolution of the energy landscape: attractors reached from random initial states

In the case of the HN, it is known that in the thermodynamic limit the memory states are the lowest energy states for $`\alpha <0.05`$, while for higher values of $`\alpha `$ spin glass states arise which have lower energies. Up to $`\alpha _c0.148`$ the memory states remain local minima of the energy even though they are not the global minima. Above $`\alpha _c`$ they cease even to be local minima and therefore become unstableAmit2 .<sup>1</sup><sup>1</sup>1For a schematic illustration of the evolution of the energy landscape, see fig. 2.18 of reference HKP .  The drop in energy from the memory states to the spin glass states at $`\alpha _c`$ is the latent heat. One way to observe the evolution of the energy landscape is to examine the attractors reached from an ensemble of random initial conditions which effectively samples the configuration space. Figures 7-9 show histograms for the energies of attractors sampled in this manner. In each case, we sampled a total of at least 200 random initial conditions with several realizations of the random patterns $`\xi ^\mu `$. In figure 7, for the HN, we can see that at low loading the attractors are separated into two groups, the retrieval states at $`E/N=0.5`$ and spurious states at a range of higher energies. The probability of retrieving a memory state from a random initial condition is high. With increasing $`\alpha `$, the spurious states move to lower in energies until they are below the memory states. At the same time, their basins of attraction take up a larger portion of the configuration space, as is apparent from the growing size of the spin glass peak in the histogram and the shrinking size of the retrieval state peak. Figure 8 shows that for the BGN with $`\gamma =2`$ the evolution is qualitatively similar. In figure 9, for $`\gamma =0.75,`$ we observe that an additional effect of high loading is to destabilize the uncondensed states. At low loading, the uncondensed states dominate the configuration space—almost all random initial conditions land on an uncondensed state, as was noted in LowLoading . The other two types of state begin to show their presence as the loading is increased, while the uncondensed states eventually disappear.

### IV.2 Retrieval of patterns from corrupted versions

To obtain information about the landscape in the vicinity of the memory states and about the shapes of their basins of attraction, we examined the fates of initial conditions which were not random but rather at specified initial overlaps with particular stored patterns. As in LowLoading , these initial configurations were generated by starting with a particular “target” pattern and flipping the signs of a specified number of randomly chosen nodes. For each value of the initial overlap $`b_{init}`$, we generated an ensemble of initial conditions for several different realizations of the random set of memory patterns. We then evolved these states until the dynamics converged, and evaluated $`b_{final}`$, the final overlap with the target pattern, for each trajectory. If $`b_{final}=1`$, this signifies that all signs which were initially flipped have been corrected and the target pattern has been retrieved perfectly. If $`b_{init}<b_{final}<1`$, then the pattern has been retrieved imperfectly. The final state is closer to the stored pattern, but not all errors have been corrected. If $`b_{init}>b_{final}`$, then the trajectory has moved farther away from the stored pattern. As discussed in section II, the ability of the network to correct sign errors depends on the competition among the signal, the crosstalk, and the local potential. The initial states currently discussed have at least a moderately large overlap with the target pattern, resulting in a signal, and random overlaps with the other memory patterns, resulting in crosstalk.

Consider now the HN. At low loading, there is very little crosstalk. The crosstalk term is typically of order $`\sqrt{p/N}`$. If the HN is set in an initial condition whose overlap with the target pattern $`\xi `$ is $`b_{init}>\sqrt{p/N}`$, then most nodes experience a net field which is aligned with the target pattern. Nodes which are initially misaligned with the target state ($`x_i=\xi _i`$) will tend to change their signs and align with $`\xi _i`$. Each node which realigns in this way increases the value of $`b`$ and thus increases the strength of the signal acting on the remaining nodes. As a result, even if in the initial state some fraction of the nodes experience a net field opposed to $`\xi _i`$, eventually the growing signal may overcome the crosstalk and correct those nodes as well. Therefore for low loading, as long as the initial state has $`b_{init}>\sqrt{p/N}`$, the probably of completely retrieving the target pattern is close to unity. As the loading ratio $`\alpha `$ increases, however, the typical crosstalk becomes stronger, and a higher signal is required to overcome the crosstalk noise. Therefore sign errors are not likely to be corrected unless the initial overlap is above a threshold, which grows higher with increasing $`\alpha `$. If the crosstalk is too large, then some signs which are initially aligned with the pattern may be flipped, and the state may move away from the target pattern instead of toward it. Each node which flips out of alignment with the target pattern reduces the size of the overlap and hence of the signal, which makes other nodes more susceptible to crosstalk-induced errors, and a cascade of errors can occur. The critical loading $`\alpha _c`$ is the level at which even a state with $`b_{init}=1`$ becomes unstable against such a cascade of errors.

Figure 10A shows a scatter plot of $`b_{final}`$ vs. $`b_{init}`$ for a HN slightly below critical loading. There is a threshold overlap for retrieval. If $`b_{init}>0.5`$ the signal is strong enough to correct errors and the majority of trajectories flow towards the target pattern. For $`b_{init}<0.5`$ the majority of trajectories move instead away from the target pattern. As $`\alpha `$ increases, the threshold value of $`b_{init}`$ for retrieval increases until at $`\alpha _c`$ it reaches 1. The plot in figure 10A is for a network with $`N=2000`$; experiments with networks of different sizes reveal that the retrieval threshold becomes sharper as $`N`$ increases.

In the BGN, on the other hand, the dynamics of error correction is more complicated due to the local potential. At strong coupling $`\gamma 1`$, the potential barriers against sign flips are less important than at weak coupling. As a result the BGN in this regime behaves in many respects like the HN. It is not surprising, then, that the scatter plot of $`b_{final}`$ vs. $`b_{init}`$ for a BGN with $`\gamma =2`$ slightly below its critical loading (figure 10B) appears qualitatively similar to the corresponding figure 10A for the HN. There is a threshold (approximately $`b_{init}=0.6`$) below which the probability of fully retrieving the target pattern drops sharply. Above this threshold, the signal is evidently strong enough to correct most sign errors. A key difference, however, is that even below this threshold the average $`b_{final}`$ is larger than $`b_{init}.`$ This means that the majority of a trajectories move toward the target pattern rather than away from it, but only some of the sign errors are corrected, not all.

In the case of the BGN with $`\gamma =0.5`$, the local potential barriers are important, and interesting dynamics results from the competition among the signal from the target pattern, the crosstalk from other patterns, and the local potential. In the $`\alpha 1`$ case LowLoading , sign errors can be corrected only if the signal is strong enough to overcome the local potential barriers, and therefore there is a threshold value of the initial overlap which is much larger than $`\sqrt{p/N}`$. For the case $`N=1000`$, $`p=5`$, for example, we find that this threshold is approximately $`b_{init}=0.5`$. For $`b_{init}>0.5`$ there is a large probability that the target pattern is retrieved perfectly, while for $`b_{init}<0.5`$ there is a large probability that the network will be stuck in an uncondensed state with $`b_{final}=b_{init}`$. In the intermediate range $`0.4<b_{init}<0.6`$, there is also a significant probability that the trajectory is attracted to an asymmetric spurious state in which $`b_{final}`$ is large but not unity and there are larger than random overlaps with one or more other memory patterns. This behaviour is illustrated by the scatter plot of figure 11A. As the loading increases (Figure 11A-D), something surprising happens: at first, the basins of attraction of the memory states *expand* slightly, contrary to what one would expect from the HN. The frontier of the uncondensed states is pushed back to lower values of $`b_{init}`$. For $`p=50`$, or $`\alpha =0.05`$ (fig. 11D) almost all states with $`b_{init}>0.1`$ undergo some motion toward the target pattern, even if the pattern is not retrieved perfectly. As $`p`$ increases further, the $`b_{final}`$ vs. $`b_{init}`$ plot preserves the approximate shape of fig. 11D.

The dynamics of retrieval and error correction for the BGN with weak coupling and high loading is evidently quite different from that of the HN. In the $`\alpha =0.05`$ case of figure 11D, $`b_{final}`$ is always $`1`$ if $`b_{int}=1`$, which means that the memory state is stable. However, it is retrieved only imperfectly when sign errors are introduced: if $`b_{init}<1`$ then $`b_{init}<b_{final}<1`$. This indicates that some nodes remain bistable and cannot be corrected. If the initial state is close to the target pattern, then the majority of errors are corrected, but that fraction decreases with decreasing $`b_{init}`$. Such a partial correction of errors does not often occur in the HN case. In the latter case, an initial condition either flows all the way to the retrieval state or moves away from it toward another attractor. The retrieval state may itself have a small number of errors due to crosstalk, but the presence of these errors does not depend on the initial state. The energy landscape of the HN in the neighborhood of a memory state apparently has the shape of a smooth basin — once the basin is entered, the trajectory usually runs without obstruction to the attractor at the bottom. For the weakly coupled BGN, on the other hand, the landscape appears to have the structure of a “funnel,” funnel1 funnel2 Karplus  i.e., a sequence of local minima at decreasing energies, with low potential barriers separating each state from the next. There is a region of configuration space which has an overall tilt toward the retrieval state, but in which there are many local minima which may trap the trajectory before it reaches the retrieval state. This is illustrated schematically in figure 12. Funnel-shaped energy landscapes were first suggested in the context of protein folding dynamics.

## V A Comment on Small Networks

So far, this paper and the companion paper LowLoading  have focused mainly on large networks of $`N=1000`$ or more. However, some applications of neural network algorithms to robotics and other areas make use of networks of only 20 or 100 nodes. Experimental studies of BGN-like chemical reactor networks ChemNetE ChemNetM  used fewer than 10 nodes. The current work on the BGN was inspired in part by results suggesting that the BGN could store many more patterns than the HN, with fewer spurious attractors IWANN . The latter results were inferred from a few selected cases using very small networks, and so we attempted to test the genericness of these results for a variety of small networks as well as for larger networks.

It should be noted that with very small networks, there are great fluctuations in properties depending on the particular set of memory patterns chosen, as it is impossible to ignore the mutual correlations among patterns. Results for the maximum storage capacity of a network are well-defined only in the thermodynamic limit $`N\mathrm{}`$, and for small $`N`$ even an HN may in particular cases be able to store more than $`0.14N`$ stable patterns. For these reasons one cannot draw strong general conclusions about storage capacity based on small networks alone, but it is nonetheless instructive to make some comparative studies of pattern stability in small networks. For several small values of $`N`$ and $`p`$, we generated random sets of stored patterns and tested their stability, using the HN and the BGN at several different values of $`\gamma `$. We counted the average fraction of memory patterns which were stable, $`f_{stable}`$.<sup>2</sup><sup>2</sup>2 In the thermodynamic limit there is an approximate permutation symmetry among the memory patterns, so that in general either all will be stable or none will be. By contrast, in small networks it is not unusual for some memorized patterns to reamain stable while others are not. In addition, we estimated the fraction of configuration space $`V_{spur}`$ occupied by spurious attractors by following the trajectories of random initial conditions inside the hypercube $`\left|x_i\right|<2.`$ These results are collected in table 1. In all cases the results were averaged over at least 100 sets of randomly generated patterns. In generating random sets of patterns, we did not eliminate cases in which two or more patterns are identical. The results show that, contrary to the selected examples discussed in Chinarov  and IWANN , the spurious attractors do *not* generically occupy much less configuration space in the BGN than in the HN (in fact they normally occupy more, especially at low values of $`\gamma `$). However, for the most part the percentage of memory patterns which are stable is larger in the BGN than in the HN as long as $`\gamma 1`$. As with larger networks, the BGN becomes most similar to the HN when $`\gamma >1,`$ while lower $`\gamma `$ leads to increased pattern stability. The increased pattern stability at low $`\gamma `$ is associated, however, with smaller basins of attraction for the memory states and therefore with a greater volume of phase space occupied by spurious states.

## VI  <br>Discussion

We have studied the properties of the BGN at high memory loading $`\alpha `$. Our results can be summarized as follows: For high values of $`\gamma `$, such as 2, there is a first-order transition similar to that of a HN. For $`\gamma =2`$, the transition occurs at a critical loading $`\alpha _c0.1`$, which is lower than the critical loading of $`\alpha 0.148`$ for the HN. As $`\gamma `$ decreases, the critical loading increases and the phase transition evidently becomes weaker and eventually disappears. For $`\gamma =1`$ the critical loading is apparently $`\alpha _c0.17`$ (higher than for the HN) and the phase transition is much less pronounced for the finite-size networks we have studied. For $`\gamma =0.5`$ there is no evidence of a phase transition at all and patterns are stable with very few errors up to $`\alpha 0.3`$.

A phenomenon that occurs in the BGN much more than in the HN is the partial retrieval of a pattern, whereby the dynamics corrects some sign errors in a pattern without correcting all of them. This is especially noticeable in the case of low $`\gamma `$ and high $`\alpha `$. This phenomenon suggests that in this case the energy landscape in the vicinity of a stored pattern has the shape of a funnel rather than a smooth basin of attraction. By a smooth basin of attraction, we mean a connected neighborhood which is sloped toward an attractor and in which there are few local minima that might obstruct a trajectory once it has begun to flow toward the attractor. By a funnel we mean a region with an overall average slope toward an attractor which however contains many other local minima in which the trajectory might become stuck before reaching the bottom. The presence of these local minima is due to the bistability of the individual nodes and the local potential barriers against spin flips. These same potential barriers are also responsible for stabilizing the patterns. They make it less likely that crosstalk noise will induce an error in a pattern, but they also can inhibit the correction of an error which is present initially. Funnel-shaped landscapes were first examined in the context of protein folding dynamicsfunnel1 funnel2 . It was suggested that at a finite temperature such a landscape would allow the protein dynamics efficiently to find the global minimum of energy in spite of the presence of numerous local minima separated by potential barriers. If indeed the landscape of the BGN at low $`\gamma `$ forms a funnel, then it is possible that the introduction of some stochastic noise (i.e., finite temperature) could improve pattern retrieval by allowing trajectories to jump over the comparatively small potential barriers into lower energy minima, just as in the protein case. The effect of stochastic noise could be a fruitful subject for further study. An interesting question is whether there is an optimum level of noise which would improve the pattern retrieval ability of the BGN with low $`\gamma `$ while at the same time maintaining the larger storage capacity. Intriguing questions remain concerning the dynamics of the BGN at low $`\gamma `$. The apparent initial expansion of the basins of attraction of the memory states with increasing loading is counterintuitive, and the patterns visible in figure 11 hint at some dynamics which is not yet fully understood.

###### Acknowledgements.

This work was supported by Materials and Manufacturing Ontario (MMO), a provincial centre of excellence. Discussions with V. Chinarov are also acknowledged.