The phase-space distribution of the M 81 satellite system

The spatial distribution of dwarf galaxies around their host galaxies is a critical test for the standard model of cosmology because it probes the dynamics of dark matter halos and is independent of the internal baryonic processes of galaxies. Comoving planes of satellites have been found around the Milky Way, the Andromeda galaxy, and the nearby CenA galaxy, which seems to be at odds with the standard model of galaxy formation. Another nearby galaxy group, with a putative ﬂattened distribution of dwarf galaxies, is the M81 group. In this paper, we present a quantitative analysis of the distribution of the M81 satellites using a Hough transform to detect linear structures. Using this method, we conﬁrm a ﬂattened distribution of the dwarf galaxies. Depending on the morphological type, we ﬁnd a minor-to-major axis ratio of the satellite distribution of 0.5 (all types) or 0.3 (dSph), which is in line with previous results for the M81 group. Comparing the orientation of this ﬂattened structure in 3D with the surrounding large-scale matter distribution, we ﬁnd a strong alignment with the local sheet and the planes of satellites around the Andromeda galaxy and CenA. Furthermore, the satellite system seems to be lopsided. Employing line-of-sight velocities for a subsample of the dwarfs, we ﬁnd no signal of corotation. Comparing the ﬂattening and motion of the M81 dwarf galaxy system with TNG50 of the IllustrisTNG suite we ﬁnd good agreement between observations and simulations, but caution that i) velocity information of half of the satellite population is still missing, ii) current velocities mainly come from dwarf irregulars clustered around NGC 3077, which may indicate an infall of a dwarf galaxy group, and iii) some of the dwarfs in our sample may be tidal dwarf galaxies. From the missing velocities, we predict that the observed frequency within IllustrisTNG may still range between 2 to 29%. Any ﬁnal conclusions about the agreement or disagreement with cosmological models needs to wait for a more complete picture of the dwarf galaxy system.


Introduction
The distribution and motion of dwarf galaxies around the Milky Way (Lynden-Bell 1976;Metz et al. 2008;Pawlowski et al. 2012;Taibi et al. 2023) and the Andromeda galaxy (Koch & Grebel 2006;McConnachie & Irwin 2006;Ibata et al. 2013) has sparked an on-going debate on whether they pose a challenge to the current Λ+Cold Dark Matter (ΛCDM) standard model of cosmology (see e.g.Kroupa et al. 2005;Libeskind et al. 2005Libeskind et al. , 2007;;Pawlowski et al. 2012;Ibata et al. 2014b;Cautun et al. 2015;Pawlowski & Kroupa 2020;Sawala et al. 2022 and the reviews by Pawlowski 2018Pawlowski , 2021)).The structure and kinematics of satellite systems pose a strong test for ΛCDM because they are driven by the gravitational dynamics on scales of hundreds of kpc, so they do not depend on internal baryonic processes (Pawlowski 2018;Müller et al. 2018a).This has pushed several teams to seek for similar structures outside the Local Group (Ibata et al. 2014a;Müller et al. 2018b;Paudel et al. 2021;Heesters et al. 2021).Clear evidence was found for our neighbor Cen A and its satellite system, showing a statistically significant correlation in phase-space (Tully et al. 2015;Müller et al. 2016Müller et al. , 2018aMüller et al. , 2021;;Kanehisa et al. 2023b).For the Sculptor group, Martínez-Delgado et al. (2021) pointed out a flattened distribution of satellites, as well as some signs of coherent motion, but the system needs additional follow-up observations of their velocities to assess the situation.For other systems in the Local Volume -a sphere of 10 Mpc around our point of viewclaims of flattened distributions have been made (Müller et al. 2017(Müller et al. , 2018b)).Intriguingly, most flattened structures seem to be aligned with the local cosmic web (Libeskind et al. 2015(Libeskind et al. , 2019)).This may point towards a common formation scenario for these structures.
A nearby galaxy group with a well studied satellite population is the M 81 group of galaxies at a distance of ≈3.7 Mpc (Ferrarese et al. 2000;Chiboucas et al. 2013).The central region of the group is characterized by three major galaxies (M81, M82 and NGC 3077) that are interacting with each other, as demonstrated by a complex HI network of filaments, tidal tails and candidate tidal dwarf galaxies (Yun et al. 1994;Chynoweth et al. 2008;Weisz et al. 2008).Due to dynamical friction, it has been argued that such compact arrangements of three galaxies are unlikely be found in a ΛCDM cosmology (Oehm et al. 2017).
The M 81 group has been a target for dwarf galaxy searches by different teams.Chiboucas et al. (2009) used the Canada France Hawaii Telescope (CFHT) telescope to survey a large field of 65 deg 2 , discovering 22 dwarf galaxy candidates, of which 14 were confirmed based on Hubble Space Telescope (HST) follow-up observations (Chiboucas et al. 2013).More recently, the group has been studied with the Hyper Suprime Cam (HSC).Okamoto et al. (2015Okamoto et al. ( , 2019) ) surveyed a field of 6.5 deg2 , discovering two additional faint dwarf galaxies.Finally, Bell et al. (2022) combined archival HSC data (seven ≈1.5deg 2 fields) and found six dwarf galaxy candidates.
In the seminal work by Chiboucas et al. (2013), the authors noted that the satellite structure seemed to be flattened, with an on-sky rms thickness of 61 kpc at the distance of M 81, when considering the dwarf spheroidals of the satellites, which is also well-aligned with the local sheet.In this work, we aim to provide a quantitative analysis of the spatial and kinematic properties of the M 81 group and a first -but not final -comparison to ΛCDM simulations of this system.In Sec. 2, we discuss the data and methods we are using, in Sec. 3 we provide a discussion on the distribution and motion of the satellite system, in Sec. 4 we compare the observations to cosmological simulations and discuss major caveats, and in Sec. 5 we provide a summary and conclusions.We use the catalog of group members from Chiboucas et al. (2013), which includes their own HST follow-up observations of the dwarf galaxy candidates found in a CFHT survey of 65 deg 2 (Chiboucas et al. 2009), as well as literature data.We consider galaxies with < −18.0 in absolute magnitudes in the rband as satellites.This effectively removes the two gravitationally dominant galaxies M 81 and M 82 from the list, but still includes the bright galaxies NGC 2976 (with M r = −18.0mag) and NGC 3077 (with M r = −17.8mag).M 81 and M 82 have stellar luminosities (derived from the K s band) of 10 10.95 M ⊙ and 10 10.59 M ⊙ (Karachentsev et al. 2013), respectively.To date, the Chiboucas et al. (2013) data represents the most complete survey of the M 81 group and its surroundings, going beyond the second turnaround radius of M 81 (which is ≈230 kpc).Because Chiboucas et al. (2013) do not compile the distances of the liter-ature data, we have taken the values from the online version of the Local Volume catalog (Karachentsev et al. 2013) 1 .A recent deep survey by Bell et al. (2022) uncovered another six ultrafaint dwarf galaxies, which are clustered around NGC 3077.We refrain from adding these to our list of satellites, because their HSC survey covered only 50 kpc around M 81, which would bias any study of the distribution of dwarf galaxies in this group.In Table 1 we provide the list of dwarf galaxies used in our study.

Data and Methods
How may the survey and its footprint bias the study of the dwarf galaxy population?The M 81 group lies in a region in the sky where galactic cirrus is dominant and this may cause trouble with assessing the distribution of the satellites.However, Chiboucas et al. (2013) made extensive tests with artificial star detection to probe their selection criteria and found that their detection of dwarf galaxies is not biased by cirrus (see their Fig.33 for the cirrus overlaid on top of their survey field).It is evident that, even in regions with strong cirrus, dwarf galaxy candidates were detected.Another issue may arise from the survey footprint.The flattened distribution found by Chiboucas et al. (2013) is aligned with the diagonal of the footprint, see Fig. 1.Because the diagonals will maximize the radial distance for which dwarfs can be found, it might prefer finding linear structures along these lines.To mitigate that, we would need to extend the survey footprint and see where the dwarf galaxy detection drops to the background.However, there seems to be a drop of dwarf galaxy detections towards the border of the survey footprint, which indicates that the dwarf galaxy population is sampled well enough towards the edges.We note though that some known galaxies like NGC 2403 or UGC4483 are generally considered to be members of the M 81 group (Karachentsev & Kashibadze 2005), but are well outside of the survey footprint, with projected separations of ∼15 deg (corresponding to roughly 1 Mpc).Here, we do not consider them as satellites of the M 81 group.
To study the flattening of the satellite distribution, we employed the Hough transformation as described in Heesters et al. (2021).In short, the Hough transform (Hough 1959(Hough , 1962) ) was originally developed as a tool to detect straight lines and other simple shapes in digital images.The principle of the method is a voting system in which all points in an image decide on the bestfit parameters from a predefined set of slope and intercept pairs.For every data point we define a range of lines that cross the point from different directions.With these lines the point votes for the best fit that should describe the overall linear distribution.Since the image space has no bounds for slope and intercept, the lines are transformed to the so-called Hough space via the relation ρ = x cos(θ) + y sin(θ). (1) Here, x and y are the coordinates of the data point in the image space.The parameter ρ is the orthogonal distance from the origin to the line passing through the data point.Finally, θ is the angle between the ρ vector and the x-axis.In this space, the pa- Fig. 2: Left: The M81 group satellite system in relative coordinates, with M81 (star) at the origin and indicating projected separations at the distance of M 81 in kpc.M 82 is shown as an empty star.The dwarf satellites of different morphologies are shown as circles and squares.We emphasize the abundance of dwarf spheroidals (dSphs) as black circles compared to other morphologies as gray squares.The best-fit line is determined via the Hough method, which finds all but four satellites (upper left corner) to be members (circled gray) of a flattened structure.The rms dimensions of this structure are indicated with black arrows.Top: the results using all morphological types of satellites, Bottom: only using the dSphs.Right: The M 81 satellite system in the Hough parameter space.The parameter ρ plotted on the y-axis is the normal distance from the fit line to the origin, while θ is the angle between this ρ vector and the x-axis in the image space.Each curve corresponds to a series of possible parameter pairs for one dwarf.The curves are z-score normalized along each axis, i.e., they have a mean of zero and a standard deviation of one.The red dot indicates the parameter pair at the center of the densest crossing region (grey circle), which is optimized to maximize the number of structure members and the flatness of the system simultaneously.The curves which do not cross the grey region stem from the declared outliers in the scatter plot left.single crossing point in this parameter space.For scattered data points, however, the curves no longer cross in a single point.For this case we identify the region with the highest over-density of crossing lines and adopt a variable search area, allowing us to probe different scales.We optimize the area such that the struc-ture flatness and the number of voting members are maximized simultaneously.This method allows for an educated estimation of the members of a potentially correlated satellite structure, i.e. a plane close to seen edge-on, in the absence of three dimensional information.The success of the method has been demon- Fig. 3: The satellite system of M 81 in Cartesian Supergalactic coordinates, shifted to the center of the satellite distribution.The dSphs are indicated in dark grey dots, the dIrr as squares in a lighter grey, and the main galaxies M 81 and M 82 as white stars.The thin black lines indicate the distance uncertainties.The top panels represent an edge-on view of the plane, the bottom panels a face-on view.In the bottom right panel the velocities are indicated as thick lines, as well as the color, with red indicating galaxies moving away from us in respect to the mean velocity of the system, and blue towards us.The color-coding is the same as in Fig 5 .strated on the M31 system (Heesters et al. 2021), for which only about half of the satellites appear to be part of a phase-space correlation (Ibata et al. 2013).Details about this fitting technique can be found in Heesters et al. (2021).

Satellite distribution
In this section we look into the spatial and kinematic distribution of the satellite system of M 81.

Planes of satellites
In Figure 2 we present the 2D distribution of the galaxies in the M 81 group and the morphological types of the satellites.A simple total least square (TLS) fit of the satellite system in 2D reveals an axis ratio of b/a = 0.79 (semi minor axis b = 91 kpc, semi major axis a = 115 kpc), but this approach does not consider the possibility that there may be outliers not belonging to a flattened structure, so artificially increasing the b/a ratio.The Hough method identifies a significantly flattened structure along the diagonal of the field.Out of the 34 satellites, 30 follow this elongated structure.Of the four satellites (HS117, d1028+70, DDO82, d1042+70) not being part, only one is a dwarf spheroidal (d1042+70).Removing these four outliers, we measure a minor-to-major axis ratio b/a = 0.50 (semi minor axis b = 61 kpc, semi major axis a = 122 kpc), which is similar to what is found for Cen A (b/a ≈0.5, Müller et al. 2016), but spatially less flattened than the Local Group planes (we note that the Local Group planes are studied in 3D due to better distance and thus 3D position accuracy).It is also consistent with the rms width estimated by Chiboucas et al. (2013).If we repeat the steps (i.e.Hough transformation and removing the outliers) for the dSph only, we get b/a = 0.34 (semi minor axis b = 42 kpc, semi major axis a = 124 kpc), with 16 out of 19 dSph belonging to the flattened structure.However, two dwarf galaxies -KK77 and F8D1 -which were previously considered in the Hough fitting, are now not found to be part of the flattened structure.This is interesting, because F8D1 has a disrupted profile exhibiting a tidal tail (Žemaitis et al. 2023) which is aligned approximately along the minor axis of the flattened structure.If we consider this shape as a tracer of the motion of the dwarf, it is moving out from the planar structure.If we consider the entire dSph sample in a TLS fit, we measure an axis ratio of b/a = 0.62 (semi minor axis b = 73 kpc, semi major axis a = 117 kpc).That the flattening is higher for the dSph population compared to the total population is in line with the finding of Chiboucas et al. (2013), who suggested there exists a flattened structure consisting of the dwarf spheroidals in the M 81 group.Fig. 4: The orientation of the poles of planes-of-satellites with respect to the surrounding large-scale matter distribution.The eigenvector ê1 corresponds to the direction of fastest collapse and is pointing towards the local void.The three lines (black, red, and blue) correspond a cone of 60 degrees off from the eigenvectors.
How is the flattened structure oriented in respect to the local universe?One way of answering this question is to examine how the flattened structure is aligned with respect to the surrounding large-scale matter distribution, as characterized by the velocity-shear tensor, specifically its principle eigenvectors (e.g Hoffman et al. 2012).The eigenvector corresponding to slowest collapse, (ê 3 ), defines the spines of cosmic filaments whereas the axis of greatest compression (ê 1 ) is associated with the normal of cosmic sheets (see Libeskind et al. 2018, and references therein).Libeskind et al. (2015Libeskind et al. ( , 2019) ) performed such an analysis on the known dwarf galaxy planes (i.e.MW, M31, Centaurus A) in 3D using a quasi-linear reconstruction of the Local Universe from CF2 (Hoffman et al. 2018;Tully et al. 2013) and found that many of their normals tend to be closely aligned with ê1 , the eigenvectors of the shear tensor corresponding to the axis of greatest compression.To test whether the M 81 satellite population is aligned with any of these vectors, we can transform the system into supergalactic coordinates (de Vaucouleurs 1956;Lahav et al. 1998).In the supergalactic coordinates the Local Sheet -including the Local Group, NGC253, Cen A, and M 81 -is aligned with the supergalactic SGX-SGY plane, which corresponds to the ê2 -ê 3 plane.The satellite system in in 3D supergalactic coordinates is presented in Fig. 3.We studied the alignment of the satellite system with a principal component analysis (pca) using the python package sklearn (Pedregosa et al. 2011).The pca is a method to find the eigenvectors of a data set.Here, we use it to estimate the direction with the smallest variance, which for a flattened distribution of data points corresponds to its normal vector.Considering all dwarfs -including the four outliers which we deem not be part of the plane-of-satellites -the normal for the planar structure in supergalactic coordinates is [−0.1314,+0.0374, +0.9906].The normal including all dwarf galaxies in the flattened structure is [−0.0899,+0.0772, +0.9930].The angle between the two normals is only 3 degrees, so they align well and can be regarded as equal.The position of the normal in respect to the shear tensor is given in Fig. 4.
The flattened satellite structure of M 81 has its normal aligned close to ê1 of the velocity-shear tensor defined by the large-scale matter distribution in the local universe.This is consistent with Chiboucas et al. (2013) who noted "an apparent flattening of the distribution of gas-deficient systems to the supergalactic equatorial plane".However, we note that employing the Hough transformation we find both an alignment considering dSphs only, and dSphs and dIrrs together.The alignment of M 81's satellite system with the Local Sheet is also for the satellite systems of M 31 and Cen A (Libeskind et al. 2015(Libeskind et al. , 2019)).However, for the best studied plane-of-satellite system, this is, around the Milky Way no such alignment can be found.This may indicate that these structures have different origins.
It is interesting to note that the four dwarf galaxies (HS117, d1028+70, DDO82, d1042+70) not being part of the main structure form some kind of elongated spur as well, which is going parallel to the main structure.It has a minor-to-major axis ratio b/a = 0.28 (semi minor axis b = 13 kpc, semi major axis a = 46 kpc).There is a gap of 222 kpc between the best-fitting plane from these four dwarfs and the main flattened stucture.This gap is well visible in 2D in Fig. 1 and in 3D in Fig. 3.

Lopsidedness
Do we find evidence for an asymmetric distribution of the satellite system?This can be studied using different metrics testing for lopsidedness in the group.Here, we use two simple approaches, these are, a) the distance of the centroid from M 81, and b) a hemisphere approach.
In the first approach, we test whether the center of the satellite distribution coincided with the dominant galaxy M 81.We define the center or centroid of the Hough fit as Here N is the of structure members and r i are the member coordinates.The centers for the two Hough fit lines (using all dwarfs and using only the dSph) are presented in Fig 2 .They are separated from each other by 23 kpc.The centroids for both samples are off from the position of M 81 by ≈40 kpc.Is this significant?To test this, we run a Monte Carlo simulation where we generated for each iteration 30 random positions uniformely drawn within the box defined by the plane-of-satellites (i.e. with b = 61 kpc and a = 122 kpc).Finding the center of the point cloud for each iteration, we estimate that an offset of 40 kpc happens in only 0.2% of the cases, in other words, is significant at the 3σ level.This may indicate that the dwarf satellite population is not symmetrically distributed around M 81, which could be a physical effect caused by the on-going interaction between M81, M82 and NGC3077.It could, also arise from an infalling group of dwarfs, see the subsection 3.3 for such hints.
Another way to test lopsidedness is splitting the distribution of the satellites in two hemipsheres.Such lopsidedness has been found within the Local Group (Conn et al. 2013), as well as statistically in galaxy groups in the SDSS survey (Libeskind et al. 2016).A study by Pawlowski et al. (2017) of these SDSS systems in ΛCDM cosmological simulations has found that this anisotropy is not in tension with observations.Kanehisa et al. (submitted) however finds the M31 satellite system's lopsidedness to be significantly in tension with ΛCDM simulation expectations.Therefore it is interesting to see whether we find this in the M 81 group or not.Taking the minor axis of the satellite distribution as a separation line, which we fix at M 81's position (the brightest galaxy in the group), we find that 11 satellite galaxies are on one side and 23 on the other.Assuming a Bernouilli distribution for a satellite being on one side or another (i.e. a fair coin flip), we estimate that finding 11 or less satellites on one side has a probability of 5.8% (note that we have to consider two cases, i.e. finding 11 or less and finding 23 or more in a Bernouilli experiment).This is not passing the 3σ threshold and we therefore do not find it to be significant.

Kinematics
Velocity measurements are available for 19 dwarf galaxies, of which 16 are members of the reported flattened structure.This is comparable to available line-of-sight data for the Andromeda galaxy and Cen A. However, these galaxy groups are more simple in terms of their dynamics, because they are both dominated by a single giant galaxy.On the other hand, the M 81 group consists of at least two main galaxies -M 81 and M 82 -as well as several massive dwarfs such as NGC 2976 and NGC 3077.This makes an analysis of the dynamics of the group not straightforward.For example, M 81 has a line-of-sight velocity of -38 km s −1 (Dickel & Rood 1978) and M 82 of 183 km s −1 (Appleton et al. 1981).The satellites' mean line-of-sight velocity is 20 km s −1 , which is offset by ≈60 km s −1 from M 81's velocity.In comparison, the difference between the velocity of Cen A and the mean of its satellite system is 1 km s −1 (Müller et al. 2021).If we take a simple approach and use the mean satellite velocity of the dwarfs that are found to be part of the structure (29 km s −1 ) as an anchor which coincides with the luminosity weighted velocity of M 81/M 82 (29 km s −1 ), we can produce a position-velocity (PV) diagram for the flattened structure, see Fig. 5.In general, we expect that a fully co-rotating system would occupy two out of the four quadrants on opposing sides, and a fully pressure-supported system all four quadrants.For the M 81 system, a clear signal of corotation is not obvious.The two opposing quadrants each populate 9 and 7 satellites, respectively.In this plot, we used the 2D center of the flattened dwarf system as found with the Hough fit as the zero point for the distance.There may be other choices to be made, e.g.either using M 81 as the center, or a weighted mean between M 81 and M 82.Using the stellar luminosities as a proxy for the latter, we find only a minor difference (11 kpc) to the position of M 81, which is negligible.Fixing the center at M 81, most dwarf galaxies would populate the two left quadrants (it would shift the vertical line to the right in Fig 5).
Assuming the flattened satellite system of M 81 is observed edge-on, we can study the line-of-sight velocities of the dwarfs from a face-on point-of-view of the system.Such a plot is presented in Fig. 3, bottom right panel.The line-of-sight-velocities will be in the plane of this panel.There are a few interesting points to note.Almost all velocities are coming from dwarfs on one side, which coincides with the position of NGC 3077.And most of these are dIrrs.Could it be that at least a few of these dwarfs were a distinct group including NGC 3077 and are falling in?Because the HI in dwarf galaxies should get stripped when accreted by a massive galaxy such as M 81 (Spekkens et al. 2014), it is possible that we are observing their (first) infall.Moreover, between M 81 and NGC 3077, there are five candidate tidal dwarf galaxies (d0959+68, BK3N, Garland, A0952+69, and Holm IX) that are embedded in the HI tidal material, so they may further complicate the dynamical interpretation of the system (see, e.g., Lelli et al. 2015 for bona-fide tidal dwarf galax-ies in other interacting systems).Out of these five candidates all but d0959+68 have velocity measurements.Two would be co-moving as expected for a co-rotating plane-of-satellites, two would be off.Tidal dwarf galaxies may offer a possible formation scenario for the creation of planes-of-satellites, with numerical fly-by models showing that up to 50% of the tidal dwarf galaxies may end up on counter-rotating orbits (Pawlowski et al. 2011).However, as long as we do not have proper motion measurements for these dwarfs, we cannot assess whether they are truly co/counter-rotating, only that they may be consistent with such a motion (see e.g.Kanehisa et al. 2023b).
With the current incomplete picture of the satellite system (i.e. the missing velocities of half the dwarf galaxies) we cannot draw firm conclusions about the dynamical of the system.But currently, there is no strong evidence for co-rotation of the satellite system.

Comparison to Illustris-TNG50
To probe the significance of M81's flattened satellite distribution, we use data from the IllustrisTNG suite (Pillepich et al. 2018) of publicly available, large-volume hydrodynamic simulations which adopts cosmological parameters from Planck (Planck Collaboration et al. 2016): ega Λ = 0.6911, ega m = 0.3089, and h = 0.6774.To best resolve the fainter dwarf population in the M81 group, we make use of the highresolution TNG50-1 run -spanning a simulation volume of length 51.7 Mpc with dark and gas particle resolutions of M DM = 4.5 × 10 5 M ⊙ and M gas = 8.5 × 10 4 M ⊙ respectively.
We define M81's simulated analogs as having a stellar mass between 5 × 10 10 M ⊙ and 2 × 10 11 M ⊙ while also requiring an associated halo mass M 200 between 5×10 11 M ⊙ and 2×10 12 M ⊙ (as motivated by M81's dynamical mass estimate of 1.17 × 10 12 M ⊙ from Oehm et al. 2017).Analogs must also be isolated from other dwarf associations -specifically, we reject all systems with a companion halo 250 − 1000 kpc away with a halo mass above 0.25M 200 .This gives us our sample of M 81 analogs.Additionally, we inspect each M81 system analog for a corresponding M82 analog, defined as a subhalo within 200 kpc of the host galaxy with a stellar mass above 0.25 M * ,M81 -thus accommodating the M81-M82 stellar mass ratio of ∼ 0.43 (Karachentsev et al. 2013).In the first step we analyze the full analog sample.In the second step, we restrict the sample by only considering systems that also feature an M82 analog and repeat our analysis.
Each M81 analog is mock-observed from 10 isotropically distributed directions with corresponding distances drawn from M81's distance modulus of µ = 27.84 ± 0.09.Since M81's satellite population is complete within a projected 230 kpc radius to a limiting magnitude of −10 in the r band, we identify all subhaloes within a projected radius of 250 kpc with a depth of 500 kpc with respect to their host galaxy.The subhaloes are ranked by their stellar mass (then dark mass once no stellar particles are available), and the 34 most massive are taken as the realization's satellite sample.Finally, distance errors are drawn from M81's observed satellites and randomly associated to the simulated satellites along the mock-observed line-of-sight.
In order to assess the prevalence of similarly flattened structures in simulations, we perform the described Hough fitting procedure on the 5530 (340 featuring an M82 analog) M81 analog systems.For this, we use the same parameter search area which was found to optimize the flattening and the number of structure members in the observed M81 satellite system.We compare both the flattening and the degree of kinematic correlation, i.e., the number of corotating satellites.The Hough fitting rejects 0 to 12 outliers (compared to four for the observations).The mean number of rejected outliers is 4 (as is the number of observed rejected outliers) with a standard deviation of 1.8.We have looked into a potential correlation between the number of outliers and the axis ratio of the resulting structures.In general we would expect that having a smaller sample (i.e. more outliers) would lead to a flatter structure (Pawlowski et al. 2017).However, there is no such correlation when we compare the distributions of b/a ratios for the different samples with outliers between 0 to 12.This is due to how we set-up the search radius during the Hough fitting on the simulated systems.We require it to be the same as for the observed system (maximizing flatness and number of objects considered in the fitting).If we would not require this, we would indeed find a correlation between the thickness and the sample size.
Because not all satellites have available velocity measurements, we assess the fraction of observed structure members that has such estimates (53%).When finding a linear substructure among outliers, the Hough fitting technique is not forced to include the same number of satellites in the simulations that are found to be part of the structure (30).As discussed before, the method selects a variable number of dwarfs, which is very similar to the number of structure members in observations.Because the structure populations are slightly variable in every simulated analog system, we select 53% of the velocities from the structure member satellites in the simulated systems.We use the satellites with the highest dark matter mass (as a proxy for their stellar mass) for this and use these to determine the fraction of corotating satellites along the major axis of a given structure.Similarly to the data analysis, we use the mean velocity of these satellites and the minor axis of their spatial distribution as anchor points for their corotation.The described comparison results in three measures of significance or p-values: the frequency of systems that match or exceed the observed axis radio, the degree of kinematic correlation, and the simultaneous occurrence of both properties.We find the corresponding p-values for the full sample and the sample featuring an M82 analog, respectively: p f lat = 0.206 (0.294) , p kin = 0.724 (0.735) and p both = 0.151 (0.212).We thus estimate that the observed structure occurs in approximately 15.1% (21.2%) of M81-like systems.The results for both samples are illustrated in Figure 6.
Since not all M81 satellites have available velocity measurements at this time, the degree of corotation in the structure members is still unknown.Inspecting the structure members we determine that the minimum number of possible corotating structure members is 15 (50%) and the maximum number is 23 (76.6%).In a second comparison between the data and the simulated systems, we consider the velocities of all structure members in simulations to determine the fraction of corotation.We study all possible outcomes in terms of corotation given the missing velocities in the data and compare all possibilities with the expectations from the simulated analog systems.This results in a p-value for every possible fraction of corotation and allows us to make a prediction on how many satellites are expected to corotate in order to be consistent with ΛCDM.The results of this comparison are illustrated in Figure 7.The fractions of simulated systems for which both the flatness and degree of corotation match or exceed the possible values with future observations range from 2.1% to 29.4%.This illustrates that the consistency with the expectations from cosmological simulations is still uncertain due to the missing velocities and therefore unknown degree of corotation in the reported structure.

Summary and conclusions
The plane-of-satellites problem has been identified as one of the biggest challenges to cosmology on small scales (Sales et al. 2022).Therefore, it is imperative to study different galaxy groups and test whether the motion and distribution of satellite galaxies follows predictions from cosmological simulations or not.In this work, we quantified previous claims of a flattened structure of dwarf galaxies in the compact M 81 group.We find that 30 out of 34 dwarf galaxies follow a flattened structure with a minor-to-major axis ratio of b/a = 0.5, which is similar to that of Cen A (Müller et al. 2019) and of early-type galaxies in the nearby universe (10-45 Mpc, Heesters et al. 2021) from the MATLAS survey (Habas et al. 2020;Poulain et al. 2021).If we consider only the dwarf spheroidals, the flattening increases to a ratio of 0.34.This is still thicker than the best studied planes around the Milky Way and the Andromeda galaxy (c/a < 0.2).However, for the Milky Way and Andromeda systems, accurate distance estimates are available, so the axis ratio can be estimated in 3D.Here, the thickness was measured in 2D.If the planar structure is significantly inclined with respect to the line of sight, the projected 2D thickness will increase.The difference between the thickness of the dSph-only and the dSph+dIrr samples may also arise from the lower numbers of tracers.The fewer tracers, the more likely it is to get a thinner spatial distribution (Pawlowski et al. 2017).
How consistent is the flattened structure with cosmological predictions?We have compared the on-sky flattening and motion of the satellites in the M 81 group with TNG50 from the IllustrisTNG suite.In contrast to the Local Group and Centaurus system, we find agreement within 2σ with ΛCDM expectations.This is an interesting result because M 81 is in a different dynamical state than the Milky Way with its quiet merger history, and Andromeda and Cen A with putative major mergers 2-3 Gyr ago (Hammer et al. 2018;Wang et al. 2020).
Considering the possible coherent motion of satellites, which is the crucial property of satellite systems to test ΛCDM, we find no correlation for the M 81 satellites.This is in contrast to the Milky Way, the M 31, and Cen A satellites, which all seem to follow a coherent motion around their hosts.Why could this be?
The M 81 group is a compact group which is currently undergoing a major merger and is therefore not in dynamical equilibrium.This is evident from strong signs of interactions between M 81, M 82, and NGC3077, especially in HI (Chynoweth et al. 2008).Even if there was a co-rotating plane-of-satellites around M 81 before, such an encounter may destroy any signs of corotation.Simulating the effect of major mergers on the creation or destruction of planes-of-satellites could help understand the observations of the M 81 group.In this respect, simulations of the major merger of Cen A hints toward a connection between the merger and the plane-of-satellites (Wang et al. 2020), but more detailed simulations are needed to reproduce the Cen A system.However, studies using TNG100 from the IllustrisTNG suite of cosmological simulations show a negligible contribution of the major merger history on the phase-space coherence of satellite galaxies (Kanehisa et al. 2023a).Under this context, we would neither expect the formation or destruction of a corotating satellite system around M 81.
Another explanation could be that we are seeing a mixing of several dwarf galaxy systems.It may be that a group of dwarf galaxies is falling in together.This is hinted at by the clustered distribution of dIrrs towards one side of M 81.In addition, at least 4 out of 19 satellites for which we have velocity information (d0952+69, BK3N, HolmIX, and Garland) are candidate tidal dwarf galaxies, being embedded in the HI tidal material between M 81 and NGC 3077.Two of those are reducing the phase-space signal.Before ruling out a coherent motion of satellites, it is imperative to sample the full dwarf galaxy population around M 81, especially dwarf spheroidal galaxies for which we can rule out a recent infall or a recent tidal formation.There are 15 more dwarf galaxies awaiting velocity measurements, and many of those are dSphs, for which we find a more significant flattening than for the whole population.By studying the possible degrees of corotation in the reported flattened structure in light of the missing velocity measurements, we can make a prediction on the expected fraction of corotating satellites such that the observed structure is consistent with its simulated analogs.The fractions of simulated analogs which match or exceed both the observed flatness and kinematic correlation range from 2.1% to 29.4%.This indicates that the consistency with expectations from state-of-the-art cosmological simulations is still uncertain until the missing 15 velocity measurements are obtained.

Fig. 1 :
Fig.1: The survey footprint fromChiboucas et al. (2009Chiboucas et al. ( , 2013)).The two main galaxies M 81 and M 82 are indicated as stars, the dwarf spheroidals as dots, the dwarf irregulars as squares.The two brightest satellites (NGC 2976 and NGC 3077) are marked.Colors indicate the observed velocities, gray stands for no velocity measurement currently available.
rameter ranges are closed with θ ∈ [−π/2, π/2] and ρ ∈ [−d, d],where d is the diagonal of the image.The votes by each data point are generated by calculating the values for ρ corresponding to a range of θ values via equation 1.This process amounts to a dotted curve in the Hough space for every data point.Each of these dots represents one vote this is, one line that passes through one of the data points.The region in the Hough space where the highest number of such lines cross represents the parameters with the highest number of votes at the end of this process, optimally describing the data set at hand.In an ideal scenario, where all data points are arranged on a perfectly straight line, there is a Fig. 5: A position-velocity diagram of the satellite system of the M 81 group.Blue denotes a satellite moving towards us, red away from us, relative to the mean group velocity (v=29 km/s).The gray symbols indicate the structure members without available line-of-sight velocities.The x-axis denotes the distance from the minor axis along the major axis.A negative sign on the distance is assigned based on their position on the left (positive x) of the minor axis.The circles highlight the dwarf spheroidals, the squares are dwarfs of all other morphologies.The two stars indicate the two main galaxies M 81 and M 82.

Fig. 6 :
Fig.6: Comparison of the M81 flattened structure parameters with M81-analogs from TNG50.On the x-axis, we plot the short over long axis ratio (b/a), and on the y-axis the number of corotating satellites.The green crosses show the parameters of the M81analogs in simulations while the pink star shows the parameter pair of the data.The top histogram shows the distribution of the axis ratios of the structures around simulated analogs, and the histogram on the right shows the corresponding number of kinematically correlated dwarfs.The pink-shaded regions illustrate the analog systems that match or exceed the parameters found in the data.The inset texts indicate the corresponding frequencies of such systems in simulations.Left: results for all simulated M81 analogs.Right: results for the sub-sample of analogs that also feature an M82 analog.

Fig. 7 :
Fig.7: Same as Figure6(right) but illustrating the possible outcomes given missing velocity measurements.The colored stars show the possible axis ratios and fractions of corotation in the data.The flatness of the structure will not change with new velocity information because distance measurements are available for all dwarfs.The fractions of corotation can vary from 0.5 (15/30) to 0.76 (23/30).For each possible scenario the fraction of simulated systems which match or exceed the degree of kinematic correlation is shown on the histogram to the right.The fractions systems show both flatness and degree of corotation as or more extreme as the data are shown in the center.