How do the dynamics of the Milky Way -- Large Magellanic Cloud system affect gamma-ray constraints on particle dark matter?

Previous studies on astrophysical dark matter (DM) constraints have all assumed that the Milky Way's (MW) DM halo can be modelled in isolation. However, recent work suggests that the MW's largest dwarf satellite, the Large Magellanic Cloud (LMC), has a mass of 10-20$\%$ that of the MW and is currently merging with our Galaxy. As a result, the DM haloes of the MW and LMC are expected to be strongly deformed. We here address and quantify the impact of the dynamical response caused by the passage of the LMC through the MW on the prospects for indirect DM searches. Utilising a set of state-of-the-art numerical simulations of the evolution of the MW-LMC system, we derive the DM distribution in both galaxies at the present time based on the Basis Function Expansion formalism. Consequently, we build $J$-factor all-sky maps of the MW-LMC system in order to study the impact of the LMC passage on gamma-ray indirect searches for thermally produced DM annihilating in the outer MW halo as well as within the LMC halo standalone. We conduct a detailed analysis of 12 years of Fermi-LAT data that incorporates various large-scale gamma-ray emission components and we quantify the systematic uncertainty associated with the imperfect knowledge of the astrophysical gamma-ray sources. We find that the dynamical response caused by the LMC passage can alter the constraints on the velocity-averaged annihilation cross section for weak scale particle DM at a level comparable to the existing observational uncertainty of the MW halo's density profile and total mass.


INTRODUCTION
The Large Magellanic Cloud (LMC) is believed to be on its first approach to the Milky Way (MW) since the early Universe (Besla et al. 2007). This is supported by many lines of evidence which show that it still hosts a massive dark matter (DM) halo, in line with expectations from abundance matching, ∼ 2×10 11 M (e.g. Behroozi et al. 2013;Moster et al. 2013). First, the nearby presence of the Small Magellanic Cloud (SMC) requires an LMC mass of ∼ 10 11 M in order to remain bound to the LMC (Kallivayalil et al. 2013). Similar analyses show that the recently discovered Magellanic satellites also require similarly high LMC masses, ∼ (1 − 2) × 10 11 M Patel et al. 2020) in order to have originally been bound to the LMC. The LMC is massive enough that it induces a strong reflex motion in the MW (Gómez et al. 2015), as evidenced by its effect on the timing argument with Andromeda and the nearby Hubble flow, which require a mass of ∼ 2.5 × 10 11 M (Peñarrubia et al. 2016). This reflex motion is also seen in the MW's stellar halo, both in its kinematics   (Belokurov et al. 2019;Conroy et al. 2021), all requiring a mass >∼ 10 11 M . Finally, its effect has been seen and characterized in many stellar streams around the MW, giving masses of (1.3 − 1.8) × 10 11 M Koposov et al. 2019;Shipp et al. 2019Shipp et al. , 2021Vasiliev et al. 2021).
Such a massive LMC halo is ∼ (10 − 20)% of the MW's mass (e.g. Wang et al. 2020), suggesting that this merger should have a substantial effect on the dark matter haloes of both galaxies. This has been studied with simulations (e.g. Laporte et al. 2018;Garavito-Camargo et al. 2019Petersen & Peñarrubia 2020) which have shown that substantial DM deformations are expected.
Several works have also explored the observational consequences of these deformations. Vasiliev et al. (2021) and Lilleengen et al. (2022) showed that these deformations can affect the Sagittarius and Orphan-Chenab stream, respectively. These effects are larger than current observational uncertainties, suggesting that in the future it will be possible to measure these deformations. Conroy et al. (2021) claimed the detection of the MW-LMC haloes' deformations through a positive correlation between giant halo stars in Gaia and WISE data and the simulations from Garavito-Camargo et al. (2019). Such deformations can also alter the velocity distribution of DM. Besla et al. (2019) and Donaldson et al. (2022) studied the impact of the MW-LMC dynamics on the local DM velocity distribution, finding an enhanced reach of direct DM detection experiments.
In this work, for the first time, we will explore another avenue for characterising these deforming haloes more directly: gamma-ray searches for DM signals. In the standard paradigm, 80% of the matter content of the Universe, i.e. the DM, can be explained by new particles beyond the ones in the Standard Model of particle physics. In particular, weak scale DM particle candidates may be thermally produced in the early universe through interactions with Standard Model particles. These so-called Weakly Interacting Massive Particles (WIMPs) are still present today, and can self-annihilate or decay into final stable products contributing to the fluxes of cosmic gamma rays and charged cosmic rays (see, e.g., the review on this topic in Cirelli (2012)). The spatial distribution of gamma rays from DM annihilation or decay in the sky depends on the distribution of DM in the target of interest, through the integral along the line of sight (l.o.s) of the DM density (squared in the case of annihilation). Therefore, any change in the expected DM spatial density, , affects the large-scale morphology of the DM signal. In this respect, a perfectly legitimate question to ask is: How does the MW-LMC dynamics affect gamma-ray searches for DM? We quantify the answer to this question by analysing all-sky data of the Fermi Large Area Telescope (Fermi-LAT) and search for DM at high Galactic latitudes. This approach follows traditionally performed searches for DM with Fermi-LAT data (Ackermann et al. 2012;Zechlin et al. 2018;Chang et al. 2018), and builds on the modelling and optimisation of the astrophysical background and foreground components of the gamma-ray high-energy sky.
In addition, we also study the impact of our DM spatial model on the constraints derived from gamma-ray observations of a region centred around the LMC, to be compared with previous Fermi-LAT analyses (Buckley et al. 2015), and the recent spatial model from Regis et al. (2021).
The main novelties of the present work are: (i) the modelling of the DM gamma-ray all-sky signal based on state-of-the-art simulations of the MW-LMC interaction (Donaldson et al. 2022;Lilleengen et al. 2022); and (ii) the quantification of the uncertainty on the DM limits issued by one of the most accurate models of the MW potential and its associated uncertainties (McMillan 2017), based on which we can properly assess the impact of the MW-LMC dynamics in gamma-ray DM searches. As a result, our limits on DM at high latitude represent the most up-to-date and robust constraints from Fermi-LAT gammaray observations. The paper is organised as follows: In Sec. 2, we describe the modelling of the DM signal, namely its spatial distribution, based on the outcome of simulations of the MW-LMC interaction, and spectrum. We dedicate Sec. 3 to explaining the details of the Fermi-LAT analysis, statistical framework and fitting procedure for the high-latitude sky and LMC regions. This section is complemented by Appendix A about the modelling of the astrophysical fore-and back-ground gamma-ray components. In Sec. 4, we present the new constraints on the particle DM parameter space from the high-latitude and LMC regions. We discuss the impact of varying the interstellar emission model, i.e. the dominant source of background modelling systematic uncertainties in Appendix D. We draw our conclusions in Sec. 6.

Simulation of the Milky Way -Large Magellanic Cloud dynamics
In order to explore the deformations of the MW and LMC, we use a suite of -body simulations run with the code . Unlike other -body codes which evaluate forces with a hierarchical tree (e.g. Appel 1985; Barnes & Hut 1986) or a particle mesh (e.g. Klypin & Shandarin 1983;White et al. 1983), this code evolves -body particles by using a basis function expansion (BFE). In particular, uses biorthogonal basis functions for the potential and density which satisfy the Poisson equation. The angular structure of each model is described by spherical harmonics, with harmonic indices ℓ and | | ℓ (ℓ = 0 is the monopole, ℓ = 1 is the dipole, ℓ = 2 is the quadrupole, and so on), while the radial structure is described by basis functions per harmonic order derived from the Poisson equation (Weinberg 1999). At each timestep, the coefficients of each basis function are estimated by summing each potential term in the expansion over the location of the particles (see eq. 5 in Petersen et al. 2022). With these coefficients, the forces can be readily computed by differentiating the potential. This technique has several advantages for our study (see Petersen et al. 2022 for more details). First, it is computationally efficient, scaling as O ∝ where is the number of particles used, allowing for more particles at a reduced computational cost. Second, the forces are less noisy than standard gravity solvers, allowing us to study the subtleties of the MW and LMC deformations without worrying about noise. Lastly, the BFE of the density allows us to quickly determine the density when modelling gamma-ray signals.
In this work, we make use of two sets of simulations of the ongoing MW-LMC merger. First, we use the simulation of Lilleengen et al. (2022). Their MW and LMC system is based on the results of Erkal et al. (2019) who measured the MW and LMC potentials with the Orphan stellar stream. In particular, the MW is initialized as a Miyamoto-Nagai disc (Miyamoto & Nagai 1975) with a mass of 6.8×10 10 M , a scale radius of 3 kpc, and a scale height of 0.28 kpc, a Hernquist bulge (Hernquist 1990) with a mass of 5 × 10 9 M and a scale radius of 0.5 kpc, and an Navarro-Frenk-White (NFW) halo (Navarro et al. 1996) with a mass of 7.92 × 10 11 M , a scale radius of 12.8 kpc, and a concentration of 15.3. The LMC is modelled as a Hernquist DM halo with a mass of 1.25 × 10 11 M and a scale radius of 14.9 kpc. These values all come from the best-fit model of Erkal et al. (2019) assuming a spherical DM halo for the MW. As a result, we dub this simulation the 'Erkal19' model. Lilleengen et al. (2022) show that the MW and LMC DM haloes experience strong deformations, most notably in the dipole of the MW and in the quadrupole of the LMC.
Second, we use the simulation suite of Donaldson et al. (2022) which uses the same technique and considers four different MW-LMC models. These simulations are built to roughly match the rotation curve and total mass constraints of the MW (Eilers et al. 2019;Eadie & Jurić 2019) and LMC (Van Marel & Kallivayalil 2014;Peñarrubia et al. 2016;Erkal et al. 2019). The MW model consists of an exponential disc with a mass of 5 × 10 10 M , a scale radius of 3 kpc, and a sech 2 scale height of 0.6 kpc, and a dark matter halo with a profile given by ( ) = 0 ( )˜− 1 (1 +˜) − ( ) with scaled radius = / where is the scale radius, and the truncation function ( ) = 0.5 (1 − erf [( − trunc )/ trunc ]) with trunc = 430 kpc and trunc = 54 kpc. One can set this profile to be either an NFW ( = 2) or Hernquist ( = 3) profile dark matter halo. We require a similar mass enclosed at 50 kpc by tuning the respective scale radii of the MW models: ,NFWMW = 15 kpc and ,HernquistMW = 44 kpc. The total mass of the NFW halo is 1 × 10 12 M , while the total mass of the Hernquist halo is 0.94 × 10 12 M . We build two LMCs, modelled as a dark matter halo only, again following the same truncated halo profile as above, with ,NFWLMC = 33.8 kpc and ,HernquistLMC = 63 kpc. The total masses are set to be 2.5 × 10 11 M for the NFW LMC and 2.35 × 10 11 M for the Hernquist LMC.
The simulation suite is constructed as a grid of four models by mixing the MW and LMC halo profiles. That is, one simulation is an NFW MW and NFW LMC, one is an NFW MW and Hernquist LMC, one is a Hernquist MW and NFW LMC, and one is a Hernquist MW and Hernquist LMC. We label this second set of simulations by referring to the halo profiles of the MW and the LMC at the stage of simulation initialisation, i.e. either an NFW or a Hernquist (HERN) profile. When analyzing the gamma-ray signals expected from these models, we take an approach inspired by Lilleengen et al. (2022) and consider the full multipole expansions as well as the monopoles for comparison. The monopole terms describe the spherical representation of the MW and LMC haloes, but changes over time as the coefficients of the monopole radial basis functions vary. Due to the relatively short travel times of gamma rays through the Milky Way (∼ 1 Myr to travel 300 kpc) compared to the timescales over which the basis function changes substantially (∼ 50 − 100 Myr, see fig. 1 of Lilleengen et al. 2022), we only consider the coefficients at the present-day. By comparing these two, we can see how much the deformations affect the predicted gamma-ray signal. We consider this difference as the most robust result of this work. Indeed, as a word of caution, we notice that, while the initial models of Lilleengen et al. (2022) and Donaldson et al. (2022) were consistent with the MW gravitational potential constraints, the full consistency with the MW potential has not been a posteriori checked for the finally deformed models. That is, the rotation curve and mass enclosed constraints originally imposed may no longer be met.
Therefore, in order to quantify the impact of the absolute value of the new constraints for gamma-ray DM searches, in addition to these two sets of models, we also consider a static MW model. For this potential, we use the results of McMillan (2017) who modelled the MW with a bulge, four discs (thin, thick, H , and H 2 ), and an NFW DM halo. McMillan (2017) fit this model to a range of data and constraints: maser data, the solar velocity, terminal velocity curves, the vertical force near the plane of the disc, and a mass constraint based on satellite kinematics (see Sec. 3 of McMillan 2017 for more details). While these constraints are primarily within the plane of the MW disc, it represents one of the most accurate models of the MW potential. To explore how the uncertainties in the MW potential affect the predicted gamma-ray signal, we sample over the posterior chains from McMillan (2017). This allows us to (a) quantify the systematic uncertainties on the high-latitude DM limits from the MW gravitational potential, and (b) properly assess the impact of the variations of the limits induced by the MW-LMC dynamics.

Dark matter-induced gamma-ray signal
In this work, we consider the gamma-ray emission resulting from pair-annihilating thermally produced DM in the MW and LMC halo. We restrict ourselves to the prompt gamma-ray component of these interactions neglecting potential secondary or tertiary contributions from particle cascades triggered by the primary annihilation products. The expected (prompt) differential gamma-ray flux dΦ /d /dΩ at the top of the Earth's atmosphere reads (see, e.g., Cirelli et al. 2011;Bringmann & Weniger 2012) where refers to the direction of the line-of-sight in Galactic coordinates and quantifies the gamma rays' energy. The DM-induced gamma-ray flux factors into two contributions under the assumption of velocity-independent annihilation cross section -the so-called -wave annihilation process. The term in the first parenthesis is commonly referred to as J -factor while the term in the second parenthesis includes and describes the particle physics model chosen for the DM candidate under study. In what follows, we provide further details about the ingredients required to compute both contributions to the DM gamma-ray signal.

J -factor all-sky maps for the Milky Way
In order to produce J -factor maps of the MW that incorporate deformations induced by the dynamics of the MW-LMC system, we take the density from the present-day snapshots of the simulations in Lilleengen et al. (2022) and Donaldson et al. (2022), and measure the square of the DM density along lines of sight on a HEALPix grid (Górski et al. 2005) with side = 64, resulting in 49,152 lines of sight. Along these lines of sight, the density is evaluated in the centre of 1 kpc-sized bins out to 100 kpc. This choice of resolution is not a limitation of the BFE, which can faithfully reproduce structure on ∼ 100 pc scales, but rather motivated by the chosen resolution for the Fermi-LAT gamma-ray data set (see Sec. 3.1). For these mock observations, the Sun is placed at a distance of 8.249 kpc from the Galactic centre (Gravity Collaboration et al. 2020) in the plane of the MW disc. As discussed in Section 2.1, we create J -factor maps for the full BFE as well as solely the monopole term for the MW standalone or the combined MW-LMC system. The outlined extraction scheme allows us to directly perform the l.o.s. integral of 2 where the l.o.s. direction is given by the central coordinates of a particular HEALPix pixel. To this end, we linearly interpolate the density values at the extracted positions and numerically integrate the result from 0 to 100 kpc.
In contrast, we derive two-dimensional all-sky J -factor maps of the static MW model of McMillan (2017) (Navarro et al. 1996;Navarro et al. 1997) adopted for the MW: (i) distance of the Sun to the Galactic centre, , (ii) scale radius of the NFW profile, (iii) virial radius 200 of the MW and (iv) the DM density at the position of the Sun. To illustrate the expected deformations induced by the MW-LMC dynamics, we provide a direct comparison in terms of the Erkal19 simulation suite in Fig. 1. The panels display the outer MW halo in the southern hemisphere of the sky at −20 • , which turns out be the optimal part of the sky to perform the Fermi-LAT gammaray analysis (see Sec. 3.4). The left panel shows the J -factor map associated to the monopole term of the BFE, which has been used to decompose the DM density distribution in the MW, whereas the central panel displays the J -factor map of the MW halo derived from the full BFE of the Erkal19 simulation. Confronting the monopole contribution with the full BFE yields the most faithful assessment of the expected deformations since the initial conditions are the same for both scenarios. To better highlight the differences, we show the relative ratio of the two quantities in the right panel.

J -factor all-sky maps for the Large Magellanic Cloud
To study the impact of the MW-LMC dynamics on the prospects of indirect searches for DM in the LMC, we repeat the approach outlined in Sec. 2.2.1 for the LMC halo, which results in all-sky maps of the LMC J -factors including the full BFE or only its monopole term for the five simulation suites. Although the DM density has been extracted at a time slice that corresponds to the current stage of the evolution of the MW-LMC system, the position of the LMC differs between the simulations. Moreover, it is not necessarily aligned with the astronomically determined centre of this galaxy, which itself is debated in the literature (Van Marel et al. 2002;Kim et al. 1998;Van Marel & Kallivayalil 2014). For definiteness, we define the LMC centre to be located at (ℓ, ) = (280.54 • , −32.51 • ) which is the favoured rotational centre derived from stellar kinematics (Van Marel et al. 2002). All HEALPix J -factor maps are rotated such that the pixel with the largest value coincides with this position thus aligning the DM halo with the stellar centre of the LMC.
To enable a comparison with the expectations from a static LMC DM halo profile, we adopt the parametrization of the NFW and Hernquist profile from Tab. 1 in Regis et al. (2021). The authors of the latter study have examined the radio emission from the centre of the LMC, and obtained constraints on thermally produced WIMP DM under the assumption of different LMC halo shapes whose parameters have been determined via a fit to rotation curve data (Kim et al. 1998), i.e. the inner parts of the LMC. We derive all-sky Jfactor maps from these two static profiles (NFW and Hernquist) with CLUMPY 1 . The thereby generated all-sky J -factor maps are by default centred on the Galactic centre. We thus rotate the CLUMPY output on the stellar centre of the LMC as before.
In full analogy to Fig. 1, we visualise the degree of deformations of the LMC DM halo profiles via the Erkal19 simulation in Fig. 2. The left panel displays the LMC J -factor map taking only into account the monopole term of the BFE whereas the central panel shows a dynamically perturbed LMC halo according to the full BFE of the Erkal19 simulation. The right panel complements both profiles by highlighting the relative differences between the selected cases. Especially the latter panel illustrates the forward wake in the northern hemisphere induced by the LMC passage.
In comparison to spherically symmetric, static LMC density profiles fitted to rotation curve data -as done, e.g., in Regis et al. (2021) -we find that the central part of the LMC is predicted less dense. We stress that the reduced central J -factor is related to the initial conditions of the simulation: the LMC models in both Lilleengen et al. (2022) and Donaldson et al. (2022) are dark matter only and do not include a stellar component for the LMC, which will affect the central dark matter distribution. The halo parameters in Regis et al. (2021) have been fixed via fits to stellar rotation curves and thus incorporate information about the small-scale behaviour around the centre of the LMC while the initial LMC haloes in both Lilleengen et al. (2022) and Donaldson et al. (2022) aim to match the enclosed LMC mass at larger radii, 9 kpc. Since the total mass is rather insensitive to the innermost part of a DM halo, differences between simulated and static LMC profiles may occur. 1 Since the virial radius of the LMC is larger than the distance of the Solar System to this object, we use the 'galactic' mode of CLUMPY by defining the LMC centre as the new reference Galactic centre.

Annihilation spectra and gamma-ray flux
In this work, we are referring to a generic thermally produced particle DM candidate at the weak scale, such as those belonging to the class of WIMPs. In addition, we assume that the WIMP candidate is a Majorana fermion ( = 1) annihilating into a single Standard Model particle species , i.e. the branching ratio for this process reads ≡ 1. We consider a single annihilation channel, namely →¯to present our results. The associated differential gamma-ray spectrum d /d , stating the expected average number of photons with energy per annihilation event, is taken from PPPC (Cirelli et al. 2011). Eventually, the DM gamma-ray flux depends on two parameters, the DM mass as well as the velocityweighted, thermally averaged annihilation cross section ann that determines the intensity of the signal, which is the parameter we ultimately want to constrain.
We stress that the main purpose of our study is to analyse the impact of deformations of the MW halo due to the passage of the LMC on indirect gamma-ray DM searches. These dynamically induced deformations are -at least for -wave annihilation processes -merely affecting the spatial morphology of the DM signal, and their impact is the same no matter what the chosen spectrum is. Thus, the choice of the exact particle DM model plays a minor role on our results and discussion. Consequently, we focus only on one annihilation channel in the main text and provide the results for the + − -channel in Appendix E.
3 Fermi-LAT DATA ANALYSIS

Data selection
The analysis is based on ∼12 years of Fermi-LAT data (Pass8 reconstruction standard) taken between the 4th August 2008 and the 3rd September 2020. The considered reconstructed gamma-ray energies range from 500 MeV to 500 GeV while we focus on those photons classified as belonging to the ULTRACLEANVETO event class and FRONT+BACK event type. We apply further cuts on the selected sample of photons via the requirement of DATA_QUAL==1 && LAT_CONFIG==1 and restricting the event zenith angle to < 90 • , which reduces the contamination of this sample by Earth limb photons. All work that requires selection, manipulation or simulation of Fermi-LAT data has been conducted via use of the Fermi Science Tools 2 (version 2.0.8).
We perform a binned log-likelihood analysis for which we bin the selected data as all-sky maps in the HEALPix format (Górski et al. 2005) with side = 64 -resulting in a mean pixel spacing of ≈ 0.9 • -as well as 30 logarithmically spaced energy bins. Due to the scarcity of photon events at the highest energies of the LAT's sensitivity range, we rebin the LAT data above 7 GeV as well as all astrophysical background and signal all-sky maps (see Sec. 3.3) into larger macro bins. Hence, the high-energy range is included within the analysis framework by creating the following two macro bins:

Statistical framework
We employ a template-based analysis to constrain DM annihilation in the MW halo, a well-known procedure for Fermi-LAT gamma-ray analyses. To this end, we seek to describe the processed LAT data map via a set of astrophysical background templates { } ∈ , which is supposed to capture as best as possible the expected different types of gamma-ray sources in and outside of the MW. To these background components we add the DM signal and evaluate whether it is preferred by the data and if such a preference is not statistically significant (see Eq. C1 in Appendix C), we set upper limits on the strength of the DM contribution. The foundation of the statistical framework to conduct this kind of analysis is the binned Poisson likelihood function where denotes the spatial pixels of the all-sky map and the energy bins of the templates. The linear combination of our set of background and signal templates µ is called the model whereas n represents the data map to which the model is fitted by maximising the value of the likelihood function. In detail, this linear combination of components is defined as follows: where, again, the index denotes the energy bins used in this analysis. Such a model definition relies on two kinds of normalisation parameters: • Background normalisation parameters • A single signal normalisation parameter, DM , which enables us to exploit both the spectral and spatial shape of the signal component (more details concerning the morphology of the signal are given in Sec. 2.2.1).
We modify the standard Poisson likelihood function in Eq. 2 in two directions: (i), we work with the logarithm of this function, thus turning the statistical inference into a series of function minimizations which are numerically more accessible with well-tested algorithms and software packages; (ii), we employ a pixel weighing scheme to incorporate the impact of instrumental systematic uncertainties. To this end, we adopt the weighted Poisson log-likelihood prescription developed by the Fermi-LAT collaboration in Abdollahi et al. (2020), which reads Per energy bin, each template pixel is assigned a weight whose value we obtain via the Fermi Science Tools calling its routines gteffbkg, gtalphabkg and gtwtsmap. These weights are essentially obtained via integration in space and energy of the provided model or Fermi-LAT data 4 . Since we aim to incorporate the effect of instrumental systematic uncertainties, we rely on a "data-driven" weight calculation approach, that is, we utilise the LAT data themselves to derive the weights. Thus, pixels are penalised in particularly bright parts of the gamma-ray sky (point-like/extended sources, diffuse emission along the MW's disc) taking into account potential sources of systematic errors like contamination of the data sample by charged cosmic-ray events, calibration of the instrument's point spread function (PSF) or spectral mismodelling of large-scale diffuse sources. Throughout the analysis, we keep the level of systematic uncertainties to 3% (for all energy bins) following the approach of the Fermi-LAT collaboration Abdollahi et al. (2020).
Hypothesis testing -the procedure to discriminate between competing, alternative descriptions of reality, here represented by our model in Eq. 3 -is implemented via the log-likelihood ratio test statistic, which in the relevant case of setting upper limits reads This construction relies on the profiled likelihood function (that is, treating the background normalisations , ∈ as nuisance parameters) as discussed in Cowan et al. (2011). Model parameters marked with· denote the best-fit values found via minimization of the log-likelihood function. The test statistic in Eq. (5) only depends on the DM normalisation. Moreover, possible values of DM smaller than the best-fit value are discarded, thus, the TS-distribution follows a half-2 -distribution with one degree of freedom (see Sec. 3.6 of Cowan et al. (2011)). Therefore, we set an 95% confidence level (C.L.) upper limit on DM where the test statistic attains a value of 2.71.

Astrophysical background model selection
The gamma-ray sky as seen by the LAT can be decomposed into a combination of a multitude of distinct contributions, which either originate in Galactic or extragalactic sources. To constrain the extended signal such as gamma rays from DM annihilation in the outer MW halo, we incorporate the same astrophysical contributions considered in Calore et al. (2022b). We summarise the components below and refer the interested reader to Appendix A for more details.
• The interstellar emission (IE) originates from cosmic-ray interactions with gas and low-energy ambient photon fields. Among different IE models and based on the findings of Calore et al. (2022b), we consider as a baseline IE model the henceforth called Lorimer I, taken from the set of realisations considered in the "1st Fermi LAT Supernova Remnant Catalog" 5 (Acero et al. 2016). Three alternative 4 The technical aspects of these routines' implementation is provided at https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/ weighted_like.pdf or Appendix B of Abdollahi et al. (2020). 5 The model files have been made public by the Fermi-LAT collaboration at: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/1st_ SNR_catalog/. models utilised in Calore et al. (2022b) -called foreground model A, B and C -are adopted from Ackermann et al. (2015). More details about the IE models are given in Appendix A.
• The isotropic diffuse gamma-ray background (IGRB) is a largescale contribution to the gamma-ray sky which is spatially isotropic and believed to originate from the superposition of many, subthreshold, sources (Fornasa & Sánchez-Conde 2015).
• The modelling of the resolved point-like and extended sources is based on a 10-year data set, i.e. the so-called 4FGL-DR2 (Abdollahi et al. 2020;Ballet et al. 2020).
• We also model other large-scale extended gamma-ray emissions from the Fermi Bubbles (FB), Loop I, the Sun and the Moon, following standard practice in LAT data analyses.
Passing from these background (and signal) models to templates containing the expected photon counts per sky direction is achieved by dedicated routines of the Fermi Science Tools. For our statistical tests, we generate "infinite statistics" or Asimov realisations (Cowan et al. 2011) of the background and signal models via gtmodel, which internally processes the full convolution of the input model files with the LAT's instrument response functions associated with the selected gamma-ray data set (c.f. 3.1). We include the LAT's energy dispersion for all background and signal components by adding the parser argument edisp_bins=-1.

Fitting procedure and region-of-interest optimization
We employ the following general analysis rundown and reasoning to statistically soundly and robustly assess the implications on DM pair-annihilation in the MW halo from Fermi-LAT data: (i) Generating a baseline gamma-ray sky model from a fit of the full set of astrophysical background templates to the all-sky data.
(ii) Including the signal component, preparing the MW halo analysis by shrinking the total region of interest (ROI) to a smaller fraction of the sky that yields a good agreement between the TS-distribution (Eq. 5) with respect to LAT data and baseline model as input data n.
(iii) Setting upper limits on the DM pair-annihilation cross section with respect to the optimised ROI and a particular scenario of signal templates.

Deriving a baseline model of the gamma-ray sky
The importance of a baseline model entirely composed of the background templates considered in this analysis lies in its utility for all future statistical inference. Such a model may be used as an alternative data map that is guaranteed to contain only known gamma-ray emitters, which is not necessarily true for the experimentally observed gamma-ray data. Hence, whenever the TS-distribution as a function of the signal normalisation (Eq. 5) shows a comparable behaviour with respect to baseline model and real LAT data, the selected sky region can reliably be described via the set of background and signal templates at hand. Since the baseline model is a combination of "infinite statistics" templates, we can draw Poisson realisations to quantify the expected statistical scatter of the eventually reported upper limits.
Since we employ the same astrophysical gamma-ray emission components as in Calore et al. (2022b), we repeat the prescription presented therein to derive the baseline model. In short, the algorithm consists of an iterative fitting scheme that splits the entire sky into three disjoint parts defined as follows: (a) high-latitude -| | > 30 • and neglecting the "patch"-region (c.f. Abdollahi et al. 2020 The normalisation constants of a particular gamma-ray emitter are only fit in the part of the sky where it mainly contributes to while these parameters are held fixed at the thereby obtained best-fit values when the other regions are addressed. The exact details about the assignment of sky regions to particular background components are provided in Sec. IV B of Calore et al. (2022b).
We run the iterative fit for 100 times to eventually derive a baseline fit to the all-sky LAT data. Besides this general scheme, there are a few technical modifications necessary to incorporate the IE models foreground model A, B and C -henceforth abbreviated as FGMA, FGMB and FGMC -into this framework as well as the wealth of sources in 4FGL-DR2, which deviates from the original recipe in Calore et al. (2022b). Treatment of IE models: All five IE models feature a single template quantifying the gamma-ray emission from inverse Compton (IC) scattering processes, which we split into three sub-templates whose boundaries coincide with the definitions of the 'high-latitude', 'outer galaxy' and 'inner galaxy' regions of the iterative fit. The same procedure is applied to the gamma-ray emission associated with the gas maps in FGMA-C. After the iterative fit, all IE-related components are multiplied by their best-fit normalisation parameters and added to form an optimised IE template that is utilised in the subsequent analysis parts. Treatment of 4FGL-DR2 sources: Combining all sources listed in 4FGL-DR2 into a single template causes the brightest sources in the template to drive the best-fit value of the template's normalisation parameter. We weaken the impact of bright sources by separating sources with an energy flux of 100 < 4·10 −10 MeV cm −2 s −1 (integrated from 100 MeV to 100 GeV) from those above this threshold. A source that surpasses this threshold is fitted individually during the all-sky fit (d) of an iteration step. We list these bright sources and some of their properties in Tab. B1 in App. B. All others are combined in a single 4FGL template. As for the IE templates, all 4FGL-DR2 templates are combined after the fit according to their best-fit values, hence creating an optimised total 4FGL template to be used in all steps that follow.

Optimising the analysis region of interest
The strategy to perform an ROI optimization before conducting any statistical inference is heavily based on the approach presented in Zechlin et al. (2018). In order to search for an ROI that yields statistically reliable upper limits on the DM annihilation cross section, we resort to the southern hemisphere -thereby circumventing most of LoopI's contribution to the gamma-ray sky -and investigate the distribution of the test statistic in Eq. 5 for both the true LAT data and the baseline model as data input vector n.
As the signal spectrum depends on the chosen DM parameters and annihilation channel, we construct a model-independent DM signal -that still exhibits the currently assumed spatial morphology of the MW halo -by exchanging the generic WIMP spectra with a power law of spectral index -2. Hence, the signal's spectrum is featureless and yields non-zero photon counts throughout the entire energy range considered in this work. The initial normalisation of the power law is somewhat arbitrary since there is no connection to a physical DM model. We hence choose the normalisation such that the expected counts in the first energy bin are maximally of order unity per pixel.
The optimisation is carried out by systematically shrinking the In addition, we mask the FBs by setting all pixels to zero whose counterpart in the FB template predicts non-zero counts. We introduce the requirement of 20 • to reduce the impact of IE along the Galactic disc. For each ROI, we scan the TS-distribution with respect to the LAT data and baseline model for TS ∈ [0, 25] as a function of DM and use the ℓ 2 -metric to quantify their mutual compatibility. We select those Galactic longitude and latitude ROI boundaries for which this metric is minimal. We stress that this optimisation step must be repeated for each of the five considered IE models.
Treatment of the optimised 4FGL-DR2 template: A special note concerns the treatment of 4FGL catalogue sources. In contrast to the iterative fit, we now mask a circular region around the central position of each source. The mask radius is energy-dependent and corresponds to the 95% containment radius of the LAT's PSF 6 for the chosen LAT event class and type. For the first three energy bins, however, we reduced the mask radius to 90% of the PSF size. The latter exception has been introduced to ensure a reasonable number of pixels to be non-zero so that a template-based analysis remains feasible.
We illustrate the TS-distribution comparison in Fig. 4 of Sec. 4.1 for a particular combination of spatial DM distribution and IE model. In what follows, we will always report the selected optimal ROI. Moreover, we test our analysis pipeline in terms of its capability to recover an injected signal in simulated data. The results of this sanity check are described and discussed in Appendix C.

Dedicated analysis of the Large Magellanic Cloud region
Since the LMC passage through the MW halo does not only induce a response in the latter but also in the LMC DM halo itself, we aim to analyse the surroundings of the LMC in a dedicated gamma-ray study. To this end, we utilise the same LAT data selected and described in Sec. 3.1 but restrict the ROI to a maximal size of 30 • × 30 • centred on the stellar position of the LMC at (ℓ, ) = (280.54 • , −32.51 • ) in agreement with the centre of the LMC J -factor maps discussed in Sec. 2.2.2. The pixel size of the binned data is set to 0.1 • × 0.1 • while the energy binning is adopted from the all-sky data set.
The general analysis and fitting strategy is completely analogous to the scheme outlined in the previous section, with the addition of emission model components of the LMC itself. A decisive characteristic of the dedicated LMC analysis is the need for an additional astrophysical background component quantifying the expected gamma-ray emission due to cosmic-ray interactions with gas and radiation fields in the LMC. To this end, we include four separate templates adopted from the set of extended templates being part of the 4FGL-DR2 catalogue. The components LMC-Galaxy, LMC-North, LMC-FarWest and LMC-30DorWest, representing the cosmic-ray induced gammaray emission of the LMC, have been derived in a previous study of the Fermi-LAT collaboration (Ackermann et al. 2016) based on a six-year data set. These models are the result of a convolution of gas column density maps of the LMC reported in Lande et al. (2012) and a data-driven approach to extract regions of significant extended gamma-ray emission associated with the LMC. Due to their data-driven nature, these models may already incorporate a contribution from an exotic gamma-ray emitter like DM if taken at face value. The authors of Buckley et al. (2015) comment on this point by analysing the correlation between the DM component and the LMC-related templates. They find a particularly sizeable correlation with the "LMC-Galaxy" template, which is centred on the stellar position of the LMC and, at the same time, the most extended among the considered ones. Such a correlation might artificially boost the constraining power of our template-based approach. We mitigate this effect by including the four LMC templates as unmasked, independent model components whose normalisation parameters are a subset of all the model's nuisance parameters that are profiled over to set upper limits on the DM annihilation cross section. Moreover, our ROI is considerably larger than the 10 • × 10 • ROI adopted in Ackermann et al. (2016) and Buckley et al. (2015). Consequently, the spatial extension of the DM component beyond the size of the four LMC templates partially brakes the degeneracy with these templates. We anticipate that the later derived upper limits that we find are indeed comparable to the results reported in Buckley et al. (2015).
The remaining differences to the scheme in Sec. 3.4 occur in the derivation of a baseline fit for a particular IE model and the optimisation of the analysis ROI. In details, these changes are: (a) All additional point-like gamma-ray sources within 3 • of the stellar position of the LMC are fitted separately by means of a template for each individual source. (b) All detected gamma-ray sources, which already were reported in the 3FGL catalogue (Acero et al. 2015) and which are at a distance of 3 • < 40 • from the stellar centre of the LMC are cast into a single 3FGL template.
( (e) The FBs as well as the gamma-ray emission from the Sun and Moon are neglected because their contribution is almost zero in the chosen ROI.
(f) All astrophysical background components are fitted in the full ROI at the same time without specifying disjoint fit regions to derive a baseline fit to the gamma-ray data in the ROI.
(g) After the baseline fit, the IE components are summed with their best-fit normalisation to an optimised IE model. The localised pointlike sources in 4FGL-DR2 are treated analogously except for the LMC-related templates, which we keep as individual templates even in the stage for setting upper limits.
(h) The ROI optimisation translates to symmetrically shrinking the width and height of the ROI up to the point where the optimal correspondence between expected TS-distribution and observed TSdistribution (with respect to the auxiliary DM signal template featuring a power law spectrum) is achieved.
(i) To explore the TS-distributions and to set upper limits, we mask all detected point-like sources in the same way as outlined before except for the positions of the LMC-related emitters.

CONSTRAINTS ON PARTICLE DARK MATTER
In this section, we present the results of the constraints on DM pair-annihilation processes in the MW and LMC haloes with the analysis framework described in Sec. 3. As mentioned in Sec. 2, the characterisation of the gravitational potential of the MW comes with a non-negligible uncertainty. The dynamical impact of the LMC passage may be regarded as a second-order effect that adds to the inherent uncertainties of the available astronomical observations. Consequently, we first explore the expected scatter of constraints on the DM parameter space induced by observational uncertainties of the MW's DM halo while in the following subsections, we shed light on the significance of including the MW's response to the LMC passage for the outcome of indirect searches for DM.
The following results for the MW outer halo have been obtained with our benchmark IE model, Lorimer I, unless stated otherwise, whereas the benchmark for the LMC analysis is FGMA. We investigate the impact of the chosen IE model on the DM constraints in Appendix D -we anticipate that limits at high latitude are mildly affected by the IE choice, and are robust in this respect.

Impact of the uncertainty of the Milky Way's gravitational potential
We utilise the assessment of the MW's gravitational potential and mass distribution in McMillan (2017) to explore the impact of their uncertainty on DM indirect searches in the outer MW DM halo. To re-iterate, the author of that study assumes a standard NFW profile (inner slope parameter = 1) to describe the MW DM halo and following this fundamental assumption derives, among others, posterior distributions for , , , and 200 to scale and normalise the NFW profile in accordance with the observational constraints. From these posterior distributions we randomly draw 200 points and generate the associated all-sky J -factor map of the MW. For each of these MW realisations, we derive Fermi-LAT upper limits on the DM pair-annihilation cross section following the scheme outlined in Sec. 3.4.
As the first step of this analysis pipeline, we search for an optimal ROI in the gamma-ray sky by successively shrinking the considered fraction of the southern hemisphere in Galactic longitude and latitude. We illustrate the statistical performance of the thus obtained optimised ROI via one particular realisation from the McMillan posterior distributions characterised by the tuple of parameters ( , , , 200 ) = (12.1 kpc, 8.3 kpc, 0.4 GeV cm −3 , 199.0 kpc) in Fig. 4. We find the best-suited analysis ROI to be defined by ℓ ∈ [−167 • , 167 • ] and ∈ [−90 • , −35 • ], which ensures a reasonable compromise between the constraining power of the remaining gamma-ray sky and the statistical robustness of the resulting upper limits. We show the Fermi-LAT data integrated over the analysis' energy range inside this optimal ROI in Fig. 3.
Both panels in Fig. 4 have been obtained from 200 Poisson realisations of the baseline fit. On one hand, the left panel of this figure shows that the TS-distribution on LAT data stays within the 68% containment band of the scatter of the baseline fit's TS-distribution. On the other hand, the width of its median's parabola deviates from its analogue on LAT data. This qualitatively indicates that the chosen ROI may contain further gamma-ray emission components which are not or not fully accounted for in the selected set of astrophysical background contributors. Hence, this ROI may feature further emission components that are degenerate with the DM signal. Consequently, the 95% C.L. upper limits fall within the 68% containment band derived from mock data for most of the scanned DM masses except for light DM below 10 GeV. Here, the constraints are slightly stronger (at the 2 level) than expected from the baseline fit. We have made sure that the chosen ROI is suitable for all of the 200 random parameter tuple ( , , , 200 ). In fact, the observed moderate fluctuation for light DM is a common feature among all of these realisations of the MW halo. From a qualitative point of view, the consistency between the TS-distribution of mock and real data in this work is similar to the benchmark scenario studied and described in Fig. 1 of Chang et al. (2018). There, the authors find the same slight deviation for light DM but they also show a more pronounced deviation for DM around the TeV scale where our TS-distributions differ only at the 1 level. However, the selected optimal ROIs in both analyses are largely disjoint since we exclude the position of the FBs. The optimised ROI in Zechlin et al. (2018) (c.f. Fig. 5 therein) shows a larger overlap with our analysis ROI. The reported accordance of the statistical expectations from their baseline fit and the corresponding performance on real data in their Fig. 4 is well in line with the results presented in this work, i.e. consistency at the 1 level for most of the tested parameter space. The comparison with both literature results corroborates that our optimal ROI provides statistically sound upper limits from DM annihilation processes in the outer MW halo.
The left panel of Fig. 5 summarises the set of constraints derived from these 200 random realisations of the MW DM halo assuming WIMP DM pair-annihilating into¯final states that eventually generate a prompt gamma-ray flux via further processes. The median 95% C.L. upper limits obtained within the optimised ROI and with respect to LAT data are denoted by a black dashed line whereas the corresponding 1 /2 containment bands are depicted as dark grey/light grey shaded bands. The impact of the observational uncertainty of the MW's gravitational potential is less than a factor of two at the 2 level across the entire range of DM masses considered in our analysis. This finding may seem astonishing at first glance. It is well known that the uncertainty of the shape of the MW's DM halo has a much larger effect on indirect searches towards the Galactic centre since the DM distribution in this region of the Galaxy can be peaked, core-like or even be completely devoid of DM (see, for example, Iocco et al. 2011Iocco et al. , 2015Iocco & Benito 2017;Karukes et al. 2019;Benito et al. 2021). However, high Galactic latitudes as inspected by us are dominated by the outer MW DM halo, that is, we probe a much larger volume of the total halo with guaranteed DM presence in order to stabilise the Galactic rotation curve of the MW. Hence, small-scale uncertainties of the MW's gravitational potential, for instance in the Galactic centre, are washed out by the fact that we investigate a large volume of the MW DM halo and probe its cumulative gravitational imprint.
Results for the + − DM annihilation channel are provided in Appendix E.

Impact of the perturbation of the Milky Way's dark matter halo caused by the Large Magellanic Cloud's passage
As discussed in Sec. 2, the passage of the LMC through the MW halo induces dynamical responses, which add to the already discussed uncertainty of the MW's gravitational potential. The respective responses in the form of wakes are particularly present in the outer MW halo trailing the LMC orbit or developing in front of its current orbital direction. Hence, this dynamical effect is supposed to be detectable in the southern hemisphere, as visualised in Fig. 1. We quantitatively examine the importance of the LMC passage via the simulations of the MW-LMC system described and discussed in Sec. 2.1. To this end, we confront in Fig. 6 the 95% C.L. upper limits on the DM pair-annihilation cross section ( −→¯) obtained from three different simulations with the previously derived uncertainty of the same limits due to observational uncertainty of the MW's gravitational potential. We distinguish constraints for the cases of taking into account solely the MW DM halo's monopole term from the BFE (black lines) and the corresponding full BFE (red lines) according to the respective simulation. As concerns the chosen simulations, we have selected two models from the Donaldson et al. (2022) suite -assuming an initial NFW (solid)/HERN (dotted) profile for both the MW and the LMC -and the Erkal19 simulation (dashed). We have checked that the optimised ROI discussed in Sec. 4.1 can also be applied to these signal DM profiles. We emphasize that the utilised J -factor maps are the sum of both the MW and LMC halo, although the additional boost due to the LMC halo is marginal.
A comparison of the black curves reveals that the different initial conditions in terms of MW halo profile are well within the 2 range of the uncertainty reported by McMillan (2017) and thus plausible representations of the physically realised halo of the MW. If we compare the results in red that incorporate the full influence of the LMC as a perturber of the MW halo with those that disregard this effect, we see that the impact on the DM constraints depends on the respective simulation. While both simulations from Donaldson et al. (2022) predict an almost negligible improvement of the upper limits, the Erkal19 simulation yields a much stronger response that induces an improvement for light DM as large as the 1 uncertainty of the MW's gravitational potential. The difference of the obtained bounds is not related to the total J -factor of the full MW halo -which is the highest for the NFW-NFW simulation -but rather correlated with the appearance of deviations from the static, spherically symmetric DM halo scenario. Local over-and under-densities induced by the dynamical response of the MW halo are most pronounced in the Erkal19 simulation, which explains the prominent effect in Fig. 6 compared to the remaining simulation suites. The difference is thus a manifestation of the initial conditions of each model simulation (i.e. assumptions about extent, profile, mass, internal structure/deformation and components) of the MW-LMC system. We discuss these properties and their interplay in Sec. 5.1.
.000 4.328 log 10 N γ  baseline fit for IE model Lorimer I in terms of the median expectation and its statistical scatter (green: 68% containment, yellow: 95% containment). The red line indicates the observed TS-distribution with respect to the select LAT data in the designated ROI. The normalisation parameter test is an auxiliary parameter without physical relevance that controls the strength of the inject DM signal following a power law with spectral index -2. (Right:) Median 95% C.L. upper limits (black) on the velocity-weighted thermally averaged DM pair-annihilation cross section depending on the assumed DM mass for the prompt gamma-ray emission from DM annihilating into¯final states and its statistical scatter (green: 68% containment, yellow: 95% containment) derived with respect to the baseline fit with Lorimer I. The red line displays the respective observed upper limits. The blue line indicates the thermal WIMP annihilation cross section for a DM particle with mass to generate the measured cosmological abundance of DM in the early universe (computed with DarkSUSY (Bringmann et al. 2018) and current Planck data (Aghanim et al. 2020)).

Indirect searches towards the Large Magellanic Cloud
We assess the constraining power of the LMC as a target for indirect DM searches and the importance of including the altered shape (and mass) of the LMC DM halo due to its gravitational interaction with the MW halo. This analysis follows the outline given in Sec. 3.5. We stress again that we select FGMA as the benchmark IE model because it shows a better performance than Lorimer I in this particular ROI. We discuss the systematic uncertainty due to the chosen IE model for this dedicated LMC analysis in Appendix D and Fig. D1 therein.
In Fig. 7, we summarise our findings regarding the analysis of the LMC's DM halo. In the left panel, we compare the constraints on the DM pair-annihilation cross section (channel: →¯) using either the full BFE (red) or only its monopole term for the LMC halo (black) of a particular simulation of the evolution of the MW-LMC system. For each of these sets of halo profiles, we optimised the square ROI centred on the LMC position in terms of its size. The optimal ROI sizes are stated in Tab. 1. We observe a mild impact of the dynamical response of the LMC on the final upper limits for the four simulations from Ref. Donaldson et al. (2022) whereas -and as we had already noticed in the case of the outer MW halo -the Erkal19 simulation suggests a more pronounced effect, which may be as large as a factor of two. We discuss the observed simulation-dependent variations in more detail in Sec. 5.2.
The right panel of Fig. 7 puts the results from the perturbed LMC haloes in the context of the current state-of-the-art in the field. We   confront a subset of the DM upper limits from the left panel (red curves) to constraints derived with standard DM haloes following the profiles from an NFW or Hernquist profile with parameters adopted from Regis et al. (2021) presented in Sec. 2.2.2 (orange lines). The latter set of upper limits hence represent a static LMC halo that does not incorporate deviations from spherical symmetry. While both types of LMC haloes agree on the strength of the constraints for light DM particles, they differ for larger DM masses. In fact, haloes from Regis et al. (2021) feature a more pronounced cusp towards the LMC's centre whereas our simulated haloes appear less peaked and smoother in general. The discussion is continued in Sec. 5.2. Representative results of the Fermi-LAT collaboration's search for DM in the LMC (Buckley et al. 2015) with five years of data are displayed with different shades of blue. We comment on the comparison to our bounds in Sec. 6.

Systematic uncertainties affecting the study of the outer Milky Way halo
In this work, we have demonstrated that accounting for the deformed dark matter haloes of the Milky Way and LMC is crucial for getting accurate cross section constraints from gamma ray searches. Indeed, Figure 6 shows that the change to the constraint from including the deformations (i.e. the difference between the black and red curves) is comparable to the uncertainty on the constraint due to uncertainties in the Milky Way's mass profile (i.e. the grey band). In addition to the effect of accounting for deformations, we also see that precise constraints also depend on how we model the Milky Way and LMC system. The reason for this is that these simulations span a wide range of Milky Way masses and concentrations which affect the strength of the dark matter deformations. In particular, the simulations from Donaldson et al. (2022) have more massive and concentrated Milky Way haloes which deform less than the simulated Milky Way halo in Lilleengen et al. (2022). For reference, we note that the model in Lilleengen et al. (2022) appears similar to models in the literature of the Milky Way-LMC interaction Rozier et al. 2022). Given this range of possibilities, we argue that the deformations can be considered as a source of systematic uncertainty on the inferred cross section until they are better characterized. We note that there are also other physical effects which have altered the Milky Way's dark matter halo and could affect the cross section constraints in similar ways. For example, the Gaia-Sausage/Enceladus (GSE, Belokurov et al. 2018;Helmi et al. 2018) merger likely brought a substantial amount of dark matter into the inner Milky Way which may still not be phase-mixed (e.g. Naidu et al. 2021;Han et al. 2022). Accounting for this dark matter would likely also lead to changes in the cross section constraints.

Systematic uncertainties affecting the study of the LMC
In Sec. 3.5 we cautioned that the use of the data-driven astrophysical templates for the LMC may artificially drive the DM bounds towards tighter contraints, which we aimed to avoid via the design of the analysis pipeline and the size of the chosen ROI. The Fermi-LAT collaboration's search for DM in the LMC (Buckley et al. 2015) with five years of data mitigated this effect in a different manner. To exemplify the results of this study, we show a selection of upper limits in the right panel of Fig. 7 displayed in three shades of blue representing different choices for the static LMC DM profile. These profiles are tuned to fit the rotation curve of the LMC as traced by stars and gas. The profile dubbed -is very similar to the initial LMC profiles in the NFW, HERN and Erkal19 simulations on which we base our work. In fact, the cyan line in Fig. 7 is comparable to the results from the dynamical LMC halo (red lines) but also to the bounds from the static DM halo profiles from Regis et al. (2021). The deviation of the latter bounds can be explained with the enlarged ROI compared to Buckley et al. (2015). Thus, we find corroborating evidence that our constructed analysis pipeline is not severely affected by a bias due to the astrophysical LMC templates.
The left panel of Fig. 7 illustrates and quantifies the level of the expected induced variation of upper limits on thermal DM when deformations of the LMC DM halo are included. The implications for the profiles of the obtained upper limits on the annihilation cross section vary between individual simulation suites. Including the LMC halo's deformations may either improve the constraints or weaken them by up to a factor of two. As pointed out in Sec. 5.1, the simulations themselves exhibit an intrinsic uncertainty regarding initial conditions and the definition of the MW/LMC morphology. This uncertainty consequently translates into a range of the J -factor maps of the LMC compatible with the simulated MW-LMC dynamics.
Regarding the obtained upper limits, however, the effect of deviations from spherical symmetry (a defining feature of the monopole of the BFE) must be understood on a case-by-case basis since the DM signal is degenerate with some of the astrophysical background gamma-ray components in the LMC ROI. Deformations of the LMC halo result in asymmetries of the J -factor maps, which may help to break (or worsen) these degeneracies. Considering, for example, the almost spatially uniform isotropic background it is clear that asymmetries in the J -factor maps greatly reduce its degeneracy with the DM template. On the flip-side, it is also conceivable that the morphology of a deformed LMC halo increases existing degeneracies as is the case in the NFW+NFW simulation from Lilleengen et al. (2022).
The degeneracies with the astrophysical background templates also explain the differences between the results for static NFW profiles and simulated LMC haloes displayed in the right panel of Fig. 7. On one hand, the simulation results exhibit less peaked DM density profiles towards the centre of the LMC than the profiles from Regis et al. (2021), thus reducing the features that may break degeneracies with background templates. On the other hand, the spectral shape of the DM signal is relevant too in order to improve the constraining power of the analysis. For light DM, for instance, the degeneracy can be broken by the information from the DM gamma-ray annihilation spectrum that shows a cutoff around the DM mass. This cutoff falls within the sweet spot of the LAT's sensitivity. Heavy DM with > 100 GeV, in contrast, features a break in the spectrum at tens of GeV where the LAT sensitivity starts to decrease. Hence, in this particular case, the spectral shape of the annihilation signal does not contribute as much to breaking the degeneracy with the background.

CONCLUSIONS
In the present work, we used state-of-the-art simulations of the MW encounters with the LMC to assess how the deformations of the MW and LMC dark matter haloes affect indirect detection for DM.
First, we focused on high Galactic latitudes and performed a search for a DM annihilation signal in twelve years of Fermi-LAT data. Since no significant signal was found (regardless of the DM spatial distribution adopted), we set 95% C.L. constraints on the DM annihilation cross section. In particular, we modelled the DM distribution in the Comparison between the constraints derived from simulated LMC haloes (red) using the full multipole expansion and similar DM halo profiles (orange) adopted from a recent study of the LMC in the radio band (Regis et al. 2021) quantifying the morphology of a static LMC. In addition, we provide the bounds derived in an independent gamma-ray analysis (Buckley et al. 2015)  Galactic halo following a recent MW mass modelling (McMillan 2017). The mass modelling of the MW comes with non-negligible uncertainties. We propagated this uncertainty on our final limits, and provided an uncertainty band which reflects it. Moreover, we verified that the systematic uncertainties from IE modelling are indeed mild as they induce a variance of the derived bounds by at most a factor of two for DM masses at the light and heavy end of the probed mass range. The optimal ROI sizes for each of the employed IE models are largely overlapping rendering a direct comparison sensible (see Appendix D for more details). The limits in Fig. 5 represent the most up-to-date and robust limits on DM at high latitudes derived with Fermi-LAT data. Our high-latitude DM limits are stronger than the results obtained by Zechlin et al. (2018) from their most constraining energy band from 1 to 2 GeV (c.f. orange line in the left panel of Fig. 5). As noted in Sec. 4.1, the authors performed their analysis in an ROI overlapping with ours but of reduced size. Thus, our enlarged ROI improves the constraining power allowing us to exclude thermal DM for masses 40 GeV, while Chang et al. (2018) report a slighter stronger bound for DM below 200 GeV and an exclusion for 70 GeV. To this end, the authors have derived an ROI that is rather disjoint with ours. It is closer to the Galactic centre and it takes into account gamma-ray data from the position of the FBs. The latter two differences can easily explain the increased constraining power compared to our study. We notice that our bounds do not account for the presence of sub-halos, which are expected to boost the annihilation signal, and therefore strengthening the limits, especially if very highly concentrated, see e.g. Delos & White (2022).
In the right panel of Fig. 5, we place our MW outer halo bounds in the context of existing constraints on the parameter space of thermal DM derived from different targets and cosmic-ray channels. The light purple line indicates gamma-ray constraints from the observation of dwarf spheroidal galaxies with space-borne and ground-based instruments. The results of this joint analysis (Abdalla et al. 2021) are the most state-of-the-art constraints from dwarf spheroidal galaxies using traditional inference techniques. While this set of exclusion limits outperforms our bounds over the entire probed mass range, they are affected by non-negligible systematic uncertainties due to both modelling of the DM distribution in these systems and the background modelling in and around dwarf spheroidal galaxies. As for the latter, the authors of Calore et al. (2018); Alvarez et al. (2020) designed an analysis that incorporates background modelling systematic effects, which weaken the upper limits on the DM annihilation cross section of about a factor of three. These limits (based only on classical dwarf spheroidal galaxies) are shown as a purple line. In this case, constraints from the outer MW halo are stronger than the gamma-ray limits from dwarf spheroidal galaxies. The variance induced by the uncertainty of the outer MW halo profile is hence less pronounced than the impact of background modelling uncertainties in dwarf spheroidal galaxy studies.
We stress that we assumed a very conservative approach in the selection of the region of interest for the analysis by requiring a strict compatibility between the statistical expectations derived from Poisson realisations of the baseline fit and the true Fermi-LAT data. This method yields reliable and robust limits but reduces the potential constraining power of the full data set and limits the accessible dynamically generated features of the MW-LMC system. As indicated by the right panel of Fig. 2, the northern sky may be indeed a more promising target to explore the signatures of the MW-LMC interaction. In this part of the sky, the LMC is expected to provoke a response in the forward direction of its orbit. However, finding a good agreement between astrophysical model and data is more challenging and needs better theoretical refinements. For example, if we relax the strict constraint of considering only data at high latitudes, we find that adding to the optimal ROI in Sec. 4 a counterpart in the northern hemisphere defined by −102 • < ℓ < 102 • and > 16 • (excluding the FBs region) yields the best accordance -although being far from perfect -between the expectations derived from the baseline fit and the true Fermi-LAT data. The improvement we can achieve in this way is as large as a factor of 3 for either a static or a simulated MW halo profile (see details in Appendix F).
In contrast, constraints (green band) from radio searches towards the LMC as found in Regis et al. (2021) show a larger intrinsic uncertainty than ours. Even in a conservative scenario, however, the bounds on thermal DM from the LMC's emission in radio light are stronger over the considered mass range. It should be noted that these results are derived assuming a static LMC DM profile while dynamical effects -as we have studied in this work -may alter the picture. The authors of Di Mauro & Winkler (2021) have derived a set of upper limits (dark blue) comparable to the radio LMC bounds considering the latest AMS-02 antiproton data release. Anti-proton constraints seem robust and only mildly affected by modelling uncertainties (DM halo profile, comsic-ray propagation). In the future, these bounds can be further strengthened with a re-calibrated prediction of secondary cosmic rays to achieve a better agreement with new data, which is further discussed in Calore et al. (2022a).
Thanks to the uncertainty band from the MW gravitational potential modelling, we quantified the significance of the MW-LMC dynamics' impact on indirect DM detection. For the set of simulations adopted in this work, we found that the MW-LMC interaction does not strongly affect the upper limits on thermal DM. The obtained variations are within a factor of 1.3, which we find in the Erkal19 simulation that generally yields the most pronounced dynamical responses.
The MW-LMC interactions also affect the mass distribution in the LMC itself -and the ensuing DM annihilation signal. Therefore, we derived bounds on DM annihilation from the region around the LMC as well, with dedicated modelling of the astrophysical backgrounds. The limits derived from the LMC ROI are largely consistent with the previous bounds derived by the Fermi-LAT collaboration for the case of a static NFW profile as discussed in Sec. 5.2. The level of systematic uncertainty of the constraints on thermal DM caused by the dynamical deformation of the LMC halo is not of the same order of magnitude as the one caused by the allowed range of static DM halo profiles that can reproduced the stellar rotation curve of the LMC as illustrated by the range of the blue-shaded lines in Fig. 7. However, we note that the initial states of the LMC in all simulations utilised in this work did not aim at bracketing the full margin of inner halo profiles consistent with stellar data. We can thus not quantify how the LMC's internal structure is affected by its passage through the MW in case of an aggressive assumption like the -parametrization in Buckley et al. (2015). Refined simulations with a wider range of inner DM halo slopes may be warranted to assess the full implications of the MW-LMC dynamics for indirect searches towards the LMC. Eventually, the effect may even alter the prospects of radio searches in the central region of the LMC which consequently relax or tighten the already stringent bounds on thermal DM reported by Regis et al. (2021).
In conclusions, we have shown that high Galactic latitudes have the potential to be the leading target for DM searches in gamma rays in the future. A crucial step in this direction will be the optimisation of interstellar emission models thanks to, for example, machine learning techniques and verification schemes, see e.g. Storm et al. (2017); Mishra-Sharma & Cranmer (2022); Caron et al. (2022) in the context of background model optimisation in the inner Galaxy. If we focus on the LMC region, a leap forward in the understanding of the astrophysical emission (and therefore a better modelling of the LMC astrophysical templates) is expected thanks to the up-coming observations of the Cherenkov Telescope Array, CTA, of this particular region as outlined in Acharya et al. (2018).
Looking forward, future work with stellar streams and other tracers will allow the community to robustly measure the dark matter deformations of the Milky Way and LMC dark matter haloes through their gravitational effects. Indeed, Shipp et al. (2021) show that the streams have their closest approaches with the LMC at different times, giving hope to the idea that the time dependence can be measured. Once these deformations are measured, they can be folded into the analysis, as we have done in this work, to derive the most accurate annihilation cross sections. On the flip side, if annihilating dark matter is detected in gamma rays then we will be able to directly measure the deforming dark matter haloes of the Milky Way and LMC.
For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.  (Kluyver et al. 2016), (Hunter 2007), (Harris et al. 2020), (Virtanen et al. 2020)

DATA AVAILABILITY
A python interface to integrate orbits and access the expansion model for the Erkal19 simulation can be found here: https://github. com/sophialilleengen/mwlmc. The code to reproduce the results of the Fermi-LAT gamma-ray analysis is provided at: https:// github.com/ceckner/MWLMCGammaRays.

APPENDIX A: ASTROPHYSICAL FORE-AND BACK-GROUND COMPONENTS
We here provide more details about the astrophysical fore-and background components used to fit the gamma-ray sky.
• interstellar emission (IE): A diffuse Galactic gamma-ray source emerging due to very-high energy, charged cosmic rays impinging on particles of the MW's interstellar medium, dust and radiation fields. The dominant processes that create gamma rays are 0 -decay, Bremsstrahlung and IC scattering. The modelling of this contribution is subject to many uncertainties so that we include five different models that aim to quantitatively characterise the IE. Two IE models (henceforth called Lorimer I and Lorimer II) are taken from the set of realisations considered in the "1st Fermi LAT Supernova Remnant Catalog" 7 Acero et al. (2016). While a detailed description of the models' preparation and their respective properties are given in the cited catalogue publication, we restrict ourselves to a brief summary of their main characteristics relevant to this work: The distribution of primary cosmic rays is linked to the distribution of pulsars as analysed and discussed in Lorimer et al. (2006). They are confined in a volume with propagation height = 10 kpc whereas the spin temperature of the interstellar medium is assumed to be = 150/1 · 10 5 K (Lorimer I/II). The interstellar medium's gas content is split in atomic hydrogen H and CO maps. The latter serves as a proxy for the distribution of molecular hydrogen H 2 . The total volume of the gas maps is decomposed into four Galactocentric annuli with respective extensions: 0-4 kpc: "ring 1", 4-8 kpc: "ring 2", 8-10 kpc: "ring 3" and 10-30 kpc: "ring 4". Since the IE from different rings contributes to different parts of the gamma-ray sky, we make use of this decomposition in our fitting strategy outlined in Sec. 3.4. We note that the properties of IE model Lorimer I are comparable to those of the official diffuse background model of the Fermi-LAT derived in connection with the 4FGL catalogue Abdollahi et al. (2020); Ballet et al. (2020).Therefore, we declare Lorimer I our benchmark IE model. The three remaining models -called foreground model A, B and C -are adopted 8 from a careful and detailed study of the diffuse extragalactic gamma-ray background conducted by the Fermi-LAT collaboration Ackermann et al. (2015). In contrast to Lorimer I and II, these IE instances exhibit the advantage to be directly prepared for the study of a large-scale emission component at high Galactic latitudes -a feature aligned with the aim of our analysis. Again, we refer to the cited publication to learn more about the exact composition of the models.
• isotropic diffuse gamma-ray background (IGRB): This largescale contribution to the gamma-ray sky is spatially isotropic, hence following the spatial structure of the LAT's exposure, and generated by the collective emission of distant extragalactic gamma-ray emitters too faint to be resolved individually. We adopt the IGRB spectrum associated with the selected LAT data's event class and type (see Sec. 3.1) provided by the Fermi Science Tools 9 . We stress that these spectra are only valid in combination with the official diffuse background model of the Fermi-LAT collaboration. Since our fit model (Eq. 3) renormalises each background component per energy bin, this restriction is, however, irrelevant for our cause.
• resolved point-like and extended gamma-ray sources: Based on a 10-year data set, the Fermi-LAT collaboration has published an extensive catalogue, 4FGL-DR2, (Abdollahi et al. 2020;Ballet et al. 2020) of all resolved and localised point-like or extended gamma-ray emitters in or outside of the MW. We include all of these sources with their respective properties reported in the catalogue in our analysis. A detailed description of the treatment of these sources in each analysis step is given in the following Sec. 3.4.
• LoopI: A large-scale, loop-like structure exhibiting a diffuse gamma-ray emission mainly concentrated in the northern hemisphere above the Galactic disc. We include in our model the spatial and spectral characterisation of LoopI given in Wolleben (2007). This particular model also features a non-vanishing contribution in the southern hemisphere of the gamma-ray sky.
• gamma-ray emission induced by the Sun and Moon: Since the LAT is in constant observation mode, both the Sun and the Moon cross the field of view of the telescope. Along their orbital trajectory they contribute a non-negligible gamma-ray emission. We make use of dedicated routines within the Fermi Science Tools 10 to derive a flux model for the gamma-ray emission of these two celestial bodies tailored towards the selected LAT data (see Johannesson & Orlando (2013)).

APPENDIX B: SELECTED BRIGHT 4FGL-DR2 SOURCES
We report in Tab. B1 all sources in 4FGL-DR2 whose energy flux (100 MeV -100 GeV) is above the threshold of 4 · 10 −10 MeV cm −2 s −1 . These sources are considered as individual templates during the iterative fit (see Sec. 3.4) designed to derive a baseline model of the gamma-ray sky solely comprised of the astrophysical background components introduced in Sec. 3.3. We state their name in 4FGL-DR2, their position in Galactic coordinates, the nominal energy flux and the source class using abbreviations from Abdollahi et al. (2020).

APPENDIX C: INJECTED DARK MATTER SIGNAL RECOVERY
In order to test the performance of the analysis pipeline outlined in Sec. 3.4 and to verify its trustworthiness, we conduct a simple signal recovery exercise: Utilising the baseline fit with IE model Lorimer I on the all-sky gamma-ray data, we inject a DM signal into the model and try to detect it with the chosen inference method, namely the loglikelihood ratio test statistic. The test statistic in Eq. 5 is not suited for this task since it is designed to constrain an alternative hypothesis when the data seems to prefer the background-only hypothesis. Thus, we modify the test statistic for the case of determining the detection significance according to (Cowan et al. 2011 Table B1. Summary table listing bright 4FGL-DR2 sources with an energy flux 100 4 · 10 −10 MeV cm −2 s −1 that are fit individually during the iterative fit procedure aimed to derive a baseline model of the gamma-ray sky. The table states the source's name in 4FGL-DR2, its position in Galactic longitude ℓ and latitude in degree, the nominal energy flux 100 and the source class using abbreviations from Abdollahi et al. (2020).
where we keep the notation established in Sec. 3.2. Succinctly put, this test statistic assumes that the preferred hypothesis contains a positive signal with a positive best-fit normalisationˆD M . The alternative hypothesis, in contrast, becomes the background-only hypothesis, that is, finding the best-fit background parameters under the assumption that DM = 0. The likelihood ratio thus quantifies the significance of the detected signal and for a given threshold one may claim a detection at the · level. We have applied this prescription to ascertain that no significant signal is present in the selected LAT data for all scenarios provided in the main text. Hence, deriving upper limits on the WIMP DM pair-annihilation cross section is justified.
For the purpose of our pipeline check, we conduct the following adapted approach: • For the monopole and full BFE MW + LMC haloes of the Erkal19 simulation, we prepare a signal template featuring a DM particle with mass = 100 GeV annihilating into¯pairs. • For discrete points in the range of annihilation cross section DM = ∈ 10 −28 , 10 −23 cm 3 s −1 , we inject the signal template with normalisation DM into the baseline fit of Lorimer I.
• Drawing 200 Poisson realisations of the baseline fit + signal gamma-ray sky, we perform a maximum likelihood fit with respect to this mock data and save the retrieved best-fit signal strengthˆD M .
• Again, using 200 Poisson realisations of the baseline fit + signal gamma-ray sky, we derive the upper limit on the DM annihilation cross section that we may set given the presence of a signal in the utilised mock data.
The results of this sanity check are displayed in Fig. C1, which shows the median/scatter of the recovered signal normalisation (red solid/shaded band) and the median/scatter of the upper limits (green solid/shaded band) as a function of the injected signal's strength. To guide the eye, we also denote the upper limit on this particular particle DM model derived from the baseline fit without injected signal as orange line; its 68% containment band is given as an orange-shaded band.
The general expectation for this kind of pipeline check is to recover the injected signal strength; the higher the signal strength the higher the confidence of recovering the signal, i.e. the smaller the observed scatter of the best-fit parameter. Moreover, for extremely small annihilation cross sections the upper limits with respect to mock data that contain the signal should asymptotically approach the corresponding upper limits with respect to the baseline fit standalone. Both of these features we can confirm with our sanity check. In addition, we find that the solid red line approaches the dashed black line for cross section values below the nominal upper limit marked by the vertical orange line. This suggests that the pipeline is capable of theoretically detecting a signal less luminous than the obtained upper limit. At the same time, the solid green line starts to deviate from the horizontal orange line when the injected signal strength can be recovered. We can hence conclude that the constructed analysis pipeline works as intended and it is suitable to perform the envisaged tasks. . Signal recovery test for WIMP particle DM of mass = 100 GeV annihilating into¯pairs and a DM halo featuring the full BFE of the Erkal19 simulation to model the spatial morphology. The signal template has been injected with a normalisation DM following the dashed black line into the baseline fit for IE model Lorimer I. The solid red line denotes the median value of the recovered best-fit signal normalisation for 200 Poisson realisations drawn from the mock data while the red band illustrates the expected scatter (68% containment) of the DM best-fit parameter. The solid green line represents the median 95% C.L. upper limit that follows from the same mock data with injected DM template; the 68% containment is shown as green band. Finally, the solid orange line is the median 95% C.L. upper limit on mock data that does not contain a signal contribution. Again, the orange band denotes the associated 1 scatter.

APPENDIX D: DARK MATTER CONSTRAINTS FOR DIFFERENT INTERSTELLAR EMISSION MODELS
The main results in Sec. 4 of our analysis of the MW-LMC system rely on the choice of a particular benchmark IE model. However, it is well known that the characterisation of the interstellar emission is not perfect and a great deal of data-and observation-driven as well as simulated models have been put forward to achieve a sufficiently good agreement between reality and theoretical understanding plus observational data. On one hand, since our study mostly focuses on the high-latitude in the southern hemisphere of the gamma-ray sky, our results are less severely affected by the IE directly related to the emission originating from the interstellar medium along the Galactic disc. On the other hand, gamma-ray emission due to IC scattering events is particularly hard to model and exhibits large uncertainties so that the impact of the IE model on our results is certainly nonnegligible. As motivated in Sec. 3.3, we aim to assess the systematic uncertainty of our results due to the modelling of the IE via alternative models that supplement the benchmark choices for either the outer MW halo study or the dedicated study of the LMC surroundings. In Fig. D1, we provide a set of 95% C.L. upper limits on the DM pair-annihilation cross section (again, with respect to the channel →¯) for different choices of the IE model.
The left panel of this figure concerns the impact of the IE on the study of the outer MW halo, where we exemplify the performance of all five IE models with the signal morphology according to the Jfactor map of the Erkal19 simulation including all terms of the BFE for both DM haloes. We checked that the choice of the simulation is not important to quantify the impact of the IE uncertainty; the other four simulations as well as the static MW halo from McMillan (2017) yield very similar results. Using a different IE model than Lorimer I as in the main text, requires us to re-perform the ROI optimisation, which yields optimal ROI sizes detailed in Tab. D2 for the respective case. As concerns the comparison of upper limits derived from real data and the baseline fit, we find consistency for almost the entire probed DM mass range at the 2 level and even better for FGMA to FGMC; an expected behaviour since these three models were initially created to facilitate searches at high latitudes. A slight fluctuation to stronger constraints for light DM with masses below 20 GeV is a common feature among all probed IE models. Overall, the numerical values of the derived constraints are comparable to the corresponding upper limits for the same simulation set in Fig. 6 but IE model Lorimer I. The apparent differences certainly partially arise because of the varying optimal ROI sizes that reduce or increase the total J -factor.
The right panel of Fig. D1 addresses the question of the IE model's impact on the study of the LMC region. As explained in Sec. 3.5, we include for this particular task the Galactic diffuse background model of the Fermi-LAT collaboration. Besides, we do not show the upper limits for FGMB, as it appears to behave exactly like FGMC, and Lorimer II, which does not provide a good fit to the data in general. In this case, we use the simulation based on a Hernquist profile for the MW and an NFW profile for the LMC as initial conditions from Ref. Donaldson et al. (2022). This example suffices to quantify the impact of the IE modelling as the remaining LMC halo models exhibit similar behaviour. Although we display the set of upper limits for each particular IE model on the same plot, the underlying optimised ROI sizes are different and all applied values are given in Tab. D1. However, the differences in those sizes are not noticeably altering the total J -factor we are probing so that a common plot is justified. As it turns out, varying the IE has a remarkably small effect on the resulting constraints on WIMP DM. The only exception being the IE model Lorimer I, which induces a deterioration of the limits by a factor of ∼ 2 compared to the other three IE model instances. Hence, the results of the dedicated LMC study are robust against variations in the IE modelling.  Table D1. Summary of the optimised ROI sizes for the dedicated study of the LMC environment by varying the underlying IE model. The optimisation has been performed under the assumption of the LMC halo shapes according to the simulation in Ref. Donaldson et al. (2022) using an Hernquist (MW) and NFW (LMC) profile as initial conditions.  Figure D1. 95% C.L. upper limits on the velocity-weighted thermally averaged DM pair-annihilation cross section depending on the assumed DM mass for the prompt gamma-ray emission from DM annihilating into¯final states. (Left: analysis of the outer MW halo) Variance of the observed upper limits due to changing the IE model. As an illustrative example, we display the signal morphology following the J-factor of the combined MW and LMC halo according to the Erkal19 simulation including the full BFE. Note that each IE model requires a different optimised ROI. The corresponding values are provided in Tab. D2. (Right: analysis of the LMC region) Impact of varying the IE model used to derive the constraints on DM. We show four distinct IE models: Lorimer I (dashed), Fermi-LAT IE model (v07) (solid), FGMA (dashed) and FGMC (dotted). Black lines represent the unperturbed LMC halo whereas red lines display the respective scenario including dynamical effects based on the halo profiles derived from the simulation based on a Hernquist profile for the MW and an NFW profile for the LMC as initial conditions from Ref. Donaldson et al. (2022). The optimal ROI sizes for each scenario are stated in Tab Table D2. Summary of the optimised ROI sizes for the five IE models used in the analysis of the outer MW halo. These optimised ROIs are derived with respect to the spatial morphology of the MW-LMC system found in the Erkal19 including the dynamical response of both DM haloes.

APPENDIX E: DARK MATTER CONSTRAINTS FOR AN ALTERNATIVE ANNIHILATION CHANNEL
In this appendix, we present our results for an alternative spectral DM model. Instead of the rather soft¯-channel, we here assume 100% annihilation into + − pairs that yield a hard annihilation spectrum. We repeat the analysis with respect to the optimised ROIs to re-derive the main results that have been shown in the main text in Fig. 6 for the MW-LMC system and in Fig. 7 for the LMC standalone. The corresponding constraints for the + − -channel are provided in Fig. E1.
As a general observation, the upper limits for DM particles with masses above a few tens of GeV are weaker compared to the¯-channel. This result is reasonable given the fact that the + − annihilation spectrum is harder, which shifts the most constraining part of the spectrum to higher energies. At energies of a few GeV, Fermi-LAT analyses tend to be statistics limited. As a consequence, their constraining power regarding hard spectra peaking above O (10) GeV is reduced.

APPENDIX F: DARK MATTER CONSTRAINTS FOR AN ENLARGED REGION OF INTEREST
In the main text of this work, we focused on determining an ROI in the southern hemisphere with optimal properties regarding the compatibility of the statistical expectations derived from our baseline fit and the true Fermi-LAT data set. On one side, the dynamics of the MW-LMC system generate anisotropies that are located in this particular fraction of the entire sky as illustrated by Figs. 1 and 2.
On the other side, such deformations -mostly associated with the LMC halo -are likewise present in the northern hemisphere. Their impact on indirect gamma-ray searches for DM in the outer MW halo remains unprobed.
To investigate the potential of an extended ROI that incorporates patches in both hemispheres of the gamma-ray sky, we fix the optimal southern ROI derived for the benchmark IE model Lorimer I and reperform the optimisation routine in the northern sky. With respect to the example of a static MW halo referenced in Fig. 4, we obtain the best match between expectation from mock data and performance on real Fermi-LAT data for an ROI defined by −102 • < ℓ < 102 • and > 16 • (shown in the left panel of Fig. F1). The quality of this accordance is, however, not as good as for the southern hemisphere standalone. We quantify the suitability of the obtained extended ROI in the right pannel of Fig. F1 by confronting the DM upper limits from true Fermi-LAT data (red) with the corresponding constraints (black) and their scatter from mock data. The region in the parameter space where we can exclude thermal DM enlarges by a factor of about two compared to Fig. 4. However, the derived bounds are not within the 2 containment band of the upper limits derived from the baseline fit in a broad range of the probed DM parameter space. This reduces the credibility of this set of upper limits and it furthermore indicates that our baseline astrophysics model is not entirely describing reality in such an extended sky region. In order to exploit the full potential of searches towards the outer MW halo it is thus essential to improve the current state-of-the-art of astrophysical models at high latitudes; in particular, the gamma-ray emission of extended diffuse components as the FBs or loop-like structures and the residual IE.
As concerns the sensitivity of such indirect searches to deformations of the MW and LMC haloes, we find ROIs of similar extensions in the northern hemisphere for the Erkal19 simulation. The statistical robustness that these sky regions provide is comparable to the case shown in Fig. F1. Translated to constraints on the DM annihilation cross section, we do not observe an enhanced difference between a signal morphology following the monopole term or the full BFE of the Erkal19 model. The effect of the MW-LMC dynamics remains at the level displayed in Fig. 6. To stress it again, the reliability of this statement is suffering from the lack of a good fit to reality. Improving the astrophysical background modelling at high latitudes may also improve the sensitivity to the dynamical response of the LMC passage through the MW. This paper has been typeset from a T E X/L A T E X file prepared by the author.