© The Authors 2024

1. Introduction

The current ΛCDM cosmological model predicts that galaxies grow in size and mass thanks to mergers with smaller systems (e.g. White & Rees 1978; White & Frenk 1991; Springel et al. 2005a; Genel et al. 2014; Planck Collaboration VI 2020). According to cosmological numerical simulations, this pattern extends down to the smallest mass scales (Diemand et al. 2008; Wheeler et al. 2015), where dwarf galaxies form by accreting even smaller systems. However, direct observational evidence for satellites around dwarf hosts is still poor (see Annibali & Tosi 2022, and references therein), due to the challenges associated with detecting very low-surface brightness companions. Additionally, dwarf satellites located in the least massive dark matter sub-halos are expected to be totally dark (Sawala et al. 2015; Read et al. 2017).

Among the population of dwarfs, it is expected from simulations that the ones living in low-density environments may have been unaffected by potential disruptive interactions with larger systems and that they may be still surrounded by a relatively large population of smaller galaxies and stellar streams (D’Onghia & Lake 2008; Deason et al. 2015; Wheeler et al. 2015; Jahn et al. 2019). Indeed, once a group of dwarfs has been accreted into a more massive host, like the Milky Way, it will be disassembled by tidal forces wiping out any evidence for coherent structure in about ∼5–6 Gyr (Deason et al. 2015). Therefore, isolated dwarfs are the most promising targets for studying the hierarchical assembling process at the lowest mass scales.

The Smallest Scale of Hierarchy (SSH; Annibali et al. 2020) survey is a recently completed long-term program at the Large Binocular Telescope (LBT) aimed at identifying merger signatures, typically visible as low surface brightness distortions or tidal stellar features, around a sample of 45 late-type dwarf galaxies in the local (D ≤ 11 Mpc) Universe. The deep optical SSH images are complemented with ancillary H I data from the literature providing information on the gas morphology and kinematics. To mention some of these surveys: the Little Things survey (Hunter et al. 2012), with the Very Large Array (VLA); the WHISP survey (Swaters & Balcells 2002), with the Westerbork Synthesis Radio Telescope; the VLA-ANGST survey (Ott et al. 2012); the FIGGS survey (Begum et al. 2008), with the Giant Metrewave Radio Telescope (GMRT). The ultimate goal of SSH is to characterize the frequency and properties of accretion events around local dwarf systems with different masses and in different environments; this is accomplished by means of state-of-the art hydrodynamical N-body simulations (e.g. Pascale et al. 2021, 2022) capable of reconstructing the dwarfs’ merger histories from the observed properties of gas and stars.

An exceptional example of a multiple-merging dwarf galaxy, observed as a pilot study case for SSH, is the extremely metal poor (∼3% solar oxygen abundance, Pustilnik et al. 2005; Izotov & Thuan 2009), gas-rich dwarf DDO 68 (Annibali et al. 2016), with a stellar mass of ∼10⁸ M_⊙ (Sacchi et al. 2016). DDO 68 resides in the nearby Lynx-Cancer Void (Pustilnik & Tepliakova 2011), and it exhibits a bright, 10 kpc-long stellar stream that crosses the galaxy from south to north, and a second, less luminous stellar stream to the northwest, detected through deep LBT data (Annibali et al. 2019). By means of hydrodynamical N-body simulations, Pascale et al. (2022) showed that the system’s structural and kinematic properties are well reproduced by an interaction between DDO 68 and two smaller satellite galaxies about 150 times and 20 times less massive than DDO 68. However, a third system discovered in VLA H I data (Cannon et al. 2014), and recently resolved into individual stars by Annibali et al. (2023), has been proposed to be interacting with DDO 68 as well, making it a unique case of a small dwarf caught in the process of accreting three smaller satellites. Systematic searches for merger signatures around large samples of dwarf galaxies, like those tackled by SSH and other similar surveys (Carlin et al. 2016; Higgs et al. 2016), may provide insights into the environmental dependence of accretion events around the smallest dark matter halos in the present-day Universe.

In this paper, we present two new compelling cases of dwarf-dwarf mergers identified through visual inspection of the SSH LBT images, namely NGC 5238 and UGC 8760, for which we managed to successfully reproduce the stellar and gaseous observed properties through hydrodynamical N-body simulations. Simulations of galaxies and mergers, despite being idealized setups, provide a complementary perspective to cosmological simulations, offering detailed insights into the processes shaping individual galaxies. For instance, employing higher spatial and time resolutions than in cosmological simulations enables a more detailed analysis of the internal dynamics and interactions within individual systems, facilitating direct comparisons with observational data for specific objects. Moreover, unlike the intricate dynamics of cosmological simulations (wherein concurrent processes complicate the identification of key factors molding a galaxy), isolating and studying individual mergers offer a clearer perspective on galaxy formation and evolution. As an example, studying the merger of galaxies with distinct properties can shed light on the origin of structures like streams or tidal tails. Here, we aim to explain, within a cosmologically motivated framework, the observed substructures in UGC 8760 and NGC 5238 as the outcome of interactions with satellite galaxies. The simulations showcase potential scenarios of interactions, offering plausible explanations based on observational evidence. However it is important to bear in mind that there is no claim of uniqueness in the proposed solutions. Indeed, while simulations of interacting galaxies allow for a more detailed exploration of specific problems, the parameter space to be explored is so vast that, along with computational limitations and potential simplifications in the simulation setup, it becomes nearly impossible to definitively exclude other possible scenarios.

NGC 5238 (also known as UGC 8565, MRK 1479, and IZw 64) is a gas-rich (M_H I #x2243; (2.8 − 2.9) × 10⁷ M_⊙, Cannon et al. 2016; Calzetti et al. 2015) dwarf irregular galaxy at the periphery of the Canes Venatici I group, at a distance of D = 4.51 ± 0.04 Mpc (Tully et al. 2009) from the Earth. According to its low tidal index Θ = −0.4 (a quantifiable metric of the density contrast from massive neighboring galaxies; Karachentsev et al. 2013), NGC 5238 is fairly isolated, more than 1 Mpc from any known massive galaxy. NGC 5238 was observed in the context of several surveys, such as the Hα 11HUGS survey (Kennicutt et al. 2008), the GALEX ultraviolet (UV) imaging survey (Lee et al. 2011), the infrared Spitzer Local Volume Legacy (LVL) survey (Dale et al. 2009), the HST optical ANGST survey (Dalcanton et al. 2009) and the HST UV LEGUS survey (Calzetti et al. 2015). With an oxygen abundance less than 20% solar (Lira et al. 2007) and a stellar mass of M_⋆ ∼ 10⁸ M_⊙ (Cook et al. 2014; Cannon et al. 2016), the galaxy is compatible with the mass-metallicity relation derived for star forming dwarfs (Berg et al. 2012). The star formation rate density (0.53 × 10⁻³ M_⊙ yr⁻¹ kpc⁻²), derived from modeling the color-magnitude diagram (CMD) of resolved stars in the galaxy, is consistent with the value inferred from the H_α flux (0.53 × 10⁻³ M_⊙ yr⁻¹ kpc⁻²; Cignoni et al. 2019). The VLA H I images acquired by Cannon et al. (2016) reveal an irregular and asymmetric morphology in the outermost detected regions, and significant velocity asymmetries are present throughout the disc. Nevertheless, Cannon et al. (2016) do not report evidence for an ongoing interaction considering the sensitivity and angular resolution of the VLA observations. The region of current star formation, as traced by the UV and Hα emission, exhibits a ∼300 pc offset with respect to H I peak, a behaviour already observed in other dwarf galaxies (e.g. Roychowdhury et al. 2009).

The dwarf irregular UGC 8760 (or DDO 183), at a distance of D = 3.2 ± 0.1 Mpc (Dalcanton et al. 2009; Tully et al. 2009, 2013), is another fairly isolated galaxy (Karachentsev et al. 2013) and resides in the outer zone of the Ophiuchus–Sagittarius–Capricornus void (Pustilnik et al. 2019). As for NGC 5238, multi-wavelength data of UGC 8760 are available from various surveys (e.g., LVL, 11HUGS, ANGST, VLA-ANGST, and FIGGS). Stellar mass estimates derived from different methods are in the range M_⋆ = (2 − 4)×10⁷ M_⊙ (Johnson et al. 2013; Cook et al. 2014), while the H I mass from FIGGS and VLA-ANGST data is M_H I ≲ 2.6 × 10⁷ M_⊙ (Begum et al. 2008; Ott et al. 2012). The H I disk extends for a diameter of ∼4 kpc (Roychowdhury et al. 2017) and it shows coherent rotation with some level of distortion in the inner regions (Begum et al. 2008; Ott et al. 2012). According to the resolved-star CMD modeling and the integrated UV flux, the galaxy has been actively forming stars with an average rate of ∼ 4 × 10⁻³ M_⊙ yr⁻¹ over the last hundred Myr (Lee et al. 2009; Johnson et al. 2013). To our knowledge, no spectroscopic metallicity estimates are available for this galaxy.

This paper is structured as follows: In Sect. 2, we present the new LBT data where signatures of merger events have been identified. Sect. 3 provides a brief description of the simulation setup, while Sect. 4 presents and discusses our results. Finally, in Sect. 5, we draw our conclusions.

2. LBT data

NGC 5238 and UGC 8760 were observed with the Large Binocular Camera (LBC) on the LBT in the g and r bands as part of the SSH survey. Figure 1 shows 8′×5.3′ (left panel) and 6.4′×4.2′ (right panel) portions of the larger 23′×23′ imaged fields of view, with the H I contours from VLA data from Cannon et al. (2014, NGC 5238) and from Begum et al. (2008, UGC 8760) superimposed. In the following, we briefly summarize the main properties of the stellar populations detected in the two galaxies and their associated substructures. A detailed and thorough analysis is presented in Sacchi et al. (2024). The deep LBT data reveal the presence of extended, low surface brightness stellar structures such as a northern “umbrella” and a southern plume in the case of NGC 5238, and an extended, disrupted-satellite feature to the northeast of UGC 8760.

Fig. 1. LBT images of NGC 5238 and UGC 8760. Left-hand panel: LBT g-band image of NGC 5238. The arrows show the positions of the northern and southern overdensities of stars while the red contours show lines of constant H I emission from VLA observations (Cannon et al. 2016) at the adopted distance. Right-hand panel: same as the left-hand panel but for UGC 8760. Here, the arrow shows the position of the northern overdensity of stars, while the VLA observations are from Begum et al. (2008). In both panels, the red contours correspond to 1, 2, 4, 8 M⊙ pc−2. North is up, east is to the left.

These stellar features are also visible in the resolved-star maps of Fig. 2, where we plot the distribution of both old red giant branch (RGB) stars and young blue bright stars. A detailed description of the point-source photometry performed on the LBT images and of the selection of RGB and young blue stars based on their location in the r versus g–r CMD are provided in Annibali et al. (2020, 2022). Here, it is sufficient to recall that RGB counts trace old stars with ages between ∼2 and 13 Gyr, while the bluest and brightest counts in the CMD correspond to stars typically younger than about a hundred Myr. These maps provide therefore a powerful tool to disentangle short-lasting, irregular structures due to active star-formation from distortions in the old stellar component, most likely of tidal origin.

Fig. 2. Spatial distribution of old (age >  2 Gyr) RGB stars and young (age ≲ 100 Myr) bright blue stars in NGC 5238 (left) and UGC 8760 (right). RGB maps are background-subtracted, and only pixels above a level of 3 times the background standard deviation (given by the vertical black segment in the colorbars) have been set to non-zero to mask density peaks due to red background galaxies. The blue points are individual young stars. H I contours at 1, 2, 4, 8 M⊙ pc−2 (same as in Fig. 1) are superimposed to the star count maps. North is up, east is to the left.

Indeed, in the case of NGC 5238, the “umbrella”-like feature detected in the deep LBT images coincides with a similar structure in the RGB count map, indicating that it cannot be accounted for by recent and offset bursts of star formation, but that it can be either the remnant of an accreted companion or a distortion of NGC 5238 of external origin. More in general, we note that the old stellar halo of NGC 5238 extends for a diameter of ∼360 arcsec, or ∼8 kpc, well beyond the outer displayed 1 M_⊙ pc⁻² H I contour. Similar to NGC 5238, the elongated northeast feature detected in UGC 8760 coincides with a clear RGB star overdensity. However, no extended stellar halo is detected in this galaxy. In both NGC 5238 and UGC 8760, stars younger than ∼100 Myr are confined to the innermost galaxy regions, producing significant stellar crowding and preventing the detection of faint RGB stars (Cignoni & Tosi 2010; Okamoto et al. 2015; Cignoni et al. 2019; Higgs et al. 2021).

3. Setting up the simulations

In this Section, we explain how we determined the relevant parameters of the galaxy models used to generate the initial conditions (ICs) of the simulations for NGC 5238, UGC 8760, and their associated satellites. We will explore the specific properties of the host galaxies and their satellite counterparts in Sects. 3.1 and 3.2, respectively. Additionally, Sect. 3.3 will provide a detailed description of our simulations setup.

3.1. Hosts

3.1.1. The galaxy models

Consistently with their classification (SABdm), NGC 5238 and UGC 8760 are modeled as dwarf disc galaxies with a dominant dark-matter halo, and thick stellar and gas discs. The dark matter halo follows the Dehnen (1993, γ = 0) profile

(1)

with M_dm and a the total dark matter mass and the halo scale length, respectively. Here r is the three dimensional distance from the halo center. Although the NFW profile is cosmologically motivated, we prefer to use a Dehnen profile because: (i) it has finite mass, being more practical to sample the halo particle phase-space positions, and (ii) we expect dwarf galaxies to be filled with cores of constant dark matter densities rather than cusps (Bullock & Boylan-Kolchin 2017, and reference therein).

The density distribution of the stellar disc is

(2)

In the above equation, R and z are the cylindrical coordinates of a reference frame whose z-axis is the disc’s symmetry axis, M_⋆ is the total stellar mass, h_⋆ is the stellar disc scale length and z_⋆ is the stellar disc scale height. The gas surface density follows the exponential profile

(3)

where h_gas is the gas disc scale length and M_gas is the total gas mass, given summing up the mass contributions of neutral hydrogen (M_H I) and helium (M_He):

(4)

We impose the vertical hydrostatic equilibrium

(5)

to compute the gas disc scale height, and

(6)

In Eqs. (5) and (6), Φ_tot is the total gravitational potential, ρ_gas is the gas density, P = ρ_gas(γ − 1)e its pressure, γ = 5/3 the adiabatic index and e its specific internal energy, which we assume to be constant throughout the galaxy.

We assume that the distribution function (DF) of all the collisionless particles depends on the specific energy E and the specific third component of the angular momentum L_z only, which makes all the mixed moments of the DF and the mean velocities in the vertical and radial directions null (Binney & Tremaine 2008). The non-null elements of the velocity ellipsoid (the vertical and radial velocity dispersions σ ≡ σ_z = σ_R and the velocity dispersion in the streaming direction σ_ϕ) are determined from the axisymmetric Jeans equations. We sample the halo particle velocities in Maxwellian approximation: in the vertical and radial directions, the velocity distributions are Gaussians with null mean and standard deviation equal to σ, while, in the streaming direction, the velocity distribution is a Gaussian with non-null mean v_ϕ and standard deviation (i.e. the halo is allowed to rotate). The stellar disc particles are sampled in Maxwellian approximation, but the streaming velocity is set by the epicyclic approximation. The only non-null gas velocity component is the azimuthal one, found imposing the stationary Euler equation (Eq. (21) of Springel et al. 2005b). The total potential is found iteratively, while the potentials of the discs are computed numerically from a discretized mass distribution following Eqs. (2) and (3) with a hierarchical multipole expansion based on a tree code (details in Springel et al. 2005b).

3.1.2. Fixing the free parameters

In the context of simulating interactions between galaxies, such as those observed in UGC 8760 and NGC 5238, careful consideration is given to setting: (i) the mentioned galaxy models, (ii) the orbital parameters defining the merger. For instance, we do not have a formal understanding of whether the properties of the galaxies post-interaction will necessarily resemble those pre-interaction. However, in the simulations, we set the latter to reproduce the former. Nonetheless, the selection is guided by logical assumptions and observational evidence. For instance, as we discuss in Sect. 3.3 and confirm in Sect. 4, the dynamics of the interactions represented by the simulations we propose is not destructive for the host galaxies. Also, the simulations lack star formation (an additional mechanism for gas consumption), stellar feedback, and, as we will see in Sect. 3.2, the baryonic contributions of the satellites to the total baryonic mass are minimal. With these considerations in mind, we adopt the strategy to select the parameters of the pre-interaction galaxies based on observations, even though observations refer to a post-interaction scenario.

To determine the parameters of the Dehnen halo Eq. (1) we link them to those of a cosmologically motivated NFW Navarro et al. (1996, hereafter NFW) model. In this way, when fixing the NFW model we fix the Dehnen profile. A NFW model is given by the density distribution

(7)

with ρ_c the present day critical density of the Universe (we assume H₀ = 67.4 km s⁻¹ Mpc⁻¹; Planck Collaboration VI 2020) and

(8)

The halo virial radius r_vir is defined as the distance from the center where the average halo density equals 200 times ρ_c, and it is linked to the concentration c by c ≡ r_vir/r_s. The halo virial mass M_vir is, then, the mass enclosed within r_vir. We determine the parameters of the Dehnen halo Eq. (1) using equations

(9)

The above relations establish, analitically, the link between the parameters of a given Dehnen profile and the parameters of a general NFW model. A similar approach has been used by Springel et al. (2005b) to link the parameters of the Hernquist (1990) model to that of a cosmologically motivated NFW profile. As in Springel et al. (2005b), we impose that the Dehnen total mass equals the NFW virial mass (first Equation of (9)), while we require that the Dehnen profile shares the same outer normalization of the Hernquist profile where both decrease as r⁻⁴. As an example, Fig. 3 shows the halo density distribution of NGC 5238 obtained with the three different halo models (NFW, Hernquist and Dehnen) linked via Eqs. (9). Here, the overall agreement among the different models can be visually observed. In the simulations we then adopt the Dehnen halo.

Fig. 3. Dark matter density distribution of NGC 5238 used in the simulations described via three different density models (different colors are associated to the three dark-matter models as indicated in the legend). The different models are linked via Eq. ((9), see text for details). The vertical dashed lines show the NFW, Hernquist and Dehnen scale radii using the same color coding, while the vertical blue solid line is the NFW virial radius.

In a NFW model, the maximum circular velocity is reached at a distance

(10)

where it assumes the value

(11)

with δ as in Eq. (8) and k = 2.162 a dimensionless constant¹. The quantities M_vir and r_s (thus, M_dm and a) have been determined imposing that v_c, max and R_vc, max approximately match the deprojected maximum circular velocity and the galactocentric distance from the H I kinematics.

Below we list the specific parameters adopted to set up the galaxies’ models.

NGC 5238. We impose v_c, max = 31.5 km s⁻¹ at R_vc, max = 3.15 kpc, larger than R_vc, max = 1.3 kpc measured by Cannon et al. (2016). This choice is motivated by the fact that the value of R_vc, max from Cannon et al. (2016) implies an unreasonably high concentration (c ≃ 40). Therefore, we prefer the value of R_vc, max = 3.15 kpc, which ensures a more cosmologically motivated concentration (c ≃ 20) while keeping a dynamical mass (Dehnen dark matter halo, stars and gas, see below) of 2.94 × 10⁸ M_⊙ within a 1.3 kpc radius, in very good agreement with a mass of 3 × 10⁸ M_⊙ from Cannon et al. (2016). To compute the deprojected maximum circular velocity, Cannon et al. (2016) derived an inclination of i = 50° ±5°, in agreement with the previous estimate by Lee et al. (2011) of i = 55°. For the stellar mass, we adopt M_⋆ = 0.89 × 10⁷ M_⊙, taken from Cook et al. (2014) but scaled for the closer distance (∼4.5 Mpc versus 5.2 Mpc) assumed in this work. This value is smaller than the 1.4 × 10⁸ M_⊙ stellar mass of Calzetti et al. (2015). We assume M_H I = 2.74 × 10⁷ M_⊙, within the errorbars as in Cannon et al. (2016, (2.8 ± 0.3)×10⁷ M_⊙), plus a correction for helium of M_He = 0.35 M_H I (e.g.; Leroy et al. 2008; Bolatto et al. 2013; Hardwick et al. 2022). The total gas mass is M_gas = 3.7 × 10⁷ M_⊙, as computed by Cannon et al. (2016). To determine h_⋆ and h_gas, we fitted both the RGB star count profile derived from the LEGUS catalogs² and the H I surface density profile from Cannon et al. (2016) with exponential disc models. Although the fit provides h_gas ≃ h_⋆ = 0.47 kpc, for our simulations we adopt a larger value of h_⋆ = h_gas = 0.52 kpc to keep the stellar disc stable against bar formation when evolved in isolation (see Sect. 3.3). For the same reason, we prefer the smaller stellar mass from Cook et al. (2014) to that provided by Calzetti et al. (2015). From a qualitative point of view, decreasing the stellar mass and/or increasing the scale radius, decreases the stellar disc surface density. This favors the stability of the stellar disc against gravity. Indeed, a razor-thin disc is locally unstable against gravity when its Toomre (1964) parameter , with ν the local epicyclic frequency and Σ_⋆ the stellar surface density.

UGC 8760. We require v_c, max = 25 km s⁻¹ at 3 kpc, consistent with the velocity peak observed in the H I velocity field map of UGC 8760 from the VLA-ANGST survey, assuming an inclination i ≃ 67° ±3° (Begum et al. 2008). The implied M_dm and a of the corresponding Dehnen halo are listed in Table 1, alongside all the other relevant parameters used to generate the galaxy model. For the stellar disc, we adopt a mass of M_⋆ = 2 × 10⁷ M_⊙ from Cook et al. (2014). The stellar disc scale length is 0.275 kpc, taken from Makarov et al. (2017) who modeled the galaxy B,V and I surface brightness emission via the same exponential model as in Eq. (2). We set M_H I = 2.64 × 10⁷ M_⊙, larger than the value from Georgiev et al. (2010, 1.86 × 10⁷ M_⊙), but consistent within the errors with the estimate from Begum et al. (2008, 2.59 × 10⁷ M_⊙). We assume that the helium mass is, again, M_He = 0.35 M_H I = 0.93 × 10⁷ M_⊙.

3.2. Satellites

As discussed in Sect. 2, the stellar populations observed within the substructures of NGC 5238 and UGC 8760 suggest that the hypothesis of their formation through an asymmetric event of star formation is implausible. On the contrary, an external origin appears more probable. The stellar masses of NGC 5238 and UGC 8760 range between 2 × 10⁷ M_⊙ and nearly 2 × 10⁸ M_⊙. According to Dooley et al. (2017), for galaxies with these masses there is a significant probability, 50% for the former and up to 80% for the latter, of hosting at least one satellite with a stellar mass of 10⁵ M_⊙ at most. This provides the basis for our approach to model the properties of the satellite galaxies in the simulations.

3.2.1. The satellite of NGC 5238

The satellite galaxy of NGC 5238 is a scaled-down version of its host: a disc galaxy with a dominant Dehnen dark matter halo, an exponential thick stellar disc, and an exponential gas disc. As previously mentioned, we expect dwarf galaxies with stellar mass similar to NGC 5238 or UGC 8760 to have, at most, one satellite with a few 10⁵ M_⊙ (Dooley et al. 2017). Thus, we fix the satellite stellar mass to M_⋆ = 2 × 10⁵ M_⊙, while h_⋆ = 0.15 kpc. The expected virial mass (from the stellar-to-halo mass relation in Read et al. 2017) and concentration (from the virial mass-to-concentration relation in Muñoz-Cuartas et al. 2011) are, respectively, M_vir = 3 × 10⁸ M_⊙ and c = 19.5, which we use, in combination with Eqs. (9) to fix M_dm and a. The gas mass is M_gas = 10 M_⋆ (Geha et al. 2006; Cimatti et al. 2020), while h_gas = 2 h_⋆.

3.2.2. The satellite of UGC 8760

In UGC 8760, the roundish shape of the northern overdensity of stars, along with the absence of H I emission at that location, motivates us to model its satellite as a gas-free, dark-matter dominated, spheroid rather than a small disk galaxy. We assume that dark-matter density of such a satellite follows the Hernquist model

(12)

In Eq. (12), M_dm and a are, respectively, the total mass of the dark matter halo and its scale radius. The stellar component is a spherically symmetric Sérsic (1968) model with surface density

(13)

with Γ the gamma function, M_⋆,sat the total stellar mass, R_{e, sat} the effective radius (the projected distance from the center that contains half of the total mass), m the Sérsic index, and b = b(m) as in Eq. (18) of Ciotti & Bertin (1999), which ensures that .

We assume that the halo and the stellar velocity distributions are isotropic, and we sample their particle positions and velocities directly from the corresponding ergodic DFs (Binney & Tremaine 2008)

(14)

where i labels the stellar (i = ⋆) or the dark matter (i = dm) DFs, while Φ_tot is the total potential (stars and dark matter). Note that Φ_⋆, sat is found numerically from the Poisson equation after computing ρ_⋆, sat from Eq. (13) via the Abel Inversion

(15)

As for NGC 5238, we may expect UGC 8760 to be orbited by one satellite only. We, thus, fixed the parameters of the satellite of UGC 8760 to M_⋆,sat = 6 × 10⁵ M_⊙, m = 0.75 and R_{e, sat} = 0.15 kpc, while M_{dm, sat} = 1.5 × 10⁸ M_⊙ and a_sat = 0.5 kpc. The chosen sizes and dark matter mass are such to mimic those of local dwarf elliptical or dwarf spheroidal galaxies (see, for instance Tucana, Carina, Andromeda XVI or Andromeda XVIII; McConnachie 2012). The selected satellite halo parameters are, in any case, consistent with those one would obtain based on cosmological considerations, using the same approach employed in NGC 5238.

Differently from NGC 5238, in the case of UGC 8760, it is possible to make a rather rudimentary estimation of its satellite stellar mass based on the luminosity of the substructure. Considering the uncertainty regarding the actual extent of the overdensity observed to the north, and after estimating the background in that region, we obtain M_r = −9.05 and M_g = −8.2 (at the adopted distance) as r-band and g-band absolute magnitudes. After converting the g/r magnitueds into V-band magnitude (Lupton 2005)³, considering mass-to-light ratios compatible with those of an old stellar populations (2–6), and marginalizing over a ∼0.5 mag error, the stellar mass can well be within (0.2 − 1.5)×10⁶ M_⊙, in agreement with the value we adopted.

3.3. Initial conditions

All the simulations are performed with the publicly available version of the hydrodynamical N-body code AREPO (Springel 2010; Weinberger et al. 2020). AREPO blends the advantages of both Lagrangian and Eulerian hydrodynamics by solving Euler equations on an unstructured Voronoi mesh that moves with the fluid flow. The gas moving mesh is integrated with a particle-mesh algorithm and an oct-tree approach (Barnes & Hut 1986) for solving the Poisson equation and calculating accelerations for both collisional and collisionless particles. We adopt a mesh refinement algorithm to refine (de-refine) the gas cells when they become more (less) massive than a target reference mass. This choice allows to sample high density regions with a larger number of gas cells.

The dark matter and stars of NGC 5238 and the dark matter of its satellite have a mass particle of m_part = 700 M_⊙, while the satellite’s stellar component is sampled with a particle mass of m_part = 350 M_⊙. Similarly, for all components of UGC 8760 (stars, dark matter, and gas), we adopt a particle mass of m_part = 350 M_⊙. For the satellite, we assume a dark matter particle mass of m_part = 350 M_⊙, while we adopt a lower mass of 117 M_⊙ for its stellar component to ensure adequate sampling. The dark matter halos are sampled out to 10 times the virial radius of the associated NFW halo (= 53a), where 98% of the total mass is contained. The softening of each component is such that the maximum force between the particles of that component should not be larger than the mean-field strength they generate. Table 2 summarizes the number of particles for each component and the softening lengths used in the simulations. Since the galaxy models of NGC 5238, UGC 8760 and their satellites have been sampled via Maxwellian and epicyclic approximations, they are not guaranteed to be strictly in equilibrium. Therefore, we let them evolve in isolation for 6 Gyr to make them shift toward equilibrium. We also ensured that none of the galaxies developed bars due to disc instabilities when evolved in isolation. This approach is similar to that used in Pascale et al. (2021, 2022). The outcomes of these equilibrium simulations are then used as ICs for the simulations where the galaxies interact.

In the simulations where we try to reproduce the interaction between hosts and satellites, the centers of mass of NGC 5238 and UGC 8760 are placed at the center of their respective simulation boxes. In the case of NGC 5238, the disc of the galaxy is inclined by 55° with respect to the positive z-axis. Its satellite is at the initial position (x, y, z) = (15 kpc, 0, 0) with initial velocity (v_x, v_y, v_z) = ( − 5 km s⁻¹, 0, 0), i.e. on a radial orbit that will cross the host. These ICs are tailored to increase the strength of the gravitational interaction between host and satellite, with the satellite crossing NGC 5238 in the north-south direction, as the direction of the observed substructures. The satellite’s disc lies on the xy-plane and rotates clockwise. As we discuss in Sect. 4, in this configuration ram pressure completely removes the gas from the satellite upon its first passage. Thus, specifying whether the satellite is a disc or a spheroid becomes less critical in NGC 5238 rather than in UGC 8760. It is only essential that, post-interaction, the satellite’s stellar component either falls below the surface brightness limit or undergoes sufficient disruption, which is easily accomplished given its stellar mass and size.

In the case of UGC 8760, its disc is inclined by 75° counterclockwise with respect to the positive z-axis, with the x-axis lying in the plane of the disc. Its satellite is positioned at (x, y, z) = (5 kpc, 0, 0), with initial velocity (v_x, v_y, v_z) = (0, 25 km s⁻¹, 0), i.e. on an orbit that would be circular if we did not account for dynamical friction. The discs of UGC 8760 and NGC 5238 would rotate counterclockwise if they belonged to the xy-plane. Note that 55 deg and 75 deg are just how much NGC 5238 and UGC 8760 discs are inclined with respect to the xy-plane in the ICs. They would correspond to the formal discs inclinations only if the y-axis is assumed as line-of-sight.

In all the simulations we set the target reference mass to 700 M_⊙. Note that due to mesh refinement, the initial total number of gas cells in the simulations reproducing the interactions is not the same as in the ICs of the simulations evolved in isolation.

4. Results

Here, we present the results of these two simulations that accurately replicate the observable characteristics of the target galaxies. The specific cases of NGC 5238 and UGC 8760 are discussed separately in Sects. 4.1 and 4.2, respectively.

The sequence of panels in Fig. 4 shows a selection of surface brightness maps from different snapshots of the simulations of NGC 5238 (top) and UGC 8760 (bottom). The left-most column shows the ICs of the simulations, the second and third columns are subsequent snapshots, corresponding to 0.67 Gyr and 0.88 Gyr for NGC 5238, and 0.39 Gyr and 0.69 Gyr for UGC 8760. The surface brightness maps in the third column are meant to reproduce the LBT observations displayed, for a comparison, in the rightmost panels. In the mock surface brightness maps, the line-of-sight is obtained by rotating the systems counterclockwise first about the x-axis by an angle α, then about the z-axis by an angle γ to align the position angles to those of the observations. For NGC 5238 (α, γ) = (165° ,110° ), for UGC 8760 (α, γ) = (175° ,120° ). These rotations correspond, respectively, to disc inclinations of i ≃ 40° and i ≃ 70° w.r.t. the line-of-sight. In the simulated surface brightness maps, we converted mass into LBT g-band luminosity assuming constant mass-to-light ratios of (M/L_g)/(M/L_bol)_⊙ = 3.2 for NGC 5238 and its satellite, and (M/L_g)/(M/L_bol)_⊙ = 3.2 and 6 for UGC 8760 and its satellite. Here, (M/L_g)/(M/L_bol)_⊙ = 3.2 is the mass-to-light ratio corresponding to a 5 Gyr old stellar population with metallicity Z = 0.01. This is a reasonable assumption for our simulated gas-rich dwarf galaxies because, according to numerous studies based on resolved-star CMD modeling (e.g.; Tolstoy et al. 2009), this type of systems is known to have experienced a continuous star formation since 10 Gyr ago (or earlier) to the present epoch, but we can exclude from our LBT photometry the presence of luminous stars younger than 1 Gyr in the extended low surface brightness galaxy outskirts that we intend to reproduce. On the other hand, (M/L_g)/(M/L_bol)_⊙ = 6, as expected for a 10 Gyr old stellar population with the same metallicity, is more appropriate for the satellite of UGC 8760 that we model as a gas-free spheroidal galaxy. Since the simulations do not model star formation, star particles cannot trace the dynamics of young stars. Therefore, by construction, the mock surface brightness maps shown in Fig. 4 cannot reproduce the brightest central parts of the galaxies where most of the contribution to light comes from young stars (see Fig. 2). In Fig. 4, we used the same pixel size (0.225 × 0.225 arcsec²) of the LBC camera and imposed a surface brightness limit of ≃30 mag arcsec⁻² as in the observations. The maps were smoothed with a Gaussian kernel of 0.45 arcsec (∼1 arcsec FWHM) dispersion to account for typical seeing during the observations. To help the visual comparison between data and simulations, the physical scales in the third and forth columns panels are the same.

Fig. 4. Comparison between simulations and observations. First row of panels: stellar surface brightness maps from the reference simulation of NGC 5238 and the corresponding LBT observation. The left panel shows the ICs of the simulation, while the middle panels are subsequent snapshots (as reported in each panel). The third panel corresponds to the time when the simulation reproduces the data, shown for comparison in the rightmost panel. In the simulated maps, the pixel size is 0.225 arcsec × 0.225 arcsec (as for the LBC camera), and the maps have been convolved with a Gaussian beam of 1 arcsec FWHM to match the typical seeing during observations. To generate the maps, we adopted constant mass-to-light ratios and we set the surface brightness limit to 30 mag arcsec−2, to match the depth of the data. Blue lines show the orbit of the satellite, obtained by joining the centers of mass of its dark matter halo computed with the shrinking sphere method (Power et al. 2003). The color gradient indicates the direction of the orbit itself: the satellite moves from where the orbit is darker to where it is lighter. The red contours are levels of iso-H I column density, specifically 2k M⊙ pc−2, with k = 0, …, 3. Second row of panels: same as the top row but for the UGC 8760 system. For details on the simulation setup see Sect. 3.3 and Tables 1 and 2. In the top row, third panel, we highlight in blue the stellar particles belonging to the satellite of NGC 5238 as a comparison (see text for details).

4.1. NGC 5238

In the NGC 5238 simulation, the satellite galaxy starts from a significantly large distance, with a low initial velocity, aimed at reproducing a purely radial interaction (see top-left panel of Fig. 4). Unlike the case of UGC 8760, these ICs are designed to increase the strength of the interaction between the two galaxies, a choice motivated by the pronounced irregularities observed to the north and south of the galaxy. After 0.67 Gyr, the satellite galaxy has entirely crossed NGC 5238 once, and this passage is indeed highly disruptive for both the satellite system and the host. On NGC 5238 the gravitational effect is such to significantly distort its disc. The satellite, instead, loses coherence to the point that its surface brightness reaches the detection threshold, making it extremely faint. To aid in identifying the satellite galaxy, we display its stellar particles in blue in Fig. 4 (top row, third panel from left).

One significant consequence of this interaction is the formation of a predominant spiral arms pattern. The formation of spiral arms induced by interactions and mergers has been extensively explored through numerical simulations (e.g. see Gerin et al. 1990; Lang et al. 2014; Peschken & Łokas 2019) and they are recognized as a key mechanism contributing to the development of spiral arms within galaxies. Generally, tidally induced spiral arms are regarded as transient phenomena, with their longevity depending upon the mass ratio of the interaction, as discussed in Oh et al. (2008) and Struck et al. (2011). In cases where the mass ratio is 1:10, such as in the reference simulation of NGC 5238, these spiral arms can persist for several billion years. In our simulated scenario, it is intriguing to note that the substructures located to the north and south of NGC 5238 are not comprised of stars originating from the satellite itself. Instead, these substructures manifest as extensions of the spiral pattern induced by the interaction. Particularly noteworthy is the southern structure, which is more faithfully reproduced in the g-band image by the interaction. The northern “umbrella”-like pattern, appears, instead, less accurately replicated being less extended in the simulation than in the observation. We do not consider this discrepancy in the replication of northern substructure a major concern. Indeed, as it can be noted from the RGB map in Fig. 2, the “umbrella” feature is less dense than the southern structure and likely more noticeable in the g-band map due to the presence of background galaxies which are removed, instead, when computing the RGB star counts map.

In the observations of NGC 5238 the spiral pattern is scarcely detectable within the central region of the galaxy. Nonetheless, it is crucial to emphasize, again, that the simulations do not incorporate star formation and its associated feedback. Consequently, the stellar particles primarily represent an old stellar population (prominently visible in the outer regions; Sect. 2) and they may not be very representative of the central regions, where the younger stellar component and the star forming gas predominantly influence the morphology of the galaxy (Cannon et al. 2016).

As anticipated, the satellite is scarcely detectable (see Fig. 4, third column, top panel), since it undergoes an almost complete disintegration after less than an orbital revolution. This results in a substantial decrease in surface brightness, to the extent that it falls almost below the surface brightness threshold (average surface brightness 29.5 arcsec²). As we already mentioned in Sect. 3.3, the distinction between whether the satellite is a disk or a spheroid also becomes of minor importance, since the observed features are extensions of the induced spiral pattern and since the satellite has been depleted of gas as a result of ram pressure.

As a comparison, in the LBT surface brightness maps of Fig. 4 (rightmost panel) we superimposed levels of H I isodensities taken from the VLA observations from Cannon et al. (2016) for NGC 5238 and from Begum et al. (2008) for UGC 8760. The same quantities have been computed from the simulations and shown in the first three columns. The isodensity levels correspond to 1, 2, 4 and 8 M_⊙ pc⁻². As noted by Cannon et al. (2016), the outer H I disc displays a few asymmetries, with its morphological major axis elongated from northeast to southwest. The high-column density H I gas exhibits a crescent-shaped morphology, with the maximum of the H I surface density offset from the central optical peak by approximately 300 pc to the northeast, and not co-spatial with the brightest UV and H_α emission. In the simulation, the gas distribution appears to be smoother than in the VLA observations by Cannon et al. (2016) and the offset is not reproduced. This discrepancy probably arises from the fact that our simulations are inherently adiabatic, thereby lacking the inclusion of star formation and associated feedback processes, preventing the faithful replication of distortions associated with these phenomena. Nevertheless, following the interaction between the satellite and the galaxy, we manage to reasonably replicate the spatial gas distribution. The overall extent of the gas distribution is, nevertheless, well reproduced although the outer H I gradient tends to be somewhat steeper in the observations than in the simulation post-interaction, as it can be appreciated looking at the different outer density of H I iso-contours in Fig. 4 (third and forth panels). The H I distribution in the simulation aligns well with recent, unpublished VLA observations, confirming a confined gas distribution that is less extended compared to the spatial distribution of stars. A detailed analysis of the H I distribution and its kinematics, based on these novel VLA observations, will be presented in an upcoming paper (Cannon et al., in prep.).

In Fig. 5 we show the H I density in velocity channels computed from the simulations together with the overall H I velocity map compared with the analogous quantity computed from the simulation. The observed H I velocity field has significant asymmetries throughout the disc, featuring an S-shaped velocity profile that has been suggested by Cannon et al. (2016) as a potential indicator of a warp in the disc. In the averaged velocity field from the simulation (Fig. 5, bottom-left panel), this S-shaped feature is fairly well replicated, especially in the central regions of the galaxy. As observed by Cannon et al. (2016), the galaxy channel maps reveal emission over approximately 40 km s⁻¹, peaking at 24 km s⁻¹ at 1.3 kpc. The average velocity field in the simulation exhibits a larger amplitude, although the emission by channels of velocity closely reproduces the galaxy’s velocity gradient and the galaxy’s emission across its entire extent. The higher amplitude of the average velocity field in the simulation must not be regarded as a limitation of the simulation, as the velocity field amplitude on the plane of the sky is heavily influenced by the line-of-sight. A small variation around the assumed line of sight results in a consistent reduction in the velocity field amplitude, but keeping the emission by channel less affected. Adopting i ≃ 40deg, as in our case, simply ensures a good compromise in terms of representing both the structural and kinematic components.

Fig. 5. Channel maps of H I gas distribution of NGC 5238 computed from the reference simulations of Sect. 4. The channels width is 3 km s−1, with the average channel velocity reported in each panel. The red curves show H I column densities of 0.5, 2 and 8 M⊙ pc−2 while, in black, the same 0.5, 2 and 8 M⊙ pc−2 isodensities but computed from the total gas distribution. The bottom panels show the overall gas-mass weighted line-of-sight velocity field in the simulation (left) and in the VLA observation (Cannon et al. 2016): blue colors and negative velocities correspond to the disc approaching arm, while red colors and positive velocities to the disc receding arm. The velocity field in the simulation extends down to column densities of 0.5 M⊙ pc−2. The panel size is 5.5 kpc × 5.5 kpc.

4.2. UGC 8760

In the case of UGC 8760, to facilitate the gradual approach of the satellite galaxy towards the host in a reasonable amount of time, we leverage dynamical friction which, according to our calculation, acts on a timescale that is, approximately, t_df ≃ 0.5 Gyr⁴. In addition to that, the truncation radius of the satellite is approximated as (Binney & Tremaine 2008), with M_sat, M_host, and d, respectively the total mass of the satellite, the mass of the host, and the distance between the two centers at the beginning of the simulation. The fact that the satellite’s r_t is larger than its size (Table 1) and way shorter than the initial distance d = 5 kpc ensures that any gravitational perturbation provided by UGC 8760 is progressively accounted for in the simulation and that it is not an artifact of the initial set up.

As the simulation starts, the satellite galaxy sinks towards the center in a time compatible with t_df: after half an orbit, taking less than 0.39 Gyr, the average distance between the host and the satellite has decreased to 2.4 kpc (first row, second panel of Fig. 4), while it is 1.6 kpc after 0.69 Gyr (first row, third panel of Fig. 4). At t = 0.69 Gyr the satellite galaxy is, in projection, approximately at the same distance as the old stellar overdensity observed in UGC 8760. In the simulation, the satellite is stripped of a considerable amount of dark matter in the outer regions. For a comparison, Fig. 6 shows the enclosed mass distributions of dark matter (top) and stars (bottom) of the satellite at the beginning of the simulation and after 0.69 Gyr. Here, the mass outside 3 kpc is significantly reduced from its initial value, with a fractional mass loss of approximately 21%. However, despite the decreased dynamical mass, the stellar component, which is very compact and concentrated at the center of the satellite, maintains a smooth and regular density distribution, remaining nearly unchanged compared to the initial one, despite its proximity to the host. In this configuration, the orbit of the satellite galaxy is almost co-planar with the plane of the sky. Therefore, the residual average line-of-sight velocity of the satellite with respect to the main galaxy is low, about 6 km s⁻¹.

Fig. 6. Top panel: dark matter mass distribution of the satellite of UGC 8760 as a function of the distance from the satellite center. The black-dashed line shows the mass distribution as in the ICs while the black solid after t = 0.69 Gyr. Bottom panel: same as in the top panel but showing the mass distribution of the stellar component on the satellite. The dashed blue line marks the position of the stellar component half-mass radius. The stellar mass distribution does not change from the ICs to 0.69 Gyr.

The galaxy H I density distribution is shown in Fig. 4 by the red contours. In the VLA observations of UGC 8760, the outer H I gas disc distribution is slightly bent, seemingly symmetric towards the northwest, with the outer distribution displaying a box-like structure. In the simulation, the UGC 8760 gas disc undergoes mild distortions in the outer regions due to the close passage of its satellite. However, after ≃0.69 Gyr, although the UGC 8760 gas distribution develops a similar outer boxy-like pattern, the bent structure is more prominently emphasized by the stellar distribution rather than the H I in the simulation. In the simulation, these stars serve to bridge the gap between UGC 8760 and the satellite. This bridge is evident in the RGB count map (Fig. 2, right panel) and in the simulations, but is not detectable in the g-band image. The reason for the absence of this star bridge in the g-band image is primarily attributed to its low surface brightness, coupled with significant contamination from background galaxies, which is absent in the mock image. RGB maps, on the other hand, can achieve deeper surface brightness levels as they enable the isolation of stars and the removal of a large portion of background galaxies.

Figure 7 shows the same channel maps of H I densities shown in Fig. 5, but for UGC 8760. The simulation effectively replicates the velocity field derived from H I observations of the galaxy, particularly in terms of the overall gradient and the position of the approaching and receding arms. However, as for NGC 5238, it is important to note that the observed velocity field exhibits a certain degree of perturbation and granularity that cannot be recreated by adiabatic simulations that do not account for stellar feedback phenomena. In this respect, the mock velocity field computed from the simulation does appear much smoother and regular than the observed one.

Fig. 7. Same as in Fig. 5 but showing the H I channel map of UGC 8760. Here, the channel width is 3.8 km s−1, the box size is 4 kpc × 4 kpc, while the observations are from Begum et al. (2008).

4.3. Additional simulations

In this section, we briefly present the outcomes resulting from the exploration of additional parameters in the ICs of the simulations of NGC 5238 and UGC 8760. These parameters were varied to investigate how they might influence the simulations’ ability to reproduce observational data. This sampling, however, covers only a small portion of the extensive parameter space, as a comprehensive exploration extends beyond the scope of this article.

In the case of NGC 5238, we have conducted simulations exploring the effects of varying the initial velocity of the satellite. Specifically, we tested velocities of v_x = −10 km s⁻¹ and v_x = −20 km s⁻¹, while keeping all the other orbital and galaxy models parameters the same. The two left-hand columns of panels of Fig. 8 show mock surface brightness maps of the NGC 5238 system from the two aforementioned simulations, corresponding to snapshots that reproduce the observations. The line of sight is approximately the same as that adopted in the mock images from the reference simulation, while the time from the beginning of the simulation is t = 0.84 Gyr and t = 0.72 Gyr, from left to right. Due to the varying initial velocities, the corresponding time also differs, being slightly greater for lower initial velocities. The maps also show H I isodensity contours in red, while the lower panels the corresponding line-of-sight velocity field. As can be readily observed, the outcome of the simulations closely resembles that of the reference simulation. Despite the decrease in initial velocity, the satellite is nonetheless completely stripped of gas due to ram pressure. The perturbations induced in both the stellar and gaseous discs of NGC 5238 remain nearly unchanged, indicating little dependence on this parameter. Therefore, as long as the accretion occurs approximately on a radial orbit, the proposed type of interaction seems to reproduce quite faithfully the observed properties of the system.

Fig. 8. First two columns of panels: stellar surface brightness maps (top) and line-of-sight velocity field (bottom) from two additional simulations of NGC 5238 discussed in Sect. 4.3. The columns, from left to right, show the simulations with initial velocities vx = −10 km s−1 and vx = −20 km s−1, respectively. The mock surface brightness and velocity field maps are obtained as in Figs. 4 and 5. Second two columns of panels: same as the first two columns, but for UGC 8760. In this case, instead of showing two different simulations, the columns show the same simulation but at different snapshots that both replicate the galaxy observed structure and kinematics. The merger mass ratio of the simulation is Mdm/Mdm, sat = 24.

As previously emphasized, these simulations do not represent the sole merger possibilities capable of reproducing the observed substructures in NGC 5238. As outlined in Sect. 3.3, the choice behind a radial interaction was motivated by the observed north-south directionality of the substructures, but it does not preclude other scenarios. For instance, it is reasonable to speculate that a flyby involving a satellite on a circular orbit, as in the UGC 8760 case, may similarly yield the disc of NGC 5238 to develop an outer spiral structure pattern ascribable to the north and south overdensities. From the parameters of the galaxy models presented in Table 1, it can be noted that, in many aspects, the UGC 8760 system is a scaled down version of the NGC 5238 system, with approximately half of the dark matter mass in host and satellite, with a stellar disc extended by half, and four time less massive. Indeed, despite the unfavorable line of sight, the spiral pattern also emerges in the UGC 8760 system, as illustrated in Fig. 4. Therefore, it is reasonable to speculate that a similar behavior may occur in the NGC 5238 system, provided appropriate scaling of the orbital ICs.

In the case of UGC 8760, the satellite’s stellar-to-dynamical mass ratio is only 0.5% within 3 kpc, effectively rendering stars as mere tracers of the total potential. Consequently, reducing the satellite’s stellar mass is unlikely to yield significant alterations in the simulation’s outcome after 0.69 Gyr, except for a decrease in luminosity, which can be compensated for by adjusting the mass-to-light ratio. As an alternative, we conducted simulations in which we maintained the same orbital parameters but varied the satellite’s dark matter mass, and vice versa. When keeping the orbital parameters fixed (i.e. the satellite’s initial position and velocity), the alteration in dark matter content was intended to explore interactions with different merger mass ratios. Specifically, in the first case we halved the dark matter mass, moving to a ratio M_dm/M_{dm, sat} = 24 instead of M_dm/M_{dm, sat} = 12 as in the reference simulation. Here, the gravitational perturbation of the satellite on the primary galaxy (UGC 8760) is significantly reduced. In the third and forth columns of panels of Fig. 8 we show two snapshots of this simulation where the observations are fairly well reproduced. Despite the substantially increased time required to reach UGC 8760, caused by the diminished mass and, thus, a less effective dynamic friction, the satellite’s stellar component remains relatively unperturbed. Moreover, the perturbations experienced by the stellar and gas discs of UGC 8760 are much less significant as well. Nevertheless, also this simulation appears suitable to reproduce the observed structure and kinematics of UGC 8760. In the second case, we tested ICs with purely radial orbits (as in the case of NGC 5238, where the satellite crosses the galaxy’s center, see Sects. 3.2.1 and 4.1), adopting two dark matter halos with varying amount of mass. Here, the gravitational interaction proves to be much more destructive to UGC 8760 and the coherence of its stellar and gas discs, potentially making such interactions far less likely. This exploration suggests that interactions featuring more circular orbits tend to foster gentler and less destructive outcome⁵, easily reproducing the observed structure of the UGC 8760 system, almost independently of the satellite dark matter mass within the range explored.

5. Conclusions

In this study, we presented the cases of NGC 5238 and UGC 8760, two dwarf galaxies observed as part of the SSH survey, exhibiting promising signs of recent interactions in deep LBT imaging. Stellar population studies based on these LBT observations revealed that the low-surface brightness and anomalous stellar substructures observed in the outskirts of the galaxies are predominantly composed of old stars, indicating limited association with sporadic and recent episodes of star formation. Through hydrodynamical N-body simulations presented in this work, we have shown that recent interactions with less massive satellite galaxies can reproduce reasonably well the structural and kinematic characteristics of both galaxies.

In both cases, the dwarf galaxies are represented by means of realistic ICs comprising of thick star and gas discs, embedded in a dominating dark matter halo, in agreement with state-of-the art measurements based on H I kinematics and stellar photometry. The simulations do not include star formation and feedback from supernovae explosions, implying that stellar particles trace only the kinematics of the old stellar populations of NGC 5238 and UGC 8760.

For the NGC 5238 system, the simulations suggest that a close interaction with a satellite galaxy whose stellar mass is 50 times less than that of the host, and with a merger mass ratio 1:10, can reproduce some of the observed structural irregularities in the galaxy outskirts. In the scenario proposed, the overdensities of old stars observed to the north and south of the galaxy are mainly composed by stars belonging to the disc of NGC 5238, distorted by the close interaction with its satellite. Similarly, our simulations reveal that interactions with a small spheroid with a merger mass ratio of up to 1:12, provide a compelling explanation also for the presence of the low-surface brightness companion situated to the north of UGC 8760. In terms of gas kinematics, our simulations fairly replicate the distribution of H I, as well as the overall velocity field in both cases. While a systematic exploration of the wide parameter space was not undertaken, the presented simulations suggest that the observed substructures emerge relatively rapidly after the simulation begins, potentially indicating a brief timescale for the interaction. Also, in the explored parameter space, in UGC 8760 we found that interactions with tangential orbit interactions produce more favorable outcomes compared to those with radial encounters, while radial interactions best reproduce the substructures of NGC 5238.

Within the SSH survey, NGC 5238 and UGC 8760 complement the well-known cases of DDO 68 (Annibali et al. 2016; Pascale et al. 2022) and NGC 3741 (Annibali et al. 2022), and will be accompanied by additional dwarf galaxies that will be introduced and analyzed in details in an upcoming paper Sacchi et al. (2024). All of these galaxies reside in low-density environments: within voids or on the periphery of well-known groups. Unlike dwarf galaxies in more densely populated regions, these isolated or low-density environments offer protection against disruptive interactions that, in a hierarchical scenario of structure formation, could otherwise occur with more massive hosts (Deason et al. 2015). Dooley et al. (2017) showed that for galaxies with a fixed stellar mass located relatively far from massive hosts like the Milky Way, the number of satellite galaxies orbiting around can be predicted using cosmological simulations. For instance, for a host with a stellar mass of 10⁷ − 10⁸ M_⊙ like NGC 5238 and UGC 8760, there is a high probability of being orbited by at least one satellite system with a stellar mass of ≃10⁵ M_⊙, with the number of expected satellites increasing as we consider even smaller masses (up to three for a stellar mass of 10⁴ M_⊙). Indeed an increasing number of low-mass systems have been discovered in the past years (Bellazzini et al. 2013; Erkal & Belokurov 2020; Carlin et al. 2021; Martínez-Delgado et al. 2022). In conclusions, environments like those surrounding NGC 5238 and UGC 8760, which lack nearby massive galaxies, offer odds of survival for small satellites significantly improved. This is precisely the type of environment that would favor the presence of satellites of satellites.

Acknowledgments

This paper is supported by the Italian Ministerial grant PRIN "Dwarf galaxies as probes of the Lambda Cold Dark Matter hierarchical paradigm at the smallest scales" 1.05.23.05.09. This paper is supported by the Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations. National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 “Italian Research Center for High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S (M4C2-19)” – Next Generation EU (NGEU).

References

All Tables

All Figures

Fig. 1. LBT images of NGC 5238 and UGC 8760. Left-hand panel: LBT g-band image of NGC 5238. The arrows show the positions of the northern and southern overdensities of stars while the red contours show lines of constant H I emission from VLA observations (Cannon et al. 2016) at the adopted distance. Right-hand panel: same as the left-hand panel but for UGC 8760. Here, the arrow shows the position of the northern overdensity of stars, while the VLA observations are from Begum et al. (2008). In both panels, the red contours correspond to 1, 2, 4, 8 M⊙ pc−2. North is up, east is to the left.

In the text

Fig. 2. Spatial distribution of old (age >  2 Gyr) RGB stars and young (age ≲ 100 Myr) bright blue stars in NGC 5238 (left) and UGC 8760 (right). RGB maps are background-subtracted, and only pixels above a level of 3 times the background standard deviation (given by the vertical black segment in the colorbars) have been set to non-zero to mask density peaks due to red background galaxies. The blue points are individual young stars. H I contours at 1, 2, 4, 8 M⊙ pc−2 (same as in Fig. 1) are superimposed to the star count maps. North is up, east is to the left.

In the text

Fig. 3. Dark matter density distribution of NGC 5238 used in the simulations described via three different density models (different colors are associated to the three dark-matter models as indicated in the legend). The different models are linked via Eq. ((9), see text for details). The vertical dashed lines show the NFW, Hernquist and Dehnen scale radii using the same color coding, while the vertical blue solid line is the NFW virial radius.

In the text

Fig. 4. Comparison between simulations and observations. First row of panels: stellar surface brightness maps from the reference simulation of NGC 5238 and the corresponding LBT observation. The left panel shows the ICs of the simulation, while the middle panels are subsequent snapshots (as reported in each panel). The third panel corresponds to the time when the simulation reproduces the data, shown for comparison in the rightmost panel. In the simulated maps, the pixel size is 0.225 arcsec × 0.225 arcsec (as for the LBC camera), and the maps have been convolved with a Gaussian beam of 1 arcsec FWHM to match the typical seeing during observations. To generate the maps, we adopted constant mass-to-light ratios and we set the surface brightness limit to 30 mag arcsec−2, to match the depth of the data. Blue lines show the orbit of the satellite, obtained by joining the centers of mass of its dark matter halo computed with the shrinking sphere method (Power et al. 2003). The color gradient indicates the direction of the orbit itself: the satellite moves from where the orbit is darker to where it is lighter. The red contours are levels of iso-H I column density, specifically 2k M⊙ pc−2, with k = 0, …, 3. Second row of panels: same as the top row but for the UGC 8760 system. For details on the simulation setup see Sect. 3.3 and Tables 1 and 2. In the top row, third panel, we highlight in blue the stellar particles belonging to the satellite of NGC 5238 as a comparison (see text for details).

In the text

Fig. 5. Channel maps of H I gas distribution of NGC 5238 computed from the reference simulations of Sect. 4. The channels width is 3 km s−1, with the average channel velocity reported in each panel. The red curves show H I column densities of 0.5, 2 and 8 M⊙ pc−2 while, in black, the same 0.5, 2 and 8 M⊙ pc−2 isodensities but computed from the total gas distribution. The bottom panels show the overall gas-mass weighted line-of-sight velocity field in the simulation (left) and in the VLA observation (Cannon et al. 2016): blue colors and negative velocities correspond to the disc approaching arm, while red colors and positive velocities to the disc receding arm. The velocity field in the simulation extends down to column densities of 0.5 M⊙ pc−2. The panel size is 5.5 kpc × 5.5 kpc.

In the text

Fig. 6. Top panel: dark matter mass distribution of the satellite of UGC 8760 as a function of the distance from the satellite center. The black-dashed line shows the mass distribution as in the ICs while the black solid after t = 0.69 Gyr. Bottom panel: same as in the top panel but showing the mass distribution of the stellar component on the satellite. The dashed blue line marks the position of the stellar component half-mass radius. The stellar mass distribution does not change from the ICs to 0.69 Gyr.

In the text

	Fig. 7. Same as in Fig. 5 but showing the H I channel map of UGC 8760. Here, the channel width is 3.8 km s−1, the box size is 4 kpc × 4 kpc, while the observations are from Begum et al. (2008).
In the text

Fig. 8. First two columns of panels: stellar surface brightness maps (top) and line-of-sight velocity field (bottom) from two additional simulations of NGC 5238 discussed in Sect. 4.3. The columns, from left to right, show the simulations with initial velocities vx = −10 km s−1 and vx = −20 km s−1, respectively. The mock surface brightness and velocity field maps are obtained as in Figs. 4 and 5. Second two columns of panels: same as the first two columns, but for UGC 8760. In this case, instead of showing two different simulations, the columns show the same simulation but at different snapshots that both replicate the galaxy observed structure and kinematics. The merger mass ratio of the simulation is Mdm/Mdm, sat = 24.

In the text