THE EVOLUTION OF DWARF GALAXY SATELLITES WITH DIFFERENT DARK MATTER DENSITY PROFILES IN THE ERISMOD SIMULATIONS. I. THE EARLY INFALLS

Dwarf galaxy satellites of the Milky Way (MW) and M31 are the main focus of "near-field cosmology," providing a plethora of useful constraints and tests for galaxy formation theories as well as for the nature of dark matter in the ΛCDM paradigm, yet their origin is still not completely understood. Nearly all the known satellites, located within the expected virial radius of the bright spirals, at distances smaller than 300 kpc, are classified as dwarf spheroidal galaxies (dSphs). They have fairly round stellar components with low or negligible rotation but relatively high velocity dispersions (Mateo 1998; McConnachie 2012). At the same time, nearly all of the dwarf galaxies located outside such radii and far from M31 are gas-rich dwarf irregular galaxies (dIrrs), which have ongoing or recent star formation (SF), flattened disky shapes, and exhibit substantial rotational support (v_rot/σ_⋆ > 1). The different spatial distribution of dwarf galaxy types is known as the morphology–density relation and reflects an analogous effect seen in galaxy clusters and rich groups on larger mass scales (Dressler 1980).

Environmental mechanisms have been often invoked to transform dIrrs into dSphs (Mateo 1998; McConnachie 2012). One of the mechanisms proposed to achieve such transformation is tidal stirring, namely the cumulative effect of repeated tidal shocks occurring near the pericenter of the eccentric orbits of satellites (Mayer 2010).

The process has been investigated thoroughly with collisionless N-body simulations adopting high-resolution models of dwarf galaxies interacting with rigid or symmetric models of hosts on orbits with eccentricities as high as those of cosmological subhalos (Mayer et al. 2001a, 2001b; Klimentowski et al. 2007; Kazantzidis et al. 2011a). It has also been investigated in combination with ram pressure mass removal using hydrodynamical simulations of dwarf–main-host interaction in which a diffuse gaseous halo consistent with observational constraints is added to the dark halo of the main host (Mayer et al. 2006; Mayer 2010).

These studies established tidal stirring as a major mode of morphological evolution of dwarf galaxies, turning them from faint rotationally supported disks into even fainter pressure-supported spheroidals, thus providing a natural explanation for the morphology–density relation ubiquitously observed in galaxy groups, including our own Local Group (Grebel et al. 1999).

Direct tidal heating as well as indirect heating stemming from bar and buckling instabilities triggered by tidal shocks have been both shown to play a role in tidal stirring, although the latter requires initial disks with mass-to-light ratios (M/L) at the upper end of those observed in dwarf galaxies (Mayer et al. 2001b; Kazantzidis et al. 2011a). These works also determined that ram pressure mass removal aided by tidal shocks can rapidly remove the gas content of satellites explaining why dSphs are devoid of gas and why they can have very high M/L.

Recently the effect of the dark matter density profile of dwarf galaxies on tidal stirring has been studied, motivated by the new findings of cosmological hydrodynamical simulations of field dwarf galaxy formation which have shown that feedback from supernovae and massive stars can turn the cuspy halos predicted by ΛCDM into halos with shallower inner density slopes (γ in [0.3, 0.7] range), depending on the star formation rate achieved in the galaxy (Governato et al. 2010, 2012; Pontzen & Governato 2012; Di Cintio et al. 2013; Shen et al. 2014; Trujillo-Gomez et al. 2015). It is important to note that there are also simulations that have cast doubts on the feedback-driven core formation mechanism (Oman et al. 2015). Surely the magnitude of the effect and its impact on galaxy and halo properties depend on the specific implementation of feedback (Oñorbe et al. 2015a).

There have also been attempts to show the existence of this effect specifically for dwarf galaxy satellites, albeit in this case, resolution is a potential issue (Brooks et al. 2013). Łokas et al. (2012) and Kazantzidis et al. (2013) have thus updated the previous tidal stirring simulations by carrying out N-body simulations in which disky dwarf galaxies have shallow dark matter profiles like those arising in the cosmological hydro-simulations, and have checked results against varying resolution. They have found that the transformation by tidal stirring is enhanced significantly when a shallower density profile is adopted (γ = 0.2 or 0.6). They have also found that since mass loss is enhanced as well, remnants of dwarf galaxies with shallow dark matter profiles can evolve into very faint objects such as ultra-faint dwarfs (UFDs; Łokas et al. 2012). The enhanced mass loss and transformation follow from the lower binding energy and more impulsive dynamical response of the system in the case of objects with a shallower halo profile.

However, these recent simulations, exactly like those of the last decade, are fairly idealized since they completely lack cosmological context. Hence, orbits are also idealized, being disconnected with the epoch of fall inside the host and with the shape of its gravitational potential, which is just assumed to be spherically symmetric. Instead, cosmological simulations show that the host is triaxial and grew substantially in mass and radius over the last 10 billion years (Diemand et al. 2007), which implies satellites can be accreted on increasingly wider orbits over time. Furthermore, in idealized simulations, subhalos can only interact with the main host and not between them. Also, for simplicity, these simulations have always explored a fixed or limited mass range for the satellites, rather than considering the wide range of masses from those of actual subhalos. It is expected that the initial mass, and also the stellar mass relative to that of the halo, ought to have an important effect on the response of the system to tidal shocks. Furthermore, other studies have shown that mergers or strong interactions, while rare, can occur among dwarf galaxies just before infall, turning them into dwarf spheroidals very quickly (Kazantzidis et al. 2011b), but leaving it completely open to identify the distinct signature of a remnant of one such merger as opposed to one from tidal stirring. Ideally, all these difficulties should be overcome in a fully cosmological hydro-simulation of an MW-sized galaxy, but with current resolutions around 100 pc and ∼10⁴ M_⊙ in the baryons, numerical results cannot be put on firm ground for objects that should be only a few hundreds of parsecs in size.

Therefore, in order to make progress beyond the idealized simulation we devise a numerical strategy that in the past has only been attempted for simulations of galaxies in clusters (Dubinski 1999; Mastropietro et al. 2005). In this approach, known as the "replacement technique," cosmological subhalos are replaced with high-resolution models of dwarf galaxies within the original cosmological volume. Models are constructed based on observational and theoretical constraints and matched as closely as possible to the total virial mass of the parent subhalo. After the replacement the simulation is rerun to the present epoch. In this way galaxy satellites are automatically placed on realistic orbits and are subject to the full tidal interaction history of the parent subhalos. Due to the computational facilities and improved speed and parallelism of the codes, the simulations can span the equivalent of more than 10 Gyr of evolution with this technique, while in previous cluster simulations only a small timespan, from z = 0.5 onward, could be covered (Mastropietro et al. 2005). In this paper we present the first replacement simulation carried out for an MW-sized halo, using the ErisDARK cosmological simulation carried out in the concordance ΛCDM model (Kuhlen et al. 2013; Pillepich et al. 2014). The simulations are carried out with state-of-the-art massively parallel N-body codes, PKDGRAV2 (Stadel 2001) and ChaNGa (Jetley et al. 2008, 2010; Menon et al. 2015; http://www-hpcc.astro.washington.edu/tools/changa.html) on the CRAYXE6 "Monte Rosa" of the Swiss National Supercomputing Centre.

The paper continues with Section 2 where the replacement technique and numerical simulations are presented in detail. Next the analysis methods and results are enumerated in Section 3. Section 4 contains an extensive discussion regarding the implications of the results for the origin of dSphs, including their relation to observational data. Finally Section 5 includes the conclusion and caveats of the presented work.

2.1. Initial Conditions and Selection Procedure

The initial conditions of the simulations of this paper are based on the ErisDark run (Pillepich et al. 2014), which is a dark matter only cosmological gravity run of 40 million particles using the ΛCDM cosmological model (H₀ = 73 km s⁻¹ Mpc⁻¹, n_s = 0.96, σ₈ = 0.76, Ω_m = 0.268 and Ω_Λ = 0.732). The cosmological volume contains a 1 Mpc scale zoom-in region with three levels of refinement that lead to a minimum mass of the dark matter particle of 1.2 · 10⁵ M_⊙ and a minimum softening length of 124 pc (Pillepich et al. 2014). In the zoom-in region an isolated MW-sized halo emerges, with a mass consistent with the lower limit for the estimated halo mass of the MW, which provides a good fit with current kinematical data on the MW stellar halo (Rashkov et al. 2013). The companion hydrodynamical simulations, Eris and ErisLE, have been shown to produce realistic spiral galaxies by z = 0 that reproduce nearly all the main properties of the disk and bulge of the MW (Guedes et al. 2011, 2013; Mayer 2012; Bird et al. 2013).

Satellite galaxies falling into the main host after the last major merger, which takes place at z ∼ 3, are considered for replacement with high-resolution dwarf galaxy models. We analyze the entire available set of outputs of the original simulation with the use of the Amiga Halo Finder (AHF; Gill et al. 2004; Knollmann & Knebe 2009) with a virial overdensity criterion that varies with redshift according to the cosmology (Bryan & Norman 1998). In this paper we particularly analyze the earliest infalling population of satellites among which there are survivors at z = 0. In the next paper we will present larger populations of dwarf galaxies which are accreted at later epochs.

We thus proceed in the following way. First of all, among subhaloes with bound masses greater than 10⁸ M_⊙ at z = 0, the satellite with the earliest infall is identified. The adopted threshold mass ensures that we are considering well resolved objects with at least 10³ particles (presumably their progenitors are even better resolved at higher redshift, before tidal mass loss begins), to minimize numerical effects on their internal dynamical evolution (Moore et al. 1996). The respective object enters the host's virial radius just before z = 2, when it is still not relaxed, being the result of a merger between two halos. We trace the satellite back to z = 2.87 when the two merging halos are still well separated, namely their virial volumes do not overlap. Therefore the epoch corresponding to z = 2.87 is set as the epoch of replacement of the original ErisDark's infalling subhalos with their high-resolution versions. Next we identify other subhalos that at this epoch are located between 1 and 4 R_vir of the host. We restrict the sample to the most massive ones, having pre-infall masses around 10⁹ M_⊙ and above. Among these we reject 60% that cannot be properly identified and replaced because they are either too elongated (they have axis ratios less than 0.65, possibly tidally disturbed by companions) or have a substantial fraction of their mass in subhalos. Four objects remain after this rejection procedure, with masses in the range 0.99–4 × 10⁹  M_⊙. Note that at z ∼ 2 this translates into satellites having masses in the range 80–300 times smaller than the host, which has a virial mass of 3.0 × 10¹¹M_⊙ and a virial radius of 55.6 kpc at the corresponding epoch (both relatively high since ErisDark has a fast assembly history with no major mergers after z = 3; Pillepich et al. 2014). We remark that the lower end of the satellite-to-host mass ratio is still in a regime in which dynamical friction is effective in eroding the orbital energy over many Gyr (Colpi et al. 1999), while this is a negligible effect in published idealized simulations that assume satellites of similar mass but within a host having a fixed mass ∼10¹² M_⊙ (Kazantzidis et al. 2013). Finally, we observe that none of the four objects has a mass above 10⁸ M_⊙ at z = 0 in the original ErisDark run. Therefore they are not part of the largest subhalos at the final epoch, which can well exceed this mass at low z (Paper II, Pillepich et al. 2014). In total we thus replace six subhalos (Figure 1), of which two are merging during infall and the other four are accreted individually. The halos' properties are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  The Virial Quantities of the Replaced Halos and Their Distances from the Center of the Host

Satellite	Mvir (109M⊙)	Rvir (kpc)	dhost (kpc)
A	3.83	13.0	140
B	3.04	12.0	75
C	1.51	9.5	82
D	0.94	8.1	157
E	0.40	6.1	165
F	0.33	5.7	159

Download table as:  ASCIITypeset image

2.2. Model Galaxy Generation

The replacement models are generated using the GalactICS code (Kuijken & Dubinski 1995; Widrow & Dubinski 2005; Widrow et al. 2008), which generates multi-component equilibrium models of galaxies. Each model galaxy is composed of two components: the dark matter halo and the stellar disk. It is assumed by construction that the initial model has a stellar disk since this stems from the assumption that disky dwarfs are the progenitors of dSphs, the core idea of tidal stirring. The code creates halos with negligible angular momenta and spherically symmetric density profiles as described by the relation:

$$\begin{eqnarray}&&{\rho }_{\mathrm{halo}}(r)=\displaystyle \frac{{2}^{2-\gamma }{\sigma }_{{\rm{s}}}^{2}}{4\pi {{Gr}}_{{\rm{s}}}^{2}}\displaystyle \frac{1}{{\left(\displaystyle \frac{r}{{r}_{{\rm{s}}}}\right)}^{\gamma }{\left(1+\displaystyle \frac{r}{{r}_{{\rm{s}}}}\right)}^{3-\gamma }}\displaystyle \frac{1}{2}\mathrm{erfc}\left(\displaystyle \frac{r-{r}_{{\rm{c}}}}{\sqrt{2}{r}_{{\rm{w}}}}\right).\end{eqnarray} \tag{ 1 }$$

In the previous equation the scale radius and velocity scale are denoted with r_s and σ_s, and γ stands for the cusp coefficient. In order to restrain the model to finite dimensions, a complementary error function with cut-off radius r_c and width r_w is used. The velocity distribution is computed by solving the Eddington equation for the prescribed density profile.

For each removed halo, two replacements are created: one with γ = 1.0 and another with γ = 0.6. The cut-off radius and width are set to be linear functions of the virial radius of the original object:

$$\begin{eqnarray}&&{r}_{{\rm{c}}}={r}_{\mathrm{vir}}-{r}_{{\rm{w}}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{r}_{{\rm{w}}}=\displaystyle \frac{\sqrt{2}{r}_{\mathrm{vir}}}{10}.\end{eqnarray} \tag{ 3 }$$

The scale radius is determined indirectly by assuming that the concentration of an average halo of mass M_vir at redshift z is given by the empirical relations found by Klypin et al. (2011):

$$\begin{eqnarray}&&c({M}_{\mathrm{vir}},z)={c}_{0}(z){\left(\displaystyle \frac{{M}_{\mathrm{vir}}}{{10}^{12}{h}^{-1}{M}_{\odot }}\right)}^{-0.075}\left(1+{\left(\displaystyle \frac{{M}_{\mathrm{vir}}}{{M}_{0}(z)}\right)}^{0.26}\right).\end{eqnarray} \tag{ 4 }$$

With the concentration and virial radius set, the scale length can be determined.

$$\begin{eqnarray}&&{r}_{{\rm{s}}}=\displaystyle \frac{{r}_{\mathrm{vir}}}{c}.\end{eqnarray} \tag{ 5 }$$

The two variants corresponding to the shallow and steep dark matter density profiles have by construct the same scale radius. In order to have the same mass enclosed in the same virial volumes with different cusp coefficients, the scale velocity is adjusted accordingly. The values of the characteristic quantities associated with dark matter halos with virial quantities presented in Table 1 are listed in Table 2.

Table 2.  Scale Quantities as Well as Concentrations of the Different Halos and the Properties of Their Associated Stellar Disks

Satellite	A	B	C	D	E	F
c	6.0	6.1	6.3	6.4	6.8	6.8
rs (kpc)	2.18	2.00	1.53	1.27	0.91	0.85
σsγ=1.0 (km s−1)	88.4	82.0	65.2	55.9	42.4	40.0
${\sigma }_{{\rm{s}}}^{\gamma =0.6}$ (km s−1)	84.7	78.6	62.4	53.4	40.4	38.0
M⋆ (106M⊙)	4.50	3.57	1.75	1.08	0.44	0.37
Rd (pc)	783	697	492	388	253	230
hd (pc)	80	80	80	80	80	80

Download table as:  ASCIITypeset image

The previously mentioned code can generate a disk component of a galaxy according to the cylindrical spatial distribution:

$$\begin{eqnarray}&&{\rho }_{\star }(r,z)={\rho }_{0}\mathrm{exp}\left(-\displaystyle \frac{r}{{R}_{{\rm{d}}}}\right){\left[\mathrm{sech}\left(\displaystyle \frac{z}{{z}_{{\rm{c}}}}\right)\right]}^{2},\end{eqnarray} \tag{ 6 }$$

and cylindrical distribution of the radial dispersion:

$$\begin{eqnarray}&&{\sigma }_{{\rm{r}}}(r)={\sigma }_{r0}{\left(\mathrm{exp}\left(-\displaystyle \frac{r}{{R}_{{\rm{d}}}}\right)\right)}^{1/2}.\end{eqnarray} \tag{ 7 }$$

The reader who is interested in more details on how the code constructs the galactic model is referred to the publication Widrow et al. (2008).

There are free parameters that have to be set using educated guesses based on available constraints on galactic structure at low mass scales. Most importantly, the mass of the stellar disk M_d for each object is 1.2 × 10⁻³ times its virial mass. Such a ratio is chosen to reflect the stellar-to-halo mass ratio, M_⋆/M_vir, inferred from the abundance matching analysis extrapolated to lower masses (Behroozi et al. 2013). Past simulations of disk galaxy formation presented in the Kaufmann et al. (2007) paper associate with objects with maximum circular velocities in the 24–53 km s⁻¹ range and scale height values in the 0.08–0.18 kpc range. Considering the previously mentioned results the values of the scale heights of all the objects are set to the minimum of 0.08 kpc. The choice of constant scale height is motivated by the expectation that at the masses of interest the disk thickness acquires minimum threshold value. The value is dictated by the balance between the thermodynamics of disk formation, such as the effect of the cosmic ionizing flux providing a minimum temperature floor of ∼10⁴ K, corresponding to a bulk velocity dispersion of ∼8–10 km s⁻¹, and the gravitational potential of the halo. For all the stellar disks introduced in the simulation the scale height was set to 0.08 kpc.

With the scale height and stellar-mass–to–total-mass ratio assigned, the scale length remains to be determined. The square of the scale length is naturally constrained to be proportional to the virial mass of the respective object. An observational study provided by Swaters et al. (2009) finds several dwarf galaxies with disk scale lengths as low as 0.33 kpc. Considering the inferred galaxy masses for the latter sample and the values assumed in previous work on tidal stirring based on various constraints (Kazantzidis et al. 2011a), the reference value of R_d = 0.4 kpc was set for an object with a virial mass of 10⁹M_⊙. Using such a reference value as a normalization, the following relation is used for the calculation of the scale length as a function of galaxy mass:

$$\begin{eqnarray}&&{R}_{{\rm{d}}}=0.4{\left(\displaystyle \frac{{M}_{\mathrm{vir}}}{{10}^{9}{M}_{\odot }}\right)}^{0.5}\;\mathrm{kpc}.\end{eqnarray} \tag{ 8 }$$

For each pair of objects the parameters of the stellar component are listed in Table 2.

From the original cosmological simulation, quantities such as the virial mass and virial radius are measured (Table 1). The virial quantities are used to determine all the parameters required for the generation of the pairs of high-resolution galaxies that will be used as replacements. Namely, the following quantities are determined: the dark matter halo quantities (σ_s, r_s, r_c, r_w) and the stellar disk quantities (M_stars, R_d). The remaining free parameters are the number of dark matter particles and star particles. For all the models created, both numbers were set to one million.

2.3. Stability of the Galaxies

Each of the 12 high-resolution-generated objects evolve in isolation for 1 Gyr prior to their introduction in the cosmological environment. The state of the objects at the end of the isolation run can be seen in Figures 2 and 3. The purpose of these simulations is to verify the stability of the models and allow the system to remove instabilities prior to replacement. Initial discreteness noise is unavoidable, and it is known to seed waves in cold disklike systems which could lead to mass redistribution (Mayer 2011). These often manifest themselves as transient spiral waves, although the method adopted by our initial condition generator combined with the low self-gravity of the light disks embedded in a massive halo should minimize amplification of the waves. Furthermore, two-body relaxation between more massive dark matter particles and lighter stellar particles may also lead to some spurious evolution (Mayer et al. 2001a, 2001b). All isolation runs were performed using the pkdgrav2 code (Stadel 2001). The softening lengths of the particles were set according to the equation:

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{dark}}=0.06\cdot {\left(\displaystyle \frac{{M}_{\mathrm{vir}}}{{10}^{9}{M}_{\odot }}\right)}^{1/3}\mathrm{kpc},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{star}}=0.02\cdot {\left(\displaystyle \frac{{M}_{\mathrm{stars}}}{4\cdot {10}^{7}{M}_{\odot }}\right)}^{1/3}\mathrm{kpc},\end{eqnarray} \tag{ 10 }$$

and are listed in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Masses and Force Resolutions of Dark Matter and Stellar Particles

Satellite	Msp (M⊙)	Mdmp (M⊙)	sp (pc)	dmp (pc)
A	4.50	3830	9.6	93.8
B	3.57	3040	8.9	86.9
C	1.75	1510	7.1	68.8
D	1.08	940	6.0	58.8
E	0.44	400	4.5	44.2
F	0.37	330	4.2	41.6

Download table as:  ASCIITypeset image

The baryonic distributions of the constructed galaxy models with a cusp coefficient equal to unity are presented in Figures 2 and 3. The previously mentioned images correspond to the objects at the end of the isolation run and at the moment of replacement. Figure 4 shows the cylindrically averaged surface density at the beginning of the isolation run and at the end. Small fluctuations are visible as small amplitude waves, which cause departure from perfect axisymmetry. Despite the existence of the local density fluctuations for both versions there are no significant changes happening on scales comparable to the scale length during the isolation run. Moreover, the analysis of the global mass distribution including the dominant dark matter component as shown by the circular velocity curves in Figure 5 along with the radial distribution of the cumulative mass present in Figure 6, do not show any appreciable changes. We conclude that the stability of the models is satisfactory.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.4. Replacement Method and Numerical Simulations

The resulting raw simulation data are analyzed using the AHF in the same manner as with the original ErisDark simulation data. From the set of halos generated by the AHF code, the six pairs of satellites are identified through the indices of bound star particles. As output of the AHF code, we use in our satellite analysis only the list of bound dark matter and star particles. There are cases when the halo finder does not detect one of the objects at a particular epoch either due to the pericenter shock or temporary interference with other large subhalos. Such events will appear as unavoidable interruption of the flow of data points on the analysis plots. As soon as the object is identified again by the AHF code, the flow of data points resumes.

The center of the satellite is defined as being the position of the stellar density peak of the respective object. From this point, a 3D half-light radius (r_1/2) can be defined as the radius of the sphere centered on the previously mentioned density peak which includes and excludes half of the mass of bound stellar particles. All the values of the half-light radius associated with each object at the available evaluation epochs are presented in Figure 7. The next step in the preparation of the satellite data for analysis is the alignment of the object. This is performed by defining a set of particles that includes both dark matter and stars inside a sphere centered on the stellar density peak and with a radius equal to two half-light radii. All the bound particles of the object, including the ones exterior to the two half-light radii spheres, are rotated such that the moment of inertia tensor, corresponding to the set of stellar particles interior to the previously mentioned sphere, is diagonalized. This process is repeated for each individual object at each epoch of measurement.

Zoom In Zoom Out Reset image size

Once the satellite is identified and analyzed there are a number of diagnostics that are routinely measured, as described in the following subsections. In addition to dynamical masses, dark matter masses, and stellar masses (see Tables 4 and 5 for the results at the starting and final time), the key diagnostics used to quantify the morphology of the model galaxies are the shape as measured by the axis ratios, the ratio between 1D velocity dispersion σ_* and rotational velocity v_rot of the stars (note that we refer here to the actual stellar rotational velocity rather than to the circular velocity), and the surface density profiles of the stars. The characteristic rotational velocity (v_rot) is measured around the half-mass radius, which is defined for the whole 3D mass distribution. For the calculation of the velocity dispersion (σ_*), the entire bound stellar population is considered.

Table 4.  Bound Stellar and Total Masses Inside the 500 pc Radius Spheres Centered on the Density Peaks of the Satellite Galaxies with Steep (Columns 2 and 3) and Shallow (Columns 4 and 5) Density Profiles at z = 2.87

Satellite	M500 (107M⊙)	${M}_{500}^{\star }({10}^{5}{M}_{\odot })$	${M}_{500}({10}^{7}{M}_{\odot })$	${M}_{500}^{\star }({10}^{5}{M}_{\odot })$
A	8.75	6.26	4.32	6.09
B	7.98	5.91	3.99	5.82
C	6.08	4.88	3.37	4.80
D	4.96	4.06	2.91	4.09
E	3.44	2.67	2.21	2.69
F	3.17	2.40	2.06	2.40

Download table as:  ASCIITypeset image

Table 5.  Characteristic Masses and Mass Ratios at Last Epoch for the Galaxies with Original γ = 1.0 (Upper Rows) and γ = 0.6 (Lower Rows)

Satellite	A	B	C	D	EF
M500 (107M⊙)	0.64	3.62	3.21	3.07	4.44
${M}_{500}^{\star }({10}^{5}{M}_{\odot })$	0.25	3.49	3.23	3.37	3.50
${M}_{500}^{\mathrm{total}}/{M}_{500}^{\star }$	256	104	99	91	127
${M}_{\mathrm{bound}}^{\mathrm{total}}/{M}_{\mathrm{bound}}^{\star }$	238	102	112	128	370
M500 (107M⊙)	0.14	0.39	0.60	0.71	2.61
${M}_{500}^{\star }({10}^{5}{M}_{\odot })$	0.06	0.72	1.42	1.75	3.42
${M}_{500}^{\mathrm{total}}/{M}_{500}^{\star }$	234	54	42	41	76
${M}_{\mathrm{bound}}^{\mathrm{total}}/{M}_{\mathrm{bound}}^{\star }$	234	39	42	47	312

Download table as:  ASCIITypeset image

Our classification scheme to decide whether or not the final remnant is a dSph is more stringent than the one used in Kazantzidis et al. (2011a, 2013). They used two criteria simultaneously, namely, that the object had to have c/a > 0.5 and v_rot/σ_* < 1. With the latter criterion, outliers from the typical structural properties of dSphs, such as the Tucana dSph, which has v_rot/σ_* ∼ 1, would be included. Here instead, we use only the criterion v_rot/σ_* < 0.5, which is indeed the case for nearly all observed dSphs and which, as we will see, also automatically guarantees that c/a > 0.5 in our simulated remnants at z = 0. Finally, for brevity we will refer to "cored models" when we talk about models with inner slopes of the dark halo γ = 0.6 and to "cuspy models" when we discuss models with a Navarro–Frenk–White (NFW) profile.

3.1. Dimensions along the Principal Axes

With the object aligned and principal axes identified, three maps of the projected stellar mass corresponding to the three pairs of principal axes can be generated.  A square grid with size 4 r_1/2 × 4 r_1/2 and a million square cells is defined perpendicular to each of the vectors defining the axes directions. For each cell a projected mass is defined that equals the mass of all bound star particles whose projection on the grid is included in the respective cell. Using the information in the projected mass maps, the dimensions of the object along the principal axes can be determined. In the process of dimensions calculations, ellipses containing cells with different levels of mass density are being fitted. The distribution of these ellipses around the center of the image allows for the calculation of the ratio of the dimensions. For each of the three planes (x, y), (y, z), (z, x) the ratio between the shortest and longest dimensions, c and a, are measured: ${\left(\frac{c}{a}\right)}_{x,y},{\left(\frac{c}{a}\right)}_{y,z},{\left(\frac{c}{a}\right)}_{z,x}$.

After the previously presented analysis is performed the first available scalar indicator of the transformation of the object can be determined:

$$\begin{eqnarray}&&\displaystyle \frac{c}{a}=\mathrm{min}\left({\left(\displaystyle \frac{c}{a}\right)}_{x,y},{\left(\displaystyle \frac{c}{a}\right)}_{y,z},{\left(\displaystyle \frac{c}{a}\right)}_{z,x}\right).\end{eqnarray} \tag{ 11 }$$

The time evolution of the previously defined quantity is presented in Figure 8.

Zoom In Zoom Out Reset image size

3.2. Velocity Dispersion and Tangential Velocity

With the directions of the principal axes determined, the velocity dispersion is computed in the direction of each axis. For the calculation of the previously mentioned quantities the population of stars contains exhaustively all gravitational bound stellar particles. The velocity dispersion of the system is:

$$\begin{eqnarray}&&{\sigma }_{\star }={\left(\displaystyle \frac{{\sigma }_{x}^{2}+{\sigma }_{y}^{2}+{\sigma }_{z}^{2}}{3}\right)}^{1/2}.\end{eqnarray} \tag{ 12 }$$

All bound stellar particles contribute to the calculation of σ_⋆.

In order to compute the rotational velocity associated with the system, three cylindrical shell volumes are defined. Each of the volumes is defined with a central axis constrained to pass through the maximum density region of the object and to be parallel to a principal axis dimension. The inner and the outer radii of the volumes are 0.95 r_1/2 and 1.05 r_1/2.

For all the stars bound to the object and included in the previously defined volume, the component of their velocity tangential to any concentric circle centered on the axis of the cylinder is considered for the average. The resulting average is the velocity v_i,j, where (i, j) ∈ (x, y), (y, z), (z_x).

With the dispersion and rotation quantities determined, the ratio:

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{rot}}}{{\sigma }_{\star }}=\displaystyle \frac{\mathrm{max}({v}_{x,y},{v}_{y,z},{v}_{z,x})}{{\sigma }_{\star }}\end{eqnarray} \tag{ 13 }$$

can be used as another indicator of the degree of transformation from a rotationally supported satellite to a dispersion-supported galaxy (Table 6). In Figure 9 the decrease of the ratio due to the effect of the tides can be observed.

Zoom In Zoom Out Reset image size

Table 6.  Characteristic Velocities and Velocity Ratio at the Last Epoch for Galaxies with Original γ = 1.0 (Upper Rows) and γ = 0.6 (Lower Rows)

Satellite	A	B	C	D	EF
σ⋆ (km s−1)	6.20	11.04	9.88	10.01	11.46
vcirc (km s−1)	9.90	18.85	17.33	16.99	20.96
σ⋆/vcirc	0.63	0.59	0.57	0.59	0.55
σ⋆ (km s−1)	⋯	3.52	4.54	4.92	8.95
${{\rm{v}}}_{\mathrm{circ}}(\mathrm{km}\;{{\rm{s}}}^{-1})$	⋯	6.0	8.0	8.1	17.8
σ⋆/vcirc	⋯	0.59	0.57	0.61	0.50

Download table as:  ASCIITypeset image

3.3. The Dependencies of r_1/2 and σ_* on the Total Bound Mass

3.4. Evolution of Satellites

3.4.1. Galaxy A

Originally the most massive of the studied objects, satellite A enters the host halo on an orbit that brings it very close to the ErisDark's center. Soon after the first pericenter passage, owing to the effect of dynamical friction, the dwarf galaxy settles on an eccentric orbit with pericenter distances smaller than 5 kpc and apocenter distances less than 50 kpc. Before the end of the simulation the cuspy version of the object completes 15 orbits, while the cored version completes 14 orbits, as can be seen from the top left panel of Figure 15. In both cases the satellite is found at z = 0 close to its maximum distance from the host. The 3D half-light radius decreases from the original 1300 pc values to 220 pc and 110 pc for the NFW version and non-NFW version, respectively. The time evolution of the half-mass radius can be inspected in Figure 7.

The tides reach the stellar disk of the cored object and unbind 10% of the stellar mass before the second pericenter is reached. Ninety percent of the stellar mass of the cored dwarf is lost prior to completion of the fourth orbit (Figure 16). In the case of the counterpart with a steep density profile, seven orbits have to be completed before losing the same amount of stars.

Based on the projected mass distribution of the satellite, it can be observed that the difference in the halo profile has no clear impact on the evolution of the shape of the satellite (Figures 10 and 11). The ratios between the shortest and longest dimension corresponding to both versions reaches the value 0.5 during the third pericenter passage and remains above this value throughout the simulation (Figure 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From a kinematical perspective, the ratios between the rotational velocity and the velocity dispersion associated with the respective half-light radius decrease to the 0.5 level differently for the two variants (Figure 9). Before the first orbit is completed the cored satellite is described by a velocity ratio of 0.5. At the same time the velocity ratio of the NFW variant decreases to values in the 1.0–1.7 range. Later during the third orbit, the value 0.5 is reached by the NFW variant as well. It should be noted that the final size of the object with γ = 0.6 is comparable to the softening length of its dark matter particles. Therefore the final stages in the evolution of the object should be regarded with caution as the response of the object to tides is affected by the limited resolution (Moore et al. 1996).

The evolutionary tracks of the descriptive parameters c/a, v_rot/σ_⋆, and M_⋆ are available in Figures 8, 9, and 16, respectively.

3.4.2. Galaxy B

Unlike the previous dwarf galaxy, the second object has a significantly wider orbit. The apocenter distance varies around the 100 kpc value with a pericenter distance less than 10 kpc. During the 11 Gyr timespan of the simulation the satellite completes six orbits, and at the epoch corresponding to z = 0 it is approaching its 7th apocenter (Figure 16). The value of half-mass radius shows large fluctuations in the case of the variant with a shallow dark matter density profile, with large variations when the object is in the proximity of the periapsis. The half-light radius of the alternative presents smaller fluctuations. Each of the versions of dwarf galaxies starts with a half-light radius of 1160 pc and ends with a value of 1010 pc in the case of the cored dwarf and 900 pc in the case of the cuspy dwarf (Figure 7). The stellar disks of the two versions start to lose significant amounts of stars during the second pericenter passage. By the end of the simulations the version with the higher central density has lost half of the initial stellar mass. In the case of the other variant this level of loss is reached before the fourth apocenter epoch. Moreover, at z = 0, around 90% of the original stellar mass is lost (Figure 16).

The velocity ratio corresponding to the NFW profile reaches the level of 1.5 at redshift 0. At the same epoch the cored dwarf reaches a value of 0.5 (Figure 8).

Spatially, the two versions of the object maintain similar ratios of the minor and major dimensions until epoch 8 Gyr. Afterward, in the case of the variant with a cusp coefficient of 0.6, the ratio grows faster. The divergence in the velocity ratio appears at the same time, namely, after the epoch corresponding to the fourth pericenter. Finally, at z = 0, c/a fluctuates around 0.8 for the non-NFW variant and around 0.5 for the NFW variant (Figure 9). When studying the projected density maps it is evident that the cupsy object maintains a two-armed spiral pattern while the cored dwarf approaches spherical symmetry (Figures 10 and 12).

Zoom In Zoom Out Reset image size

3.4.3. Galaxy C

The third most massive object is originally on a trajectory with a pericenter distance of 40 kpc and an apocenter distance of 110 kpc. Soon after one orbit is completed it decays to a trajectory closer to the host. The following orbits are characterized by pericenter distances in the range 10–20 kpc and apocenter distances in the 90–100 kpc range. Like the previous satellite, before the last epoch it completes six orbits (Figure 15).

Kinematically, the characteristics of the two versions diverge after epoch 7.5 Gyr which roughly corresponds to the third apocenter passage (Figure 8). The indicator of the spheroidal transformation c/a clearly reflects the large difference in the transformation level (Figure 9). While the cuspy object maintains a ratio around 1.5, the cored object reaches values below 0.2. The shape indicators of the two variants diverge earlier during the third pericenter passage and eventually converge to 0.7 values at z = 0, roughly the time of the last apocenter passage.

A visual inspection of the density map discloses a more complex case (Figures 10 and 12). The shape of the stellar distribution in the inner regions resembles, case by case, the distribution of the satellite B pair. At the same time, the outer low density and spatially more extended regions are very different. Particularly they are more isotropically distributed in the case of object C than in the case of object B for both versions of the respective object.

The stellar masses of the two versions of satellite C remain approximately constant until the 7 Gyr epoch. From this moment on the stellar mass decreases from 1.7 × 10⁶  M_⊙ to 1.0 × 10⁶ M_⊙ in the case of the cuspy dwarf and to 4 × 10⁵ M_⊙ in its cored counterpart. The original value of the half-light radius for both versions was 810 pc, while the final values are 680 pc and 650 pc in the cuspy and cored dwarf, respectively (Figure 8). Fluctuations similar to ones observed with object B are present, although with smaller amplitudes. Moreover, they start to be visible at later epochs than in the case of object B.

3.4.4. Galaxy D

The object D is bound to ErisDark on an orbit with a pericenter distance of about 20 kpc and an apocenter distance in the 70–80 kpc range. During the simulation timeframe, the dwarf galaxy passes through seven pericenters, and at z = 0 it is approaching the eighth (Figure 15). The tides of the host galaxy start to significantly affect the stellar component just before the 7 Gyr epoch. From this epoch onward the cuspy object loses 3 × 10⁵  M_⊙ mass. During the same period its cored counterpart loses 8 × 10⁵  M_⊙ (Figure 16). Similarly to the previous cases, the stellar component in the cored model begins to lose substantial mass when the onset of large fluctuations in the values of the c/a ratio is also observed, reflecting the more impulsive response to tidal shocks relative to the cuspy case (Figure 9). The fluctuations are particularly enhanced for the cored version. Close to the present epoch c/a reaches ∼0.7 for both variants of the model.

Similarly, the velocity ratio starts to diverge at the 7 Gyr epoch and converge to similar values of 0.5 at z = 0 (Figure 8). Shapewise the final state of the two versions appear to be similar (Figures 11 and 13). Moreover, the evolution of the half-light radius is similar for both cases, showing small fluctuations and a small decrease from the original 640 pc values to 570 pc and 580 pc for the cuspy and cored object, respectively (Figure 7).

Zoom In Zoom Out Reset image size

3.4.5. Galaxy EF

The dwarf galaxy EF, which is the result of a merger between two objects comparable in size, evolves on a wide orbit with the pericenter distance varying from 30 to 50 kpc and apocenter distances with values as high as 160 kpc (Figure 15). After infall the satellite completes three orbits. After the merger, the stellar mass of the object remains virtually constant for both variants (Figure 16). The ratio v_rot/σ_⋆ stabilizes around 0.8 for the object resulting from the merging of two cuspy dwarfs. For the other variant, v_rot/σ_⋆ continues to decrease until it reaches values close to 0.1 at the 13 Gyr epoch (Figure 9). In the case of the shape ratio c/a the object resulting from the merger of the cuspy subhalos has on average a higher value; after 12 Gyr the average c/a is around 0.8 for the cuspy pair and 0.7 for the cored pair (Figure 8).

The original values of the half-mass radius prior to the merging of the two objects were 420 pc for object E and 390 pc for object F. Immediately after the merger the radius of the resulting objects is 550 pc and during the remaining simulation time it increases to 610 pc, with no significant differences between the cuspy and core variants (Figure 7). This highlights a clear divergence from the behavior of the tidally stirred dwarfs, in which the half-light radius always decreases.

In the current work the evolution of five pairs of satellites entering a MW-sized halo has been presented. Each pair consisted of two objects differentiated by a single parameter, namely the cusp coefficient of the dark matter halo. In each pair, each of the objects was composed of a stellar disk component and a dark matter halo. The virial mass of the objects spanned one order of magnitude from 3.3 × 10⁸ M_⊙ to 3.83 × 10⁹  M_⊙. The simulation commences at z = 2.87 when all dwarf galaxies are in the proximity of their future host halo, yet outside its virial volume. The original distances to the host center ranged from 75 to 165 kpc. Since the last major merger experienced by the host halo takes place before z = 3, the described population of objects is representative of the earliest infall objects and thus the most tidally affected population. On average, galaxies with later infall epochs will be affected by the tides to a lesser extent.

All of the objects with cored density profiles match our definition of dSphs before z = 0. At the same time only two out of four objects with cuspy density profiles satisfy the criterion. The more efficient transformation of cored models is expected as their response to tides is more impulsive owing to their lower gravitational binding profile for comparable mass (Kazantzidis et al. 2013). In particular the binding energy decreases much faster toward the center, allowing the tides to disturb the stars found closer to the center of the satellite (Łokas et al. 2012).

Figures 17–20 show the time evolution of various structural properties of all models, highlighting how the cored models generally better fit the range of values that classical dSphs and UFDs span for such properties (McConnachie 2012). We caution, though, that there is a lack of objects with small half-light radii, which may reflect the relatively small sample of initial conditions adopted here (see next section).

The half-light radius, velocity, and shape ratios have significant fluctuations on orbital timescales. Their fluctuations appear to be in phase with the orbital position and grow in amplitude as more stellar mass is unbound from the respective satellite. In general the objects with an original γ = 0.6 show significantly greater amplitudes in the fluctuations than the ones with γ = 1.0, again reflecting the different binding energy profiles which place cored models in the more impulsive regime for the response to tides. The minimum of the fluctuations in half-light radius is reached when the object is near its periapsis, while the maxima is reached a few hundred million years afterward.

In general the satellites with γ = 0.6 are also transformed into dwarf spheroidals faster than their counterparts. The final dark-matter–to–stellar mass ratios are significantly higher in the cuspy remnants (see Table 5), yet the final values in cored remnants are still high enough to match those of classical dSphs and UFDs. Assuming a stellar M/L_* ∼ 3–4, appropriate for an old stellar population (Table 5) would yield values in the range M/L = 90–1000 for the cored remnants. Instead, velocity dispersions and dynamical masses of cuspy remnants indeed appear to exceed observed values in models B and D because of the very high central dark matter density. Note this result is not coincidental; rather, it reflects the fact that while the core remnants have radii of less than 1 kpc, the initial dark-matter–to–stellar mass ratio was very high even inside the central kiloparsec in order to match the constraints from abundance matching with the chosen density profiles of dark matter.

As can be observed in the first panel of Figure 22, objects from the two classes of cored and cuspy dwarf galaxies react differently to mass loss. The dependency of the half-light radius on mass loss is stronger in NFW objects than in non-NFW objects. Cored subhalos that have lost the same amount of mass as their corresponding cuspy subhalos harbor dwarf galaxies with larger radii than the ones contained in the latter. This also explains why the half-light radii of the two subhalo versions remain almost identical over many gigayears (Figure 7). In contrast to the simulations with King models presented in Peñarrubia et al. (2008), the quantities measured in the ErisMod simulations present significant fluctuations on small mass intervals. This is to be expected since in our simulations the objects do not have spherical symmetry. The lack of symmetry allows for dependencies on other parameters such as the orientation of the stellar disk with respect to the orbital plane. Moreover, we use a cosmological N-body environment that contains potentially perturbing substructure. Nevertheless the general trends observed for the King models remain valid for the presented simulations.

All the tidally stirred remnants are found at $R\lt 100\;{\rm{kpc}}$ from the host at z = 0. Within such a distance from the MW and M31 all the galaxies in the mass range of the remnants (stellar masses 10⁴–10⁶${M}_{\odot }$) are dSphs (or UFDs, which structurally appear similar to classical dSphs, as shown by McConnachie 2012). Based on such observations and the presented results we state that unless the progenitors of present-day satellites formed with higher initial pressure support than modeled here, for example, due to the chaotic high-z accretion environment and the effect of SN feedback at earlier times in their own progenitors, it seems natural to assert that they had cored dark matter distributions at their epochs of infall.

Even if the progenitor had high original pressure support, the existence of a steep dark matter density profile will still be questionable. The effects of the supernovae before infall would not be limited to the move toward a more pressure supported stellar system. It would also force the departure from the NFW distribution by displacing matter from the central high-density regions toward the periphery (Governato et al. 2010). Therefore, for the dwarf galaxies with masses relevant to the ErisMod simulations a more pressure-supported progenitor would also imply a progenitor with a shallower dark matter distribution.

Note that idealized simulations already hinted at such a more efficient transformation of the cored models, still they never found such a strong effect because the mass ratio adopted between the satellite and the main host was at least an order of magnitude higher than that used in this work. The reason is that here the host halo is relatively low in mass (closer to ∼10¹¹M_⊙ than to 10¹² M_⊙ when the most massive satellites fall in at z ∼ 2), an effect that only a fully cosmological simulation can account for. The case for cored progenitors is suggested by the notion that dwarf galaxy satellites, which had the highest star formation rates as required to reach a stellar mass in the range of 10⁶–10⁷ M_⊙ already at z > 2, are those in which feedback was most effective to turn cusps into cores (Governato et al. 2010, 2012; Shen et al. 2014; Oñorbe et al. 2015b). These are the objects in which the SF efficiency, defined as the ratio between stellar mass assembled and dark halo mass, was particularly high (Madau et al. 2014). A high SF efficiency in classical dSphs versus similarly low mass field dIrrs is indeed seen in SF histories based on hi-res color–magnitude diagrams (e.g., Skillman et al. 2014; A. Aparicio et al. 2016, in preparation).

Interestingly, gas-poor dwarfs with disk features have been detected in clusters (Lisker et al. 2006) and have been explained as harrassed galaxies that have accreted onto the cluster core only recently. Therefore they have not been able to transform fully into dwarf spheroidals. It is tempting to speculate that such a population would be present in clusters even if the progenitors have dark matter cores before infall since a much larger fraction of galaxies accrete at z < 1 in clusters as opposed to z > 1 for satellites into galaxy-sized halos. This is strongly suggested by the evolution of some of our models, which also in the cored variant transform fully into dSphs only after nearly 10 billion years (models B and D in Figure 9).

Moreover, in order for the velocity ratio of a tidally affected object to reach the 0.5 threshold for a spheroidal it must lose more than 25% of its original stellar mass (Figure 17). The total fractional mass loss was 1/3 at a minimum (in the cuspy version of model satellite D). This occurs for both cuspy and core variants and implies that a tidal stream has to form if the object is transformed into a dSph. The total mass of the stream will thus be comparable to the mass of the associated satellite. Detecting tidal features or tails might be difficult because, as our simulations show, tails and streams show up at a surface density that is between 10³ and 10⁴ times lower than the central surface density of the object (Figures 14, 21), which corresponds to 7–10 mag lower relative to the central surface brightness. The lowest values apply to cuspy models in which stellar removal is less efficient (Figures 14 and 21). This confirms earlier results of Mayer et al. (2002). The transition to a marginally bound tidal tail and fully unbound stream following the orbital path is marked by the flattening of the profile at large radii, typically occurring at 3–6 kpc, or 3–6 times the half-light radius (Figures 14 and 21), confirming past work with a variety of simulation techniques and models (e.g., Johnston et al. 1999; Mayer et al. 2002). This means that a very wide field of view would be needed to unequivocally detect the tidal debris surrounding the dwarf in addition to very deep photometry. The actual shape of the stellar profile also varies with projection and over time (Figure 14), but the basic facts just highlighted hold. We also notice that in some cases the change of slope of the surface density profile between the inner bound and outer unbound stellar component is more evident in the cuspy variant, probably because in the cored models tidal mass loss occurs all the way to the very center, leading to more continuous removal of stars as a function of radius.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of the observable structural parameters for the simulated objects shown in Figures 18, 19, and Figure 20 also highlights that cored models can evolve enough in stellar mass and radius to end up in the region populated by UFDs, confirming a previous claim of Łokas et al. (2012) obtained with idealized simulations. In general the objects evolve from the proximity of the locus populated by dwarf irregulars toward the locus of UFD spheroidals by passing through the classical dwarf spheroidal locus. Except for the satellite A pair all galaxies can be found in the locus of classical dwarf spheroidals and the transition region between the faintest irregular dwarf galaxies and dwarf spheroidals. Rather than being necessarily reionization fossils, as it has been claimed in the literature (Ricotti 2010), some UFDs might be just more extreme outcomes of tidal stirring. This scenario has the attractive aspect of allowing us to explain naturally why there is a substantial continuum between the properties of classical dSphs and those of UFDs. Such a population of UFDs would have cored profiles, while by contrast a population of pristine UFDs formed by gas collapse in mini-halos would likely maintain its original cuspy profile as feedback from SF would have a negligible effect (Peñarrubia et al. 2012).

Despite the low final stellar mass such a population would yield no tension with the proposed core formation mechanism via feedback (Peñarrubia et al. 2012) because it would be the remains of a population of much more massive dwarf galaxies that were severely affected by the tides of the host. Indeed among our remnants the UFD-like object originates from one of the most massive infalling dwarfs whose initial luminosity would be comparable to Fornax, exactly the mass scale at which the effect of feedback should be important (Governato et al. 2012; Shen et al. 2014). Finding UFDs with cored profiles is thus a potential test of our scenario. UFDs forming from an early infall population of satellites would also be consistent with the notion that these galaxies typically lack an intermediate or young population (Brown et al. 2012) as SF would be truncated soon after infall as gas is rapidly removed by tides and ram pressure when the cosmic ionizing background is still at its peak amplitude (Mayer et al. 2007).

The smallest two objects replaced merge while orbiting the host galaxy at significantly larger pericenter and apocenter distances in comparison to the other satellites. The effect of the tides is minimal with insignificant stellar mass loss and small effects in the shape and velocity ratios for the merged objects with γ = 1.0. In the case of the γ = 0.6 mergers the velocity ratio further decreases due to the gravitational interaction with the host. If one considers the condition v_rot/σ_⋆ < 0.5 as the definition of a dwarf spheroidal, then the merging pair with γ = 0.6 clearly reaches this state at the end. At the same time the merger remnant of galaxies with γ = 1.0 fails the velocity criterion, although shapewise the latter appears to be closer to a spheroid than the former. The formation of a significantly round object albeit with rotation still dynamically important cannot be caused by tidal stirring only since in that case the transformation of the initial disky shape is correlated with the redistribution of angular momentum due to tidally induced instabilities and tidal heating.

Hence finding a dwarf with photometric properties analogous to those of classical dSphs but with a significant rotational velocity might be the signature of dwarf–dwarf mergers. Indeed low luminosity ellipticals, whose kinematics are reproduced well by mergers of disks with intermediate masses, do retain significant rotation (Cappellari et al. 2011). One such rare object might indeed be the Tucana dSph, which is on a very wide trajectory around the primary and exhibits appreciable rotation (Fraternali et al. 2009).

Figure 14 also shows that the change of slope in the case of the merger remnant occurs much more gradually, with a clear flattening not appearing before 6–7 kpc from the center. This is not surprising as in this case the morphological evolution is not accompanied by significant tidal mass removal. The different shape of the profile relative to tidally stirred remnants is very interesting because it hints at a further possible way to distinguish between a dSph formed by tidal stirring and one formed via dwarf galaxies mergers. The presence of a swift change from one exponential to another should be the signature of a tidally affected stellar system.

The surface stellar mass measurements can be used as proxies for estimating the surface brightness of the simulated objects. We assume that the mass-to-light ratios of the objects at the final epoch are in the 2–4 M_⊙/L_⊙ range. At a distance of 5 kpc from the center of the dwarf galaxy, the surface brightness values are in the range 34–39 mag arcsec⁻² if ${(M/L)}_{\star }=4\times {M}_{\odot }/{L}_{\odot }$ and 33–38 mag arcsec⁻² if ${(M/L)}_{\star }=2\times {M}_{\odot }/{L}_{\odot }$, respectively. At a distance of 8 kpc the values are 35–38 mag arcsec⁻² for ${(M/L)}_{\star }=4\times {M}_{\odot }/{L}_{\odot }$ and 34–37 mag arcsec⁻² for ${(M/L)}_{\star }=4\times {M}_{\odot }/{L}_{\odot }$, respectively. From the panels of Figure 21 it can be inferred that different orientations of the object with respect to the observer can lead to variations of the surface brightness of about 2 mag arcsec⁻² for both measurement at 5 and 8 kpc distances from the satellites' center. Moreover, the measurements at different epochs suggest that according to the orbital position the surface brightness can change by 3 mag arcsec⁻² at 5 kpc from center and 5 mag arcsec⁻² at 8 kpc, respectively. The low values of the surface brightness presented above can explain why tidal stellar streams associated with low luminosity dwarf satellites present deep inside our galaxy's virial volume have not been detected to date. For most of the tidal dwarf spheroidals, the transitions from one exponential to another take place at distances smaller than 3 kpc from the center (Figure 14). At these points of sudden slope change the surface brightness values are in the range 32–35 mag arcsec⁻² for ${(M/L)}_{\star }=2\times {M}_{\odot }/{L}_{\odot }$ and 33–36 mag arcsec⁻² for ${(M/L)}_{\star }=4\times {M}_{\odot }/{L}_{\odot }$, respectively. Since these transitional regions of a tidally steered satellite are in general brighter than the nearby tidal stream regions, they should be easier to detect.

To the authors' knowledge the current work addresses for the first time the importance of the dark matter density profile for the dynamical evolution and morphological transformation of disky satellites into spheroidals in a cosmological context. Therefore the work bridges the gap between detailed simulations of dwarf–host interaction and fully cosmological hydrodynamical simulations of MW-sized galaxies which at the moment do not have the resolution to reliably address this subject. This is an important step forward since the assembly history of the host, the infall orbits of the satellites, and the structure of the potential well of the host, which all concur to determine the effect of tidal shocks, are all self-consistently taken into account. As a by-product, dwarf–dwarf interactions are also taken into account, and we witnessed a particular example of an interaction that gave rise to a merger.

An important overall conclusion of this work is that tidal stirring is confirmed to be an effective mechanism for transforming disky field dwarfs into objects with properties resembling dwarf spheroidals and even UFDs. Most importantly, it now appears clear that the morphological transformation induced by tidal stirring is a natural and inevitable consequence of hierarchical structure formation. The transformation is enhanced when the satellites have shallow dark matter density profiles, as in our cored models, confirming previous results based on idealized simulations of individual dwarf–host interaction (Kazantzidis et al. 2013). In the latter case, both tidal mass loss of the stellar component as well as tidal heating and internal redistribution are augmented as the response of the galaxy is impulsive down to its central regions. As a consequence remnants can become as faint as UFDs.

We also found that major mergers of satellite galaxies with shallow density profiles can effectively evolve into spheroidals on wide orbits. Since the merger happens before the satellites are accreted by the host, when they have lower velocities, it naturally leads to a dwarf spheroidal orbiting far from the center of the host. Especially in the cuspy case the remnant retains significant rotation despite achieving a spheroidal shape, which would not happen in tidally stirred remnants. This is a plausible scenario for the origin of distant dSphs such as Tucana and Cetus, and confirms a previous claim based on idealized simulations (Kazantzidis et al. 2011a).

Another important conclusion of this work is the identification of two signatures that could separate dwarf spheroidals formed due to a merger on large orbits from the ones formed due to the strong tides generated by the host galaxy. The first signature is the presence of a tidal stellar stream with a total mass comparable to the associated satellite. The second signature can be found in the stellar density/surface brightness profile of the object. It should appear as a fast transition on lengths of a few hundreds of parsecs from the steep central exponential to a significantly shallower one. The objects that are not the source of such a stream or do not have a fast departure from the central brightness exponential decay most probably did not become spheroidal due to the gravity of their current hosts. They were probably formed through the merger of comparable galaxies.

There are of course important caveats in the simulations, which primarily stem from the fact that they are carried out in the dissipationless approximation, namely without treating hydrodynamics, SF, and feedback. The stars are effectively eternal. They do not undergo birth or evolution into stellar remnants. Another caveat is the absence of the baryonic disk of the primary. This would increase significantly the central density of the host, within 10–20 kpc, enhancing tidal mass loss and tidal stirring for satellites with orbits having pericenters smaller than ∼10–20 kpc (D'Onghia et al. 2010). Given the orbits of the satellites analyzed in this paper, this effect would play a role at least for some of them (e.g., model A and model B, which have the smallest pericenters). The neglected effect should enhance the tidal transformation, but may lead to complete tidal disruption of objects that would have otherwise ended up as UFDs, or produce UFDs out of cuspy progenitors (Peñarrubia et al. 2010). Another caveat is the lack of a gaseous component in the disky dwarf progenitors, but previous work has shown that this would be readily removed by ram pressure in the galactic diffuse corona over 1–3 Gyr for the mass range considered here (v_max < 45 km s⁻¹), especially at z > 1, when the cosmic photoionizing background maintains the ISM warm and loosely bound, thus suppressing further SF even for the gas that has not been removed yet (Mayer 2010).

The model galaxies used in the simulations were designed with the goal of creating plausible models of field dwarf galaxies at z = 2–3, motivated by cosmological hydro-simulations (e.g., Shen et al. 2014) as well as by general considerations based on modern constraints on galaxy formation theory, such as the existence of the stellar-to-halo mass relation inferred by abundance matching at varying redshift. Therefore such models are not meant to reproduce the properties of observed present-day field dwarf galaxies. Two noticeable differences between observed dIrr and non-perturbed simulated galaxies are the tendencies of the numerically built objects to be more dark matter dominated and have a comparatively large scale length relative to observed field dIrrs (Figures 19 and 20). Both disparities are understandable in terms of different star formation history. Present-day field dwarf galaxies have probably spent all of their lifetimes in a relatively benign environment, without major external influences. This relative isolation would allow for continued growth of the stellar population while the inner region of the dark matter halo, which is what is measured by kinematical observations, would grow very little from high to low z (Diemer & Kravtsov 2014). Conversely, at z > 2–3, dwarf galaxies except those with virial masses signficantly above 10⁹  M_⊙ would have had their SF suppressed until z ∼ 1–2 due to the feedback of the cosmic ionizing UV flux on cold gas, which quenches the formation of the cold star-forming molecular phase (Gnedin & Kaurov 2014; Shen et al. 2014). For this reason the dIrr galaxies presented in Figures 19 and 20 with comparable half-mass radii and dynamical masses are found to have higher stellar masses than the model galaxies.

Furthermore, our modeling approach relies on the assumption that the progenitors of dwarf spheroidal satellites started out their lives as disky dwarf galaxies. Some cosmological simulations show that this is a well justified assumption for dwarf galaxies forming in halos in the range 10⁹–10¹⁰ M_⊙ (Governato et al. 2010; Teyssier et al. 2013; Shen et al. 2014), while at lower masses, field dwarfs can have a rather thick and turbulent disk, with a ratio ${v}_{\mathrm{rot}}/{\sigma }_{\star }\sim 1$ (S. Shen et al. 2016, in preparation), which appears consistent with the kinematics of the faintest dIrrs in the Local Group (Mateo 1998; McConnachie 2012). In contrast to this, some recent simulations that use stronger stellar and SN feedback find low values for ${v}_{\mathrm{rot}}/{\sigma }_{\star }$ to be common among newly formed dIrrs (Wheeler et al. 2015).

A thicker, kinematically hotter disk could have a different response to tidal shocks. On one hand, the triggering of non-axisymmetric instabilities, such as bars and spiral arms, would be weaker due to the reduced self-gravity and on the other hand, the disk-binding energy will be lower, so the direct effect of tidal heating could be stronger. Kazantzidis et al. (2011a) have considered thicker disks in the wide range of initial conditions they explored in their tidal stirring simulations, finding small differences in the transformation as the initial disk thickness was varied by a factor of three, with only a slight tendency of thicker disks to transition faster to a nearly isotropic configuration. However, future cosmological hydrodynamical simulations with enough resolution to study the morphological evolution of satellites with the same level of detail as done here will be required to assess the impact of disk thickness and kinematics, of both gas and stars, on the efficiency of the transformation.

Finally, here we have considered a relatively small sample of satellites as we focused on the earliest infalling objects with a sizable mass, large enough to host luminous dwarfs and be sufficiently well resolved. The advantage of the small sample is that we could study in detail the different dynamical evolution between cuspy and cored variants of the same object. In addition, since these early infallers are exposed to the tides of the primary for the maximum amount of time possible we can safely assume that the results provide an upper limit on how deeply the structure of the dwarfs can be modified by tidal shocks in cored versus cuspy variants. As a result, we inferred that the population of satellites of the MW should be dominated by remnants of cored objects otherwise we would see a population of dwarfs with intermediate properties between dwarf irregulars and dwarf spheroidals at ∼100–200 kpc galactocentric distances, which is not seen in either the MW or the M31 halo. In a forthcoming paper we will study a much larger sample of satellites which includes objects that are accreted at low redshifts. The size of the sample that we will present will be of the same order of the number of luminous satellites known so far to lie within the putative virial radius of the MW, thus allowing us to make statistically meaningful statements on the evolution of the entire population of satellites.

The first author is grateful to Annalisa Pillepich, Davide Fiacconi, Michela Mapelli, Doug Potter, Stelios Kazantzidis, and Rok Ro$\check{s}$kar for their help. Funding for the work was provided by the Institutes for Theoretical Physics and Computational Science at the University of Zurich.