The Evolution of Magellanic-like Galaxy Pairs and the Production of Magellanic Stream Analogs in Simulations with Tides, Ram Pressure, and Stellar Feedback

Though observations and measurements suggest that the MS is composed of gas extracted from the MCs, the nature of the actual process is not clear. Ram pressure stripping caused by the MCs' motion through the hot gaseous halo of the MW could explain the origin of the TA but could not explain the origin of the LA, which, as its name indicates, is a leading structure. Tidal interaction with the MW (and between the MCs) can explain the origin of the MS (Connors et al. 2006; Diaz & Bekki 2012; Venzmer et al. 2012; Yozin & Bekki 2014; Pardy et al. 2018). These tidal models could also explain the irregular morphologies of the MCs and the presence of a bridge between them (Růžička et al. 2010; Guglielmo et al. 2014; Mackey et al. 2018). They can also produce both the LA and TA with a sufficiently large density and extent, and with the correct orientation. All models that successfully reproduce these structures require at least one and possibly two passages of the MCs near the MW. D'Onghia & Fox (2016) point out that tidal models based on a single passage fail to reproduce the filamentary structure of the MS.

While models that only include tidal interaction can produce a trailing and a leading structure under the right conditions, they may fail to reproduce several observed properties. Measurements of velocity gradients along the stream indicate that gas extracted from the MCs slows down as it moves away from them, the opposite of what might be expected from tides alone (Venzmer et al. 2012). Also, the chemical abundances of the LA match those of the SMC, while the kinematics suggest that the LA originated from the LMC, something difficult to explain from tides alone (Fox et al. 2018; Tepper-García et al. 2019).

Consequently, attention has been shifted to models that include a gaseous halo and the effects of ram pressure (Hammer et al. 2015; Tepper-García et al. 2015, 2019; Wang et al. 2019). Ram pressure could either reshape the MS once the gas has been extracted from the MCs by tidal forces, as suggested by Venzmer et al. (2012), or it could contribute to the actual extraction of the gas. These models are significantly more complex than purely tidal models, because they depend heavily on the detailed structure of the MW halo, and not merely on the MW mass. The strength of ram pressure, and its effectiveness in stripping gas from dwarf galaxies (DGs), depends on the orbital velocities of the dwarfs and the density of the gaseous halo inside which they are moving.

The orbital velocities depend on the total mass of the host galaxy, which can be determined fairly accurately through dynamical measurements. The density, or total mass, of the gaseous halo can be estimated by measuring absorption in spectra of distant QSOs, and the results have been a matter of debate. Prochaska et al. (2017) conducted a survey of 13 nearby L_* galaxies (the "COS-Halos" sample) and found cumulative halo gas masses of order  ∼8 × 10¹⁰ M_⊙ at a radius of 100 kpc and  ∼ 10¹⁰ M_⊙ within 50 kpc. This is considerably higher than an earlier analysis of lines of sight through the MW alone (Miller & Bregman 2013), which found a halo gas mass of 3.8 × 10¹⁰ M_⊙ at a radius of 200 kpc. A subsequent analysis (Bregman et al. 2018) of the COS-Halos sample, the "Stock-Bowen" samples, and additional MW lines of sight agrees with the lower gas mass for L_* and MW halos, although with a flatter profile than that found by Miller & Bregman (2013).

While tides, ram pressure, or a combination of both can possibly strip gas located in the outer regions of the MCs, these processes might not be powerful enough to extract gas located deep inside the potential wells of the MCs. Salem et al. (2015) found, in their simulations, that ram pressure stripping could account for only a few percent of the mass of the MS. By using a denser MW halo, the simulations of Tepper-García et al. (2015) extracted more gas from the MCs, but this gas then dissolves into the halo. Also, if ram pressure is sufficiently strong to form the TA, it might affect the morphology of the LA or even prevent its formation. Essentially, in all such models, ram pressure must be strong enough to extract the gas from the MCs, yet not so strong as to destroy the MS once the gas has been extracted, and these two objectives seem difficult to reconcile.

If ram pressure is kept at a sufficiently low level to allow the formation of the LA and to prevent the evaporation of the MS into the halo, it might then need some "assistance," to help extract the gas from the MCs in the first place. This could result from internal processes within the MCs, or from their mutual interaction. Besla et al. (2010, 2012) have shown that the influence of the MCs on each other may be more effective than the MW's influence. Here the one-armed spiral and warped stellar bar of the LMC are produced by a close encounter with the SMC, and the MS is produced by the LMC tidally stripping the SMC. Diaz & Bekki (2012) similarly found that recent close encounters between the MCs could play a critical role by ejecting the gas that forms the MS.

Also, several authors have suggested that galactic outflows resulting from stellar formation and feedback could push the gas out of the potential wells. This can produce a "Magellanic Corona" (Lucchini et al. 2020), which can then be stripped by tides or ram pressure (Fox et al. 2018; Pardy et al. 2018). Simulations have shown that supernova (SN) driven outflows can be quite efficient in removing gas from DGs (Murakami & Babul 1999), and this process is the likely explanation for the stellar-mass/halo-mass relation at low masses (e.g., Rodríguez-Puebla et al. 2015). With the inclusion of galactic outflows in the models, the ram pressure need not be as high, and this might help reproduce the morphology and kinematics of the LA and prevent the dissolution of the MS into the halo.

In this present study, we are not attempting to reproduce the Magellanic system specifically. Instead, we focus on the broader issue of the evolution of Magellanic-like systems composed of two interacting DGs orbiting in the vicinity of a massive host galaxy. Dwarf companions are commonly found in numerical simulations of hierarchical structure formation in a ΛCDM universe (Sales et al. 2013; Besla et al. 2018; Kallivayalil et al. 2018). Though few double-dwarf systems are known at present, some have indeed been observed. Besla et al. (2018) studied the distribution of low-redshift (z < 0.0252) dwarfs in the Sloan Digital Sky Survey and found that roughly 1 out of 25 dwarfs has a companion of comparable mass. They limited their study to dwarfs that are isolated from more massive galaxies, but while the galaxy mass distribution may be different in richer environments (e.g., Faltenbacher et al. 2010), we should still expect to find dwarf pairs.

Our objective is to understand the role played by star formation and feedback, galactic outflows, tides, and ram pressure in the evolution of Magellanic-like systems, as well as how the interplay between these various processes determines the morphology, kinematics, and chemical properties of the resulting structures. In an earlier paper (Williamson & Martel 2018, hereafter WM18), we presented a study of the evolution of a single DG orbiting within the gaseous halo of a massive host galaxy. We found that the combined effects of tides and ram pressure shape the morphology and determine the metallicity of the outflows but have little effect on the final metallicity of the DG itself. In this paper, we extend this early work to the case of two interacting DGs orbiting within the halo of a common massive host galaxy. This is the highest-resolution numerical study of the evolution of Magellanic-like systems that includes the combined effects of galactic outflows, tides, and ram pressure. The remainder of this paper is organized as follows: In Section 2, we describe the numerical algorithm and the initial conditions used for the simulations. Results are presented in Section 3, and their implications are discussed in Section 4. Our conclusions are presented in Section 5.

Our simulations are performed using the simulation code unmodified from WM18, but with the initial conditions altered to include two DGs, and with other run-time parameters modified. Fuller descriptions of the algorithms are given in our previous papers (Williamson et al. 2016a, 2016b; WM18), but here we provide a brief summary of the code, directly quoting from WM18.

We use a version of the GCD+ smoothed particle hydrodynamics (SPH) code (Kawata & Gibson 2003; Barnes et al. 2012; Kawata et al. 2013, 2014). This code includes a stochastic star formation model that relaxes the single stellar population assumption, allowing different star particles to represent stars of different masses. Star particles return energy and metals to the interstellar medium through SNe and stellar winds. We do not include stellar photoionization and radiation pressure and may somewhat overpredict the star formation rates (SFRs) for our simulated galaxies (Hopkins et al. 2014).

Chemical reactions (e.g., formation of H₂) are not explicitly calculated during the simulation (unlike, e.g., Yozin & Bekki 2014), but cooling rates are pretabulated with CLOUDY using a full chemical network. The cooling rates depend on the metallicity, density, and temperature of the gas.

Smoothing lengths are calculated dynamically through an iterative method, so that each particle has  ≈58 neighboring particles. We apply a minimum smoothing length of 2 pc, which means that particles in very dense regions have  >58 neighbors. We set the Plummer-equivalent force softening length to 2 pc.

The metal content of particles is tracked throughout the simulation, and a subgrid diffusion model allows metals to spread between particles, tracking eight separate species. Our version of GCD+ includes modified algorithms for metal deposition and diffusion, as described in Williamson et al. (2016a), and a dynamic background potential to represent the varying tidal forces on a satellite galaxy moving through a host galaxy potential, as described in Williamson et al. (2016b).

The center of mass of the DG system is stationary in our simulation frame, and the host potential is applied as a tidal force, with the position of the host halo directly integrated. To avoid the computational expense of modeling the entire host gaseous halo, we treat it using a "moving box" method, where halo gas particles are generated or deleted on the boundaries of a cubic box centered on the center of mass of the DG system. This allows halo gas particles to have the same mass resolution as the DG particles, preventing a source of numerical error. The temperature and density of particles entering the box through the leading edge are calculated from an analytic profile based on the location of the dwarf system. Gas particles located on the side edges are frozen to maintain (rough) hydrostatic equilibrium in the host halo. This method is described in further detail in WM18, and a similar method was introduced in Nichols et al. (2015). Here we have increased the box width from 160 to 320 kpc to capture both dwarfs and the gas stream.

2.2.1. Host Model

We use the same MW-like host galaxy model and parameters as Run B of WM18. The host galaxy consists of a dark matter halo of mass 10¹² M_⊙ with an NFW profile (Navarro et al. 1997) with a concentration parameter c = 12. It also contains a gaseous component with a β-profile, as used in Salem et al. (2015), with the temperature profile given in Makino et al. (1998). In our simulations, the DGs are never close to the core of the host galaxy, and so we approximate the β-profile using the large r limit as

$$\begin{eqnarray}&&n(r)={n}_{0}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-3\beta }.\end{eqnarray} \tag{ 1 }$$

The parameters for the MW β-profile are taken from Miller & Bregman (2013), which is on the lower end of MW gas halo densities. In WM18 we performed simulations with both an MW gas halo and a gas halo with densities scaled down by a factor of 10. In both models, both an LA and a TA were produced from a single Magellanic-like dwarf. To reduce computational expense while maintaining consistent baryon particle masses between both the dwarfs and the MW halo, most of the simulations in this work use the downscaled halo model, setting n₀ = 0.046 cm⁻³, r_c = 0.35 kpc, and β = 0.71. We do, however, run a single test (Run G) at the higher mass, with n₀ = 0.46 cm⁻³, but a shorter simulation time. We do not change the particle mass for Run G—the halo is modeled with more gas particles than in runs A–F and is therefore more computationally expensive. We mostly confine our analysis to the downscaled simulations, Runs A–F, and use Run G as a convergence check. In all these runs, the initial halo metallicity is 10⁻² Z_⊙. We also perform simulations (Runs H–K) without a host gas halo, to further constrain the effects of ram pressure. Two simulations (Runs J and K) are also performed without the host gravity field, as a basis to demonstrate the effects of the host galaxy tides.

2.2.2. DG Models

2.2.3. Orbits

Proper-motion studies suggest that the MCs have had a close encounter with each other in the past (Zivick et al. 2018) and may be approaching the MW for the first time (Kallivayalil et al. 2006a, 2006b; Besla et al. 2007; Piatek et al. 2008). To determine MC-like initial conditions for the DGs' orbit around the host halo, we take the LMC and SMC position and velocities of Zivick et al. (2018) and Kallivayalil et al. (2013) and integrate this motion forward and backward by 5000 Myr, as shown in Figure 1. The exact details of this motion will depend on the MW potential and MC proper and radial velocities, but here we are interested in indicative initial conditions and not a precise model of the exact history of the MCs. We use these as the basis for two sets of dwarf orbits around the host.

Zoom In Zoom Out Reset image size

We base one set of ICs on the predicted position and velocity of the MCs about 1500–2000 Myr in the past, where the dwarfs are initially  ∼220 kpc from the MW center. We slightly tweak these initial conditions from our predicted past coordinates, to increase the pericenter of the DGs' orbit. Our predictions showed a pericenter of  <10 kpc, which would heavily disrupt the system and not produce MC-like galaxies. This setup is used for Runs A, B, and H. We plot the initial coordinates for Run A in Figure 1 for comparison with our predicted MC past orbits.

The second set of ICs is a circular orbit with a radius of  ∼110 kpc, which allows us to investigate the effects of tides and ram pressure over an extended period of time, in an idealized environment. In all simulations, the host galaxy is at rest at the origin. This setup is used for Runs C–G and I. We plot the initial coordinates for Run D in Figure 1 for comparison with our predicted MC past orbits.

As our DGs do not have a mass ratio similar to that of the MCs, we are not able to directly take our dwarf–dwarf orbital properties from the MCs. Instead, we produced a series of low-resolution test runs to explore which dwarf–dwarf orbits produced MC-like interactions. For simplicity, the orbits of the dwarfs around each other are coplanar with their orbit around the host galaxy, which we define as the x–y plane, i.e., initially v_z = 0 and z = 0 for all galaxies. Both dwarfs rotate in the same direction, inclined 90° from the orbital planes. We vary the dwarf–dwarf relative tangential velocity and relative position in these test runs, while keeping their initial separation constant at Δr = 20 kpc and initial relative radial velocities at zero. We discarded runs where the dwarfs immediately merged, or quickly escaped from each other without interaction. It appears that the specific present-day configuration of the MCs, as a pair of strongly interacting DGs, under strong influence of a host galaxy, may be a rare situation. As the strong interactions appear to be short-lived before the dwarfs escape each other or merge, some fine-tuning is required to ensure that the timing is correct. From our test runs, we found that the resulting "interesting" simulations had dwarf–dwarf relative velocities of  ±30 km s⁻¹ and  ± 50 km s⁻¹, where a positive or negative velocity indicates that the dwarf's orbit (the orbit of the dwarfs around each other) is prograde or retrograde, respectively, relative to the dwarfs' orbit around the host galaxy. We selected these initial conditions to rerun the simulations at full resolution.

The initial coordinates of both galaxies in all full-resolution runs are summarized in Table 1, as well as the pericenter radius reached by each DG in each simulation. Column (1) gives the name of each run. Columns (2)–(5) give the initial position and velocity of the center of mass of the dwarf system. Columns (6)–(9) give the position and velocities of the dwarfs relative to each other. Columns (10) and (11) give the pericenters for each dwarf. Column (12) gives the gas density n₀ used for the density profile of the host galaxy gas halo (Equation (1)). The pericenters are directly obtained from the simulations and therefore include any hydrodynamic effects on the galaxy orbits.

Table 1.  Parameters of the Simulations

Run	xCOM,i	yCOM,i	vx,COM,i	vy,COM,i	Δxi	Δyi	Δvx,i	Δvy,i	rper,1	rper,2	n0
	(kpc)	(kpc)	(km s−1)	(km s−1)	(kpc)	(kpc)	(km s−1)	(km s−1)	(kpc)	(kpc)	(cm−3)
A	0	−220	−60	120	0	20	50	0	61	49	0.046
B	0	−230	−85	120	0	20	50	0	89	82	0.046
C	0	−110	−190	0	0	20	−30	0	96	96	0.046
D	0	−110	−190	0	0	20	−50	0	87	79	0.046
E	−100	−10	0	190	0	20	30	0	100	97	0.046
F	−100	−10	0	190	0	20	50	0	102	64	0.046
G	0	−110	−190	0	0	20	−50	0	86	79	0.460
H	0	−230	−85	120	0	20	50	0	89	82	0
I	0	−110	−190	0	0	20	−50	0	87	79	0
J	−	−	−	−	0	20	50	0	−	−	0
K	−	−	−	−	0	20	−50	0	−	−	0

Note. Note that z = v_z = 0 for all galaxies initially. Subscripts i indicate initial values. Pericenters are directly measured from the simulation. Runs J and K do not show COM distances from host halo center, as these are performed in isolation.

Download table as:  ASCIITypeset image

To summarize, Runs A, B, and H represent eccentric orbits around the host galaxy, with the same relative velocity between the dwarfs. In all three runs, the dwarf orbits are prograde with respect to their orbit around the host galaxy and have the same speed—we found that these eccentric orbits were particularly sensitive to parameter variations that would cause the dwarfs to either merge or escape each other almost immediately. Runs C–F and I represent five runs with circular orbits around the host, with the dwarf's orbit being either prograde or retrograde, and with two different relative velocities between the dwarfs. As stated in Section 2.2.1, Run G is similar to Run D but with a higher host halo gas density, while the mass resolution of Run G is the same as in the other runs—the run has more particles and is therefore more computationally expensive. Runs H–K are performed without a host gas halo, and Runs J and K are performed without the host gravity field—i.e., these are simulations of isolated dwarf pairs. Run H is equivalent to Run B without a gas halo. Run I is essentially Run D without a gas halo. Run J has the same dwarf–dwarf configuration as Runs A, B, F, and H, but in isolation. Similarly, Run K has the same dwarf–dwarf configuration as Runs D, G, and I.

We define the center of each DG as the local minimum of the gravitational potential, using the pykdgrav Python library to recalculate the potential. As we note in Section 3.3, due to disturbances from interactions and strong feedback, this center may not correspond to other centers of interest, such as the metallicity slope center.

A visualization of the evolution of the system for Run D is shown in Figure 4. Results for Runs B–E are similar. In both dwarfs, star formation proceeds rapidly with strong feedback, producing SN bubbles at t = 500 Myr, as seen in the first four panels of the top row. The stellar distribution is irregular at early times as new stars are formed stochastically from the feedback-disturbed gas (see Kawata et al. 2014), though a hint of spiral arm structure may be visible. The feedback expels gas, which forms a halo around each dwarf, as seen in the rightmost panel. Tidal forces between the dwarfs drive some of this gas into a bridge between the dwarfs. At this early stage, many of the characteristic properties of the MC system are already formed.

Zoom In Zoom Out Reset image size

By t = 1000 Myr, tidal interaction between the dwarfs becomes more dominant, and tidal arms can be seen in the stellar distribution. The outflows have also been shaped by the tides of the host into a stream. By t = 1500 Myr, the gas density inside the dwarfs has been reduced, due to expulsion of gas and consumption by star formation. Fewer SN bubbles are visible. Tidal stirring also causes the dwarf disks to thicken. Finally, by t = 2000 Myr, the dwarfs have merged. This remnant is a boxy elliptical galaxy, with long gas streams.

A visualization of the evolution of the system for Run A is shown in Figure 5. As noted above, Runs A and F do not result in a merger. Here we still see anisotropy in the stellar distribution from stochastic star formation, but tidal tails from the interaction between the dwarfs are less apparent. The dwarfs do still form an extended gas stream, and a gas bridge between them. After their close encounter, both dwarfs settle back into fairly regular disks, but inclined relative to their initial orientation. This inclination is caused by the tidal field of the host galaxy. For further comparison, we show in Figure 6 the final stage of the evolution, at t = 2000 Myr, of a single-dwarf system. This is "Run B" from WM18. Again, the gas and stars are disturbed primarily by feedback, which produces an outflow that is shaped by the tidal field of the host. The dwarf itself settles into a fairly regular disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These results suggest that most of the irregularity in the stellar distributions of the dwarfs is caused by the dwarfs themselves—through feedback and mutual interactions—and not by the host galaxy. The thickening of the dwarf disks is also caused by interactions between the dwarfs, and eventually by their merger. This is in agreement with the classic picture of the morphological origin of dwarf ellipticals from tidal stirring (Moore et al. 1996; Mayer et al. 2001; Smith et al. 2010). The host galaxy has a gentler effect on the dwarf morphologies and only dominates the shape of the outflowing stream.

The top row of Figure 7 and the bottom row of Figure 3 show the SFR for both dwarfs (blue and orange lines) and the sum of the two (green lines) for Runs A–F and H–K. In Runs B–E, the dwarfs eventually merge, and from that point only the green line is plotted. The SFR rises to a peak at t ∼ 1200 Myr, driven by internal gravitational collapse and by tidal stirring. Star formation then slowly shuts down as the gas is consumed or expelled.

Zoom In Zoom Out Reset image size

The peak of star formation occurs at roughly the same time in all runs, even in the absence of an external tidal field, suggesting that internal processes dominate star formation—star formation increases as instabilities form, and then slows down as the fuel has been consumed. The single-dwarf models presented in WM18 all had a similar star formation peak at t ∼ 1200 Myr, as do Runs J and K in the absence of an external tidal field. Dwarf–dwarf interactions have a significant secondary effect. Runs B–E, which have close encounters that end in a merger, have a higher peak of star formation than Runs A and F. In particular, Run D shows two star formation peaks, triggered by two close encounters separated by  ∼1000 Myr. In these runs, the final merger results in a starburst, during which a significant fraction of the remaining gas is converted to stars. Ram pressure and host tides seem to have little effect on star formation. The two dwarfs in Run F move apart as they follow different trajectories, but they do not have significantly different SFRs. Runs A and B have similar orbits around the host galaxy but show significantly different SFRs owing to interactions between the dwarfs. Runs J and K, in the absence of a tidal field, show similar (low) SFRs to Run F, as the dwarfs interact gently over a long period of time rather than having strong star formation triggering encounters. Overall, the picture is that star formation is primarily driven by internal processes, secondarily affected by close encounters with the other DG, and not noticeably affected by the host galaxy, except in how the host galaxy affects the dwarf–dwarf orbits.

We track the masses of stars and of total and "warm/cool" gas for Runs A–F in the second row of Figure 7. Here "warm/cool" gas is defined as gas with temperatures below 2 × 10⁴ K, equivalent to the "neutral gas" definition in Wang et al. (2019). We divide the gas and stars into the DG population, defined as gas and stars within 10 pc of the center of either dwarf, and "expelled" gas and stars as gas and stars located outside that radius, but that originated from the dwarfs, not from the host galaxy. The mass of ejected stars is quite low, and this (orange) line is not visible in many plots. The (pink) warm/gas line is also not often visible, because almost all of the expelled gas is in this warm/cool phase and the pink and brown lines almost always coincide in the figure—the vast majority of hot (≥2 × 10⁴ K) gas outside of the dwarfs is that which originated from the host galaxy halo.

The strongest effect on the gas mass of the dwarfs is star formation. The gas in the dwarfs is consumed by star formation and converted into stars. Feedback pushes out only a fraction of the warm/cool gas into the stream, while the rest is consumed by star formation. Further gas is expelled when galaxies have a close encounter or merge, both due to the disruption of the encounter, and due to the resulting burst of star formation. This is more clearly shown in the third row of Figure 7, where we plot the fraction of gas within the dwarfs and expelled from them. Generally, the expelled gas is comparable to or greater than the retained gas by the end of the simulations, consistent with observations of the MS (Fox et al. 2014; Barger et al. 2017; Richter et al. 2017). Although only a fraction of the initial gas mass is expelled, the consumption of gas by star formation causes the expelled and retained masses to be similar. This suggests that it does not require extremely effective stripping or outflows for the gas mass of the MS to be similar to that of the MCs—modest outflows plus gas depletion by ongoing star formation can achieve this instead.

The final ratio between the expelled gas and retained gas masses can be considered an approximate proxy for the strength of interactions between the dwarfs. In Run F, the dwarfs only pass by each other early on, the star formation peak is low, and most of the gas is retained in the dwarfs. In Run A, the galaxies have a longer and closer encounter, and the retained gas fraction is smaller. In Runs B, D, and E, the dwarfs undergo at least one close encounter before merging, and most of the remaining gas ends up in the stream. Of course, this is not the only factor—ongoing star formation will also continue to expel and consume gas. Despite the strong interaction and sharp star formation peak in Run C, the final expelled gas fraction is lower than that of the similar Runs D and E. This is because Runs D and E have stronger ongoing star formation in the final 1000 Myr of the simulation.

We can further illustrate the relative importance of the final starburst and instabilities triggered by close encounters by comparing Run E with Runs B–D. In Run E, the encounter is not as "close" as in Runs B–D, with a minimum separation of 4.5 kpc. The resulting effect is weaker, and as a result, the amount of gas remaining before the final merger takes place is larger, resulting in a much stronger starburst. The peak SFR in Run E is larger than in Run C, and more than twice as large as in Runs B and D, resulting in a sharp increase in stellar mass at t = 1000 Myr. We note that the competing effects of instabilities and mergers conspire to produce a comparable stellar mass by the end for the simulations.

Both star formation and interactions between the dwarfs are of similar importance for driving outflows, while host tides and ram pressure have only a small effect (beyond the host perturbing the dwarf orbits and modulating their interactions somewhat). Dwarf–dwarf interactions are effectively counted twice here—the interactions can directly expel material through tidal forces, and indirectly by triggering bursts of star formation.

The "effective yield" of star formation is plotted in the fourth row of Figure 7. We define this as

$$\begin{eqnarray}&&{\eta }_{\mathrm{eff}}(t)=\displaystyle \frac{{Z}_{g}(t)-{Z}_{g}(0)}{{M}_{\mathrm{SF}}(t)},\end{eqnarray} \tag{ 2 }$$

where Z_g(t) is the total mass of metals (of all tracked species) in the gas phase within a 10 kpc distance from the center of each galaxy, M_SF(t) is the total mass of stars formed since the beginning of the simulation, and η_eff(t) is the effective yield. This value measures the production of metals from star formation, as well as the ejection of metals through feedback. Note that the very large yields at very early times are not representative, but are caused by the initial distribution of ages and masses of stars—stars present in the initial conditions at t = 0 produce metal enrichment before any new stars have formed.

After about 500 Myr, the yield in all simulations reaches a steady state of η_eff ∼ 0.02. Here, star formation is ongoing, but outflows are not yet strong enough to carry many metals out of the galaxy. At this point, the rate of increase of metal mass is simply proportional to the SFR. The effective yield therefore becomes close to the "true" yield we might expect in a closed box, which observations have shown to be ∼10⁻² (Garnett 2002; Tremonti et al. 2004; see also the sample yields for GCD+ in Figure 2 of Kawata & Gibson 2003).

Rapid star formation then produces a metal-rich outflow at t ∼ 1000 Myr, causing the effective yield to drop. This effect is much more significant in Runs B–E, where mergers trigger strong bursts of star formation and feedback. After this, continual star formation causes the effective yield to slowly recover back toward the steady state. In Runs A and F, the weak outflows caused by weak star formation cause only a small drop in effective yield at this time, and the effective yield remains closer to the true yield.

Figure 8 shows the disk-plane radial profiles of the abundance ratios [Fe/H] and [O/Fe] of each dwarf, or their merger remnant, at various epochs. Note that in the runs resulting in mergers (B, C, D, and E), the system is not fully relaxed by the end of the simulation. These profiles are centered on the minimum of the potential well, which does not exactly correspond with the location of the metallicity peak, resulting in the little "glitches" seen in the central parts of the profiles at late times.

Zoom In Zoom Out Reset image size

In all runs, star formation is more intense in the center of the dwarfs than in the outer regions, resulting in more metal production in the center and negative [Fe/H] slopes. As star formation continues, metallicities continue to increase in the center. The final peak metallicity in the dwarf centers generally relates to the strength of star formation. Runs B–D produce a higher peak metallicity than Runs A and F. This star formation is triggered by mergers, and the runs without mergers show smaller high-metallicity cores. The peak metallicities in the final snapshot are quite high, which may relate to our feedback method, as we discuss below.

The metallicity slopes at t = 2500 Myr are quite steep in Runs A and F but are much shallower in the merger remnants of Runs B–E. Here, the merger has disturbed and mixed the gas, producing a flatter metallicity slope—although metals are still somewhat centrally concentrated. Although we do include metal diffusion in the simulations, the dominant effect is the direct advection of metal-rich gas through flows driven by feedback and tidal stirring. Some star formation occurs in outer regions, but the main form of enrichment beyond the central few parsecs is from bulk transport of metal from the inner region.

The evolution of the [O/Fe] profiles relative to the [Fe/H] profiles is primarily determined by the time delay between the onset of Type II SNe (which produce oxygen and other α-elements) and Type Ia SNe (which produce most of the iron). The values of [O/Fe] are initially high because Type II SNe explode almost immediately, the lifetime of their progenitors being very short on the scale of the simulations. Eventually, Type Ia SNe explode, resulting in an increase in iron abundance and a corresponding drop in [O/Fe]. Generally, we see a positive slope develop, due to ongoing star formation producing α-enriched gas in the dwarf centers. The observed MCs appear to have flatter slopes, and our models may be underestimate mixing from feedback and turbulence (see Section 4).

3.4.1. Production and Evolution

Figures 9 and 10 show the evolution of the gas column density for Runs A–F and Runs H–K, respectively. We use a two-variable color scheme to illustrate the origin of the gas, with Σ₁ and Σ₂ being the column density of gas extracted from the first and second dwarf, respectively. The contribution of gas from the host halo to the stream is very small and is not plotted here.

Zoom In Zoom Out Reset image size

The outflows from the two dwarfs produce halos that are initially unmixed, except in the region where the halos directly overlap—the halos quickly become larger than the interdwarf distance. As the dwarfs orbit around the host galaxy, their gas halos are stretched by tidal forces into streams. In the absence of a host galaxy (Runs J and K, Figure 10), the outflows produce a largely spherical halo instead. Both arms of the streams tend to contain gas from both dwarfs. Although largely well mixed, typically one tail of the stream consists of somewhat more gas from one dwarf than the other. The dominant dwarf for each tail is simply the dwarf that is on that side of the system—as our dwarfs have the same mass, no one dwarf is preferentially stripping the other. This effect is particularly strong where mixing is weak—in Runs A and B, which have elliptical orbits, and in Run F, where the dwarfs escape. Some mixing does occur—the tails in Runs A and B are less mixed at t = 1000 Myr than at t = 2000 Myr—but the difference is still visible by the end of the simulation, where we can clearly see which arm came from which dwarf.

Zoom In Zoom Out Reset image size

The difference in the level of mixing between the simulations is due to how rapidly the dwarfs orbit each other. In all runs, the dwarfs are initially separated from each other by the same distance, and with two different relative speeds (30 and 50 km s⁻¹). However, the dwarf–dwarf orbit can be modified by the tidal influence of the host galaxy, with a strength depending on the orbit of the dwarfs through the halo, and this influence tends to dominate over small variations in the relative velocities of the dwarfs. In Runs C–E, the clouds orbit each other rapidly, producing a mixed outflow stream (see also Figure 2). In Runs A and B, the tidal force of the host increases the period of the dwarf–dwarf orbits, and there is less mixing. Similarly, in Run F, the relative orientation of the dwarf–dwarf orbit allows tidal forces to separate the dwarfs from each other, and mixing is reduced.

All simulations produce two steams stretching in opposite directions from the dwarfs, a leading one and a trailing one, similar to the LA and TA found in the Magellanic system. However, these labels are somewhat arbitrary and ambiguous. For example, in Run A, an LA and a TA are formed, but the dwarfs overtake the LA in their orbit. This causes the LA to become a TA. Meanwhile, the TA swings around to become an LA. The two arms effectively swap places. To avoid ambiguity, we define the LAs and TAs according to the velocities and morphology of the system at the snapshot of interest—i.e., the "LA" is the arm currently closest to the dwarfs' direction of motion, regardless of its historical origin.

In Run F, where the dwarfs end up moving apart at late times, we note the persistence of the bridge between them, which stretches over a distance of 200 kpc. This further demonstrates how tidal forces dominate over ram pressure in these simulations—the bridge is stretched by tides and is not significantly stripped by ram pressure.

3.4.2. Effects of Ram Pressure on Evolution

Ram pressure does have some effect, particularly in Runs A and B, where the dwarfs plunge the deepest into the host halo. We should expect ram pressure to push gas in the opposite direction of the dwarfs' motion (yellow arrows in Figures 9 and 10). By t = 1500 Myr and t = 2000 Myr, the LAs of the streams do appear to be bent backward. This confinement of the LA is also visible (but less dominant) in some others runs. Runs H and I, which lack any ram pressure, show similar stream structures to Runs B and D, although the streams of Runs H and I are more extended and diffuse. Here, the effect of the host gas halo is not to produce streams through ram pressure but to constrain the outflows into narrow streams and to somewhat limit the distance they propagate to. Interestingly, in some runs there are perhaps indications of a two-tailed TA, where one tail is produced by tides and the other is produced by ram pressure. This is clearest in Run D at t = 1500 Myr, where the fork of the TA is directly opposite the dwarfs' direction of motion and appears to be influenced by ram pressure.

We performed one additional run, Run G, at a higher host halo gas density, closer to that of the MW as measured by Bregman et al. (2018). This run uses the same orbits as Run D, the only difference being the mass of the host galaxy (see Table 1). A visualization of the evolution of the system for this run is shown in Figure 11. To limit computational cost, we stopped this simulation at t = 1700 Myr.

Zoom In Zoom Out Reset image size

Comparing Figures 4 and 11, there are no systematic differences in the internal evolution of the dwarfs. Both show similar disturbances, with only stochastic differences. We also find that SFRs are very similar.

The key differences are in the large-scale evolution of the streams. Here we again find that pressure confinement has a stronger effect than ram pressure stripping. An LA and a TA are still both present and largely point in the same direction as in Run D. Ram pressure stripping has a small effect, in that these arms are somewhat curved "backward" at the ends. But the more significant effect is that the arms do not extend quite as far beyond the dwarfs. Increasing the gas density by a factor of 10 has approximately halved the extent of the arms. This is caused by the pressure of the gas in the halo resisting the outflows and tidal forces that stretch out the MS.

We emphasize that pressure confinement constrains both the LA and TA. This differs from ram pressure stripping, which enhances the TA but constrains (or prevents) any LA.

3.4.3. Properties

Figure 12 shows the time evolution of the gas metallicity, for Runs A–F. The upper limit of the scale corresponds to solar abundance. As the two dwarfs are initially identical, we should not expect that metallicity can uniquely trace a stream to its source. Instead, the streams generally have low metallicities. They largely trace their origin to the gas in the low-metallicity outer regions of the dwarfs (Section 3.3) that has been entrained by outflows. As the streams take time to propagate, and the DGs are continuing to produce metals through star formation, the streams are sampling the lower metallicity of an earlier era in evolution of the DGs. However, a strong burst of star formation can also produce a high-metallicity outflow that adds an enriched component to the stream. The resulting streams are a combination of low-metallicity and high-metallicity gas, which does not exactly match the metallicity of either DG. This highlights the difficulty in attempting to trace the origin of the MS by comparing abundances to those of the MCs.

Zoom In Zoom Out Reset image size

The metallicity gradients of the simulated streams are also not uniform, even in sign. These gradients depend on the nature of the outflow. When the outflow is gentle (see Run D and Run E at t = 2000 Myr), the metallicity gradient is negative, as enriched gas mixes with the metal-poor medium as it flows. When the outflow is dramatic (see Run B at t = 2000 Myr), the metallicity gradient is positive. Here, an enriched front plows through the medium, leaving behind a lower-density wake (see Figures 9 and 10) that mixes with the metal-poor medium. However, in general, our simulated MS shows lower metallicities than the DGs, consistent with observations (Fox et al. 2018).

We also examine the kinematics of these systems. We calculate the line-of-sight velocity (v_r) and position in Magellanic coordinates (Nidever et al. 2008), for all gas particles that originated in the dwarfs, at t = 1500 Myr. We define our Magellanic coordinates by setting L_MS = 0° and B_MS = 0° at the center of mass of the dwarf system, where L_MS is the azimuthal position along the orbital plane as seen from the halo center and B_MS is defined along the orthogonal great circle. We plot the L_MS − v_r distribution in Figure 13. Here we have reordered our runs so that models with similar orbits can be more clearly compared. Naïvely, we might expect that tides and ram pressure would produce opposite kinematic signatures—tides would show acceleration away from the dwarfs, while ram pressure would show deceleration toward the host galaxy velocity in the TA but no LA at all. In these simulations that combine tides and ram pressure, we often find multiple kinematic components, with the exact structure and gradient depending on the orbital configuration.

Zoom In Zoom Out Reset image size

Runs A, B, and H have strong radial velocity gradients. Here, the dwarfs have passed periapsis with the host galaxy and are on nearly radial orbits out of the system. Run A has a stronger velocity gradient than Run B because the small difference in initial velocities means that the dwarfs in Run A have traveled farther on their orbit and therefore are on a more radial trajectory. Run H has the same orbit as Run B but excludes ram pressure, providing a basis to demonstrate the effects of ram pressure. In this case, tides alone (Run H) produce a simple linear velocity gradient along the streams and through the merged DG. Ram pressure (in Run B) has little effect on the TA but truncates the LA and brings its velocity down to that of the halo, inverting the velocity gradient seen from tides alone.

The remaining runs have nearly circular orbits, and therefore small radial velocities and radial velocity gradients. Runs D, G, and I have the same initial orbital configuration, but Run D has moderate halo gas density, Run G has strong halo gas density, and Run I has no halo gas. The difference between Runs D and I is similar to the difference between Runs B and H—ram pressure causes the simple linear velocity gradient in the tides-only run to be flattened in the LA and causes the LA to be somewhat truncated. Run G shows truncation of both the LA and TA. Here, pressure confinement is again clearly dominating over ram pressure stripping.

Runs C, E, and F show a clear S shape in their velocities. The inner gradient is produced by the rotation of the dwarf–dwarf system, and the outer gradient by tides and host gas pressure (as in Runs D, G, and I). In these three runs, the dwarfs have retrograde orbits at their merger or closest encounter, and the gas expelled during the interaction has an inverted velocity gradient. While the dwarfs in Run C initially have a prograde orbit, the influence of the host tides causes the dwarfs to merge in a retrograde sense. This can be seen in the top row of Figure 2. The velocity gradient from the dwarf–dwarf rotation is less visible in Runs D, G, and I, where the opposite process happens. There, the dwarf orbits are initially in a retrograde orbit, but the host tides cause the dwarfs to merge and eject gas in a prograde fashion.

We can compare our kinematics plot with the position–velocity diagram of Nidever et al. (2010), which we have also plotted in Figure 13. Here the rotation of the Magellanic system is visible, as is a long linear velocity gradient in the TA that flattens at large distances, which we have shown is consistent with a tidal origin for the TA, somewhat truncated by host halo pressure. The LA is short and has a flat velocity gradient. From our simulations, we interpret this as a tidally stripped tail, decelerated by halo pressure.

We plot the gas column density in Magellanic coordinates centered on the center of mass of the dwarfs, as seen from the center of the host, in Figure 14. This is the total surface density of all gas species but can be compared qualitatively with the H i maps of Putman et al. (2003) and Nidever et al. (2010). We include the H i map of Nidever et al. (2010) in Figure 14, for comparison. In Runs A, B, and H, the apparent sizes of the streams are considerably shortened by projection effects, because the DGs are on nearly radial orbits and produce nearly radial streams. This foreshortening is particularly noticeable in Run A, where the arms are most radial. While we do produce extended streams, the projected size of these streams is particularly parameter dependent, and (as stated above) we do not attempt to fine-tune our models to produce the exact geometry of the Magellanic system.

Zoom In Zoom Out Reset image size

In general, the projected column density maps show that most runs produce LAs truncated by halo pressure and a somewhat more extended TA. Runs H and I, without ram pressure, show LAs and TAs of similar extent. This further supports that tidal stripping with pressure confinement can explain the general extent and morphology of the MS.