N -body simulations of the Magellanic stream

A suite of high-resolution N -body simulations of the Magellanic Clouds–Milky Way system are presented and compared directly with newly available data from the H I Parkes All-Sky Survey (HIPASS). We show that the interaction between Small Magellanic Clouds (SMC) and Large Magellanic Clouds results in both a spatial and kinematical bifurcation of both the stream and the leading arm. The spatial bifurcation of the stream is readily apparent in the HIPASS data, and the kinematical bifurcation is also tentatively identiﬁed. This bifurcation provides strong support for the tidal disruption origin for the Magellanic stream. A ﬁducial model for the Magellanic Clouds (MCs) is presented upon completion of an extensive parameter survey of the potential orbital conﬁgurations of the MCs and the viable initial boundary conditions for the disc of the SMC. The impact of the choice of these critical parameters upon the ﬁnal conﬁgurations of the stream and leading arm is detailed.


I N T RO D U C T I O N
The progressive collapse and merging associated with hierarchical clustering within a cold dark matter cosmology, while dominated by activity at early epochs (redshifts z ∼ 2-5; e.g. Murali et al. 2002), continues to the present-day and is readily observable even in the local Universe. In our own Milky Way (MW), ongoing satellite disruption and accretion events include the Sagittarius dwarf (Ibata, Gilmore & Irwin 1994), the Canis Major dwarf (Martin et al. 2005, and references therein), and perhaps the most spectacular of all, the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC). The disruption and accretion of the Magellanic Clouds (MCs) is perhaps best appreciated through the nearly circum-Galactic polar ring of gas -the Magellanic stream -emanating from the clouds (Mathewson, Cleary & Murray 1974;Putman et al. 1998).
Two primary, competing, scenarios have been postulated to explain the origin of the Magellanic stream (MS): (i) ram-pressure stripping of LMC and SMC gas due to the motion of the clouds through the tenuous coronal gas in the Galactic halo (Moore & Davis 1994;Mastropietro et al. 2005) -this 'drag' scenario faces difficulty in explaining the leading arm feature (LAF) observed by Putman et al. (1998); (ii) tidal disruption of the SMC (Murai & Fujimoto 1980). Because of its ability to simultaneously produce both trailing and leading streams of gas, the 'tidal' model has gathered considerable support in the literature (Gardiner & Noguchi E-mail: tconnors@astro.swin.edu.au 1996;Gardiner 1999;Yoshizawa & Noguchi 2003, hereafter, GN96, G99 and YN03, respectively).
Constraints upon theoretical models of the MS have improved dramatically with the release of the H I Parkes All-Sky Survey (HIPASS; Barnes et al. 2001) data set. The extent and fine-scale structure of the MS can now be appreciated to a level not previously possible, in particular the unequivocal evidence for the existence of the LAF (Putman et al. 1998). Motivated by HIPASS, we have initiated a programme of high-resolution N-body modelling of the MS aimed at de-constructing the temporal evolution of the MC-MW interaction. Here, we use ∼30 times higher resolution than previous models such as GN96, G99 and YN03, 1 giving us a H I mass resolution of ∼5600 M and a H I flux resolution of 7 Jy km s −1 . In addition, we construct a detailed H I map, to compare with the HIPASS data directly and quantitatively. The combination of a higher-resolution simulation and the simulated H I map enables us to argue about more detailed features in the MS and the LAF. As a result, this paper shows more convincingly that the observed LAF and MS can be produced by the leading and trailing streams of the SMC, induced by tidal interactions with the MW and LMC. Our fiducial model is found after an extensive parameter survey of the potential orbital configurations of the MCs and the different initial condition of the SMC. Based on the parameter survey, we also demonstrate how the final features of the MS and the LAF are sensitive to the initial configuration of the SMC. Our new high-resolution simulations also reveal that the tidal interactions create spatial and kinematical bifurcation in the MS and LAF. We present this prediction based on our simulations, however we also demonstrate that the existing HIPASS data show such bifurcations.
The work described here (Paper I) is the first in a series of papers aimed at providing a definitive model for the MS and the LAF. Here, we concentrate solely upon the effects of gravity, describing the SMC disc with collisionless particles, and ignore the baryon physics (including hydrodynamics, star formation, energy feedback and chemical enrichment). In Paper II, the effects of baryonic physics will be detailed.
The layout of this paper is as follows. In Section 2, we provide a description of the suite of simulations generated to date. In Section 3, we show the results from our best model, and compare them with the observational data. In Section 4, we demonstrate how the final features of the MS and LAF depend upon the orbits of the MCs and the initial properties of the SMC disc. Finally, in Section 5, the discussion and future directions for our work are presented.

N U M E R I C A L S I M U L AT I O N
The framework upon which our simulations are based parallels that described by GN96 and YN03. The MW and LMC are taken to be fixed potentials, while the SMC is treated as an ensemble of self-gravitating particles, in recognition of the fact that the MS is thought to originate from the tidal disruption of the SMC disc (e.g. Gardiner, Sawa & Fujimoto 1994, hereafter GSF94;GN96;Maddison, Kawata & Gibson 2002;YN03; but see also Mastropietro et al. 2005). The orbits of the MW and LMC with respect to the SMC are pre-calculated, the procedure for which is outlined in Section 2.1 and GSF94. The initial boundary conditions for the 'live' SMC model is described in Section 2.2, and its evolution explored in Section 2.3.

The orbits of the LMC and SMC
We assume a spherically symmetric potential for the MW halo. The MCs sample a volume of the halo where the gravitational potential is insensitive to both the expected central cusp and the outer regions where the density profile may be steeper than an isothermal profile. We hence assume a constant rotational velocity of V c = 220 km s −1 within the MW halo. Thus, the potential is described by with a mass enclosed within r kpc of M G (< r ) = 5.6 × 10 11 V c 220 km s −1 2 r 50 kpc M .
We assume that there is little disc contribution to the potential at the typical Galactocentric distances of the LMC and SMC (∼50 kpc). The haloes of both the MW and LMC are assumed to be invariant for the duration of the simulation (2.5 Gyr -i.e. a look-back time corresponding to redshift z ∼ 0.2). Plummer potentials are adopted for both the LMC and SMC, i.e.
where r L,S are the positions of the clouds relative to the MW centre, and K L,S are the core radii, set to 3 and 2 kpc for the LMC and SMC, respectively. In the fiducial model, we assume a constant mass of M L = 2 × 10 10 M for the LMC, and M S = 3 × 10 9 M for the SMC.
Adopting these parameters, we next backwards integrate the orbits of the SMC and LMC, using as boundary conditions the current tangential and radial velocities and positions of the clouds listed in Table 1 (see also Murai & Fujimoto 1980;GSF94;GN96;Lin, Jones & Klemola 1995;Bekki et al. 2004). The values we adopted for our work were chosen to match those of GN96, and are consistent with the extant literature.
The equations of motion of the clouds about the stationary MW arë and where the potentials φ L , φ S and φ G refer to the LMC and SMC and the Galaxy, respectively (Murai & Fujimoto 1980). Since we do not model the MW as live particles, dynamical friction is modelled as per GSF94. F L and F S are the dampening force between the Galaxy and each of the clouds: where ln ( ) is the Coulomb logarithm. Although equation (6) is a simple analytical formulation for dynamical friction that fails to model accurately the orbit of the clouds for more than 5 Gyr, according to recent studies (Hashimoto, Funato & Makino 2003;Just & Peñarrubia 2005), it gives a good approximation to the predictions of full N-body simulations over the ∼2.5 Gyr period we focus on in this paper. We note also, the strong effect the Coulomb logarithm has on the orbits. Hashimoto et al. (2003) suggest that the ln ( ) value of 3 advocated by Binney & Tremaine (1987) for the LMC is too large, causing the system to evolve too fast. Thus, for this paper, we set ln ( ) = 1 to enable direct comparison with the previous works.  Our fiducial model is initiated at a look-back time of 2.5 Gyr when the clouds are at apo-Galacticon (with the tidal forces between all three bodies being minimized). The subsequent orbit appears in Fig. 1, and the distance between the clouds and Galaxy is shown in Fig. 2.

SMC model
The initial dynamical configuration for the SMC was constructed using GALACTICS (Kuijken & Dubinski 1995). We employed a bulgeless equilibrium model with a truncated exponential disc, and a dark matter King-profile halo. The relevant GALACTICS parameters for the fiducial SMC model appear in Table 2 (a detailed explanation of the individual parameters is provided by Kuijken & Dubinski 1995). The disc possesses an exponential profile with scale radius  d = 3.5 kpc, smoothly truncated beyond R t = 7 kpc (with 95 per cent of both the disc and halo masses being within r d,95 ∼ r h,95 ∼ 7 kpc), to give a disc with total radial extent (where the face-on surface density of the disc reaches 1 M pc −2 ) of 7.5 kpc (compared to the halo with radial extent ∼14 kpc). The rotation curve peaks at ∼2 kpc with a velocity of ∼45 km s −1 , and turns over to become approximately constant, giving a total SMC mass of 3 × 10 9 M , with the disc mass being 1.5 × 10 9 M . The central velocity dispersion of the disc was chosen to be near the current H I velocity dispersion of the observed SMC (Fig. 8), and the dark matter halo velocity dispersion is similar to the value of ∼25 km s −1 obtained from observations of the stellar halo carbon stars and planetary nebulae (Dopita et al. 1985;Hardy, Suntzeff & Azzopardi 1989;Hatzidimitriou et al. 1997). The Toomre Q-parameter at the disc half-mass radius is Q = 1.4. Our fiducial SMC model assumes a somewhat different scalelength and total extent for the SMC disc when compared with earlier studies (GSF94; GN96; YN03). We should stress that one cannot simply take as initial conditions, observed parameters of the SMC at the present day, since the SMC has (obviously) been significantly disturbed by interactions with the MW and LMC; we instead use initial parameters consistent with other relatively isolated dwarf disc galaxies (and integrate forward to ensure the final characteristics is consistent with that observed today). Dwarfs with H I mass 1.3-2.0 × 10 9 M (comparable to the SMC) have H I disc scalelengths ranging from 1.6 to 4.4 kpc, and smaller optical (i.e. stellar) scalelengths of 0.9-2.7 kpc (Swaters et al. 2002). These same galaxies have H I discs whose face-on corrected H I density drops to 1 M pc −2 at a radius between 10 and 12 kpc.
The reasons for the larger truncation radius compared to previous work are twofold. First, we suggest the tidal radius of the SMC disc is larger than the 5 kpc proposed by GN96. Several Gyr ago, the MW was likely to have been slightly less massive (secular halo growth), the SMC slightly more massive (less tidal disruption), and the perigalactic distance of the SMC was larger. In the estimated orbit of the SMC, at a look-back time of 2.5 Gyr (Fig. 2)

Disc
Halo 2.47 × 10 9 3.5 7.0 0.35 1.1 35 −5.72 × 10 3 29 1.0 0.1 0.96 distance r p was 60 kpc, and the apogalactic distance r a was 120 kpc. However, the MW mass enclosed within the larger SMC orbit has then increased. The tidal radius of the SMC is (Faber & Lin 1983) where the eccentricity e = (1 − r 2 p /r 2 a ) 0.5 , and thus r t = 6.3 kpc when r p = 60 kpc. We also suggest that it is not unreasonable that the disc was initially somewhat larger than the tidal radius when the SMC started to experience tidal stripping.
We performed stability tests on all initial SMC models, in the absence of any external potential (the equilibrium run), to ensure the initial models were indeed in equilibrium. Using the GCD+ parallel tree N-body code described by Kawata & Gibson (2003), we encountered only a minimal degree of disc heating and newly introduced spiral structure (generally with two symmetric arms), after 2.5 Gyr. The SMC models proper were then evolved by the same code for 2.5 Gyr further (the interaction run).

Interaction simulations
We adopt an SMC-centric non-inertial (but non-rotating) coordinate system for our interaction simulations, one in which the SMC disc lies on the xy plane. The orbits of the LMC and MW from Section 2.1 are translated and rotated to this coordinate system, and the SMC model constructed in Section 2.2 evolved within this new coordinate system under the influence of the now 'orbiting' MW and LMC. There are two degrees of freedom for the current inclination angle of the SMC, both currently unknown. We define θ and φ following fig. 1 of GN96, survey the full range of both, and find that the final features of the simulated MS are sensitive to this angle. We adopt in the fiducial model an angle of (θ , φ) = (45 • , 210 • ), as discussed in Section 4. This choice leads to a trailing tidal stream with an orientation consistent with the observed MS, and a leading arm with shape qualitatively similar to the observed LAF. This angle is mildly different from that adopted in GN96, (θ, φ) = (45 • , 230 • ), but we found that the new value leads to a better match to the HIPASS data set (data which were not available to GN96).
The particles in the SMC are assumed to be collisionless and their dynamical evolution calculated with GCD+. The acceleration applied to the ith particle is described (GN96) bÿ where is the softening length and m j is the mass of the jth particle. The position of r is measured in the SMC-centric coordinate frame. The first term is the self-gravity of the SMC particles; the second and third terms are the forces on the particle resulting from the Galaxy and LMC, and can be derived from the respective potentials in equations (1) and (3): The final two terms arise from needing to correct for the integration of the equations of motion in a non-inertial reference frame centred on the SMC. We use 200 000 disc and 200 000 halo particles to describe the SMC. This corresponds to a resolution ∼30 times greater than that employed by GN96, G99 and YN03. Such high-resolution allows us to examine features of the MS, LAF and SMC, in a manner not previously possible, since smaller fractional differences in particle density become statistically significant. We adopt approximately equal masses for the halo and disc particles of 7.6 × 10 3 and 7.4 × 10 3 M , respectively, and employ softening lengths d,h of ∼65 pc.

F I D U C I A L M O D E L
In this section, we show the results of our fiducial model, the basic parameters for which are listed in Section 2 and Table 3. Since we are primarily interested in the formation of the MS and LAF, which we assume are both stripped H I gas components from the SMC disc, our analysis focuses on the particles which were initially within the original SMC disc. In what follows, our simulation products are compared closely with the empirical HIPASS data set.

Evolution
In Fig. 3, we present a series of snapshots in the Galactocentric coordinate system, with the orbit of the MCs overlayed. Consistent with earlier models (e.g. GN96; YN03), when the SMC experienced a close encounter with the MW and the SMC ∼ 1.5 Gyr ago, the edge of the disc of the SMC began to be drawn out, which formed the tidal features that later became the LAF and MS. By T = −1.0 Gyr, at the subsequent apo-Galacticon, the LAF becomes more prominent, whilst the MS was still under development. By T = −0.3 Gyr, much of the initially stripped material in the leading tidal arm had been pulled back into the intercloud region (ICR), and the material still in the LAF was brought within 3 kpc of the solar circle. By this time, the MS and LAF morphology resemble that seen today. The next encounter with the MW and SMC at T = −0.2 Gyr caused little obvious consequence to either the MS or LAF. It did, however, cause the dispersion of much of the material that had been within the SMC. Much of this material either ended up in the ICR (although most of the ICR material was already in an ICR structure before this event), or contributed to the large spread in radial extent of the SMC. At the current time, the ICR extends radially from ∼30 to ∼80 kpc, and the SMC from ∼45 to ∼60 kpc.  (1), (2) --0.05-0.65 × 10 9 (9), (10), (11) Initial H I mass (M ) 1.5 × 10 9 1.1-2.0 × 10 9 0.5-1.8 × 10 9 a Stellar Mass within SMC (M ) --0.58-1.8 × 10 9 (12), (13) Total Mass within SMC (M ) --0.9-2.4 × 10 9 (14), (15), (16), (17) Initial total mass (M ) 3 × 10 9 3.0-3.6 × 10 9 1.1-3.6 × 10 9 b Initial radius of H I disc (r d,95 ; kpc) 7 4-7 10-12 c  We also ran the simulation 'forward' in time to T = +0.25 Gyr and found that the ICR undergoes mass-loss with material being dragged out into two tidal tails, separated radially and kinematically from the main MS and LAF features. This is consistent with what is expected by Brüns et al. (2005). By T = +0.25 Gyr, the material in the ICR has been dragged 55 kpc towards the MS in the plane of the sky, giving an angular separation of 35 • between the SMC and tip of the ICR.

H I column density distribution
In this section, we compare the present-day H I column density distribution between the empirical HIPASS data set (Putman et al. 1998) and our simulation results. Fig. 4 displays the H I column density map on a zenith equal area (ZEA) projection. Here, we arbitrarily define the regions corresponding to the MS, LAF, ICR and SMC as shown in the right-hand panel of Fig. 4. We converted the observed 21-cm flux to H I column density, after Barnes et al. (2001). To construct H I column density maps from the simulation results, we assume that the disc particles in the SMC are purely gaseous, and the H I mass fraction is 0.76. The column densities within the SMC and ICR region will be somewhat overestimated, as we neglect currently any associated stellar and ionized components. On the other hand, Fig. 5 demonstrates that most particles in the MS and the LAF at T = 0 originate in the outer edge of the initial SMC disc, where there was likely a lack of stars. To date, searches for a putative stellar component to the MS have proven unsuccessful (e.g. Philip 1976a,b;Recillas-Cruz 1982;Brueck & Hawkins 1983). Regardless, the H I fraction of 0.76 remains technically an upper limit, and thus our predictions should considered as upper limits. For these reasons the comparisons between the simulation and the HIPASS data focus mainly upon the properties of the MS and LAF. Fig. 4 shows that in our fiducial model the gross features of the observed MS are reproduced, and the LAF appears as a consequence of tidal interactions. This confirms previous studies, such as GN96 and YN03, which suggested that the MS and LAF features originate from gas stripped from the SMC disc by the tidal interaction with the LMC and MW.
The left-hand panel of Fig. 7, itself a new representation of the HIPASS data set, demonstrates that the observed leading arm extends above the Galactic plane to latitude b ∼ 30 • . While the full extent of the LAF above the plane remains a matter of debate (e.g. Brüns et al. 2005), the metallicity of its gas is consistent with an SMC origin (Lu et al. 1998;Gibson et al. 2000), supporting the tidal disruption scenario for the MS (and LAF). Furthermore, the observed H I cloud distribution seems to show 'a kink' near (l, b) = (310 • , 0 • ), where two further components of the LAF are seen north of the Galactic plane (labelled LAF II and III in Brüns et al. 2005). Although the exact position of the kink is inconsistent with the empirical data, our simulation does naturally predict its existence.
The tail of the observed MS shows spatial bifurcation near (l, b) = (300 • , −70 • ) and (l, b) = (80 • , −55 • ) in Fig. 4, with the two components forming an apparent twisting double helix-like structure (Putman et al. 2003). This bifurcation is not apparent in previous studies such as GN96 and YN03; higher-resolution simulations enable us to study such subtle features. Since this might be further evidence of the tidal interaction between the LMC and SMC, we return to this issue later.
An advantage of the present work is that the H I data available to constrain the models are significantly improved beyond that of GN96 or YN03. As in those previous studies, while the gross features of the observed MS and LAF are reproduced by our fiducial model, there remain subtle discrepancies between simulations and data. While the simulated MS is both broader and more extended than that observed (and hence the mean column density is somewhat lower than that encountered in the HIPASS data set), the derived H I gas mass (assuming a heliocentric distance of 57 kpc) from HIPASS is within 50 per cent of that of our simulated fiducial model (3.5 × 10 8 and 2.4 × 10 8 M , respectively). Fig. 6 Fig. 4). This demonstrates that the distance to the simulated MS is not constant, instead increasing across its length. Hence negating the distance ambiguity by obtaining a total flux gives us a fairer comparison than a total mass. Doing so results in an observed MS total H I flux of 4.8 × 10 5 Jy km s −1 , while the fiducial simulated MS has a total H I flux of 2.3 × 10 5 Jy km s −1 , still within a factor of ∼2 of the HIPASS data.
Conversely, the simulated LAF has a predicted H I mass and associated flux (7.3 × 10 7 M and 7.9 × 10 4 Jy km s −1 , respectively) both factors of ∼2 greater than that inferred from the HIPASS data set (3.5 × 10 7 M and 4.6 × 10 4 Jy km s −1 , respectively). 2 One clear difference between observation and simulation is that of the geometry of the LAF, in particular that of the projected deflection angle between the LAF and a Great Circle aligned with the MS proper (Fig. 7). In addition, the simulated LAF extends above the Galactic plane beyond that observed. We will discuss possible solutions to these apparent problems in Section 5. Fig. 8 shows the first and second moment maps for both our fiducial simulation and that derived from the HIPASS data set. In order to remove the large velocity gradient along the MS (which acts to obscure fine kinematical details within the maps), the first moment map is shown as the distribution of the velocity of v GSR − v MS , where v MS is the mean trend of velocity across the observed MS and LAF, and is defined in terms of two fitted Fourier components, v MS = 86 sin(θ MC + 2. • 5) − 92 sin(2(θ MC − 26. • 3)). The second moment map is the distribution of the velocity dispersion of H I gas.  In the first moment map, the observed velocity trend along our simulated MS is shown to be globally consistent with that of the empirical data. Fig. 9 displays v LSR against the Magellanic longitude, θ MC , and demonstrates perhaps more clearly that the mean velocity of the simulated MS is consistent with that observed. In contrast, the line-of-sight velocity of the simulated LAF is significantly larger than that observed, although it does follow the general trend of decreasing velocity at θ MC 0.

H I kinematics
The second moment maps seen in the lower panels in Fig. 8 indicate that the velocity dispersion of the simulated LAF is roughly consistent with, although slightly higher, than that observed. However, the velocity dispersion of the simulated MS is somewhat greater than that inferred from the HIPASS data set. This might be due to the neglect of gas dissipation via radiative cooling in the current suite of simulations, as will be examined further in Paper II.
We end by drawing attention to the evidence of a bifurcation in v LSR (right-hand panel of Fig. 9) within the MS and LAF. The next section discusses this bifurcation in more detail.

Spatial and velocity bifurcation
Our fiducial model shows a bifurcation in the MS, as seen in the ZEA spatial distribution of the H I column density of Fig. 4. Within the simulation, this bifurcation occurs both radially and tangentially. Fig. 4 displays the tangential bifurcation, and shows that there are two stream components that appear to follow a twisting topology governed by the orbits of the MCs -both the orbits and the   (Figs 3 and  6). At the tail of the MS, there is a third component, well separated from the other two streams -extending from 170 to 220 kpc. It is visible only through its heliocentric distance from the Sun, and coincides with the (l, b) position of the tip of the main streams. Also apparent in Figs 6 and 7, are the radial and tangential bifurcations of the LAF, respectively.
The right-hand panel of Fig. 9 demonstrates that the bifurcation of the MS appears also in v LSR . Plotted in Fig. 9 are lines showing the history (and future) of the MCs v LSR . Only one of the bifurcated components follows the orbit of the clouds (primarily the LMC), while the other possesses a higher velocity. Interestingly, there is a second velocity component at the position of the SMC (as well as at the head of the MS and the LAF). We remind the reader that Mathewson & Ford (1984) observed two velocity components within the SMC itself, with a separation of ∼50 km s −1 . While the separation of our two velocity components is much larger, it might indicate that the two observed velocity components are caused by a similar process.
Snapshots of the simulation, similar to those seen in Fig. 3, hint at the origin of these bifurcations. Prior to the first major peri-Galacticon at T = −1.5 Gyr, an encounter with the LMC at T = −2.2 Gyr drew the particles from the tip of the SMC disc closest to the LMC (and furthest from the MW), which resulted in the particles that eventually came to reside in the most distant MS component, at a Galactocentric distance of 170-220 kpc. The MS particles did not become 'distinct' from the SMC proper until T = −1.5 Gyr; at T = −1.05 Gyr, the MS then received an impulse from an LMC encounter, which caused the spatial bifurcation of the MS. The MS was then given a 'kick' by the LMC at the subsequent apo-Galacticon at T = −0.55 Gyr. This encounter at T = −0.55 Gyr resulted in the MS being broken into two kinematic components, resulting in the apparent velocity bifurcation.
On the other hand, a portion of the LAF comes from the particles on the same side of the disc, but closer to the SMC centre at the T = −2.2 Gyr encounter. The opposite side of the edge of the disc mostly consisted of particles that end up at the current time in the ICR. At T = −0.9 Gyr, the LMC passed through the LAF, splitting it into two bifurcated radial components.
Our models seem to create naturally both spatial and velocity bifurcations, via tidal interaction between the LMC and SMC. If this is the case, any observed bifurcation features would be strong evidence supporting the tidal formation scenario for the MS. Observationally, the left-hand panel of Fig. 4 shows that there is spatial bifurcation in the observed H I distribution. However, from this data alone, it is difficult to ascertain from the v LSR distribution of Fig. 9 whether there is an observed velocity bifurcation. Nevertheless, both Fig. 9 of this work, and Brüns et al. (2005), with data of higher velocity resolution, show evidence for a bifurcation in v LSR along the MS, with the two components observed in this work being separated by approximately 100 km s −1 at θ MC ∼ 100 • , and two components being visible in fig. 3 of Brüns et al. (2005) between the interface region and Galactic plane. There is also a hint of bifurcation (with separation of ∼100 km s −1 ) in the LAF at θ MC ∼ 45 • . Unfortunately, the current quality of the observational data is not enough to lead to a firm conclusion. Higher quality data cubes with improved velocity resolution will provide critical information concerning this apparent bifurcated structure.

PA R A M E T E R D E P E N D E N C E S
As explained in Section 2, our simulations involve several parameters which are not well constrained by current observations. In this section, we briefly demonstrate how the final configurations of the MS and LAF are sensitive to these free parameters. We varied parameters over a wide range of parameter space in our survey, to obtain our fiducial model. The parameter survey was performed at a lower resolution (25 000 disc and 25 000 halo particles) than the final fiducial model. We have confirmed that in the fiducial model the results of the lower-resolution simulation are consistent with the higher-resolution simulation, although, as presented above, the choice of a higher resolution enables us to discuss more detailed features. For example, in Fig. 10, we show the H I column density map for the lower-resolution model with identical parameters as our fiducial model. The distribution of the H I column density is roughly consistent with that shown in Fig. 4, although the higher-resolution models affords an improved examination of the finer-scale structures intrinsic to the simulated streams. The bifurcation alluded to in the previous section is an example of such a feature, and is perhaps best appreciated through an inspection of Fig. 11, where the flux 'excess' of the high-resolution model with respect to the lowerresolution model is presented in the v LSR versus θ MC plane, and the bifurcation becomes readily apparent. Most of the parameter space surveyed is summarized in Table 3. The scaleheight of the SMC disc, the ratio of SMC disc mass (H I mass in Table 3) to SMC halo (DM) mass, the velocity dispersion of the SMC disc, and the initial and final total mass of the SMC (varying the mass of the SMC for the purposes of the orbit calculation) were found not to be important to the evolution of the system. Both the velocity dispersion, and the ratio of H I to DM mass of the disc, simply scale the velocity dispersion and quantity of H I found in the final streams and clouds in a linear fashion.
In the parameter survey, the orbits of the MW, LMC and SMC were derived assuming the mass of the SMC is constant. However, we found that the SMC bound mass decreases from 3 × 10 9 to 1.5 × 10 9 M approximately linearly with time between T = −2.25 Gyr (near the first interaction between the LMC and SMC), and T = 0. Such mass-loss may affect the orbit and the final features of the MS and LAF (e.g. Zhao 2004;Knebe et al. 2005). Thus, we explore the effects of decreasing the mass of the SMC linearly with time from 3.0 × 10 9 to 1.5 × 10 9 M . We found that there is little change to the orbits, since the LMC, rather than the less massive SMC, primarily determines the orbit of both bodies in equations (4) and (5). As a result, we confirmed that the final features of the MS and the LAF are also not influenced heavily by the time evolution of the SMC mass. Thus, in all other simulations in the parameter survey, we use the orbits predicted with no evolution of the SMC mass. We now highlight the influence of the most important input parameters.

The initial scalelength of the SMC disc
We found that the final H I distributions are sensitive to the scalelength of the initial SMC disc. Fig. 12 shows the column density map of a model with a smaller scalelength (1.4 kpc) for the initial SMC disc. This reduced scalelength results in a lower total H I flux for the MS, as the initially more concentrated SMC mass distribution results in less material being stripped from the disc. Since the fiducial model does not have a total H I flux high enough to match perfectly the observations (Section 3.2), we conclude that reducing the SMC scalelength is not appropriate. It is also worth noting that Figure 12. H I column density map for models with a smaller scalelength (1.4 kpc) for the initial SMC disc. the smaller scalelength model leads to a less significant bifurcation in the simulated MS.
To quantify the difference between the models, we have measured the H I masses and total fluxes within the regions delineated as the MS in Fig. 4. As a result, the low-resolution simulation of the fiducial model has an H I mass in the MS of 2.4 × 10 8 M and total H I flux of 2.5 × 10 5 Jy km s −1 . The low-resolution simulation has a very similar mass to that of the high-resolution case (see Section 3.2). On the other hand, the small scalelength model leads to a MS H I mass and total flux of 1.8 × 10 8 M and 1.9 × 10 5 Jy km s −1 , respectively, which is significantly smaller than that of the fiducial model. Since the fiducial model has a mass which is somewhat lower than that inferred from the HIPASS data set, we again conclude that the models derived with the reduced initial scalelength for the SMC perform worse than the fiducial selection.
The H I survey of late-type dwarf galaxies by Swaters et al. (2002) suggested that the range of scalelengths of the gas disc for galaxies which have a similar H I mass to that of the SMC is 1.5-4.5 kpc. The scalelength of our fiducial model (3.5 kpc) appears reasonable, and thus the progenitor of the SMC is considered to have had a large H I disc before the SMC fell towards the MW.
We note that choosing a smaller scale radius is similar to setting a small truncation radius, in that the SMC material is concentrated more strongly towards the centre of the SMC, where it is more difficult to tidally strip. We obtain similar results to the above when we choose smaller SMC truncation radii, such that the stream is substantially retarded when the truncation radius is reduced from 7 kpc to 4 kpc. In this situation, the MS H I mass and total flux are reduced to 7.2 × 10 7 M and 7.0 × 10 4 Jy km s −1 . The LAF mass is scaled in a similar manner in both cases.

The inclination angle of the SMC
As mentioned in Section 2.3, the current inclination angle of the SMC is unknown, and the initial disc angle θ and φ are free parameters in our simulations. Fig. 13 shows the H I column density distribution for the models with (θ , φ) = (30 • , 210 • ) and (45 • , 230 • ), to demonstrate how these angles affect the final distribution of the MS and LAF. We remind the reader that we use (θ, φ) = (45 • , 210 • ) in the fiducial model, and therefore small differences of only 20 • in the initial inclination angle have quite a marked effect on the details of the final distribution. If accurate observations of the current incli- nation angle of the SMC were to be made, it would provide a strong constraint on any putative model of the formation of the MS.

The mass of the LMC
Another unknown parameter is that of the mass of the LMC itself. At the time of the GN96 study, the mass of the LMC was believed to be ∼2 × 10 10 M (e.g. Schommer et al. 1992). However, recently, some authors claim the lower mass of the LMC (e.g. van der Marel et al. 2002 who suggested that the mass of the LMC within R < 8.9 kpc is about 9 × 10 9 M ). Motivated by such claims, we also ran models with varying orbital parameters, and in particular, a lower LMC mass. Fig. 14 shows the result of a model with an LMC mass of 1.5 × 10 10 M . Even such small differences in the LMC mass cause a large change in the orbits of the LMC and SMC. Since the MS follows the past orbit of the MCs, the angle of the MS is radically different between the model with small LMC mass, and both the fiducial model and the observed MS. This result suggests that such a low-mass LMC is unlikely, if the MS is the result of tidal interactions.

S U M M A RY A N D D I S C U S S I O N S
We have carried out high-resolution N-body simulations of the history of the SMC disturbed by tidal interactions with the MW and LMC. We have surveyed most of the possible parameter space for the SMC orbit and the properties of the initial SMC disc, and found the best model. The increased numerical resolution of ∼7 Jy km s −1 per particle for the H I flux, ∼30 times higher than the previous studies (GN96; G99; YN03), made it possible for us to examine the detailed features of both leading and trailing tidal streams. Taking advantage of this higher resolution, we for the first time made a direct and quantitative comparison of the simulation results with newly available high-quality observational data from HIPASS. We convolved the HIPASS data set with the identical software tools used to analyse the simulated data sets, and compared the results in identical manners. Such comparisons confirm the conclusions of previous studies -that the existence of the LAF and MS can be explained by the leading and trailing streams of the SMC created by tidal interactions with the MW and LMC. However, our quantitative comparison have revealed extant problems with the models (some minor; some more significant). We found that even in our best model  (1) the shape of the MS is too extended in both width and length (Section 3.2); (2) the total H I flux of the MS is too low, and thus the mass of the MS is too low (Section 3.2); (3) the angle the LAF emanates from the MCs is not entirely consistent with observations (Section 3.2); (4) and nor is its total H I flux consistent with observations (Section 3.2); (5) the velocity dispersion of the MS is too high (Section 3.3) and (6) the line of sight velocity of the LAF is too high (Section 3.3).
These problems suggest that additional physics may be required to explain the observed properties of the MS and LAF, quantitatively. Obvious physics which this study excludes are gas physics, such as hydrodynamics and dissipation by radiative cooling. Simulations in Maddison et al. (2002) indicate that dissipation causes the MS to become narrower, which might lead to a reduced velocity dispersion. Ram pressure (Moore & Davis 1994;Mastropietro et al. 2005) is another physical process which our numerical model does not currently take into account. Ram pressure (or drag -see G99) is expected to shorten the leading arm and increase the gas density in the MS, which should help to solve the deficiencies of the model at recreating the LAF and MS with their correct shapes and densities. If it were to also bring the MS significantly closer than the 57 kpc fixed for the observational data within this paper, then the mass of the observed stream and total flux of the modelled stream would be increased, bringing them back into agreement. However, drag might lead to more problems for the length of the MS. Another possibly important mechanism is supernovae feedback which might aid in ejecting gas from the SMC and/or LMC, and help to increase the gas density in the MS and LAF. Finally, we note that the LAF passed very close to the centre of the MW 0.2 Gyr ago, and since we modelled the MW potential as spherically symmetric with a constant rotational velocity of 220 km s −1 , any deviations from this (such as the unknown contribution from the MW disc; e.g. Fich & Tremaine 1991) would affect our modelled LAF in particular. We are in the midst of introducing these physical processes into our simulations, and will report upon their respective effects in Paper II.
Another benefit of the higher-resolution simulations performed here is the identification of the bifurcation of the MS and LAF both spatially and kinematically. Our simulations predict that if the MS is created by a tidal interaction with the LMC, the bifurcation would appear both in the H I column density map and the line-of-sight velocity versus Magellanic longitude plane. Current observations are consistent with the existence of a spatial bifurcation in the H I column density map (left-hand panel of Fig. 4). In the velocity versus Magellanic longitude plane, it is difficult to make a firm conclusion from the current HIPASS data, although it is interesting that the SMC itself is found to consist of at least two velocity components (Mathewson & Ford 1984), perhaps caused by the same tidal disruption processes that form the bifurcation in our models. Observations with high velocity resolution and sensitivity may be able to test our prediction of bifurcation within the MS. If confirmed, it would provide strong evidence that the MS was created by tidal interactions.

AC K N OW L E D G M E N T S
We would like to thank Masafumi Noguchi who kindly provided his code for the orbit calculation, Mary Putman for the carefully flagged MS portion of the HIPASS HVC cubes, and the anonymous referee for many useful suggestions to improve the manuscript. TWC would like to thank Jeremy Bailin, Chris Power, Chris Thom, Virginia Kilborn, Lister Staveley-Smith and Geraint Lewis for helpful discussions, Stuart Gill for his PVIEW visualization tool source code during development, and Christian Brüns for an advance copy of his paper before publication. The support of the Australian Research Council and the Australian and Victorian Partnership for Advanced Computing, the latter through its Expertise Grant Programme, is gratefully acknowledged. DK acknowledges the financial support of the JSPS, through Postdoctoral Fellowship for research abroad. In this work we make use of data from HIPASS, which is operated by the Australia Telescope National Facility, CSIRO.