text:
# AD Leo from X-Rays to Radio: Are Flares Responsible for the Heating of Stellar Coronae?

## 1. Introduction

While large, episodic stellar flares heat coronal plasma rather efficiently up to 100 MK for minutes to hours, it is the potentially large number of small, and even ‘undetected’ flares (“microflares”, “nanoflares”, e.g., Parker 1988) that have recently attracted the attention of both solar and stellar research. Solar observations show that the flare occurrence rate is distributed in energy as a power law,

$$\frac{\mathrm{d}N}{\mathrm{d}E}E^\alpha $$
(1)
(e.g., Lin et al. 1984). The value of $`\alpha `$ determines the importance of the low- or high-energetic tail of the distribution: The integration of the energy over the energy distribution (1)

$$L_\mathrm{X}_{\mathrm{E}_{\mathrm{min}}}^{\mathrm{E}_{\mathrm{max}}}\frac{\mathrm{d}N}{\mathrm{d}E}EdE$$
(2)
with $`\alpha 2`$ produces arbitrarily large total emission rates if $`E_{\mathrm{min}}0`$ (microflares, nanoflares). There is ample evidence that the ensemble of flares play a fundamental role in the heating of (quiescent) coronae of magnetically active stars:

* Active stars emit strong gyrosynchrotron radio emission during quiescence: Evidence for accelerated (MeV) electrons as in solar flares (Güdel 1994).
* This non-thermal emission correlates with the overall X-ray radiation the same way as solar flares do (Güdel & Benz 1993; Benz & Güdel 1994).
* The optical U band flare frequency correlates with the quiescent X-ray stellar luminosity (Doyle & Butler 1985).
* Transition region lines show broadening probably related to explosive events (Wood et al. 1996).
* Small flare events with energies of the order of $`10^{27}10^{28}`$ ergs have become observable with the Hubble Space Telescope in cool M dwarfs, proving their ubiquity in active stellar atmospheres (Robinson et al. 1995).
* The structuring of X-ray emitting active coronae is reminiscent of the thermal structure of flares (Güdel et al. 1997) and may be explained by the superposition of a distribution of statistical flares (Kopp & Poletto 1993; Güdel 1997).
* Solar observations now show the importance of micro- and nanoflares in coronal energy release: $`\alpha =2.32.6`$ for microflares in the quiet corona (Krucker & Benz 1998).

Studies of statistical flare distributions have been rare in the stellar context, due to the paucity of relevant data sets. Collura et al. (1988) found a power-law index $`\alpha =1.52`$ from EXOSAT observations of dMe stars, while Osten & Brown (1999) report $`\alpha =1.6`$ for a sample of RS CVn binary systems observed with EUVE. In a series of papers (Audard, Güdel, & Guinan 1999; Audard et al. 2000; Güdel et al. 2000a) we have been investigating systematically the role of statistical flares in the overall coronal heating of active stars. From a large sample of EUVE observations, we found that

* relatively steep ($`\alpha =22.5`$) power laws dominate the flare rate distributions, indicating that small flares are important in these coronae.
* Statistical heating by flares leads to dominant coronal temperatures in agreement with the measurements.

## 2. New Observations of AD Leo: Flares over and over again

An extremely long EUVE observation of AD Leo was approved during EUVE’s cycle 7 between April 2 and May 16, 1999 (with a few minor time gaps) comprising 900 ksec of on-source exposure time. The DS instrument and the three spectrometers (SW, MW, LW) were used. Between May 1 and May 15, we obtained a total of 270 ks of exposure time with BeppoSAX, covering about 8 days within this interval. We obtained data from the LECS (0.1–10 keV), the MECS (2–10 keV), and the PDS (15–400 keV) instruments. On April 29, a 10 hr integration was carried out with the VLA. We used the 2 cm, 3.6 cm and 6 cm bands with the array in its D configuration. Two optical photometry observatories (Villanova University and Crimean Astrophysical Observatory) observed AD Leo in the U, B, V, R, and I bands during several nights in April and May.

AD Leo is a dM4.5e star with a rotation period of 2.7 d (according to our \[EFG\] measurements, $`P_{\mathrm{rot}}`$ appears to be 1.7 d only). An overview of the EUVE and BeppoSAX observations is shown in Figures 1$``$2 (also Güdel et al. 2000b). The EUVE light curve is extremely variable. Most of the fluctuations visible in Fig. 1 are real (the error bars being smaller than the visible fluctuations).

## 3. Measuring the Importance of Flare Heating

The light curves of AD Leo at hand are ideal for further investigation of statistical flares. We address two questions: (i) What is the statistical distribution of the visible flares in energy? (Given the lack of detailed spectroscopy for each flare, we approximate the total energy by the total number of counts detected, times a constant count-to-energy conversion factor.) (ii) Can the complete emission, including the apparently quiescent radiation, be explained by a superposition of statistical flares with a distribution that is compatible with (i)? We use two approaches to address these issues.

### 3.1. Kolmogorov-Smirnov Statistical Tests of the Flux Distribution

Here, we investigate the distribution of count rate values in the light curve. We proceed as follows:

We determine the average flare profile through autocorrelation analysis from the observations.

We simulate light curves that are composed of statistical flares. These are distributed in energy according to a power law (1). We use 5480 bins, which is about ten times the number of bins in the EUVE light curve (if binned to one data point per satellite orbit), i.e., we simulate ten statistical realisations of our observation.

The simulation is renormalized to the observation. To this end, we determine the cumulative distribution of flux values both for the observation and for the simulations (i.e., number of bins with a flux exceeding a given flux - see Figures 3$``$4). We normalize the model flux to the observed flux by adjusting the middle portion of the cumulative distribution, which lies above the noise level but at fluxes that are frequently attained.

Statistical noise corresponding to the observations is added.

We then perform a Kolmogorov-Smirnov (KS) statistical test between simulation and observation by comparing the number distribution of flux values in the available bins. Option: The dominant ‘quiescent’ part which is difficult to model can be subtracted beforehand to avoid problems with overlapping flare wings, rotational modulation effects, some true long-term variability etc.

Figures 3$``$4 show best-fit realisations for different $`\alpha `$ values. The top panel shows the data used, the second row refers to the optimum $`\alpha `$ found, while the other panels show cases for too large and too small $`\alpha `$. In each case, the solution with the highest confidence was identified by varying the number of flares in the simulation.

The left panels show the first 500 bins of the simulated light curves, while the right panels show the cumulative distributions of count rates attained for the observation (black, solid) and the simulation (red, dashed). We do not use the lowest 10% of the distribution since those flux levels may be influenced by slowly varying emission in the ‘quiescent’ emission. The vertical differences between the simulated and the observed distribution are indicated dotted (blue) around the zero level: deviations to positive values indicate a locally too soft simulated distribution, negative deviations indicate locally too hard simulated distributions. The maximum deviation is indicated by a blue vertical bar. The (green) dashed horizontal line indicates the lower threshold for identifying the maximum vertical difference.

Since the last segment in the EUVE DS light curve suffered from incursion into the DS dead spot and from elevated particle radiation, we omitted that segment in our analysis. We plot in Fig. 5 the best results for different $`\alpha `$ values, for three cases: (i) Modeling of the EUVE DS emission above the ‘quiescent’ level (A); (ii) modeling of the complete EUVE DS emission (B); (iii) modeling the complete BeppoSAX LECS light curve (C). Acceptable models ($`>`$90% confidence) are found for $`\alpha =2.22.5`$.

### 3.2. Effect of Weak Flares

The high time-resolution of the EUV events ($`1ms`$) also affords the possibility of detecting weaker flares in the light-curve (see Kashyap et al. 2000). By comparing the distribution of photon arrival-time differences in the data stream with that expected from a model of flare distributions with specific power-law indices, we can account for the weak (but numerous) flares that contribute to the emission. Specifically, we

* adopt a power-law flare-distribution model (Equation 1) described by the power-law index $`\alpha `$ and the normalization $`N_f`$, and generate a high-time-resolution light-curve via Monte-Carlo simulation;
* add a constant component $`N_b`$ to the flare light-curve to account for true steady emission, flare emission too weak to be distinguished from steady emission, and background;
* obtain a set of event times from the model light curve;
* compare the observed distribution of arrival-time differences with that derived from the generated model event-times and compute a test statistic similar to $`\chi _r^2`$; and
* carry out the comparison over a grid of parameter values ($`\alpha `$, $`N_f`$, $`N_b`$) to find the best-fit and confidence range.

We have run the above algorithm on the first half of segment II of the AD Leo EUVE observation (see Figure 1) and find that $`\alpha =2.1\pm 0.05`$, a smaller value than found above, but yet $`>2`$. These data are clearly dominated by the large flare; an analysis of longer segments, which is expected to reduce the bias due to this flare and provide a more realistic assessment of the value of $`\alpha `$, is in progress.

## 4. Discussion and Conclusions

The analysis performed so far on EUVE light curves (Audard et al. 1999; 2000; also Kashyap et al. 2000) and the new results presented here suggest a distribution of flare energies according to a power law with a steep index: $`\alpha 23`$. If the flare rate distribution continues down to levels of average solar flares, then the complete stellar corona could be heated solely by the energy released in statistical flares.

The results suggest values of $`\alpha 2.12.5`$ independent of whether only the detected flares are investigated, or whether the complete emission is simulated. This suggests that (i) the extrapolation of the flare distribution to undetected flare energy levels could add large amounts of emission, and (ii) that the complete statistical distribution of count rate levels is indeed compatible with such a distribution. The two points require a lower threshold for the power-law flare distribution to confine the emission to the observed level. This has implicitly been taken into account as our simulated power-law distributions contain flares only with energies above a pre-set lower threshold. The value of the latter is determined by the normalization to the observation.

Flares could then also explain why the temperatures of active stars become increasingly hotter with ‘increasing activity’: The higher flare rate keeps more (high-density) plasma at hot temperatures, and thus the hotter plasma increasingly dominates the X-ray emission (Güdel 1997; Audard et al. 2000).

#### Acknowledgments.

We thank the CEA/EUVE, VLA, and BeppoSAX staff for their efforts to coordinate these observations. M. A. acknowledges support from the Swiss National Science Foundation, grants 2100-049343 and 2000-058827. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. SRON is financially supported by the Netherlands Organization for Scientific Research (NWO).

## References

Audard, M., Güdel, M., & Guinan, E. F. 1999, ApJ, 513, L53

Audard, M., Güdel, M., Drake, J. J., & Kashyap, V. L. 2000, ApJ, 541, 396

Antonucci, E., et al. 1984, ApJ, 413, 786

Benz, A.O., & Güdel, M. 1994, A&A, 285, 621

Collura, A., Pasquini, L., & Schmitt, J. H. M. M. 1988, A&A, 205, 197

Doyle, J. G., & Butler, C. J. 1985, Nature, 313, 378

Güdel, M. 1994, ApJS, 90, 743

Güdel, M. 1997, ApJ, 480, L121

Güdel, M., Benz, A. O. 1993, ApJ, 405, L63

Güdel, M., Guinan, E. F., & Skinner, S. L. 1997, ApJ, 483, 947

Güdel, M., et al. 2000a, ApJ, in preparation

Güdel, M., Audard, M., Guinan, E. F., Mewe, R., Drake, J. J., & Alekseev, I. Y. 2000b, in Eleventh Workshop on Cool Stars, Stellar Systems, and the Sun. ed. R. J. Garc a López, R. Rebolo, & M. R. Zapatero Osorio (San Francisco: ASP), in press

Kashyap, V. L., Drake, J. J., Audard, M., & Güdel, M. 2000, Bull. AAS 196, #13.05

Kopp, R. A., & Poletto, G. 1993, ApJ, 418, 496

Krucker, S., & Benz, A. O. 1998, ApJ, 501, L213

Lin, R. P., Schwartz, R. A., Kane, S. R., Pelling, R. M., & Hurley, K. 1984, ApJ, 283, 421

Osten, R. A., & Brown, A. 1999, ApJ, 515, 746

Parker, E. N. 1988, ApJ, 330, 474

Robinson, R. D., Carpenter, K. G., Percival, J. W., & Bookbinder, J. A. 1995, ApJ, 451, 795

Wood, B. E., Harper, G. M., Linsky, J. L., & Dempsey, R. C. 1996, ApJ, 458, 761


format: text

text:
# Dynamics of directed graphs: the world-wide Web

## I Introduction

Understanding the evolution of complex networks  remains a challenging problem of modern statistical physics. Modeling the dynamic structure of these networks has significant practical applications, for instance in the design of search and control algorithms and in predicting emergence of new phenomena in information networks . An evolving network consists of increasing number of nodes, which are connected by links according to certain rules specific to each network. Although the spatial distribution of nodes in a network is random, a different structure can be observed when the distribution of links is considered. In this respect two major groups  are: (1) random networks, with the number of links attached to a node fluctuating around an average value, and (2) networks with hierarchical structure of links, manifested in the power-law decay of the distributions of node ranks. The structure of connections has immediate impact on the accessibility of nodes and on functional stability of the network. For instance, the hierarchically connected networks are shown  to be robust to large failure rates, however, they are vulnerable to removal of a few ’key nodes’. The inhomogeneous structure of connections has been found in metabolic cycles , in various social networks , including the science citation network , and the Internet  and the world-wide Web .

The world-wide Web is not a separate physical networks. Rather it is a subset of the Internet that uses a specific protocol (hypertext transfer protocol) and directed links to access a variety of data, representing an example of distributed client/server computing . Based on these properties a hierarchical structure of links emerges that is peculiar to the world-wide Web. Recently large-scale measurements in the world-wide Web that involved $`10^8`$ nodes reveal  that: (a) Both the distributions of incoming links and distribution of outgoing links follow a power-law behavior with different exponents; (b) The distribution of size of connected component—number of linked nodes in the web crawls—also exhibits scaling behavior; and (c) The Web shows a very intricate structure of connections, with a large strongly connected component in the center.

Recently a great deal of effort was devoted to understand emergence of the scale-free structure of links in various complex networks . An important ingredient of the dynamics—preferential attachment  was shown to lead to the power-law distribution for incoming links (the outgoing links are fixed by the rules of the model ). The model of preferential attachments for incoming links has been generalized to include initial attractiveness of nodes in Ref. , where the authors proved analytically the existence of the scaling region for large evolution times. Recently it was demonstrated in Ref.  how a giant connected component can occur in the model of random directed graphs assuming arbitrary and statistically independent distributions of incoming and outgoing links. The predictions of the random graph theory, however, show substantial quantitative differences compared to the properties of the real world-wide Web .

In this work we simulate growth and response of a directed graph with the dynamic rules that are motivated by prominent features of the world-wide Web. These dynamic rules yield the statistically correlated distributions of outgoing and incoming links. Another important feature of the model is that the links between pairs of nodes are not fixed in time, but may vary on the time scale of the network’s evolution (updates of links). This property, that may be found in some metabolic cycles, is not shared with many other social networks, whose physical links are either fixed or vary on much slower time scale. We argue that the (at present high) frequency of updates of the outgoing links, which is peculiar to the agents in the Web, is essential for the observed inhomogeneous structure of connections. We find a consistent agreement with the recent empirical data on the real world-wide Web both in the structure of links and in the response of the network.

The paper is organized as follows: In Sec. II we introduce the growth rules of the network and determine the emergent rank distributions of outgoing and incoming links. In Sec. III we study the temporal fractal structures in the growth phase of the network and determine the corresponding distributions of the first-return links. Sec. IV is devoted to investigation of the connected components on the fully grown network and fractal properties of noise. A short summary and the discussion of the results is given in Sec. V.

## II Growth model and rank distributions

We suggest the following simplified set of rules to incorporate basic features pertinent to the real world-wide Web: (i) Directed nature of linking. The world-wide Web represents a directed graph with nodes corresponding to Web pages, and arcs corresponding to hyperlinks between pages . (j) Growth and rearrangements at unique time scale. At each time unit $`t`$ a new node $`i=t`$ is added to the network (growth) and a number $`M(t)`$ of new links are distributed among the nodes (update of links) following two rules specified below. A fraction $`f_0(t)\alpha M(t)`$ of new links are outgoing links from the new added node $`i=t`$, whereas the rest $`f_1(t)(1\alpha )M(t)`$ are the updated links at other nodes in the network. Hence, the relevant parameter in the model is the ratio of updated and added links at each time step, i.e., $`\beta f_1(t)/f_0(t)=(1\alpha )/\alpha `$, which is independent of the actual number of links $`M(t)`$. Thus, in the model $`M(t)`$ represents a net increment of the number of outgoing links at time step $`t`$. Variations in $`M(t)`$ are, in principle, not restricted, being caused by adding new connections between old nodes (resulting in the increase of the number of links), and/or by removing some earlier links (causing decrease of the number of links). We assume that variations in $`M(t)`$ are such that we can define an average value $`M\overline{M(t)}`$, which can be considered as a constant in first approximation. In practice, the number of nodes and the number of links in the network increases with time, so that reasonable values for the average $`M`$ are positive. For consistency, we keep $`M=1`$ throughout this work. Different values of $`M`$ do not affect the universal properties in the scaling region (i.e., for large evolution times). Taking much larger $`M`$ values, however, considerably increases computation time. (k) Preferential update (k1) and preferential attachment (k2). These properties represent a paradigm of social behavior, for instance preferential attachments are driven by “popularity” of a node, as discussed in earlier models of evolving social networks. In addition, here we have that not all nodes are getting updated at every time step, rather only a few nodes update at a time. Moreover, some of the nodes update outgoing links more frequently than others. This features can be formulated in terms of probabilities as follows: Where do the links come from? We assume that apart from the new node, updating at time $`t`$ occurs with larger probability at most active nodes (preferential activity), i.e., an outgoing link from the node $`ki`$ appears with the probability

$$Prob1=\frac{\alpha M+q_{out}(k,t)}{(1+\alpha )Mi}.$$
(1)
A new link goes to the node $`n`$ (whose age is $`tn`$) with the probability

$$Prob2=\frac{\alpha M+q_{in}(n,t)}{(1+\alpha )Mi},$$
(2)
i.e., a “popular” node attracts even more links. Here $`q_{out}(n,t)`$ and $`q_{in}(n,t)`$ are the dynamically varying number of outgoing and incoming links, respectively, at the node $`n`$ at current time $`t`$. It is assumed that at the time of addition of a node $`i`$ to the network $`q_{out}(i,i)=q_{in}(i,i)=0`$. Note that the second rule in Eq. (2) is formally equivalent to the model of preferential attachment for incoming links . However, in our model the role of the parameter $`\alpha `$ in this equation is precisely determined through the rule in Eq. (1), which regulates outgoing links: (1) For $`\alpha 1/(1+\beta )<1`$ a fraction $`(1\alpha )M`$ of outgoing links in Eq. (1) refers to updates at earlier nodes; (2) For $`\alpha =1`$ (equivalent to $`\beta =0`$) no updates are allowed. New links originate only from the new added node; (3) Formally one can apply Eqs. (1-2) for $`\alpha `$ strictly larger than unity, however, in this region $`\beta `$ becomes negative, suggesting that the number of links in the network is not conserved. Therefore only the incoming links (2) can be counted correctly when $`\alpha `$ is in the interval $`\mathrm{}>\alpha >1`$; (4) In the mathematical limit $`\alpha \mathrm{}`$, or in practical cases when $`\alpha \overline{q}/M`$, the effects of the dynamical quantities $`q_{out}(k,t)`$ and $`q_{in}(n,t)`$ become negligible, and we recover the case of fully random network.

In what follows we would like to demonstrate that the features (i), (j), (k1) and (k2) explained above are essential for the behavior observed in the Web. Earlier models  concentrated only on the preferential attachment rule (k2), while largely neglecting properties (i), (j) and (k1). The most complete account of the preferential attachment model is given in Ref. . On the other side, the random graph model  considered the directed nature of network (i), while neglecting the details of the linking rules (j), (k1) and (k2).

Using the rules in Eqs. (1)-(2) we grow a large network of $`N=10^6`$ nodes. We measure the ranking of nodes in this network by the distributions of outgoing links and incoming links, $`P(q_{out})`$ and $`P(q_{in})`$, respectively, which are shown in Fig. 1 for parameter $`\beta =3`$, corresponding to $`\alpha =0.25`$ in Eqs. (1-2). For comparison we also simulate the distributions of outgoing and incoming links in the case of fully random directed network. For a range of values of their arguments the cumulative distributions for finite $`\beta `$ values exhibit a power-law behavior according to

$$P(q_{out})q_{out}^{(\tau _{out}1)};P(q_{in})q_{in}^{(\tau _{in}1)}.$$
(3)
By increasing the parameter $`\alpha `$ in Eqs. (1)-(2) in the range $`(0,1)`$, corresponding to decrease of the relative fraction of updated links $`\beta `$ in the interval $`(\mathrm{},0)`$, the slopes of the distributions increase. The scaling exponents $`\tau _{in}`$ and $`\tau _{out}`$ given in the inset to Fig. 1 are parametrized by $`\alpha `$. Our numerical results for $`\tau _{in}`$ agree with the exact solution  $`\tau _{in}=2+\alpha `$ for the incoming links. For the distribution of outgoing links no analytical results are available. Our numerical results in the inset to Fig. 1 can be well approximated with the linear dependence $`\tau _{out}2+3\alpha `$. Therefore, the difference between the two exponents $`\mathrm{\Delta }\tau _{out}\tau _{in}2\alpha `$ persists for any finite $`0<\alpha <1`$. Comparing these results with the available empirical data we estimate the parameter $`\alpha =0.22\pm 0.1`$. Here we used the following values as the experimental estimates for the exponents $`\tau _{in}=2.16\pm 0.1`$, obtained as the average value from the data in  and last reference in , and $`\tau _{out}=2.62\pm 0.1`$ from data in , by fitting only the straight part of the curves while avoiding the noisy tails and the curvature at small $`q`$. It is clear that these values are not definitive. More careful fits of the data and further measurements are necessary in order to reduce the experimental error bars. Hence, the approximative value of the control parameter $`\beta `$ in the current state of the Web is $`\beta 3`$, i.e., in the average three updated links come to one added link at each evolution step. We now discuss the properties of the network for this particular value of the control parameter $`\beta `$.

## III Dynamic fractal structures

In the dynamical systems with large-scale organization the occurrence of algebraically decaying distributions can be linked to the formation of certain fractal structures  in the system. In the case of complex evolving networks, however, the existence of such structures is less clear, due to the spatial randomness of the graph. To show that the power-law decay of the distributions of node ranks, shown in Fig. 1, is associated with certain fractal structures on the network we examine the temporal pattern of linking between nodes. A part of the pattern is shown in Fig. 2 (top panel). Clustering in the patterns indicate frequently active nodes, whereas voids of various sizes indicate lesser activity at corresponding nodes. Visually heterogeneous picture suggests that both pattern of origins of links and pattern of targets have fractal structure. The fractality of these patterns can be quantified, for instance, by computing the distribution of time intervals between two successive linking at a given node (return time). The distribution of time intervals $`\mathrm{\Delta }t`$ between two successive linking to a selected node in the network $`P(\mathrm{\Delta }t)`$ measured for the network of $`10,000`$ nodes and the parameter $`\beta =3`$ is given in Fig. 2 (lower panel). The algebraic decay of the distribution is a signature of the underlying fractal structure of the pattern. The distribution of time intervals between two successive links from a selected node also shows a power-law tail for large $`\mathrm{\Delta }t`$. At small intervals $`\mathrm{\Delta }t`$ the pattern is nearly random. The two distributions coincide, indicating mutual correlations in the patterns, for time intervals $`\mathrm{\Delta }t>100`$ (that are accessible for large evolution times $`t`$). To emphasize the relation to Fig. 1 we recall that the number of links of either kind at a node accumulates with time. Hence, the regions of large time intervals between successive linking at a node, $`\mathrm{\Delta }t`$ in Fig. 2, correspond to small number of links at that node, $`q`$ in Fig. 1, and vice versa.

## IV Dynamic critical states

The scaling behavior of the rank distributions can be examined in terms of the potential dynamic critical states by measuring the statistics of triggered avalanches  on the network. We consider an analog of the avalanche size in the networks: the size of a cluster of nodes which are physically accessible along directed links starting from a random node in the network. A network is grown using the rules in Eqs. (1)-(2). We then select a random node and make a list of nodes which are connected by outgoing links from that node. In the next step we make a new list following outgoing links from the nodes on the previous list, keeping only the nodes which are different in the two lists, and so on. The process ends up when no new nodes can be reached—the list is empty. The number of different lists before an empty list occurs can be recognized as the depth of the connected component. Technically, we preserve ranks $`q_{in}`$ and $`q_{out}`$ for each node from the growth phase and, in the case of small networks, the exact physical links between the nodes. In large networks only ranks are preserved and links are searched using the rule in Eq. (1). For a large ensemble of networks the results are expected to be statistically the same. In this way the size of a cluster of linked nodes represents a response of the system on the random excitation. The process is reminiscent of the invasion percolation on a random graph, where a new edge $`j`$, which already does not belong to the graph, is invaded along an outgoing hyperlink $`ij`$ from the node $`i`$ if the node rank $`q_{out}(i)`$ exceeds unity. The probability to add the edge $`j`$ to the graph is given by the preference rule in Eq. (1). In a more familiar case of the invasion percolation in physical wetting on the lattice  that probability is given by the least resistance rule. To our knowledge the problem of invasion percolation on random graphs has not been studied so far .

The differential distribution of size of connected clusters of nodes is given in Fig. 3 for different values of the parameter $`\beta `$. The distribution of depth of the same clusters is given in the inset to Fig. 3. In the case $`\beta =3`$ we find that small clusters follow a power-law distribution with the exponent $`\tau _s=2.79\pm 0.12`$, comparable with the measured  value $`\tau _s=2.52`$. (Note that in the case of Internet the corresponding exponent is estimated  as $`\tau _s^I=1.90`$.) In Fig. 3 the occurrence of the peak in the distributions at large clusters $`S_010^2`$ indicates that a large number of the avalanches have approximately same size. This is a signature of the existence of a large subset of strongly connected nodes (a giant component): Once a ‘crawl’ enters the subset of strongly connected nodes, it explores it entirely. By taking a larger number of nodes $`N`$ the position of the peak moves towards larger values approximately as $`S_0N/7`$.

Before discussing the distributions in Fig. 3, we demonstrate that the network with the hierarchical structure of links is characterized by a fractal noise. In the search of a connected cluster described above we examine detailed variation of the number of nodes added to the cluster at each step of investigation. We obtain the noisy signal shown Fig. 4. Properties of the signal vary with the control parameter of the dynamics $`\beta `$. Considering the number of steps as a total elapsed time of investigation, the Fourier spectrum of the signal is shown in the top panel in Fig. 4. It shows a region of correlated behavior with power-law decay between the upper cut off at high frequencies and lower cut off (due to finite size of the network). In addition, the spectrum in the case of supercritical network exhibits a peak at a characteristic frequency $`f_0`$, which is absent in the case of the critical parameter $`\beta =\beta _C(N)`$ (cf. Fig. 4).

## V Discussion and conclusions

In the theory of critical states , appearance of the peak in the distribution of avalanches, such as the distributions shown in Fig. 3, indicates that the system is supercritical. In the critical state the absence of any characteristic scale is manifested in the purely algebraic decay of the distribution of cluster size until a cut-off, which depends on the size of the network $`N`$. On the other hand, the form of distribution in Fig. 3 suggests that a critical point exists that can be reached by varying a relevant parameter of the dynamics . In fact, by decreasing the parameter $`\beta `$ we see in Fig. 3 that the slope $`\tau _s`$ decreases and the peak eventually disappears at a critical value  of the parameter $`\beta =\beta _c0.081`$. Beyond the critical value of the control parameter the scaling behavior of the size and depth of connected components is entirely lost. In the critical state the scaling exponents are (cf. inset to Fig. 1): $`\tau _{in}2.925`$, $`\tau _s=\tau _d=1`$. The distribution of outgoing links shows very sharp decay, the exponent $`\tau _{out}`$ is difficult to measure. It should be stressed that the absence of the peak in the critical state suggests that no giant component can be formed. Rather, the network consists here of many small groups of well interlinked nodes and percolating directed links between these groups.

It should be stressed that the special limit of our model $`\beta =0`$ (i.e., $`\alpha =1`$), corresponding to the original model of preferential attraction of incoming links proposed by Barabasi, Albert and Jeong , belongs to an entirely different class of scaling behavior as regards the distribution of size and depth of connected clusters: these distributions are flat between lower and upper cutoffs (see Fig. 3). The only distribution which enjoys a power law in this limit is the distribution of incoming links, which has the exact  exponent $`\tau _{in}=3`$ (see also inset to Fig. 1). In this case ($`\beta =0`$) updates at already existing nodes are no longer possible. This feature—freezing of the outgoing links—as well as quantitative disagreement in the exponent $`\gamma `$ for distribution of incoming links, makes the model of Ref.  inappropriate as a model of the world-wide Web dynamics. A realistic network with frozen outgoing links is the network of scientific citations , where physical links, corresponding to the cited references in already published papers, remain fixed in time.

In conclusion, we have demonstrated that certain salient features of the dynamics of the world-wide Web require more careful modeling, compared to models of a generic complex evolving network with preferential attraction of links  and models of random graphs with rewiring . The dynamic structure and functioning of the world-wide Web is deeply rooted in the activity of the agents who are creating the outgoing links. The hierarchical structure of outgoing links, which is documented by measurements in the real world-wide Web , is related to the bias activity of the agents in the rule (1) of our model. A random selection of the active agent (see curve (c) in Fig. 1) fails to describe the distribution of outgoing links in the real Web. Temporal variation of the outgoing links inside the network, i.e., updates of links, has twofold consequences on the global structure of the network: (i) When the updates of the outgoing links are allowed (i.e., when $`\beta `$ strictly larger than zero), the structure of incoming links qualitatively changes, the exponent $`\tau _{in}(\beta )<3`$. The distribution of incoming links in our model coincides with the analytical results of Dorogovtsev et al.  in the scaling region. In addition, the distribution of size and depth of connected clusters becomes hierarchical when $`\beta >0`$, in agreement with results found in the real Web, which is not the case when no updates are allowed ($`\beta =0`$). (ii) Varying the frequency of updates $`\beta `$ implies changes in the Web structure, both in the outgoing and incoming links. This has immediate impact on the accessibility of nodes. The correlations between the outgoing and incoming links suggests that the local structure of the network is qualitatively different compared to the case without updates. Here we did not study in detail the probability that two nodes are linked (clustering coefficient). Rather we studied the physical properties of the network which make the background of the observed behavior: fractal temporal patterns and number of nodes that join the cluster at one time unit.

By comparison of the simulated results and the data obtained in the real world-wide Web we obtained a systematic agreement both for the distributions of outgoing and incoming links and for the size of connected components when a single control parameter $`\beta `$ is fixed to $`\beta =3`$ within error bars of the data. To our knowledge no previous model claimed to describe the world-wide Web achieved such degree of consistency. This makes us believe that the present simplified model takes into account some basic features of the dynamics of Web. It would be interesting to estimate $`\beta `$ by directly measuring the average number of updated links in the world-wide Web relative to the number of added links originating from each added node.

The structure of links shown by the distributions in Figs. 1 and 2 is closely related to the evolution of the number of links at a given node and to other local properties of the network. A detailed study of these properties within the present model requires additional work that remains to be done in the future. Here we expect, based on the analytical results for the incoming links of Ref. , that in our model the average number of links $`<q_{in}(i,t)>`$ at a node $`i`$ will decay in time $`ti`$ as $`<q(i,t)>(i/t)^{\gamma _{in}}`$, where the exponent $`\gamma _{in}`$ is given by the exact scaling relation  $`\gamma _{in}=1/(\tau _{in}1)`$, leading to $`\gamma _{in}=1/(1+\alpha )0.86`$. Assuming that the density of the outgoing links at a given node also exhibits scaling behavior, as Figs. 1-2 suggest, we can predict a slower decay for the average number of outgoing links at a given node. The expected exponent is $`\gamma _{out}1/(1+3\alpha )0.6`$.

The simulation results that we reported here suggest that the emergent structure of links in the world-wide Web is strongly related to the updating policy of the agents: who updates and how often. It remains to understand the potentially more intricate reverse effect: how the amount of information currently stored in the Web influences the conduct of the agents, with implications on self-tuning of the control parameter $`\beta `$. In our model the structure of the network at $`\beta >\beta _c(N)`$ appears to be supercritical, where $`\beta _c(N)<0.1`$ in a large network. Therefore, it is plausible to expect that the evolutionary selected values of $`\beta `$ will be much smaller than the currently estimated value $`\beta 3`$, if the scenario of self-tuning of the control parameter is active in the real Web. In this respect the current state of the real Web can be regarded as a transient rather than a stationary state.

Finally, our results support the conclusion that in the directed graphs two growth rules are necessary to describe the dynamics of the outgoing and incoming links, respectively. In the case of the world-wide Web the statistically correlated distributions of outgoing and incoming links appear as a fundamental feature of the evolution of the Web. This conclusion may serve as a starting point for the future modeling of the real world-wide Web in terms of the master equations.

###### Acknowledgements.

This work was supported by the Ministry of Science and Technology of the Republic of Slovenia. I am grateful to Deepak Dhar, Vyatcheslav Priezzhev and Álvaro Corral for their response and for the constructive criticism on this work.


text:
# Close stellar encounters with planetesimal discs: The dynamics of asymmetry in the 𝛽 Pictoris system

## 1 Introduction

The starlight of $`\beta `$ Pic scatters off dust grains in its edge-on circumstellar disc (Smith & Terrile 1984). Optical imaging reveals that the disc exhibits asymmetrical structure with respect an ideal disc on projected radial scales from $``$50$`\mathrm{au}`$ (Burrows et al. 1995, Mouillet et al. 1997) to $``$1000$`\mathrm{au}`$ (Smith & Terrile 1987, Kalas & Jewitt 1995). Kalas & Jewitt (1995) presented sensitive $`R`$-band images that showed the northeast (NE) extension could be detected out to $`790`$$`\mathrm{au}`$, whereas the southwest (SW) extension was detectable out to just $`650`$$`\mathrm{au}`$, corresponding to a $``$20 per cent length asymmetry. More recently, the longer NE midplane extension has been discovered to carry several brightness enhancements from $``$500 to $`800`$$`\mathrm{au}`$ (Kalas et al. 2000; hereinafter KLSS), that have been attributed to a series of rings within the disc’s midplane. Similar features in the SW extension were not detected and so eccentric rings were hypothesised.

KLSS presented a dynamical model in which the stellar fly-by scenario (Kalas & Jewitt 1995) for generating the disc asymmetries was applied to the $`\beta `$ Pic system. Their finding was that the asymmetry types present in the $`\beta `$ Pic disc have analogues in the dynamical response of a particle disc to a perturbing stellar fly-by encounter. Furthermore, they demonstrated that these features could occur simultaneously with the formation of eccentric circumstellar rings. A simple analysis of the perturbed disc model was found to be in good quantitative agreement with the observations of both $`\beta `$ Pic’s disc asymmetries and the midplane features; for particular values of the viewing angle and epoch, and fly-by parameters. In this paper we present new observations and measurements of the length asymmetry of $`\beta `$ Pic’s scattered-light disc, and then go on to simulate the dynamics of a stellar flyby encounter with an initially symmetrical disc of particles.

Similar numerical modelling to that presented here and in KLSS has been considered previously in the context of angular momentum transfer between the perturber’s orbit and a protostellar disc (Clarke & Pringle 1993, Hall, Clarke & Pringle 1996), the dynamics of the Edgeworth-Kuiper belt (Ida, Larwood & Burkert 2000), and the tidal interaction of spiral galaxies (Toomre & Toomre 1972). Here we shall consider a wider range of encounter velocities (from bound orbits through to hyperbolic trajectories) and stellar masses (from M dwarfs to A stars, when scaled to the mass of $`\beta `$ Pic) than has been studied before in any context. We also include a discussion of the orbital dynamics of captured particles with respect to the perturber. This paper can be considered as a companion to KLSS since here we also supplement their finding of circumstellar ring formation in similar simulations with an explanation of the dynamical origin of the features.

As in other works (Whitmire, Matese & Tomley 1988, Mouillet et al. 1997, KLSS) we assume that the simulation particles’ distribution represents the distribution of an underlying disc of planetesimals, which are parent bodies whose infrequent and therefore dynamically inconsequential collisions supply the dust scattering surface that is observed in the real system. A more complete model should include consideration of grain removal and generation processes in determining the appearence of the perturbed disc (reviewed by Backman & Paresce 1993). However, we defer treatment of those issues to future work, and focus here on the first-stage problem of the dynamics. This stage is important in deducing the mass and orbital parameters of the postulated stellar perturber, for which it is possible to perform star catalogue searches (e.g. Kalas, Deltorn & Larwood 2000), and we proceed to do this by examining the length asymmetry, and other measures of the disc response, as a function of stellar encounter parameters. However, we choose to examine the length asymmetry as the principal diagnostic measure of the disc response since it is the most distinct such feature both in the observations and the dynamical models. In our future work we shall also consider the subtler asymmetries in greater detail.

In Section $`2`$ we present deep optical coronagraphic images of the disc. In Section $`3`$ we describe the results of our numerical investigations of the dynamics of tidally disrupted particle discs, and provide general discussion of the induced asymmetries and particle stripping. In Section $`4`$ we discuss our findings and summarise the conclusions.

## 2 New observational results

New optical coronagraphic observations of the environs of $`\beta `$ Pic were obtained with the University of Hawaii $`2.2`$ m telescope on 2000 January 31. A key goal was to maximize sensitivity to the outer regions of the disc. The coronagraph used has 12 interchangeable, circular, focal-plane occulting masks, from which we selected a mask $`12.5`$ arcsec in diameter. This particular mask size is approximately two times larger than the minimum size allowable by the seeing conditions, and it partially obscures the image out to $`120`$$`\mathrm{au}`$ projected radius about $`\beta `$ Pic. However, it also minimizes the scattered light in the field that dominates the background noise. Additionally, the occulting mask increases the efficiency of the observations by maximizing the integration time permitted before CCD saturation occurs just beyond its edge.

We obtained 122 integrations of 30 seconds each through an $`R`$-band filter ($`\lambda `$0.67 $`\mu `$m). Re-imaging optics provided a pixel scale of $`0.407`$ arcsec pixel-1 on a $`2048\times 2048`$ CCD. The field was masked by the coronagraphic focal plane assembly which has a circular clear aperture with diameter $`336`$ arcsec. The airmass of observation ($`>`$3.0) combined with high winds, produced image quality varying from 1.0 to 1.8 arcsec. For final data analysis, we selected 77 frames with the image quality characterized by stellar full-width at half-maximum (FWHM) $`1.2`$ arcsec.

Data were bias-subtracted, flat-fielded, and sky-subtracted using the median sky value in the corners of the field for each frame. Comparison stars were also observed for the purpose of subtracting the $`\beta `$ Pic PSF. However, this subtraction step adds noise and so we chose to remove the background using an artificial point-spread function (PSF). The artificial PSF is a figure of rotation for a seventh-order polynomial that is a least squares fit to $`\beta `$ Pic’s PSF sampled in radial directions perpendicular to the disc midplane.

Fig. $`1`$ shows the final image of $`\beta `$ Pic with effective integration time 2310 seconds. The 3$`\sigma `$ background noise corresponds to a sensitivity of $`R=+25.8`$ mag. The midplane of the SW extension is difficult to trace visually as a distinct linear feature beyond $``$55 arcsec ($`1062`$$`\mathrm{au}`$), whereas the NE midplane is traceable at up to $``$95 arcsec ($`1835`$$`\mathrm{au}`$) radius from $`\beta `$ Pic. By these measurements the length asymmetry is approximately $`70`$ per cent.

A key problem in comparing the lengths of the two extensions is that the morphology of the SW extension departs significantly from a linear feature, becoming vertically distended and amorphous at large radii. To compensate for this difference in vertical structure we rotate the image such that the disk midplane lies along pixel rows and bin the data in the vertical direction. We do this by summing the light at each pixel location with the light from its two vertically-adjacent pixels. This vertical binning procedure smoothes the data and equalizes the NE and SW FWHM measurements when taking vertical cuts at given projected radii along the midplane extensions. Smoothed contours of the resulting image are shown in Fig. $`2`$. Here, the 25.0 mag arcec-2 contour traces the NE extension to 95 arcsec radius, and the SW disk extension is traced to $``$75 arcsec ($`1450`$$`\mathrm{au}`$). The length asymmetry by this measurement is therefore $``$25 per cent, which is fractionally larger than the length asymmetry measured by Kalas & Jewitt (1995) at about a half of the radial scale probed here. The asymmetric vertical morphology evident in the Fig. $`2`$ contours also confirms the *butterfly asymmetry* (see Section $`3.3`$) reported by Kalas & Jewitt (1995).

### 2.1 Summary

The new data reported above show that the $`\beta `$ Pic disk has a significantly greater radial extent ($`1835`$$`\mathrm{au}`$ for the NE disk extension; $`1450`$$`\mathrm{au}`$ for the SW) than has been previously measured. Smith & Terrile (1987) reported the maximum extent of $`\beta `$ Pic’s disc at $`1300`$$`\mathrm{au}`$ (i.e 67 arcsec) for the NE midplane at optical wavelengths. However, this was known as a sensitivity-limited value rather than the true physical extent of the disc. Here too we are unable to define the outer disc extents above the level of the background noise.

The existence of a disc miplane out to $``$103$`\mathrm{au}`$ radius, as well as the $``$102$`\mathrm{au}`$ vertical thickness of the SW extension (see Fig. $`4`$ in Kalas & Jewitt 1995), is unexpected given our current understanding of planetesimal formation as applied to the Solar System (reviewed by Weidenschilling & Cuzzi 1993). Thus we must conclude that planetesimals were created at much smaller radii than they are presently inferred to exist in $`\beta `$ Pic, and subsequently suffered redistribution. Moreover this redistribution must allow for the existence of a lopsided and vertically flared disc. Global redistribution of this kind is a general outcome of a close stellar fly-by encounter (KLSS). Below we shall present numerical investigations of the dynamics involved.

## 3 Dynamical models

In the fly-by model presented in KLSS a transient system of circumstellar eccentric rings was found to form in the perturbed disc of particles, as shown in Fig. $`3`$. Fly-by parameters, and viewing angle and epoch, were chosen to give the most satisfactory fit possible to the positions of the midplane brightness peaks found in the observations, and to the large-scale disc asymmetries. The model presented in KLSS will serve here as a reference case with which we compare the model behaviour as we vary the parameters about its particular values. In that model the computational length unit scaled to $``$270$`\mathrm{au}`$ with the initial disc radius ($`2`$ units) and distance of pericentre ($`2.6`$ units) for the perturber being $``$540$`\mathrm{au}`$ and $``$700$`\mathrm{au}`$ respectively. Since $`G`$ and the primary mass were set to unity in the computations, the computational time unit scaled to $``$530$`\mathrm{yr}`$ (i.e. the period at unit radius divided by 2$`\pi `$). Their system of computational units is used here in similar numerical simulations, however we point out that the scalings mentioned above are only used for reference, in order to maintain generality in the presentation of our results.

The new numerical calculations presented below employ a fifth-order accurate Runge-Kutta-Fehlberg integrator with a variable time-step size. While not as accurate as the Bulirsch-Stoer method (e.g. as used by Hall et al. 1996), this method allows for the fast integration of much larger numbers of particles ($``$104-5) with moderate accuracy. Time-step size control is implemented as described in Press et al. (1992), by specifying a maximum relative error $`ϵ`$ in the velocity magnitude. The value of $`ϵ`$ was taken to be as small as practicable. Unless stated otherwise we take $`ϵ`$10-9. Thus, in a typical run of $``$105 time-steps, the net relative error in the energy and angular momentum of any given particle is $``$10-4. In practice the net relative error for most of the particles will be substantially smaller than this since the time-step size is controlled by a minority of particles that fall close to either of the stellar masses. To prevent too large a reduction in the time-step size owing to this effect we remove particles from the calculation if they approach either star to within a radius $`R_\mathrm{p}`$ of the primary or a radius $`R_\mathrm{s}`$ of the secondary. Unless stated otherwise we take: $`R_\mathrm{p}=0.03`$ and $`R_\mathrm{s}=0.01`$. The former value was chosen as the minimum required to ensure that the time-step size would not be controlled by close approaches of particles to the primary rather than by close approaches to the secondary, for all model parameters considered here. The secondary was integrated as if it were a test particle in motion about a fixed primary of mass equal to the total binary mass. The gravitational forces on particles are then prescribed as a function of position at each time-level with the addition of a fictitious force to compensate for the non-inertial origin (see below). As mentioned above, the test particles are collisionless and their distribution is taken as a tracer of the observed dust presumed to be produced by infrequent collisions in the real system.

The origin of time is taken as the instant of pericentre passage, hence computational time-levels may be negative or positive. Throughout we shall use $`10^4`$ test particles initialised in Keplerian circular coplanar (zero inclination) orbits, placed at random in the radial interval $`0.5`$$`2`$ such that the surface density profile decreases according to a power law of radius to the three-halves. The perturber has inital position and velocity components appropriate to the chosen orbital inclination, eccentricity, and mass. The perturber is initialised at a distance of ten times the initial outer disc radius (i.e. at $`20`$ units). We shall consider a single pericentre distance of $`q=2.6`$ units, thus the initial disc size is $``$0.8$`q`$. The length and mass units, and therefore the time unit, used in our calculations can be scaled so as to apply to any system within the range of dimensionless parameters that we cover. Thus, given our choice of relative disc size, fixing $`q`$ does not affect the generality of our models. Also, since the initial disc is axisymmetric, the longitude of periastron for the perturber’s orbit is arbitrary. In all our particle position plots the perturber orbit is orientated as in Fig. $`3`$, except for differences in orbital inclination.

### 3.1 Circumstellar ring formation

In order to identify the dynamical origin of the ring features simulated by KLSS we consider the system response to different components of the tidal potential (an approach also used recently by Pfalzner & Kley 2000, for gas discs). The total potential is given by $`\mathrm{\Psi }=\mathrm{\Psi }_\mathrm{p}+\mathrm{\Psi }_\mathrm{s}`$, where the point mass potential at position $`𝐫`$ due to a primary of mass $`M_\mathrm{p}`$, to which we have fixed the origin, is:

$$\mathrm{\Psi }_\mathrm{p}=\frac{GM_\mathrm{p}}{𝐫}.$$
(1)
The tidal potential, due to the secondary with mass $`M_\mathrm{s}`$ and instantaneous position $`𝐫_\mathrm{s}`$, is given by:

$$\mathrm{\Psi }_\mathrm{s}=\frac{GM_\mathrm{s}}{𝐫𝐫_\mathrm{s}}+\frac{GM_\mathrm{s}𝐫𝐫_\mathrm{s}}{|𝐫_\mathrm{s}|^3}.$$
(2)
The first term on the right-hand side of (2) is the direct potential due to the secondary star; the second term is the indirect potential due to the acceleration of the origin. In a stellar fly-by this acceleration corresponds to the reflex motion of the primary star due to a gravitational tug from the secondary star.

In Fig. $`4`$ we give particle position plots for fly-by simulations in which: we use only the direct term in the tidal potential, use only the indirect term in the tidal potential, use the full tidal potential. In the case that employs only the direct tidal term the disc is pulled towards the perturber as it approaches pericentre. In the case using only the indirect tidal term the disc is repelled by comparison as the perturber approaches pericentre (in the inertial frame this effect corresponds to acceleration of the primary towards the secondary). In these first two cases eccentric nested rings are observed to form, although they appear as stronger features in the case where only the indirect tidal term is used. Since the effectiveness of the direct tidal interaction depends on the disc being able to closely approach the perturber, the role of the indirect tidal term when the full tidal potential is used is towards weakening the ring features that result from the direct tidal interaction, and to cause approximate pericentric anti-alignment of the two sets of rings. This scheme in the alignment of the particle orbits gives a length asymmetry between opposing sides of the disc when viewed edge-on, with the strongest ring features along the longest extension.

In the above the features that appear to be ring-like are actually formed from distended tightly-wound one-armed kinematic spiral patterns. Different levels of advancement in the winding process can be seen at the third and fourth time-levels depicted in Fig. $`4`$. This is more clearly illustrated in Fig. $`5`$ where we plot particle positions according to the *true anomaly* (i.e. the angle between the instantaneous position vector and the line of pericentre) for the last frame given in Fig. $`4`$ and for a time much later after pericentre passage, for the model using only the indirect term in the tidal potential. \[Although we have now identified the origin of the ring features, since over much of the azimuth the appearence of the features resembles nested eccentric rings, we shall continue with our original usage of that description.\]

Winding of the kinematic spiral occurs as particles move around their orbits with local orbital frequency $`\mathrm{\Omega }`$. Thus the time-scale for the generation of a ring after pericentre passage is $`\mathrm{\Omega }^1`$. The lack of significant width to the spiral in Fig. $`5`$ as compared to Fig. $`4`$, indicates that the width of the rings is mainly due to the spread in the particle eccentricities. The rings will appear to propagate outwards and dissolve due to the growth of incoherence in the phase relationships of particles. Consider two neighbouring particles at radii or semimajor axes $`a_\mathrm{o}`$ and $`a_\mathrm{i}`$, where $`a_\mathrm{o}>a_\mathrm{i}`$. The decoherence time-scale $`t_\mathrm{d}`$ for the particle pair is given by:

$$t_\mathrm{d}^1=\frac{\mathrm{\Omega }(a_\mathrm{i})}{2\pi }\left[1\left(\frac{a_\mathrm{i}}{a_\mathrm{o}}\right)^{3/2}\right].$$
(3)
In $`\beta `$ Pic, and the simulations, the width of the midplane brightness features is $``$5–10 per cent of their radii (KLSS). For dissolution of the ring features we require the tight winding of a spiral down to the scale of a ring width, thus for a ring at radius $`r`$:

$$t_\mathrm{d}10^5\left(\frac{r}{500\mathrm{au}}\right)^{3/2}\mathrm{yr}.$$
(4)
This agrees with the simulations of KLSS in which it is found that the ring features dissolve in about ten local orbital periods. KLSS also found that viewing the perturbed model disc at approximately $`10^5`$$`\mathrm{yr}`$ after pericentre passage gave the most satisfactory fit to the observed disc, which is consistent with the innermost ring feature being detected at $``$500$`\mathrm{au}`$.

### 3.2 Length asymmetry in coplanar models

As discussed above, a length asymmetry in the disc is readily obtained from a close fly-by encounter. The length asymmetry in $`\beta `$ Pic’s disc is also the most easily measured of its asymmetric features. Below we shall continue to numerically investigate the form and size of the length asymmetry as a function of fly-by parameters.

#### 3.2.1 Parabolic versus hyperbolic trajectories

A secondary star of orbital eccentricity $`e`$ with respect to the primary star has relative velocity magnitude $`V_\mathrm{q}`$ at pericentre distance $`q`$ given by:

$$V_\mathrm{q}^2=(1+e)\frac{G(M_\mathrm{p}+M_\mathrm{s})}{q}.$$
(5)
Hence for the scaling used by KLSS, eccentricities of $`1`$, $`2`$, $`5`$, and $`10`$ correspond to pericentre velocities of $``$2.5, $`3.1`$, $`4.3`$, and $`5.9`$km s-1 respectively. In this Section we shall consider coplanar encounters of uniform pericentre distance $`2.6`$, a perturber mass of $`0.3`$, and these orbital eccentricities. The corresponding trajectories and the initial disc size used in these simulations are shown in Fig. $`6`$.

A useful diagnostic measure of the length asymmetry induced in the disc is given by plotting the final apocentre distance against the longitude of pericentre for each perturbed particle orbit. By *final* we shall mean that the secondary mass has passed through pericentre and has returned to a sufficiently large radius that its affect on the positions and velocities of the particles bound to the primary is negligible. In our simulations we integrated the system for a total time of $`120`$ units, which means that the final perturber distance was always much greater than the initial distance, corresponding to a time of at least $`80`$ units after pericentre passage. This total run-time was also generally required in order for the numbers of particles either bound or unbound to each of the two stellar masses to stabilise (in practice this meant to approach a constant value with less than about $`10`$ per cent fluctuations; cf. Hall et al. 1996). This information is summarised in Table $`1`$, for all the $`e1`$ models presented in this paper.

In Fig. $`7`$ we compare the disc response for encounters of different eccentricity. In all cases there are two oppositely directed extensions. The extension near to a longitude of $`2`$ contains the strong ring features, and the other extension (near longitude $`+1`$) contains the weak ring features, as discussed above. Respectively we shall refer to these as extension A and extension B. Thus extension A forms as a response to the indirect component of the tidal potential, and extension B forms as a response to the direct component of the tidal potential.

As the encounter velocity (equivalently, the orbital eccentricity) increases the peak of extension A diminishes in amplitude, while the peak of extension B increases its amplitude. Further examination of the extension B results illuminates a key difference between the parabolic ($`e=1`$) and hyperbolic ($`e>1`$) orbit cases. In the former, extension B is compact, but in the latter cases all show a much longer extension B. In Fig. $`8`$ we show the final apocentre distance against the longitude of pericentre for the parabolic case in which only the direct part of the tidal potential is used. As we might have expected extension A is not in evidence (being due to the indirect part of the tidal potential), and extension B is similar in nature to that found in the hyperbolic orbit cases. Thus we conclude that in the full tidal potential case the indirect component acts on the disc prior to pericentre passage such that the disc response due to the direct component near pericentre is modified as compared with a calculation that uses only the direct term (also discussed in Section $`3.1`$). Hence the more impulsive (i.e. hyperbolic) encounters yield a response that is similar to that found for the parabolic case that uses the direct term only, being manifest in the form of extension B. It follows that higher velocity (i.e. more eccentric) perturbers will have extension B forms more like that found in Fig. $`8`$, and smaller amplitudes for extension A. This can be seen in Fig. $`7`$.

In Fig. $`9`$ we give the final eccentricities of the disc particles for different eccentricities. In the parabolic orbit case, extension A particles correspond to the wing-shaped distribution having eccentricities approximately in the range $`0.2`$$`0.6`$. Superposed with this are the extension B particles that have generally very high eccentricites. The change in the relative lengths of extensions A and B with increasing perturber eccentricity is also apparant. Generally, the excited eccentricities of the disc particles are lower for larger perturber eccentricities, which can be explained from the simplistic impulse model for tidal interaction (Lin & Papaloizou 1979): higher relative velocity results in weaker coupling between the particles and the perturber.

#### 3.2.2 Variation of the perturber mass

In Fig. $`10`$ we show the final distance of apocentre against the final longitude of pericentre for parabolic encounters with different masses for the perturber. We consider the mass ratio $`\mu =M_\mathrm{s}/M_\mathrm{p}`$. Increasing the perturber mass from $`\mu =0.1`$ to $`\mu =1.0`$, we find that extension A increases in amplitude and extension B diminishes in amplitude. The former effect is due to the larger perturber mass inducing a larger amplitude reflex motion of the primary star; the latter effect evidently results from the more massive perturber being able to strip more particles from the disc in the direct interaction that sets up extension B (see Section 3.4).

The lowest mass model with $`\mu =0.1`$ has an extension B that has a similar extended form to the high eccentricity cases discussed above. This is due to the low mass of the perturber giving only a small reflex motion of the primary, and so a relatively weak interaction prior to pericentre passage, that necessarily results in a more impulsive response. Only the $`\mu =0.3`$ model gives a compact extension B and extension A. This case has in terms of a model for the length asymmetry in the disc of $`\beta `$ Pic the most appropriate amplitude ratio for the extensions A and B of those coplanar models tried here. The ratio of lengths of apocentre for this model is approximately $`9/6=1.5`$, the length ratios of the other models clearly being very different from unity. The length ratio between extensions A and B is given in Table $`1`$, for all our models with $`e1`$. Where the ratio is clearly much different from unity we simply state either $`1`$ or $`1`$, since the interpretation of the length of an extension as defined by a small number of particles at large distances is not clear.

We note that the size of extension A can be reduced by increasing the perturber’s eccentricity, so it may be possible that a more massive ($`\mu >0.3`$) and eccentric orbit ($`e1`$) perturber could fit the observational results. However, since eccentricities larger than approximately unity do not display a compact extension B, we favour the perturber mass originally considered by KLSS (i.e. $`\mu 0.3`$, for $`e1`$).

### 3.3 Inclined orbits

In general, if the perturber’s orbit and the initial disc plane are not coplanar, then the tidal interaction occurring between the disc particles and the perturber can excite vertical motions in the particle orbits. Here we shall consider this effect in terms of the inclination $`i`$ of the perturber’s orbit to the initial disc plane. This angle is generated by rotating the perturber’s orbit about the *latus rectum* (i.e. the line perpendicular to the line of pericentre, that also passes through the origin), which then becomes the *line of nodes* (e.g. see Fig. $`3`$).

KLSS found that (for $`i=30`$°) extension B became vertically flared, with much weaker inclinations in extension A, as can be seen in Fig. $`3`$. From the preceeding discussion we deduce that the flaring results from the stronger direct tidal interaction, and the smaller inclinations generated in extension A result from the weaker indirect tidal interaction. The resulting morphology resembles the butterfly asymmetry (Kalas & Jewitt 1995, KLSS) in which the vertical height of the imaged disc is greatest north of the midplane in the southwest extension (see Figs. $`1`$ and $`2`$), and south of the midplane in the northeast extension, with the overall greatest vertical height in the former.

#### 3.3.1 Length asymmetry

In Fig. $`11`$ we give the final distance of apocentre against the longitude of pericentre for parabolic encounters with $`\mu =0.3`$, and $`i=30`$° and $`i=60`$°. In the lower inclination case the disc response is very similar to the coplanar encounter presented above, except that extension A is shorter by about $`20`$ per cent. In the higher inclination case the response resembles that of a coplanar hyperbolic encounter, having a very short extension A and a very long extension B. As before, this is consistent with the notion that the disc interaction prior to pericentre passage is a key difference between hyperbolic and parabolic encounters. In the $`i=60`$° case the high inclination takes the perturber orbit further from the disc, making the pericentre interaction more impulsive. This shows that we require an orbit with not too high an inclination, if we wish to preserve the compact nature of extension B. Additionally, the other angle possible, which corresponds to rotation of the orbit about the line of pericentre, need not be considered since for low values of $`i`$ that angle takes the orbit yet further from the disc. In other words, we only consider an argument of pericentre of $`90`$°.

In the higher inclination case the length ratio of extension A to extension B is clearly much less than unity. In the lower inclination case extension B is compact, but the peak is less distinct than for the coplanar case. We measure its height by taking the intersection of the tangents to the sides of the envelope of the extension particles near its base. This yields a height of approximately $`6`$, which is similar to the coplanar model. Thus we deduce a ratio $`7/61.2`$. This is the best fitting value, in terms of the length asymmetry, of those models presented in this paper. Scaling according to the KLSS model the lengths of the extensions are $`1890`$$`\mathrm{au}`$ and $`1620`$$`\mathrm{au}`$, which is in good agreement with the new determination of the size of the disc presented in Section $`2`$.

The exact numerical value of the length ratio determined for the $`\beta `$ Pic disc could be recovered in our models by taking a slightly smaller orbital inclination than $`30`$°. That would not be productive at this stage, however, since here we are not concerned with the details of translating our results into a scattered light disc model, which could introduce a small correction factor. However, for reference, a model with $`i=15`$°was tried and was found to give a length ratio of $`1.4`$. Thus we expect the appropriate inclination to lie between $`15`$° and $`30`$°, i.e. closer to the lower inclination model above than to coplanar.

#### 3.3.2 Flaring and warps

In Fig. $`12`$ we show the final eccentricity and inclination distributions for the disc particles in the two cases mentioned above. The amplitude of the response is lower for higher values of $`i`$, for the reasons indicated above. The eccentricity response resembles that of a coplanar hyperbolic orbit, as expected from the preceeding discussion. The inclinations associated with extension A are much smaller than for those associated with extension B. However, the extension A inclinations are higher in the $`i=60`$° model, owing to input from the direct interaction because of increased impulsiveness of the interaction occuring at the ascending node in that case. Thus, by symmetry of the orbit with respect to the initial disc, we expect the inclinations in each extension to become comparable when orthogonality is approached by increasing $`i`$. This provides a further need for a low inclination encounter as the disc of $`\beta `$ Pic is strongly lopsided with respect to vertical flaring (Kalas & Jewitt 1995).

For these two cases we also give the final pericentre distance distributions in Fig. $`13`$. In the low inclination case we see that particles from extension B can plunge to pericentre distances of $`0.2`$. KLSS proposed that these particles could explain the existence of an inner ’warp’ in the scattered light dust disc of $`\beta `$ Pic at $``$50$`\mathrm{au}`$ (Burrows et al. 1995, Mouillet et al. 1997). We can see from Fig. $`13`$, that for the higher inclination encounter, extension B particles have significantly higher perturbed pericentre distances (down to about $`0.7`$; or $``$200$`\mathrm{au}`$). This is a direct consequence of the lower eccentricities produced at comparable semimajor axes in extension B for the two cases. Thus the minimum pericentre distance is minimised in coplanar cases. Correspondingly, larger minimum pericentre distances are also charateristic of hyperbolic encounters. Fig. $`14`$ gives the final eccentricity and pericentre distributions for a hyperbolic encounter with $`e=2`$ and $`i=30`$°. If the explanation of the inner warp suggested by KLSS is confirmed then this also supports the need for a low inclination and relative velocity encounter.

### 3.4 Planetesimal capture

As can be seen in Fig. $`4`$ the close fly-by of a disc of particles by a perturbing star can result in the very close approach of some of the disc particles to the perturber. Since the unperturbed disc does not extend all the way out to pericentre, this occurs by particles streaming through the vicinity of the instantaneous Lagrange Point between the two stars, about the time of closest approach. In the case that particles become captured by the perturber the capture efficiency is generally low (see also Hall et al. 1996). Table $`1`$ gives the percentage of initial disc particles that are captured by the perturber for the models discussed above, and from it we can see that the lower eccentricity and inclination encounters are required, in order to effect capture at all.

For the model parameters used by KLSS, the percentage of captured particles is $``$10 per cent of the initial particle number, but note that if the initial disc were allowed to extend to pericentre this figure could be fractionally larger (Hall et al. 1996). In the case of $`\beta `$ Pic, considering the maximal parent body disc of Backman & Paresce (1993), the captured material could therefore represent up to about $`10^3`$ earth masses of planetesimals which would subsequently orbit the perturber and produce a scattered light signature similar to Vega-like sources. This would not however uniquely identify a perturber candidate since approximately 15 per cent of main sequence stars of all types are known to have such discs (Backman & Gillett 1987, Aumann 1988), but given the encounter parameters we find here, it seems likely that a perturber candidate should have this signature.

None the less we may find crude approximations for the orbital properties of the captured planetesimals. We shall make an analogy with the theory of mass transfer in semi-detached binaries, in which the captured matter orbits the secondary with the specific angular momentum of relative motion at the Lagrange Point (see Frank, King & Raine 1992). Here we deduce that if a particle in the outermost parts of the disc is captured close to the perturber’s pericentre, then it will orbit the perturber with specific angular momentum:

$$\left[\sqrt{\frac{G(1+e)(M_\mathrm{p}+M_\mathrm{s})}{q}}\sqrt{\frac{GM_\mathrm{p}}{R_\mathrm{o}}}\right]b,$$
(6)
where $`R_\mathrm{o}`$ is the initial outer radius of the disc and $`b`$ is the separation of the secondary from the Lagrange Point at the instant of pericentre passage, approximated by:

$$\frac{b}{q}=\frac{1}{2}+0.227\mathrm{log}(\mu ).$$
(7)
In the above we have taken the circular velocity of a particle lying on the binary axis relative to the perturber at pericentre for the input particle velocity. Although the instantaneous Lagrange Point is located well within the envelope of the initial disc, the rapid eccentricity changes in the disc prior to particle transfer cause compression of the disc such that the outer disc particles pass through the vicinity of the Lagrange Point with approximately their initial velocities (e.g. see Fig. $`4`$).

With respect to the perturber, the specific angular momentum of the captured particle is: $`\sqrt{GM_\mathrm{s}a^{}(1e^2)}`$, where the characteristic semimajor axis and eccentricity with respect to the secondary are respectively $`a^{}`$ and $`e^{}`$. Hence we have:

$`a^{}(1e^2)`$ $`=`$ $`q(1+e){\displaystyle \frac{(1+\mu )}{\mu }}\left({\displaystyle \frac{b}{q}}\right)^2`$ (8)
$`\times \left[1\sqrt{{\displaystyle \frac{q/R_\mathrm{o}}{(1+e)(1+\mu )}}}\right]^2.`$
Thus for a model with $`e=1`$, $`\mu =0.3`$, and $`q=2.6`$, we find $`a^{}(1e^2)0.3`$. This equation does not apply when the relative velocity of the perturber and disc particles is sufficiently high that particle transfer is energetically disallowed. As shown in Table $`1`$ we find that particles are captured when $`e2`$ but not when $`e5`$.

The above is in fair agreement with the simulations. In Fig. $`15`$ we give the positions, eccentricities, and semimajor axes of captured particles relative to the perturber for the above model. We find the bound particles have a median value $`a^{}0.70`$, from which we calculate $`e^{}0.76`$. The actual median value is $`e^{}0.85`$ for the particles, which is $`10`$ per cent larger than predicted. This is likely due to the continued interaction of the primary and the captured particles as the perturber receeds from pericentre. The particle transfer does not occur in a well defined stream which results in a large spread in semimajor axes and eccentricities. However, the captured particles appear to manifest as an asymmetrical disc about the perturber (upper panel of Fig. $`15`$), and show clustering in $`e^{}`$ and $`a^{}`$ (lower panel of Fig. $`15`$), indicating some level of coherence in the transfer (seen in Fig. $`4`$). Scaling according to KLSS, we find that $`90`$% of the disc particles are located within $`500`$$`\mathrm{au}`$ of the perturber. For the KLSS model (i.e. the previous model augmented with an orbital inclination of $`30`$° for the perturber), the median eccentricity and semimajor axis are each found to be $`0.7`$, consistent with the higher relative disc velocity expected for the inclined encounter.

The luminosity of an M0V dwarf star (of mass $`0.5`$ M, corresponding to $`\mu =0.3`$) is $``$102 times smaller than for $`\beta `$ Pic. If the amount of scattering surface in the disc is proportional to the square of the number of parent bodies (e.g. Backman & Paresce 1993), then as an upper limit the brightness of the perturber’s disc is expected to be $``$104 times fainter than that around $`\beta `$ Pic, and $``$102 times fainter relative to the central star. Thus it is plausible that the perturber may demonstrate a measurable scattered light disc as a result of collisional erosion amongst the captured planetesimals at the distance to $`\beta `$ Pic.

### 3.5 Close encounters with a bound perturber

Since the innermost midplane brightness feature identified by KLSS is measured to be at $``$500$`\mathrm{au}`$ the age of the rings in the $`\beta `$ Pic disc should be $``$105$`\mathrm{yr}`$, as calculated above. Although this figure could be revised up or down according to future observations, we shall take this to mean that the epoch of the fly-by encounter was about $`10^5`$$`\mathrm{yr}`$ in the past. Use of the *Hipparcos* survey data reveals that this is also the minimum time-scale we must look back from the present in order to find a perturber fly-by at $`10^5`$$`\mathrm{au}`$ (Kalas et al. 2000). The eighteen most closely approaching stars identified by Kalas et al. (2000) have a mean encounter velocity $``$40$`\mathrm{km}\mathrm{s}^1`$, which is considerably higher than that deduced from our dynamical calculations. Thus we are prompted to consider companions that are or were physically associated with $`\beta `$ Pic.

Repeated close encounters of the bound star with the disc may destroy the asymmetrical disc structure and ring features set up by the first encounter. Since for $`e<1`$ the semimajor axis of the perturber can be written $`q/(1e)`$, by demanding the orbital period to be greater than some value $`P_{\mathrm{min}}`$, we deduce:

$$e>1q\left[\frac{4\pi ^2}{G(M_\mathrm{p}+M_\mathrm{s})P_{\mathrm{min}}^2}\right]^{1/3}.$$
(9)
Thus if we take $`P_{\mathrm{min}}=10^5`$$`\mathrm{yr}`$ we infer that the perturbed binary eccentricity needs to be greater than $`0.75`$, in order that we may view the ring features at their inferred age.

In Fig. $`6`$ we show a trajectory with $`e=0.75`$. The similarity of such high eccentricty bound orbits to the parabolic trajectory suggests that the disc response should be similar. In Fig. $`16`$ we give the apocentre distance versus longitude of pericentre plots for a coplanar encounter simulation with $`e=0.8`$, $`\mu =0.3`$, and $`q=2.6`$. The perturber starts at apocentre and we give results for the first and second completed orbital periods. Clearly the length asymmetry is completely disrupted after the second encounter. In Fig. $`17`$ we give the corresponding particle positions, which show that the ring structure is also destroyed after the second encounter. Thus as suspected only a single fly-by can be allowed.

Possible scenarios in which a single pericentre passage could have occurred in the recent past are: 1) the perturber was originally a bound companion to $`\beta `$ Pic with a sufficiently large pericentre distance that no significant interaction with the disc occurred over its lifetime, and sufficently large semimajor axis and eccentricity that it could be perturbed into the required trajectory by the close fly-by of a massive field star $``$105$`\mathrm{yr}`$ in the past; 2) $`\beta `$ Pic was a hierarchical multiple system that was dynamically unstable on a time-scale comparable to its present age, resulting in the ejection of one or more stars and the development of a disc-intercepting orbit for another, or an ejected star passed by $`\beta `$ Pic on its way out of the system.

In the first scenario the three-body interaction should have occurred at least $`10^5`$$`\mathrm{yr}`$ ago. The further in the past we put this event, the more eccentric the resultant binary orbit must be in order to delay closest approach until $`10^5`$$`\mathrm{yr}`$ ago (if the perturbed star has an initially out-going trajectory). Although perhaps more attractive than the lower probability close binary encounter model used by KLSS, the problem of low relative velocity is shifted from the star-disc interaction to the interaction between the second and third bodies. In respect of the second scenario we point out that Weinberger et al. (2000) report an age of $`5`$$`\mathrm{Myr}`$ for the triple system HD 141569. This consists of a Vega-like BVe primary with a $`500`$$`\mathrm{au}`$ disc and two M stars each at $``$1000$`\mathrm{au}`$ projected separation. Weinberger et al. (2000) deduce that the relative separations involved may make the system unstable, giving the possibility that HD 141569 represents an earlier stage in the evolution of the $`\beta `$ Pic system (with an age of $`20`$$`\mathrm{Myr}`$). Studies of the dynamics involved suggest that young stars with multiple companions may not lose some of them via dynamical instability until main sequence ages (Eggleton & Kiseleva 1995), which is supported by observations of high binary frequency amongst young stars compared to field stars (Mathieu 1994).

## 4 Discussion

Here we have considered the stellar fly-by hypothesis that may account for the asymmetrical structure in the scattered light dust disc of $`\beta `$ Pic. We investigated the dynamics of the perturbed planetesimal disc in an attempt to learn more about the possible parameters of a fly-by encounter. Additionally, we have described the origin of transient circumstellar eccentric ring structures as being a general outcome of an encounter, due to the reflex motion of the primary star as the perturber passes through closest approach. The apparent ring system is actually an eccentric tightly-wound one-armed spiral pattern that gradually disappears owing to phase-mixing in the particles’ orbital motion. In inclined-orbit simulations we find that the vertical disc response also has analogues in the scattered light disc of $`\beta `$ Pic. We summarise the coplanar results in Fig. $`18`$, and the inclined-orbit results in Fig. $`19`$.

For the perturber we conclude that the mass should not be very different from $`0.5`$ M (corresponding to spectral type M0V). In this case we favour a low inclination prograde encounter with relative velocity of a few km s-1, which makes a bound perturber seem the most natural choice. One problem with this is that asymmetries would be destroyed in the second pericentre passage. Thus, given the deduced youth of the asymmetries, the perturber in this scenario is required to make a single close approach in the recent past (i.e. $``$105$`\mathrm{yr}`$ ago). This can be accounted for by making the perturber-$`\beta `$ Pic system an initial wide binary that is perturbed by a massive passing field star $``$105$`\mathrm{yr}`$ in the past, such that the perturber is sent into a close approach orbit about $`\beta `$ Pic. A more probable solution is to arrange for $`\beta `$ Pic to have originally been a hierachical multiple system (cf. HD 141569; Weinberger et al. 2000) that became dynamically unstable on the age of the star, and with the result being that one component was sent into a disc-intercepting orbit. In any case the perturber should have captured planetesimals from the primary’s disc and may exhibit a detectable scattered light disc as a result. A survey of a more extended field around $`\beta `$ Pic than has been considered previously is in progress in order to test our hypotheses.

The dynamical response of the models suggests that the $`\beta `$ Pic disc might be separable into different components corresponding to groupings of perturbed particle orbits. For example, we find that outwards of $``$250$`\mathrm{au}`$ the SW extension may not have a distinct midplane owing to pumping of particle inclinations. The positions close to pericentre of the moderately eccentric NE midplane particles could account for the measurement of a SW midplane out to $``$650$`\mathrm{au}`$, in superposition with the more diffuse vertically and radially extended component. When the optical image is vertically compressed (Fig. $`2`$) we are able to measure the SW extension to much larger radii ($`1450`$$`\mathrm{au}`$). In our model the corresponding (very eccentric and inclined) particle orbits also intersect the disc near to their pericentres (down to $``$50$`\mathrm{au}`$) possibly accounting for the inferred inner warp of the disc, which is aligned with the outer flared envelope in the SW. In summary, there are three recognizable particle groupings that can be related to the morphology of the $`\beta `$ Pic disc. These are: the highly eccentric and inclined particles that reach apocentre in the SW (extension B particles), the moderately eccentric and weakly inclined particles that reach apocentre in the NE (extension A particles), and the relatively unperturbed particles inside $``$200$`\mathrm{au}`$ radius (Figs. $`12`$ and $`13`$).

## Acknowledgments

We are grateful to the University of Hawaii for supporting the observations. This work was funded in part by a grant from NASA held by D. Jewitt. We also thank J. Gradie, B. Zuckerman, and E. Becklin for the use of their coronagraph. The clarity of this paper was improved by comments from S. Ida. P. Kalas acknowledges financial support and hospitality given by the Astronomy Unit, Queen Mary & Westfield College, during his visit there.