text:
# Traceroute sampling makes random graphs appear to have power law degree distributions

## I Introduction

The Internet and the networks it facilitates — including the Web and email networks — are the largest artificial complex networks in existence, and understanding their structural and dynamic properties is important if we wish to understand social and technological networks in general. Moreover, efforts to design novel dynamic protocols for communication and fault tolerance are well served by knowing these properties.

One structural property of particular interest is the degree distribution at the router level of the Internet. This distribution has been inferred Govindan ; Faloutsos ; Rocketfuel ; CAIDA ; Oregon ; LookingGlass ; Opte  both by sampling traceroutes, i.e., the paths chosen by Internet routers, which approximate shortest paths in the network, and by taking “snapshots” of BGP (Border Gateway Protocol) routing tables BGP . These methods have been criticized as being noisy and imperfect Amini ; Chen . However, Lakhina et al. Lakhina  recently argued that such methods have a more fundamental flaw. Due to the fact that these methods use a small number of sources for the inference, they were able to show that only a small fraction of the edges of a graph are “visible” in such a sample. Moreover, the set of visible edges is biased in such a way that Erdős-Rényi random graphs gnp , whose underlying distribution is Poisson, will appear to have a power law degree distribution. While no one would argue that the Internet is a purely random graph, this certainly calls into question the standard measurements of power law or “scale-free” degree distributions on the Internet, and reopens the problem of characterizing Internet topology.

In this paper we explain this bias phenomenon analytically. Specifically, by modeling the growth of a breadth-first spanning tree with differential equations, we show that sampling shortest paths from a single source in an Erdős-Rényi random graph gives rise to a power law degree distribution of the form $`P(k)1/k`$, up to a cutoff $`kc`$ where $`c`$ is the average degree of the underlying graph. While sampling traceroutes from a single source is rather limited, Barford Barford  provides empirical evidence that, on the Internet, merging shortest paths from several sources leads to only marginally improved surveys of Internet topology. In the Conclusions we discuss generalizing our approach to sampling from multiple sources.

Finally, even if the Internet has a power-law degree distribution, the exponent may be rather different from the one observed in traceroute samples. In future work we plan to extend our approach to graphs with arbitrary degree distributions, to study the relationship between the observed exponent and the underlying one.

## II Internet Spanning Trees

Most mapping projects have inferred the Internet’s degree distribution by implicitly building a map of the network from the union of a large number of traceroutes from a single source, or from a small number of sources. Assuming that traceroutes are shortest paths, sampling from a single source is equivalent to building a spanning tree. We therefore model this sampling method by modeling the growth of a spanning tree on a graph.

There are several ways one might build a spanning tree. We will consider a family of methods, in which at each step, every vertex in the graph is labeled reached, pending, or unknown. The pending vertices are the leaves of the current tree, the reached vertices are those in its interior, and the unknown vertices are those not yet connected. To model traceroutes from a single source, we initialize the process by labeling the source vertex pending, and all other vertices unknown.

Then the growth of the spanning tree is given by the following pseudocode.

> | while | there are pending vertices: |
> | --- | --- |
> |  | choose a pending vertex $`v`$ |
> |  | label $`v`$ reached |
> |  | for | every unknown neighbor $`u`$ of $`v`$, |
> |  |  | label $`u`$ pending. |

The type of spanning tree is determined by how we choose which pending vertex $`v`$ we will use to extend the tree. To model shortest paths, we store the pending vertices in a queue, take them in FIFO (first-in, first-out) order, and build a breadth-first tree; if we like we can break ties randomly between vertices of the same age in the queue, which is equivalent to adding a small noise term to the length of each edge as in Lakhina . If we store the pending vertices on a stack and take them in LIFO (last-in, first-out) order, we build a depth-first tree. Finally, we can choose from among the pending vertices uniformly at random, giving a “random-first” tree.

Surprisingly, while these three processes build different trees, and traverse them in different orders, we will see in the next section that they all yield the same degree distribution when $`n`$ is large. To illustrate this, Fig. 1 shows empirical degree distributions of breadth-first, depth-first, and random-first spanning trees for a random graph $`G(n,p=c/n)`$ where $`n=10^5`$ and $`c=100`$. The three degree distributions are indistinguishable; further, they are well-matched by the analytic results given in the next section, and obey a power law with exponent $`1`$ for degrees less than $`c`$. For comparison, we also show the Poisson degree distribution of the underlying graph.

## III Analytic Results

In this section we show analytically that building spanning trees in Erdős-Rényi random graphs $`G(n,p=c/n)`$, using any of the processes described above, gives rise to an apparent power law degree distribution $`P(k)1/k`$ for $`kc`$. We focus here on the case where the average degree $`c`$ is large, but constant with respect to $`n`$; we believe our results also hold if $`c`$ is a moderately growing function of $`n`$, such as $`\mathrm{log}n`$ or $`n^ϵ`$ for small $`ϵ`$, but it seems more difficult to make our analysis rigorous in that case.

We will model the progress of the while loop described above. Let $`S(T)`$ and $`U(T)`$ denote the number of pending and unknown vertices at step $`T`$ respectively. The expected changes in these variables at each step are (where $`\mathrm{E}[]`$ denotes the expectation)

$`\mathrm{E}[U(T+1)U(T)]`$ $`=`$ $`pU(T)`$
$`\mathrm{E}[S(T+1)S(T)]`$ $`=`$ $`pU(T)1`$ (1)
Here the $`pU(T)`$ terms come from the fact that a given unknown vertex $`u`$ is connected to the chosen pending vertex $`v`$ with probability $`p`$, in which case we change its label from unknown to pending; the $`1`$ term comes from the fact that we also change $`v`$’s label from pending to reached. Moreover, these equations apply no matter how we choose $`v`$; whether $`v`$ is the “oldest” vertex (breadth-first), the “youngest” one (depth-first), or a random one (random-first). By the principle of deferred decisions, the events that $`v`$ is connected to each unknown vertex $`u`$ are independent and occur with probability $`p`$. Our experiments do indeed show that these three processes result in the same degree distribution.

Writing $`t=T/n`$, $`s(t)=S(tn)/n`$ and $`u(tn)=U(t)/n`$, the difference equations (1) become the following system of differential equations,

$`{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}`$ $`=`$ $`cu`$
$`{\displaystyle \frac{\mathrm{d}s}{\mathrm{d}t}}`$ $`=`$ $`cu1`$ (2)
With the initial conditions $`u(0)=1`$ and $`s(0)=0`$, the solution to (2) is

$$u(t)=\mathrm{e}^{ct},s(t)=1t\mathrm{e}^{ct}.$$
(3)
The algorithm ends at the smallest positive root $`t_f`$ of $`s(t)=0`$; using Lambert’s function $`W`$, defined as $`W(x)=y`$ where $`y\mathrm{e}^y=x`$, we can write

$$t_f=1+\frac{1}{c}W(c\mathrm{e}^c).$$
(4)
Note that $`t_f`$ is the fraction of vertices which are reached at the end of the process, and this is simply the size of the giant component of $`G(n,c/n)`$.

Now, we wish to calculate the degree distribution $`P(k)`$ of this tree. The degree of each vertex $`v`$ is the number of its previously unknown neighbors, plus one for the edge by which it became attached (except for the root). Now, if $`v`$ is chosen at time $`t`$, in the limit $`n\mathrm{}`$ the probability it has $`k`$ unknown neighbors is given by the Poisson distribution with mean $`m=cu(t)`$,

$$\mathrm{Poisson}(m,k)=\frac{\mathrm{e}^mm^k}{k!}.$$

Averaging over all the vertices in the tree gives

$$P(k+1)=\frac{1}{t_f}_0^{t_f}dt\mathrm{Poisson}(cu(t),k).$$

It is helpful to change the variable of integration to $`m`$. Since $`m=c\mathrm{e}^{ct}`$ we have $`\mathrm{d}m=cm\mathrm{d}t`$, and

$`P(k+1)`$ $`=`$ $`{\displaystyle \frac{1}{t_f}}{\displaystyle _{c(1t_f)}^c}dm{\displaystyle \frac{\mathrm{Poisson}(m,k)}{cm}}`$ (5)
$``$ $`{\displaystyle _{c\mathrm{e}^c}^c}dm{\displaystyle \frac{\mathrm{Poisson}(m,k)}{cm}}`$
$`=`$ $`{\displaystyle \frac{1}{ck!}}{\displaystyle _{c\mathrm{e}^c}^c}dm\mathrm{e}^mm^{k1}.`$
Here in the second line we use the fact that $`t_f1\mathrm{e}^c`$ when $`c`$ is large (i.e., the giant component encompasses almost all of the graph).

The integral in (5) is given by the difference between two incomplete Gamma functions. However, since the integrand is peaked at $`m=k1`$ and falls off exponentially for larger $`m`$, for $`kc`$ it coincides almost exactly with the full Gamma function $`\mathrm{\Gamma }(k)`$. Specifically, for any $`c>0`$ we have

$$_0^{c\mathrm{e}^c}dm\mathrm{e}^mm^{k1}<c\mathrm{e}^c$$

and, if $`k1=c(1ϵ)`$ for $`ϵ>0`$, then

$`{\displaystyle _c^{\mathrm{}}}dm\mathrm{e}^mm^{k1}`$ $`=`$ $`\mathrm{e}^cc^{k1}{\displaystyle _0^{\mathrm{}}}dx\mathrm{e}^x(1+x/c)^{k1}`$
$`<`$ $`\mathrm{e}^cc^{k1}{\displaystyle _0^{\mathrm{}}}dx\mathrm{e}^x\mathrm{e}^{x(k1)/c}`$
$`=`$ $`{\displaystyle \frac{\mathrm{e}^cc^{k1}}{ϵ}}<{\displaystyle \frac{\mathrm{e}^{(k1)}(k1)^{k1}}{ϵ}}`$
$``$ $`{\displaystyle \frac{\mathrm{\Gamma }(k)}{ϵ\sqrt{2\pi (k1)}}}`$
This is $`o(\mathrm{\Gamma }(k))`$ if $`ϵ1/\sqrt{k}`$, i.e., if $`k<cc^\alpha `$ for some $`\alpha >1/2`$. In that case we have

$$P(k+1)=(1o(1))\frac{\mathrm{\Gamma }(k)}{ck!}\frac{1}{ck}$$
(6)
giving a power law of exponent $`1`$ up to $`kc`$.

Although we omit some technical details, this derivation can be made mathematically rigorous using results of Wormald Wormald , who showed that under fairly generic conditions, the state of discrete stochastic processes like this one is well-modeled by the corresponding rescaled differential equations. Specifically it can be shown that if we condition on the initial source vertex being in the giant component, then with high probability, for all $`t`$ such that $`0<t<t_f`$, $`U(tn)=u(t)n+o(n)`$ and $`S(tn)=s(t)n+o(n)`$. It follows that with high probability our calculations give the correct degree distribution of the spanning tree within $`o(1)`$.

## IV Discussion and Conclusions

Lakhina et al. Lakhina , argued that sampling traceroutes from a small number of sources profoundly underestimates the degrees of vertices far from the source, and that this effect can cause graphs to appear to have a power law degree distribution even when their underlying distribution is Poisson. In this work, we modeled this sampling process as the construction of a spanning tree, and showed analytically that for sparse random graphs of large average degree, the apparent degree distribution does indeed obey a power law of the form $`P(k)k^1`$ for $`k`$ below the average degree. This illustrates the danger of concluding the existence of a power-law from data over too small a range of degrees, and, more specifically, the danger of sampling traceroutes from just a few sources. While our analytic results hold for traceroutes from a single source, we conjecture that if we use any constant number of sources and take the union of the resultant spanning trees, the observed distribution will remain subject to the same sampling bias, and one will again observe a power-law degree distribution.

Certainly this exponent $`1`$ is not the one observed for the Internet. For instance, Faloutsos, Faloutsos and Faloutsos Faloutsos  observed degree distributions $`P(k)k^{2.5}`$ at the router level, and $`k^{2.2}`$ for out-degrees at the inter-domain (BGP) level. While the real degree distribution of the Internet may be a power law, one possibility is that the true exponent is rather different from the observed one; this possibility was recently explored by Peterman and de los Rios paolo  (they also study another mechanism for observing apparent power laws in the degree distribution, which samples the original graph via a probabilistic pruning strategy and gives an observed exponent of $`2`$ for random graphs). In future work we will extend our differential equation model to random graphs with power law degree distributions to explore this possibility.

Although our analytic model of spanning trees illustrates one possible mechanism for apparent power-law degree distributions, another possibility is that decisions (such as those made by the border-gateway protocol) made by actual Internet routers interact in complex ways with the the underlying link-level topology. It may well be that routers only use a small fraction of the edges in the network, and that the edges they actually use give rise to an effective degree distribution with roughly a power-law form. This possibility implies that even if the real degree distribution of the Internet is something very different from a power law, it may not matter as these extra links are not utilized in normal routing decisions.

Another difference of unknown, but perhaps small, significance is that on the Internet, routing algorithms are intended to optimize for end-to-end traffic flow, and routers in the “interior” of the Internet such as those on the backbone are not themselves destinations or sources. Unlike the empirical studies cited earlier, we selected sources uniformly at random for our spanning trees.

The story of the topology of the Internet is far from over and will undoubtedly remain a topic of great interest for many years to come. However, knowing decisively that current methods are fundamentally biased will serve to push the state-of-the-art forward toward more robust methods of characterization. We see the following three lines of inquiry as enticing:

First, having analytically shown that a Poisson degree distributions can lead to apparent power laws, we naturally wish to generalize our approach to random graphs with arbitrary degree distributions, characterize the family of distributions which generate apparent power laws, and understand the relationship between observed and underlying exponents for power law distributions. Having this generalization may allow us to make firm and useful claims about the topology of the Internet.

Secondly, in what manner is the traceroute-sampled degree distribution of a random graph dependent on the length of the paths taken between a source and destination? If we simulate a “crawl” by using random walks rather than short paths, how would the observed degree distribution change?

Finally, if we use build spanning trees from $`m`$ sources and take their unions, how does the observed degree distribution vary with $`m`$? If we are correct that any constant $`m`$ leads to the same kind of bias, how does $`m`$ need to grow with the network to obtain an accurate sample?

###### Acknowledgements.

The authors are grateful to David Kempe, Mark Newman, Michel Morvan, Paolo De Los Rios, Mark Crovella and Tracy Conrad for helpful conversations, and to American Airlines for providing a conducive working environment. This work was funded by NSF grant PHY-0200909.


format: text

text:
# Driving spiral arms in the circumstellar disks of HD 100546 and HD 141569A

## 1 Introduction

Recent high angular resolution studies of nearby young stars systems have revealed structure in the disks surrounding these stars. For example, gaps and holes in disks have been observed via direct or coronagraph imaging in systems such as HR 4796A and HD 141569A (Jayawardhana et al., 1998; Augereau et al., 1999a, b; Weinberger et al., 1999), and intricate spiral structure has been revealed in the circumstellar dusty disks surrounding HD 100546 and HD 141569A (Grady et al., 2001; Clampin et al., 2003). These two disks are very large (outer radius a few hundred AU) and the spiral structure is suspected to be caused by a perturber external to the system, such as a recent stellar encounter in the case of HD 100546 or its companion binary, HD 141569B and C, in the case of HD 141569A (Grady et al., 2001; Clampin et al., 2003; Augereau & Papaloizou, 2004).

While HD 100546 is relatively isolated, most stars are born in embedded stellar clusters (Lada & Lada, 2003). HD 141569A is in a low density region near a molecular cloud and has 3 nearby young stars in its vicinity (Weinberger et al., 2000). During the time that a star is part of its birth cluster, large circumstellar disks can be perturbed by close encounters with stars (e.g., see Adams & Laughlin 2001). Many stars are part of binary systems. Consequently the nearby stellar environment may be a key factor influencing the extent and evolution of circumstellar disks. By studying the structure observed in the largest of these disks, we hope to better understand the effect of external perturbers on these disks. We can also probe for the dependence of the spiral structure on the properties of the disk itself, thus allowing us to learn about the formation and evolution of extra-solar planetary systems.

HD 100546 is a nearby southern Herbig Be star (KR Mus; B9.5 Ve; distance $`d=103\pm 6`$pc) with an age estimated $`t10`$Myrs (van den Ancker et al., 1997). HD 141569A (B9.5 Ve) is similar with $`d=99\pm 10`$pc, age $`t5`$Myrs (Weinberger et al., 2000). Infrared spectroscopy with the Infrared Space Observatory revealed solid state emission features in HD 100546 similar to those seen in the spectrum of comet Hale-Bopp (Malfait et al., 1998; Bouwman et al., 2003). While many studies of circumstellar disks have been carried out with spectra (Hu et al., 1989; Malfait et al., 1998; Li & Lunine, 2003), recent imaging studies have been able to resolve the disks out to the large distance of 3–6” or 300–600AU (Augereau et al., 1999b, 2001; Grady et al., 2001; Pantin et al., 2000; Weinberger et al., 2000, 1999). The HST/STIS and ACS studies, in particular, have revealed the presence of intricate structure in these disks, including 2-armed spiral structure at a radius of a few arcseconds or a few hundred AU from the star. The opening angle is large enough that spiral structure is clearly defined in the STIS and ACS images. For HD 100546, the amplitude of the spiral structure is probably low. From the surface brightness profiles shown by Grady et al. (2001) we estimate $`A=\frac{dS}{S}0.15`$ where $`S`$ is the surface brightness at optical wavelengths and $`dS`$ represents the size of azimuthal variations in $`S`$. Because of the difficulty removing the residual diffracted and scattered light from the star (as discussed by Grady et al. 2001), this amplitude is uncertain. For HD 141569A, the spiral structure is high amplitude/contrast in the outermost regions, and lower amplitude at smaller radii (Clampin et al., 2003). We note that because of the different instrument characteristics and filter bandpasses, care should be used in comparing the STIS images (HD 100546; Grady et al. 2001, HD141569A; Mouillet et al. 2001) with ACS images (HD 141569A; Clampin et al. 2003). While HD 141569A’s disk edge is sharply defined, HD 100546 has an envelope that extends out to a radius larger than 1000 AU (Grady et al., 2001).

As suggested by Grady et al. (2001); Clampin et al. (2003), massive external perturbers could be responsible for exciting the observed spiral structure in these circumstellar disks. Alternative possibilities include gravitational instability (e.g., Fukagawa et al. 2004), and an internal perturber such as a planet or brown dwarf (e.g., Mouillet et al. 2001). We first consider the possibility that the disk could be gravitational unstable to the formation of spiral density waves. The Toomre $`Q`$ parameter for a gaseous disk $`Q\frac{c_s\kappa }{\pi G\mathrm{\Sigma }}`$, where $`\kappa `$, the epicyclic frequency at radius $`r`$ is equivalent to the angular rotation rate $`\mathrm{\Omega }`$ for a Keplerian disk. Here $`\mathrm{\Sigma }`$ is the gas density per unit area, and $`c_s`$ is the sound speed (Binney & Tremaine, 1987). When $`Q`$ drops to near 1, the disk becomes unstable to the growth of spiral density waves. Hydrostatic equilibrium can be used to relate the disk thickness to its sound speed, $`\frac{h}{r}\frac{c_s}{v}`$, where the disk aspect ratio is $`h/r`$, $`h`$ is the vertical scale height of the disk and $`v`$ is the circular velocity. Assuming a Keplerian disk, and approximating the mass of the disk $`M_d\pi \mathrm{\Sigma }r^2`$, we find $`Q\frac{M_{}}{M_d}\frac{h}{r}`$, where $`M_{}`$ is the mass of the star. From this we see that if the disk has a low density, a very thin and cold disk is required for $`Q`$ to drop below one or two, allowing the growth of spiral density waves in the disk. Gas densities are expected to be low in circumstellar disks. For example, when $`\frac{M_{}}{M_d}100`$, only an extremely thin disk, $`h/r0.01`$ would be unstable, allowing $`Q1`$. Because this aspect ratio is extreme, we must consider alternative possibilities to account for the spiral structure.

The dispersion relation for tightly wound hydrodynamic waves in a two-dimensional gaseous disk suggests that an internal perturber (such as a planet) would produce spiral arms that are more tightly wound with increasing radius from the star (e.g., Lin & Papaloizou 1993). However, the spiral arms in the disks of HD 100546 and HD 141569A are more open (less tightly wound) with increasing radius, suggesting that the force responsible is increasing with radius. This would be true if the spiral structure was caused by an external tidal force or caused by density waves excited from an external massive body. Consequently we consider here massive perturbers that are located exterior to the location of the spiral arms.

When a star has a close encounter with a planetesimal disk, the tidal force from the star can cause an inclination change in the outer disk (Larwood & Kalas, 2001), scatter objects in an extra-solar Oort cloud analog (Kalas et al., 2001), and cause short lived spiral structure in a gaseous disk (Pfalzner, 2003; Clarke & Pringle, 1993; Boffin et al., 1998; Korycansky & Papaloizou, 1995). Alternatively, in the context of stellar binaries, previous studies have shown that if a secondary star is exterior to a disk, spiral density waves are driven into the disk (Artymowicz & Lubow, 1994; Savonije et al., 1994). Most of the binary/disk simulations of Savonije et al. (1994) showed prominent 2-armed structure similar in morphology to the outermost parts of the HD 100546 system. Tidal forces from a binary can also cause the outermost part of disk to be lopsided, as shown in the case of HD 141569A (Augereau & Papaloizou, 2004).

In this paper we consider these two scenarios that could account for the spiral structure in HD 100546 and HD 141569A. We investigate the feasibility of these scenarios by carring out hydrodynamic simulations. With these simulations we explore the types of encounters can produce 2 armed spiral structure similar to that observed in the disks of HD 100546 and HD 141569A. We aim to produce numerical simulations that exhibit two outer open spiral arms and roughly match the spiral amplitudes in the disks of HD 100546 and HD 141569A. In Section 2 we describe our simulations. In Section 3 we investigate the possibility that an object bound to the star, external to the disk, could tidally excite spiral structure and drive density waves into the disk. In Section 4 we explore the possibility that a recent nearby stellar encounter causes the observed spiral structure. In Section 5 we consider observational constraints on the possible external perturber responsible for the spiral structure in the disk of HD 100546. A summary and discussion follows.

## 2 Simulation of spiral structure caused by external perturbers

To determine the likely cause of the spiral structure in the outer parts of these circumstellar disks, we have carried out two series of simulations. Because we have little information about the vertical geometry and structure of these disks, we have restricted this initial study to 2-dimensional hydrodynamic simulations, confining our study to the plane containing the disks.

Both planetesimal and gaseous disks can be approximated as fluid disks that are described with a sound speed and a kinematic viscosity. For planetesimal disks the sound speed depends on the particle velocity dispersion, but for gaseous disks it is expected to depend on turbulent motions. Based on the images by Grady et al. (2001) we estimate that fine features in HD 100546’s disk have widths $`0.3^{\prime \prime }`$ at a radius from the star of $`3^{\prime \prime }`$. in HD 141569A’s disk exhibits features with widths $`0.6^{\prime \prime }`$ at a radius from the star of $`4^{\prime \prime }`$. The width of these features suggests that structure exists in these disks at a scale height $`h`$ with vertical aspect ratio $`h/r0.15`$. However, because we may only be seeing an illuminated surface (e.g., Dullemond & Dominik 2004), these features may not be directly related to density variations in the disk. Nevertheless if the features observed are interpreted to be in the disk, the disk is not likely to be extremely thick $`h/r>0.2`$. Consequently we assume here that the dynamics of these disks may be approximated with 2-dimensional simulations.

Scattered light optical and near-infrared imaging and thermal infrared observations show that the disks of HD 100546 and HD 141569A are dusty (Malfait et al., 1998; Bouwman et al., 2003; Li & Lunine, 2003; Augereau & Papaloizou, 2004; Weinberger et al., 1999). Both disks are strongly depleted of gas; HCO+ was not detected in HD 100546, Wilner et al. (2003), though cold CO has been detected at low levels in HD141569A Zuckerman et al. (1995). Gas has been detected in the inner few tens of AU in HD 100546 from UV absorption lines of atomic species and molecular hydrogen (Deleuil et al., 2004; Lecavelier des Etangs et al., 2003), and in CO and H3+ from HD141569A Brittain & Rettig (2002). The small dust grains detected in infrared emission from HD141569A need to be replenished on timescales shorter than the age of the system Li & Lunine (2003); Augereau & Papaloizou (2004), consequently HD 141569A is most likely a young debris disk. The vertical optical depth of the HD 141569A’s disk in visible wavelengths at 200 AU is low, $`\tau 0.024`$ (Li & Lunine, 2003). HD 100546’s disk intercepts about 50% of the stellar light and emits prominantly in the mid infrared (Bouwman et al., 2003). The low level of hot dust emitting in the near-infrared, and the presence of crystalline silicate emission features suggest that the disk is depleted at around 10AU and flared outside this radius, and that planets may have formed in the inner region (Bouwman et al., 2003). The disk is probably optically thick in visible wavelengths at 10-30 AU but may not be at larger radii (Bouwman et al., 2003). There is no direct evidence based on the modeling of the infrared spectrum that HD100546 is a debris (rather than gaseous) disk, though Grady et al. (1997) suggested that planetesimals exist in the disk based on spectroscopic variability that could be interpreted in terms of stargrazing comets.

Though the composition of these disks is not precisely known, the possible presence of planetesimals in these disks suggests that a collisionless particle simulations (e.g., Pfalzner 2003; Augereau & Papaloizou 2004) might be a better match to their physical properties. Here we are carrying out hydrodynamic simulations rather than collisionless particle simulations. The primary advantage of the grid based hydrodynamic simulations is that structure can be seen across a fairly large dynamic range in density. This allows us to probe structure at disk edges and resolve small differences in density due to spiral structure. To resolve low intensity spiral structure, particularly at the edge of a disk (such as is seen in HD141569A) or low amplitude structure (such as is seen in HD100546) a large number of simulated collisionless particles would be required. For younger systems such as HD100546, which has a high optical depth, collisions could be frequent enough that particle motions are damped. Due to shocks and viscous effects, hydrodynamic simulations will damp non-circular motion. While neither purely collisionless particle simulations or hydrodynamic simulations are perfect physical models for debris disks, because we can better resolve faint or low amplitude structure, our simulations are complementary to those carried out using simulations comprised of collisionless particles (e.g., Pfalzner 2003; Augereau & Papaloizou 2004).

Both hydrodynamic and particle simulated disks support spiral structure, however there are differences in the dispersion relation for spiral density waves. High frequency spiral density waves (similar to sound waves) are carried by hydrodynamic disks but are damped in particle disks. In both cases the effective sound speed can be estimated using an estimate for the disk thickness and vertical hydrostatic equilibrium.

Our simulations are carried out using the hydrodynamical code developed by Masset (2002); Masset & Papaloizou (2003). This code is an Eulerian polar grid code with a staggered mesh and an artificial second order viscous pressure to stabilize the shocks (see also Stone & Norman 1992). The hydrocode allows tidal interaction between one or more massive bodies (planets) and a 2D non-self-gravitating gaseous disk, and is endowed with a fast advection algorithm that removes the average azimuthal velocity for the Courant timestep limit (Masset, 2000). The simulations are performed in the non-inertial non-rotating frame centered on the primary star (similar to a heliocentric frame). A logarithmic grid in radius was adopted.

The code is scaled so that the unit length is the outer edge of the disk, the unit mass is that of the central star. Time is given in units of one orbital period for a particle in a circular orbit at the outer edge of the disk. For both HD 100546 and HD 141569A with a disk of radius $`300`$AU and stellar mass of $`M_{}2.5M_{}`$, this period is $`P_{outer}3300`$ years.

The disk aspect ratio is defined as $`h/r`$ where $`h`$ is the vertical scale height of the disk. The disk aspect ratio is assumed to be uniform and constant for each simulation. The sound speed of the gas is set from the disk aspect ratio. From hydrostatic equilibrium, we have $`h/rc_s/v_c`$ where $`v_c`$ is the velocity of an object in a circular orbit and $`c_s`$ is the sound speed or velocity dispersion. The simulations are begun with a disk with a smooth and azimuthally symmetric density profile $`\mathrm{\Sigma }(r)\left(\frac{r}{R_{min}}\right)^1`$ where $`R_{min}`$ is the inner disk radius. Because the tidal force from the perturber can pull gas outward, we set the outer boundary to be open, allowing mass to flow outward. The grid inner boundary only allows gas to escape so that the disk material may be accreted on to the central star.

Both particle and gaseous disks exhibit viscosity. In gaseous disks the viscosity is usual related to the $`\alpha `$ parameter, whereas in particle disks, an effective viscosity is dependent on the opacity of the disk, (e.g., Murray & Dermott 1999). In our simulations the disk viscosity is parametrized with the Reynolds number $`r^2\mathrm{\Omega }/\nu `$, where $`r`$ is the radius, $`\nu `$ is the kinematic viscosity, and $`\mathrm{\Omega }`$ the Keplerian angular rotation rate of a particle in a circular orbit of radius $`r`$. In these simulations, the Reynolds number is assumed to be constant with radius. The Reynolds number determines the timescale in rotation periods over which viscosity can cause accretion. The ages of these systems 5–10 million years is only a thousand times the orbital period at the outer disk edges. Because of the short number of orbital periods comprising the lifetime of these disks, viscosity probably does not strongly affect the recent dynamics. We have chosen a Reynolds number for our simulations ($`=10^6`$) which is high enough that viscous effects such as accretion should not be important during the simulations.

We ran two sets of simulations, one set for eccentric bound objects (planets or stars), and one set for stellar flybys or nearby stellar encounters. The eccentric object simulations were run with a radial and azimuthal resolution of $`N_r=150`$ and $`N_\theta =450`$ and the planet or star was begun initially at apoapse. The flybys were run with a resolution of $`N_r=200`$ and $`N_\theta =600`$, and the star was initially begun at a radius of 1.8 times $`R_{max}`$ on a parabolic orbit. Table 1 shows the characteristics of the simulations with bound planets (simulations #A–#L) or stars (simulations #M–#P), and Table 2 lists those involving stellar flybys (simulations #R–#U). Parameters varied are the mass ratio between the perturber and central star, $`q`$ , the disk aspect ratio, $`h/r`$, the eccentricity of the perturber, $`e`$, when it is bound, and the radius from the central star of the perturber’s closest approach, $`R_a`$, when the perturber is on a parabolic encounter. All perturbers were assumed to be on prograde orbits.

## 3  Bound external perturbers

We consider bound external planets and bound external stars separately. Bound low mass objects require many orbits to move material in the disk. The main mechanism for doing this is the excitation of spiral density waves at Lindblad and corotation resonances (e.g., Lin & Papaloizou 1993). However, an external star exerts such a strong tidal force that its Roche lobe can reach beyond the location of most of these resonances. Because HD 141569A has the binary HD 141569B,C in its vicinity (Weinberger et al., 2000) we can consider the possibility that the binary is bound to HD 141569A. However, HD 100546 has no known companions. If it has a bound companion, it is probably low mass. Consequently, we divide this section into two pieces, that investigating the influence of bound low mass companions (appropriate for HD 100546) and that that investigating the influence of higher mass stellar companions (appropriate for HD 141569A).

### 3.1 Bound external planets

In Figure 1 we show the effect of different mass planets exterior to the disk for a planet eccentricity of $`e=0.25`$ (simulations #C, #F, and #I). We see that prominent spiral structure is driven into these disks, similar to those exhibited in previous simulations of stellar binary/disk interactions (e.g., Savonije et al. 1994). For low eccentricity, moderate mass planets, the dominant resonance driving spiral density waves is the 2:1 or $`m=2`$ Lindblad resonance, exciting in the disk two strong spiral arms. The location of the spiral arms is near the location of this Lindblad resonance resonance or at 0.76 times the semi-major axis of the external perturber in the case of a Keplerian disk. The dispersion relation for tightly wound hydrodynamic waves in a two-dimensional gaseous disk is $`\left[m(\mathrm{\Omega }\mathrm{\Omega }_p)\right]^2=\kappa ^2+k^2c^2`$ where $`\kappa (r)`$ is the epicyclic frequency, $`\mathrm{\Omega }_p`$ the angular rotation rate of the spiral pattern, $`m`$ is the number of arms, $`c`$ is the sound speed and $`k`$ the wavevector of the spiral arms. This dispersion relation ignores the self-gravity and vertical structure of the disk. As expected from the dispersion relation, the spiral structure is more open (has a higher pitch angle or lower value of $`k`$) nearing the Lindblad resonance.

While two dominant spiral arms are seen in the simulations shown in Fig. 1, these two arms are not exactly symmetric. In other words, one spiral arm rotated about 180 degrees does not lie exactly on top of the other. As explored by Henry et al. (2003) in the case of M51’s spiral arms, the interplay between different density waves provides one possible simple explanation for such a phenomenon. In these simulations we see a feature extending toward the location of the planet that rotates with the planet, and so is probably associated with the planet’s co-rotation resonance. The asymmetries in the spiral structure are likely to be caused by the time dependent interplay between the $`m=2`$ and other density waves (e.g., $`m=1,m=3`$) excited by the planet. Deviations from perfect bisymmetric symmetry observed in the spiral structure of HD 100546 and HD 141569A could be explained by the interplay between spiral density waves.

When the external planet is on an eccentric orbit, because additional Fourier components exist in the perturbing potential, waves are driven at corotation resonances as well as at Lindblad resonances (Goldreich & Tremaine, 1979). Consequently we would expect stronger $`m=2`$ spiral structure in the simulations containing an eccentric object (e.g., Goldreich & Sari 2003). In Figure 2 we show simulations for different eccentricity planets (simulations #E, #F and #G) but the same semi-major axis. Contrary to our initial expectation, the simulations containing higher eccentricity perturbers have lower amplitude spiral structure. Examination of the simulations shows that the eccentric objects more effectively truncate the disk, in part because the perturbing object passes deeper into the disk. Nevertheless an $`m=2`$ spiral density wave is still driven into the disk, though at lower amplitude. When the perturbing object has eccentricity above 0.23, the object passes within its 2:1 Lindblad resonance, the location where we expect two armed spiral density waves are excited. Since these simulations are shown at only $`t6P_{outer}`$, we infer that that gas in the outer disk is quickly cleared away by the planet. In Figure 3 we show edge profiles corresponding to the simulations shown in Figures 1 and 2. Except for the planet in a circular orbit, the edge of the disk is within the location of the 2:1 Lindblad resonance (at $`r=0.76`$). We see from Fig. 3 that the position of the disk edge is more strongly dependent on the planet eccentricity (and the location of periapse) than on the planet mass.

One way to account for the lower amplitude spiral structure exhibited in the simulations containing higher eccentricity objects is with a model where the wave is driven in low density region at the location of the resonance (e.g., Varnière et al. 2004). Because the waves are excited in a region of lower density, the resulting amplitude of the waves is low even though there is more than one resonance responsible for driving the waves. Because the eccentric object can more effectively truncate the disk, the amplitude of the $`m=2`$ waves can actually be reduced rather than increased. In Figure 2 we also note that the spiral structure in the simulations with higher eccentricity objects is more tightly wound. This is expected from the dispersion relation which implies that the waves become more tightly wound as the distance from the resonance increases.

Figure 1, 2 show images at time $`t6P_{outer}`$ after the beginning of the simulation. This time period corresponds to about $`2\times 10^4`$ years and is significantly less than the age of the star. There are two effects that reduce the density of the outer part of the disk, viscous accretion which takes place on the viscous timescale, and the waves themselves which carry angular momentum. Because these simulations have disk Reynolds numbers of $`=10^6`$, we do not expect to see significant accretion that takes place on a viscous timescale ($`10^6`$ times the orbital period) in our simulations which are only run for 100-1000 orbital periods. However, the timescale for the spiral density waves to transfer material away from the planet is much shorter (Savonije et al., 1994; Varnière et al., 2004). From Figure 2 we can see that the simulations with eccentric objects are particularly effective at driving inward gas flow in the disk edge, whereas the simulations with zero eccentricity or low eccentricity objects accrete gas much more slowly. The disks with low eccentricity external objects maintain their strong $`m=2`$ spiral density waves to $`t>30P_{outer}`$ (of order $`10^5`$ years), whereas those with eccentric perturbers don’t.

From the dispersion relation of tightly wound arms, we expect that thicker disks with a higher sound speed would exhibit more open spiral structure (as previously seen in the simulations of Savonije et al. 1994). Figure 4 shows how the disk aspect ratio (which sets the sound speed) affects the morphology. This figure shows the same simulations at different orbital phases; the leftmost panels show the perturber at periapse and the right most panels show the perturber at apoapse. We see from this figure that the spiral density waves excited do tend to be more tightly wound when the sound speed is lower and the disk thinner. We also find that in the larger scale height disks, the slope of the disk edge is shallower, suggesting that the resonant excitation of the spiral density waves arises from a larger region in the disk (Artymowicz & Lubow, 1994; Lin & Papaloizou, 1993).

From Figure 1 we see that low mass (Jupiter sized), low eccentricity bound perturbers can excite two armed moderate amplitude spiral arms. This morphology is similar to that exhibited by the disk of HD 100546. Because of the short timescale (a few rotation periods at the outer edge of the disk) that gas is cleared to within the 2:1 Lindblad resonance, strong observed open 2-armed spiral structure as seen in the disk of HD 100546 is unlikely to be excited by a highly eccentric ($`e0.20`$) bound planetary perturber. On a timescale of a few rotation periods, the edge of the disk moves sufficiently away from the perturber that strong open two-armed structure is no longer evident.

We now compare the morphology of the simulations to that seen HD 100546. The disk of HD 100546 exhibits 2 arms that are more open at larger radii than at smaller radii. The simulated disks shown in Fig.1 also show spiral structure that is more open at larger radii and more tightly wound at smaller radii. The simulations also show more tightly wound structure at smaller radius that is not detected in the image of HD 100546. This could be in part because of the increased level of light from the central star which would make it difficult to detect fine structure in the disk. Alternatively, the disk of HD 100546 may be comprised of planetesimals; a particle disk would damp tightly wound waves. For objects on circular orbits, we find that the amplitude of the waves drops below about 15%, (as observed for HD 100546), when the mass ratio $`q0.005`$. If a bound perturber is responsible for the spiral structure in the disk of HD 100546, our simulations suggest that the perturber must have mass above 10 Jupiter masses. If the disk is thicker a higher mass perturber would be required. We will use these limits in Section 4 when we discuss observational constraints on these scenarios.

### 3.2 Bound external stars

For HD 141569A, the likely perturber is the companion binary HD 141569B,C, consisting of an M2 and an M4 star, located $`9^{\prime \prime }`$ away from HD 141569A (Weinberger et al., 2000; Augereau & Papaloizou, 2004). We estimate a mass ratio for the binary compared to HD 141569A of $`q0.2`$ based on the low mass stellar models by Baraffe et al. (2002) and effective temperature scale by (Luhman et al., 2003), consistent with $`q=0.2`$ adopted by Augereau & Papaloizou (2004). The amplitude of the outer spiral structure (at radii greater than 300 AU) in the disk of HD 141569A is much higher than that of HD 100546 and the nearby binary provides a likely candidate responsible for driving it. However $`q=0.2`$ is well above the mass ratios considered for the previous simulations. To explore the way that spirals structure could be excited in the disk of HD 141569A, we ran a set of simulations with higher mass external bound perturber with $`q=0.2`$, (see Table 1).

In Figure 5 we show simulations for different binary eccentricities (#M, #N and #O). Simulations #M, #N are shown approximately 5.5 binary orbital rotation periods (BC about A) after the start of the simulation and simulation #O is shown approximately 4.5 orbital rotation periods after the start of the simulation. It is computationaly expensive for us to run finely gridded hydro simulations for significantly longer time periods. The disk is truncated by the binary on its first few passages, but each subsequent close encounter excites spiral structure. The size of the disk is consistent with truncation at the Roche radius of the binary when the binary is at periapse. Consequently the disk outer edge is well within the location of the binary’s 2:1 Lindblad resonance. We see in the simulations that following the first few passages the disk outer edge is much higher density than the initial disk density. The binary has not only truncated the disk, and induced a sharp outer edge but has also compressed it. The high edge density, could be seen as a ring of material, and the lack of gas interior to the ring (e.g., Clampin et al. 2003) might be mistakenly interpreted as evidence that material has been cleared out by another process such as an interior planet.

Observations of HD141569A suggest that the disk is cleared between 150 and 200 AU (Mouillet et al., 2001). Clearing at these inner radii is probably not caused by the perturbations from the binary HD141569B,C, and so is probably induced by internal processes, such as planets (see discussion by Mouillet et al. 2001). The inner boundary of our simulation is sufficiently distant from the truncated edge that it should not have affected the outer disk structure. Spiral waves driven in these simulations are damped before they reach the inner regions of the simulated disk. So the presence of the clearing at 200 AU should not be able to strongly affect the morphology observed in the outer disk at 300–450 AU.

In Figure 5 we show panels preceding and following the binary’s periapse. Before periapse we see the disk is nearly axisymmetric and is quiescent. It lacks strong spiral structure. However at periapse the tidal force from the binary begins to pull two arms from the disk. Since one arm is directly excited by the binary and the other excited by the recoil of HD 141569A (Pfalzner, 2003), there is an asymmetry between the spiral arms. The disk is lopsided when the binary is near periapse. The asymmetric spiral structure, the lopsided disk, and the spiral arm pointing toward the perturber are features that are exhibited by the disk of HD 141569A as seen in scattered light (Augereau et al., 1999b; Clampin et al., 2003). We find that outer spiral arms are less strongly excited by a binary on a more highly eccentric orbit. This is likely because the binary on a higher eccentricity orbit is moving somewhat faster at periapse than one at a lower eccentricity, resulting in a weaker impulse.

We investigate the role of disk thickness by displaying in Figure 6 two simulations (#N and #P) that differ only in the disk aspect ratio. The outer spiral arms excited by the encounter with the binary are more open in the thicker disk, however the amplitude of the features is lower, particularly in the inner region. Also the outer arms point directly toward the binary, rather than leading the binary as is seen in HD 141569A.

We found the closest correspondence between the observed morphology and that seen in the simulations with binary eccentricity $`e0.2`$, a thin disk ($`h/r0.04`$) and binary approaching periapse (simulation #R). This simulation is shown in figure 7 along with the observed morphology pointed out by (Clampin et al., 2003). We support the scenario proposed by Augereau & Papaloizou (2004) in which the disk structure is influenced by perturbations from the binary. Our simulation is successful at producing open and high amplitude two armed spiral structure at the outer edge of a sharp disk, similar to that observed in the HST images at radii greater than 300 AU. It successfully accounts for the disk truncation. The simulation also exhibits some assymetries that are similar to those observed. The simulated northern spiral arm is not 180 degrees away from the southern spiral arm, but shifted or rotated westward so it is closer to the binary than it would be if it were an image of the southern spiral arm rotated by 180. The simulated disk edge is lopsided, with the northern edge thicker and more distant from the star than the southern edge. However, the observed disk edge is brighter on the southwestern side than on the northeastern side and this brightness difference is not seen in our simulation. The northern spiral arm is further east and the southern arm is further west in the simulations than the observed ones.

While we have attempted to match the outermost structure of the disk, Augereau & Papaloizou (2004) concentrated primarily on the structure of the disk edge. Our simulation accounts for the outer spiral structure but is less good at matching the assymetry or lopsided appearance of the bright rim of the disk. This assymetry was well modelled with the secular model of Augereau & Papaloizou (2004) though their simulations did not exhibit the outer spiral arms that are evident in our simulations. Augereau & Papaloizou (2004) and our simulation also differ in the predicted location of the binary periapse. Our simulation has binary periapse near the current location of the binary HD141569B,C, whereas Augereau & Papaloizou (2004) suggest that the periapse is at position angle $`130^{}`$ or on the southwestern side of the disk.

To account for the outer spiral arms we have assumed that the binary is in a prograde orbit. A prograde encounter allows open high amplitude spiral arms to be drawn out of the disk when the binary is near periapse. To truncate the disk, a much closer periapse would be required with a retrograde orbit for the binary. Retrograde encounters are also much less effective at pulling out spiral arms near closest approach. The model of Augereau & Papaloizou (2004) did not specify whether the binary was on a prograde or retrograde orbit. Periapse is on the southwestern side of the disk in their model, so the binary is either post-periapse and retrograde or pre-periapse and on a prograde orbit. Since the phenomena they illustrate and describe in the disk was caused by secular perturbations, the direction of the binary orbit may not be important in their model. Our simulations suggest that much of the outer structure in the disk is transient and excited specifically by the current close encounter. We find that two armed open spiral arms are best excited just past periapse for a prograde binary orbit. This presents a problem for the model of Augereau & Papaloizou (2004) since neither possibility (post-periapse and retrograde or pre-periapse and prograde) would be capable of exciting the outer spiral structure. Modifications of the secular model of Augereau & Papaloizou (2004) will be required to account for the outer spiral structure. We note that the outer spiral arm connecting the HD141569A disk to the stellar companions was only barely detected in the images presented by Clampin et al. (2003), whereas that on the north was detected by both Clampin et al. (2003) and Mouillet et al. (2001). If future observations fail to confirm the presence of the southern arm then a different scenario or very different orbital parameters would be required to account for only the nothern arm.

Our simulation does exhibit some lopsidedness in the disk, even though the binary eccentricity we adopted was not as high as that adopted by Augereau & Papaloizou (2004). The simulations of Augereau & Papaloizou (2004) found that a highly eccentric binary was required to account for the disk asymmetry (Augereau & Papaloizou, 2004) on long timescales. We found that a binary on a more eccentric orbit excited slightly more assymetry in the disk near periapse. However, the disk tended to be offset toward the perturber rather than toward the west, as observed Mouillet et al. (2001). In our simulations, the assymetry damps as the binary recedes following periapse. The hydro simulations explored here are unlikely to maintain the elliptical orbits exhibited by the particle simulations of (Augereau & Papaloizou, 2004). It’s possible that a particle simulation would provide a better match to the physical state of the disk on longer timescalse than our hydro simulation. Our simulations exhibit weak spiral waves driven inside 200 AU that are not observed in the observed images of HD 141569A. Because they are tightly wound, they are similar to sound waves, and would not be present if the disk is comprised of planetesimals.

While our simulation provides a reasonable match to the observed outer morphology we found previously that lower eccentricity binaries and colder disks exhibit higher amplitude and finer structure in the disk. This degeneracy makes it difficult for us to finely discriminate between models and so estimate the eccentricity of the binary’s orbit. While we find somewhat better correspondence between the simulations and observed outer spiral arm morphology with a lower eccentricity binary and relatively thin disk, we cannot exclude the possibility that the binary could have a higher eccentricity, as suggested by Augereau & Papaloizou (2004).

The chosen simulation, #R, shown Figure 7, predicts the outermost 2 spiral arms, matches reasonably well the position of the binary, accounts for the steep disk edge profile, and predicts an asymmetry in the outermost two arms that is seen. However we fail to exactly match the shape of the outermost arms and the structure in the disk edge. Our previous figures showed that this issue cannot be resolved with a binary with higher eccentricity or with a thicker disk because these simulations tend to have less well defined inner disk structure, though because of the degeneracy between binary eccentricity and disk thickness we cannot tightly constrain both quantities. We note that our simulations were carried out in 2D, and it is likely that the binary is not exactly in the same plane as HD 141569A’s disk. The worst discrepancies between our simulations and the observed disk would be in the outer regions where the tidal force from the binary on the disk is largest. By carrying out 3D simulations, it may be possible to more exactly match the morphology and so more tightly constrain both the thickness of the disk and eccentricity of the binary orbit from the disk’s morphology. We have assumed here that a hydrodynamic simulation provides a good approximation to the disk of HD 141569A. However, if this disk is a debris disk, particle simulations that include collisions might provide a better physical model.

Assuming that the binary HD141569B,C is at periapse and an eccentricity of 0.2, the period of the orbit is a few thousand years. Based on an age estimate of 5Myrs for this system, we estimate that a few thousand orbits could have occurred since the birth of the system. Augereau & Papaloizou (2004) inferred a high orbit eccentricity based on the timescale for the maintainance of a lopsided disk. If the orbit of the binary is stable, then many periapse encounters could have taken place, particularly if the binary is not on a highly eccentric orbit. While a few periapse encounters are required to truncate the disk, we expect that repeated encounters could more strongly truncate the disk. As the disk decreases in size, the binary would pull out fainter and fainter arms from the edge of HD141569A’s disk. Nevertheless we have proposed here that the outer spiral arms observed in the disk by Clampin et al. (2003) are a result of tidal excitation from the binary HD141569B,C. We note that the wide separation of HD141569B,C from HD141569A implies that this system may eventually be disrupted altogether by the passage of field stars. The probability of such a field star perturbing the HD141569A,B,C system can be estimated as $`PnvA\mathrm{\Delta }t`$ where $`n`$ is the number density of stars ($`10\mathrm{p}\mathrm{c}^3`$), $`v`$ is the field star velocity dispersion or $`10`$km/s, $`A`$ is the cross section of the binary orbit ($`\pi (1000AU)^2`$) and $`\mathrm{\Delta }t`$ is the timescale (similar to the age of the system). Based on these quantities, we estimate that the HD141569A,B,C system has a probability of $`5\%`$ of a field star encounter during the lifetime of the system. If the binary is on a more highly eccentric orbit then the probability would be higher. One interesting possiblity is that the disk of HD141569A has only recently been truncated by the binary orbit. If the binary is on an eccentric orbit, its orbit could have been perturbed by a nearby field or birth cluster star.

### 3.3 Excitation of spiral structure by flybys

In this section we consider stellar encounters as possible causes for spiral structure in circumstellar disks. This is relevant for HD 100546 which does not have any known nearby companions (Grady et al., 2001).

In Figure 8 we show simulations for parabolic stellar flybys as a function of time after the closest approach. Fig 8 shows simulations for stellar mass ratios of $`q=0.3,0.1`$ and $`0.05`$, simulations #R, #S and #T respectively. We ran these simulations so that the star passed well within the disk to minimize the effect of the outer boundary. Nevertheless at late times, the disk is clearly affected by shocks associated with the outer boundary. Our code is complementary to previous studies because it is good at showing spiral structure, but not as accurate as N-body or SPH codes at modeling elongated tidal arms (as simulated by Pfalzner 2003; Boffin et al. 1998). Initially the gas disk in our simulations exhibits a global two-armed spiral. The size of the secondary arm appears initially to depend on the mass ratio, though at later times the morphology is less assymetric or lop-sided. Pfalzner (2003) showed that the secondary arm is dependent upon the recoil of the central star and our simulations show the same dependence. Since HD 100546 exhibits two-arms, the perturber (if a flyby) probably did not have a small mass ratio. From our simulations we estimate $`q0.1`$, corresponding to $`M0.2M_{}`$, since lower mass perturbers fail to excite two strong spiral arms.

After 4 rotation periods (at the outer disk edge) little strong spiral structure remains. This short timescale is consistent with the results of previous simulations (Pfalzner, 2003), and allows us to place limits on the time since the closest approach. The spiral structure HD 100546 is likely to have been excited less than one rotation period (at the outer disk edge) ago, which is less than 3000 years ago based on the rotation period at the edge of HD 100546’s disk. For an object traveling at 10 km/s this would correspond to a distance on the sky of 1 arcminute.

Since the pitch angle and wave travel time depends on the disk sound speed, we ran a comparison simulation with $`h/r=0.1`$ (simulation #U) that is shown in Figure 9. We see that the amplitude of the spiral density waves is reduced in the thicker disk, even though they are more open. Because the amplitude of the spiral structure is lower, a higher mass perturber that encountered HD 100546 more recently would be required to account for the observed spiral structure.

In Figure 10 we show a comparison between an encounter and the observed morphology of HD 100546. The simulation is shown about 3000 years ($`t=1P_{outer}`$) after the encounter for a mass ratio $`q=0.3`$ and with disk aspect ratio $`h/r=0.04`$ (simulation #R). The simulation exhibits an asymmetry in the spiral structure similar to that observed, though the easternmost spiral arm is not as open as the observed one. This figure demonstrates that it is easy (and so not particularly meaningful) to simulate spiral structure similar to that observed in the disk of HD 100546. We must identify the actual object responsible for exciting the spiral structure to constrain this scenario further.

## 4 Observational limits on HD 100546’s perturbing object

As discussed by (Grady et al., 2001; Augereau et al., 2001) the STIS, NICMOS and ADONIS observations strongly constrain the types of nearby objects that could be associated with HD 100546. The nearby point sources detected in the STIS images are all too blue to be associated with the HD 100546 system, and so are likely to be background stellar sources (Grady et al., 2001). This puts strong limits on both of our scenarios. If a flyby caused the spiral structure, then the object responsible must be outside the field of view of the STIS image ($`20^{\prime \prime }`$ away from HD 100546). Alternatively if a low mass bound object causes the spiral structure, then the mass required to excite the spiral structure must be low enough so that it would not have been detected in the existing images. Our simulations suggest that a bound mass at least as large as 0.01 times the mass of the star is needed to match the amplitude of the spiral structure. This would correspond to a mass of $`0.02M_{}`$ or $`20M_J`$ (see Section 3.1).

The division between bound perturber and flyby is important. A bound perturber will make more than one encounter over a timescale of a few thousand years. Since a first encounter disrupts the disk, the remaining disk may not exhibit strong spiral structure during subsequent approaches. Since a flyby disrupts a previously quiescent disk, strong outer spiral structure can be excited, but it would be transient, only lasting a few thousand years. Afterward the dirupted disk could extend well past the pericenter of the encounter.

The lack of detected associated objects (point sources) in the NICMOS and STIS images can be used to estimate an upper limit on the mass of possible objects near HD 100546. The models of Baraffe et al. (2002) predict that 5$`M_J`$ and 10$`M_J`$ mass objects that are 10 Myr old have absolute H magnitudes of 15.3 and 11.8 respectively. These correspond to H magnitudes of 20.3 and 16.8 for objects at the distance of HD 100546. The NICMOS image of (Augereau et al., 2001) detected a star of H magnitude 16.1 at a distance of 4” from HD 100546, suggesting that objects somewhat fainter could have been detected nearer HD 100546. The disk brightness drops to 17.5 mag/arcsecond<sup>2</sup> in H band at a radius of 3” from the star (Augereau et al., 2001). Considering the NICMOS point spread function, a star of magnitude 19-20 could have been detected against the disk at a radius of 3” from the star. If we adopt $`10M_J`$ ($`0.01M_{}`$) as a limiting mass for a coeval object that could have been detected in the STIS and NICMOS imaging to a radius of 3”, then we can limit the possible bound companions of HD100546 within $`300AU`$ to a mass ratio $`q0.005`$. The perturber mass limit suggested by the simulations, $`M>0.02M_{}`$ is in contradiction to the limits placed on possible companions by the lack of associated objects detected in HST/STIS and NICMOS observations. We conclude that a bound low mass companion is not a good explanation for the observed spiral structure in the disk of HD 100546.

We now consider observational constraints on objects that could have previously encountered HD 100546. USNO-B1.0 is an all-sky catalog complete to a visual magnitude of 21 that presents positions accurate to 0.2” and proper motions measured across 50 years of photographic sky surveys (Monet et al., 2003). In the vicinity of HD 100546 (within 2 arcminutes) the USNO-B catalog (using surveys during 1976, 1984 and 1991) yields one nearby object with a similar proper motion to HD 100546, USNO-B1.0 ID #0197-0270043, that is 54” arcseconds away from HD 100546. This star with $`M_R=17.4`$ does not coincide with any 2MASS object, and it is quite blue ($`BR0`$) suggesting that it is a white dwarf rather than a low mass main sequence star. Consequently it is unlikely to be the same age as HD 100546 and so is unlikely to be associated with the star. Inspection of the USNO-B digitized images suggests that the proper motion is uncertain due to two nearby stars, one which does have a 2MASS counterpart.

We also searched the USNO-B1.0 catalog for objects that might not be associated with HD 100546 but could have intersected the path of this star on the sky. We find no moving candidates that are red enough to be at the same distance as HD 100546 from us and that could have approached (within 15,000 AU) the path of HD 100546 in the past 9000 years. Both the 2MASS and USNO-B catalog contain numerous red faint sources within a few arcminutes of HD 100546 that could either be candidates for low mass stellar bound companions or low mass stars that could have encountered the star in the past few thousand years. We failed to find an extremely red object that was seen in 2MASS (limiting magnitude $`K15`$) and was not seen in the USNO-B surveys (limiting magnitude $`V21`$). The 2MASS all sky survey point source catalog contains 23 stars with near-infrared color $`JK>0.9`$ within 2’ of HD 100546. These are candidate low mass main sequence stars at the distance of HD 100546 lacking measured proper motions from the USNO-B catalog. Our failure to find a catalog candidate for the object responsible suggests that it must be a low mass object. Stars later than M6 at the distance of HD 100546 (100 pc) would not be visible in the USNO-B surveys. Brown dwarfs at the distance of HD 100546 would not be seen in the 2MASS survey.

The low proper motion of HD 100546 with $`\mu `$RA $`=38`$mas/yr and $`\mu `$DEC $`=2`$mas/yr makes it difficult to pinpoint faint stars at the same proper motion using ground based data across only a 15 year baseline (USNO-B1.0 is based on surveys in 1976, 1984 and 1991 in this region), during which time HD 100546 would have moved by only $`0.6\mathrm{"}`$.

Stellar collisions or flybys in the field are unlikely (e.g. Kalas et al. 2001). Assuming a stellar density of 10 stars pc<sup>-3</sup> and a stellar velocity dispersion of 10km s<sup>-1</sup> typical of the Milky Way density in midplane, and a cross section of $`\pi (600\mathrm{A}\mathrm{U})^2`$ for HD 100546’s disk, we would estimate a probability of $`0.03`$ that a low mass star could have approached within 600 AU of HD 100546 in the past 10 Myrs. The probability is only $`3\times 10^6`$ that a field star approached HD 100546 in the past 10,000 years, the upper limit on the time since the encounter we have inferred based on the the morphology of the spiral structure. However, the birth cluster of HD 100546 may have provided a denser stellar environment, increasing the probability of a stellar collision (Kalas et al., 2001; Adams & Laughlin, 2001). As pointed out by Vieira et al. (1999), HD 100546 is located just to the south east of the dark filament DC296.2-7.9. Through their study of the extinction in this cloud Vieira et al. (1999); Corradi et al. (1997) find this cloud is likely to be just past HD 100546, 110–130 pc away, and so could be associated with HD 100546. The statistics of star colors in the region suggest a moderate mininum extinction of $`A_V0.5`$ is provided by the cloud (Vieira et al., 1999). From inspection of 2MASS images near HD 100546 we find no obvious stellar cluster associated with the dark filament. Unless a region of higher stellar density is found within a few arcminutes of HD 100546, the encounter would have been an exceedingly improbable and exceptional event.

## 5 Summary and Discussion

In this paper we have used 2D hydrodynamical simulations of disks perturbed externally by stars or planets to investigate possible scenarios which can account for outer spiral structure observed in the circumstellar dusty disks of HD 100546 and HD 141569A. We consider two types of external perturbers, those in orbit about the central star containing the disk and those involving a recent stellar encounter or flyby. Both scenarios produce morphology similar to that observed, moderately open 2-armed spiral structure.

Our simulations exhibit some trends. Thicker disks tend to have lower amplitude and more openly wound spiral density waves. Eccentric planetary perturbers more efficiently clear gas away from the perturber. At earlier times in the simulations, eccentric planetary perturbers excite strong and open spiral structure, however after a few rotation periods, the edge of the disk moves inward, and the waves excited become more tightly wound and amplitude drops. The two armed structure excited tends to be asymmetric, a phenomena we attribute to the interplay between $`m=2`$ and additional spiral density waves at different angular rotation rates that are also excited by the perturber.

When the bound object is higher (stellar) mass, the disk is truncated at the Roche lobe of the companion. Each time the companion approaches periapse, tidal arms can be pulled out from the disk and spiral density waves driven into the disk. We see from the simulations that the disk is not only truncated by the binary but also can be compressed, causing a high gas density at the disk edge resulting in the appearance of a gas ring.

We find that the outermost two spiral arms of HD 141569A’s disk are qualitatively reproduced by a perturbation from its nearby binary companion HD 141569B,C, confirming the scenario proposed by Augereau & Papaloizou (2004). However, we find closer morphological correspondence with thinner disks ($`h/r0.04`$), moderately $`e0.2`$, but not highly eccentric orbits for the external binary with the binary located near the periapse of its orbit. Because a higher binary eccentricity causes a similar effect on the morphology as a thicker disk, these two parameters are redundant and we cannot exclude the possibility that the binary eccentricity is higher. Our simulation accounts for only some of the asymmetries in the inner region (200-300 AU). It matches reasonably well the position of the binary, and predicts an asymmetry in the outermost two arms that is seen in the disk of HD 141569A. However, our simulations also are not as lopsided as the observed disk and do not exhibit the high density on the southwestern side of the disk that was seen in the observations by Clampin et al. (2003). Our simulations do not do as good a job as the secular collisionless particle simulations of Augereau & Papaloizou (2004) at accounting for the lopsided shape of the disk edge, though those simulations do not exhibit outer spiral arms. Our simulations also fail to exactly match the shape of the outermost spiral features. We suspect that a 3D simulation would be required to do this. Our simulations also exhibit an enhanced density at the disk edge (similar to that observed), and episodes of high accretion in the inner disk following the binary’s close passage. In a planetesimal dominated disk (rather than a hydro one), these episodes could correspond to episodes of increased collisions and dust production (Kenyon & Bromley, 2002).

The current uncertainty in the composition of these disks hampers our ability to simulate them. These disks are suspected to contain planetesimals, but collisions are also expected to be important over the lifetime of the disks. Here, we have opted to carry out 2D hydrodynamic simulations as a first step. In future 3D particle simulation which include collisions would provide a better physical match to the suspected properties of debris disks. The simulations of Augereau & Papaloizou (2004) investigate long timescale secular perturbations and suggest that much of the morphology of this disk has been caused by repeated encounters with the binary HD 141569B,C. In contrast our simulations suggest that some of the structure (the outer spiral arms) in HD141569A’s disk has been caused by a recent encounter with the binary. From attempting to match the orientation of the the tidally excited spiral arms, we infer that the binary HD141569B,C is near periapse and that the binary is on a prograde orbit. Our estimated location periapse differs from that predicted by Augereau & Papaloizou (2004). Future numerical work could attempt to account for both the outer spiral arms (as done here) as well as the disk edge structure (as better done by Augereau & Papaloizou 2004).

In the case of HD 100546, a low mass planetary perturber in a circular orbit might cause low amplitude two-armed spiral structure that is more open at larger radii than smaller radii. However, to match the observed spiral amplitude, $`15\%`$ , the bound object must have a mass greater than $`M>0.02M_{}`$, in contradiction to the limits placed by the lack of detection of possible associated objects in the NICMOS and STIS images of Grady et al. (2001); Augereau et al. (2001), ($`M0.01M_{}`$ if the object is the same age as HD 100546). A recent stellar encounter might also be responsible for the spiral structure. Our simulations suggest that an star of mass about $`0.1M_{}`$ passing within 600 AU of HD 100546 within a few thousand years ago could account for the spiral structure. The timescale since the encounter is long enough that the star could be outside the field of view of the STIS and NICMOS images. However, we have failed to identify a candidate from the USNO-B catalog. The probability that a field star was responsible for the collision in the past few thousand years is very low, $`10^6`$. If the star responsible for the encounter originated from a region (not yet identified) with a higher stellar density located within a few arcminutes of HD 100546 then the probability of the encounter might be higher. Here we have discussed three scenarios which could account for the observed spiral structure in the disk of HD 100546. Only an extremely cold and thin disk could be unstable to the formation of spiral density waves. For the reasons discussed above, bound perturbers also provide unlikely explanations. This has left only the possibility that a recent stellar encounter caused the spiral structure, however this scenario is highly improbable and the star responsible has not yet been identified. This implies that the observed spiral structure currently lacks a good explanation. A recent studies of the infrared spectral energy distribution of HD100546 suggests that it’s outer edge is flared (Dullemond & Dominik, 2004; Meeus et al., 2001). Consequently structure observed in the visible bands could be due variations in the vertical structure. Future scenarios for HD100546’s disk might consider flared disk models that exhibit variations in both vertical and azimuthal structure. Promising possibilities include transient spiral structure excited by planetesimal collisions (e.g., Kenyon & Bromley 2004).

The division between bound perturbers and flybys is important. While both excite spiral arms, a single flyby will not truncate the disk. However an eccentric bound perturber, because it makes multiple close passages is likely to truncate the disk. Notably HD 141569A is an example of the second case but previous studies have proposed that HD 100546 could be an example of the first case or a flyby. HD 100546 has a large disk extending out to 500AU, past the location of the spiral structure. HD 100546 also has an envelope that extends out to 1000AU (Grady et al., 2001). We expect that the process of disk truncation depends on the velocity of the perturber. If the perturber is moving fast, then an impulse approximation is appropriate and disk particles nearer the point of closest approach are more strongly affected than those on the opposite side, particularly if the object is low mass (Pfalzner, 2003). In contrast, the repeated encounters of the binary HD 141569B,C enforce a sharp outer disk boundary.

We note that the simulations presented here suffered from various restrictions. They were 2-dimensional making it impossible to investigate the 3-dimensional structure of these disks. They suffered from artifacts caused by the boundary conditions. Future studies could be carried out in 3 dimensions. Larger grids could also be used, lessening the artifacts caused by the boundary. 3D simulations may provide a better match to the outer morphology of HD 141569A’s disk, and make it possible to more tightly constrain its disk thickness and companion’s orbit. In the future we also hope to simulate the effect of bound perturbers and flybys on disks containing planets, in search of explanations for the eccentric planets proposed to explain the structure of Vega’s, HR 4796A’s and Epsilon Eridani’s disks (Wilner et al., 2002; Wyatt et al., 1999; Quillen & Thorndike, 2002). The simulations carried out here were done using gaseous disks, however these disks may contain massive planetesimals. While our flyby hydro simulations exhibit similar morphology to the collisionless simulations of Pfalzner (2003), future work which investigates the dynamics of disks containing both planetesimals that can collide and gas may be able to place constraints on the constituents of these disks from their morphology. Higher angular resolution, and deeper imaging spanning additional wavelengths will also provide more constraints on the dynamical processes affecting these disks.

Since most stars are born in stellar clusters and are binaries, encounters may strongly influence forming planetary systems. We find here that the binary companion of HD 141569A has likely truncated and compressed its disk. Future close encounters will also excite spiral density waves in the disk. This type of excitation, taking place on timescales of thousands of years, could contribute to periodic episodes of high accretion in gaseous disks or high dust production in planetesimal disks. We note that flybys excite structure for only short periods of time, so future studies can search for and identify objects capable of perturbing young stellar disks. Once the object responsible for exiting the spiral structure is identified, the nature of the encounter can be explored more quantitatively, and we expect that many of the free parameters affecting our simulations (such as orbit, disk thickness, boundary conditions and disk composition) can be better constrained.

We thank Mike Jura, Dan Watson, Pawel Artymowicz and Carol Grady for helpful discussions and comments. We thank F. Masset for his help and for providing the code. Support for this work was provided by NSF grants AST-0406823, AST-9702484, AST-0098442, DOE grant DE-FG02-00ER54600, and the Laboratory for Laser Energetics. This material is based upon work supported by the National Aeronautics and Space Administration under Grant No. NNG04GM12G issued through the Origins of Solar Systems Program and grant NAG5-8428, This research was supported in part by the National Science Foundation to the KITP under Grant No. PHY99-07949.


text:
# Ejection of bi-helical fields from the sun

## 1. Introduction

Dynamo action is possible both for helical and nonhelical turbulence \[see, e.g., Meneguzzi et al. (1980) for early numerical work\]. However, the reason one is particularly interested in helical dynamos is that there is a well-known mechanism (the $`\alpha \mathrm{\Omega }`$ mechanism) driving cyclic behavior together with latitudinal migration of a large scale field. Both cyclic behavior and latitudinal migration are important features of dynamos in the sun and many late-type stars.

The $`\alpha `$-effect is formally introduced via the longitudinally averaged induction equation. This effect characterizes the conversion of toroidal mean field into poloidal via a sequence of events, all of which produce similarly oriented loops that are tilted in a clockwise sense against the toroidal direction in the northern hemisphere, and anticlockwise in the southern. This effect is the result of the combined action of Coriolis force and nonuniformity of the turbulence, as embodied by the formula (e.g. Rüdiger & Kitchatinov 1993)

$$\alpha _0\tau ^2\overline{𝐮^2}𝛀\mathbf{}\mathrm{ln}(\rho \overline{𝐮^2}).$$
(1)
This equation, which has been obtained using different approaches (see Moffatt 1980), reflects the fact that the turbulent velocity field has attained kinetic helicity. The subscript 0 indicates that the nonlinear feedback is not included.

## 2. Connection with magnetic helicity

It is first of all important to realize that practically no net magnetic helicity can be generated in the sun. What is possible, however, is a segregation of the magnetic field into its positively and negatively helical constituents. Such a segregation can occur either in wavenumber space or in real space. The latter can be accomplished by differential rotation acting on field lines crossing the equator, while the former can be the result of the $`\alpha `$-effect. This will now be discussed in more detail.

Consider the rise of a toroidal flux tube lifted upward either by thermal or magnetic buoyancy. As it rises, it gets tilted by the Coriolis force and, as a consequence, it must develop some internal twist. This is already quite clear from a simple sketch (Fig. 1), where the flux tube is depicted as a ribbon, so one can more easily trace the induced twist.

Mathematically, magnetic helicity is the sum of writhe and twist helicities, and their sum must stay nearly unchanged (magnetic helicity conservation). In the example shown in Fig. 1, the overall structure of the tube follows the right hand rule, corresponding to positive writhe helicity. At the same time, the interval twist of the tube follows the left hand rule, corresponding to negative twist helicity. On the southern hemisphere both signs would be reversed.

## 3. Writhe and twist as driver and quencher in dynamo theory

We have already eluded to the fact that the tilted magnetic field from the $`\alpha `$-effect contributes directly (via many systematic and similar events) to the large scale poloidal field, $`\overline{𝐁}_{\mathrm{pol}}`$. Since $`\alpha `$ is positive, it should give a positive contribution to the current helicity of the mean field, i.e. $`\overline{𝐉}\overline{𝐁}>0`$. We also know from helicity spectra of magnetic flux tube experiments that the negative internal twist contributes at scales smaller than the positive writhe (Blackman & Brandenburg 2003). To put this into numbers, we expect the scales of the large scale poloidal field to be on the order of around 300 Mm, which corresponds to the latitudinal extent of the toroidal flux belts of around 20-30. The typical scale associated with the twist is expected to be smaller, perhaps 30-300 Mm, but the dividing line between small and large scales may not be very clear. (Obviously, what is small scale to a dynamo theorist may be large scale to a solar observer!)

The internal twist is of tremendous importance in dynamo theory. As explained in the previous section, if the writhe helicity is composed of large scale field, the twist helicity must be composed of small scale field (where ‘small’ could be anywhere between 30-300 Mm!). We may therefore identify the helicity from the internal twist with the current helicity of the small scale field, $`\overline{𝐣𝐛}`$. Here, lower characters denote the small scale field, i.e. we have decomposed the field as $`𝐁=\overline{𝐁}+𝐛`$ and likewise the current as $`𝐉=\overline{𝐉}+𝐣`$.

The significance of the $`\overline{𝐣𝐛}`$ term is that it gives an extra contribution to the $`\alpha `$-effect (Pouquet et al. 1976),

$$\alpha =\alpha _0+\frac{1}{3}\tau \overline{𝐣𝐛}/\rho ,$$
(2)
but it acts such as to suppress or quench the dynamo ($`\alpha _0`$ and $`\overline{𝐣𝐛}`$ tend to have opposite signs). Indeed, the $`\overline{𝐣𝐛}`$ term constitutes the main nonlinearity of $`\alpha `$-effect dynamos. There are other nonlinearities, such as a direct suppression of terms in equation (1), but such nonlinearities are never catastrophic in the sense that they do not depend on the magnetic Reynolds number111 By contrast, equation (2) can lead to a catastrophic quenching formula: use Keinigs’ (1983) formula for the saturated state, $`\alpha =\eta \overline{𝐣𝐛}/\overline{𝐁}^2`$ (which follows from the magnetic helicity equation), eliminate $`\overline{𝐣𝐛}`$ from (2), to get catastrophic quenching: $`\alpha =\alpha _0/(1+R_\mathrm{m}\overline{𝐁}^2/B_{\mathrm{eq}}^2)`$. Here we have used $`R_\mathrm{m}=\eta _\mathrm{t}/\eta `$, with $`\eta _\mathrm{t}=\frac{1}{3}\tau \overline{𝐮^2}`$ and defined $`B_{\mathrm{eq}}=(\mu _0\rho \overline{𝐮^2})^{1/2}`$ to eliminate $`\tau `$; see Eqs (18) and (40) of Blackman & Brandenburg (2002) for a more general formulation..

## 4. Why coronal mass ejections might be good for the dynamo

Active regions and coronal mass ejections are the main contributors of magnetic helicity flux from the solar surface (e.g. Démoulin et al. 2002). It might not be appropriate to associate them with small scale losses only. Instead, in view of the combined presence of large and small scale fields in a single field structure (Fig. 1), such losses occur probably simultaneously at large and small scale fields. In this section we argue that this leads to an optimal scenario for dynamos.

Phenomenologically, surface losses of large and small scale helical fields may be described by diffusion terms. This concept has been tested against simulations in the case where losses occur almost entirely at large scales (Brandenburg & Dobler 2001). The results are twofold: the saturation time is formally decreased, but only at the expense of lowering the saturation field strength.

The rest of a possible success story is still speculation, but is based on physical reasoning and a numerical experiment: if losses occur preferentially on small scales, then the saturation amplitude of the large scale field is increased. This is because $`\overline{𝐣𝐛}`$, the ‘quencher’ in dynamo theory, is constantly being removed. This is confirmed by simulation where small scale magnetic field is artificially removed (Fig. 2). Even in the absence of such losses, the saturation value of the large scale field may be sufficient for the sun, but the cycle period may be too long. The primary importance of small scale losses may actually be to maintain the observed cycle period. This needs further investigation.

## 5. Conclusions

An important missing link in the story outlined above is, in our view, a simulation of a coronal mass ejection and an analysis in terms of the magnetic helicity budget. At the same time, it will be necessary to improve mean field theory to take helicity fluxes properly into account. For attempts in that direction we refer to recent work by Kleeorin and collaborators (2003).

## References

Blackman, E. G., & Brandenburg, A. 2002, ApJ, 579, 359

Blackman, E. G., & Brandenburg, A. 2003, ApJ, 584, L99

Brandenburg, A. 2001, ApJ, 550, 824

Brandenburg, A., & Dobler, W. 2001, A&A, 369, 329

Brandenburg, A., Dobler, W., & Subramanian, K. 2002, AN, 323, 99

Démoulin, P., Mandrini, C. H., van Driel-Gesztelyi, L. 2002, ApJ, 382, 650

Keinigs, R. K. 1983, Phys. Fluids, 26, 2558

Kleeorin, N. I, Kuzanyan, K., Moss, D., et al. 2003, A&A, 409, 1097

Meneguzzi, M., Frisch, U., & Pouquet, A. 1981, PRL, 47, 1060

Moffatt, H. K. 1978, Magnetic Field Generation (Cambridge University Press)

Pouquet, A., Frisch, U., & Léorat, J. 1976, JFM, 77, 321

Rüdiger, G., & Kitchatinov, L. L. 1993, A&A, 269, 581