# The morphology and structure of stellar populations in the Fornax dwarf spheroidal galaxy from Dark Energy Survey Data

M.Y. Wang, T. de Boer, A. Pieres, T.S. Li, A. Drlica-Wagner, S.E. Koposov, A.K. Vivas, A.B. Pace, B. Santiago, A.R. Walker, et al.

► To cite this version:

M.Y. Wang, T. de Boer, A. Pieres, T.S. Li, A. Drlica-Wagner, et al.. The morphology and structure of stellar populations in the Fornax dwarf spheroidal galaxy from Dark Energy Survey Data. The Astrophysical Journal, 2019, 881, pp.118. 10.3847/1538-4357/ab31a9. hal-01885538

HAL Id: hal-01885538

https://hal.science/hal-01885538v1

Submitted on 9 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Morphology and Structure of Stellar Populations in the Fornax Dwarf Spheroidal Galaxy from Dark Energy Survey Data

M. Y. Wang<sup>1,2</sup>, T. de Boer<sup>3,4</sup>, A. Pieres<sup>5,6</sup>, T. S. Li<sup>7,8</sup>, A. Drlica-Wagner<sup>7</sup>, S. E. Koposov<sup>1,4</sup>, A. K. Vivas<sup>9</sup>, A. B. Pace<sup>2</sup>, B. Santiago<sup>6,10</sup>, A. R. Walker<sup>9</sup>, D. L. Tucker<sup>7</sup>, L. Strigari<sup>2</sup>, J. L. Marshall<sup>2</sup>, B. Yanny<sup>7</sup>, D. L. DePoy<sup>2</sup>, K. Bechtol<sup>11</sup>, A. Roodman<sup>12,13</sup>, T. M. C. Abbott<sup>9</sup>, F. B. Abdalla<sup>14,15</sup>, S. Allam<sup>7</sup>, J. Annis<sup>7</sup>, S. Avila<sup>16</sup>, E. Bertin<sup>17,18</sup>, D. Brooks<sup>14</sup>, D. L. Burke<sup>12,13</sup>, A. Carnero Rosell<sup>5,6</sup>, M. Carrasco Kind<sup>19,20</sup>, C. E. Cunha<sup>5</sup>, C. B. D'Andrea<sup>21</sup>, L. N. da Costa<sup>5,6</sup>, J. De Vicente<sup>22</sup>, S. Desai<sup>23</sup>, T. F. Eifler<sup>24,25</sup>, J. Estrada<sup>7</sup>, B. Flaugher<sup>7</sup>, J. Frieman<sup>7,8</sup>, J. García-Bellido<sup>26</sup>, D. W. Gerdes<sup>27,28</sup>, D. Gruen<sup>12,13</sup>, R. A. Gruendl<sup>19,20</sup>, G. Gutierrez<sup>7</sup>, D. L. Hollowood<sup>29</sup>, K. Honscheid<sup>30,31</sup>, D. J. James<sup>32</sup>, K. Kuehn<sup>33</sup>, N. Kuropatkin<sup>7</sup>, O. Lahav<sup>14</sup>, M. A. G. Maia<sup>5,6</sup>, R. Miquel<sup>34,35</sup>, E. Sanchez<sup>22</sup>, V. Scarpine<sup>7</sup>, I. Sevilla-Noarbe<sup>22</sup>, M. Smith<sup>36</sup>, R. C. Smith<sup>9</sup>, F. Sobreira<sup>37,6</sup>, E. Suchyta<sup>38</sup>, M. E. C. Swanson<sup>20</sup>, and G. Tarle<sup>28</sup>

(DES Collaboration)

<sup>1</sup> Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

<sup>2</sup> George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA

<sup>3</sup> Department of Physics, University of Surrey, Guildford GU2 7XH, UK

<sup>4</sup> Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

<sup>5</sup> Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil

<sup>6</sup> Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil

<sup>7</sup> Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA

<sup>8</sup> Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

<sup>9</sup> Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

<sup>10</sup> Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS - 91501-970, Brazil

<sup>11</sup> LSST, 933 North Cherry Avenue, Tucson, AZ 85721, USA

<sup>12</sup> Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, CA 94305, USA

<sup>13</sup> SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

<sup>14</sup> Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

<sup>15</sup> Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown, 6140, South Africa

<sup>16</sup> Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

<sup>17</sup> CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>18</sup> Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>19</sup> Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

<sup>20</sup> National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA

<sup>21</sup> Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

<sup>22</sup> Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

<sup>23</sup> Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

<sup>24</sup> Department of Astronomy/Steward Observatory, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

<sup>25</sup> Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA

<sup>26</sup> Instituto de Física Teorica UAM/CSIC, Universidad Autonoma de Madrid, E-28049 Madrid, Spain

<sup>27</sup> Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

<sup>28</sup> Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

<sup>29</sup> Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

<sup>30</sup> Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

<sup>31</sup> Department of Physics, The Ohio State University, Columbus, OH 43210, USA

<sup>32</sup> Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA

<sup>33</sup> Australian Astronomical Observatory, North Ryde, NSW 2113, Australia

<sup>34</sup> Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

<sup>35</sup> Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona) Spain

<sup>36</sup> School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

<sup>37</sup> Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil

<sup>38</sup> Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

Received 2018 September 21; revised 2019 July 9; accepted 2019 July 10; published 2019 August 21

## Abstract

Using deep wide-field photometry 3 yr data (Y3) from the Dark Energy Survey (DES), we present a panoramic study of the Fornax dwarf spheroidal galaxy. The data presented here—a small subset of the full survey—uniformly cover a region of 25 deg² centered on the galaxy to a depth of g ~ 23.5. We use these data to study the structural properties of Fornax, overall stellar population, and its member stars in different evolutionary phases. We also search for possible signs of tidal disturbance. Fornax is found to be significantly more spatially extended than what early studies suggested. No statistically significant distortions or signs of tidal disturbances were found down to a surface brightness limit of ~32.1 mag arcsec⁻². However, there are hints of shell-like features located ~20'–40' from the center of Fornax that may be stellar debris from past merger events. We also find that intermediate-age and young main-sequence populations show different orientation at the galaxy center and have many substructures. The deep DES Y3 data allow us to characterize the age of those young stellar substructures with great accuracy, both those previously known and those newly identified as possible overdensities in this work, on the basis of their
color–magnitude diagram morphology. We find that the youngest overdensities are all found on the eastern side of Fornax, where the Fornax field population itself is slightly younger than in the west. In summary, the high-quality DES Y3 data reveal that Fornax has many rich structures and provide insights into its complex formation history.

**Key words:** galaxies: dwarf – galaxies: individual (Fornax) – Local Group

# 1. Introduction

According to the standard $\Lambda$ cold dark matter cosmology, the universe as we see it today is built via a scale-free process in which matter collapses under gravity and iteratively forms larger and larger structures. Dark matter dominates this hierarchical structure formation process. Dwarf galaxies, the earliest formed structures and the most dark-matter-dominated objects known, are therefore crucial laboratories for studying galaxy formation and evolution processes (e.g., Kauffmann et al. 1993). Considerable efforts have been invested in the study of dwarf galaxies in the Local Group (LG), as summarized recently in McConnachie (2012).

Even though they have small physical sizes, the proximity of these galaxies results in large angular extents on the sky, of approximately a few degrees. Hence, the study of these systems has particularly flourished since the advent of wide-area imagers and multiobject spectrographs, which has allowed us to determine the properties of their resolved stellar component in great detail out to their very low surface brightness outskirts. This applies in particular to the “brightest” of the Milky Way (MW) dwarf spheroidal (dSph) galaxies, those with luminosities $-15 \lesssim M_V \lesssim -8$, for which one can gather wide-area photometric data sets with large statistics to reach well below the oldest main-sequence turn-off (MSTO) and spectroscopic data sets to provide velocities and metallicities of several hundreds of stars per galaxy. It is only through a joint analysis of the entire stellar sample that one can obtain a complete picture of the morphology and populations of dwarf galaxies.

The Fornax dSph is one of the more intensively studied dSphs of the MW. First discovered by Shapley (1938), it is one of the most luminous LG dwarf galaxies, second only to the Sagittarius dSph in terms of luminosity ($M_V = -13.4 \pm 0.3$) among the MW classical dSphs. Fornax has a total (dynamical) mass of $1.6 \times 10^8 M_\odot$ at a distance of $138 \pm 8$ kpc ($m-M)_V = 20.84 \pm 0.04$; Walker et al. 2006; Rizzi et al. 2007; Łokas 2009; Pietrzyński et al. 2009). Fornax is one of the few LG dwarf galaxies that has experienced a composite, extended star formation history covering ancient (>10 Gyr), intermediate-age (2–8 Gyr), and extremely young (100 Myr) stars (e.g., Gallart et al. 2005; de Boer et al. 2012; del Pino et al. 2013; Weisz et al. 2014). Spectroscopic studies of individual stars have found Fornax to have a broad metallicity distribution with a dominant metal-rich ([Fe/H] $\approx -0.9$ dex) component (e.g., Mateo et al. 1991; Tolstoy et al. 2001; Shetrone et al. 2003; Pont et al. 2004; Battaglia et al. 2006; Letarte et al. 2010; Suda et al. 2017). Studies of the kinematics have further shown the presence of at least three kinematically distinct populations (Battaglia et al. 2006; Amorisco & Evans 2012), indicating that it has undergone a complex formation process. Finally, Fornax also contains several globular clusters (Shapley 1938; Hodge 1961a), some of which are among the most metal-poor clusters found in the LG (Buonanno et al. 1998; Letarte et al. 2006; Larsen et al. 2012; de Boer et al. 2013).

Fornax can be easily resolved into individual stars down to the base of the red giant branch (RGB) with relatively short

exposures on 4 m class telescopes. However, given its substantial size on the sky ($r_{tidal} = 1^\circ 18 \pm 0^\circ 07$ or $2.85 \pm 0.16$ kpc), a detailed study of its stellar content requires wide-field imaging capabilities. Wide-field studies of Fornax made with photographic plates by Hodge (1961b) revealed that the ellipticity of Fornax isophotes increases with radial distance. The stellar density has also been shown to be asymmetric, with a peak density offset from the central position (de Vaucouleurs & Ables 1968; Hodge & Smith 1974). Resolved stellar studies of the stellar content within the tidal radius have shown that Fornax contains a radial population gradient, with younger, more metal-rich stars found preferentially toward the center (Battaglia et al. 2006; Coleman & de Jong 2008; de Boer et al. 2012).

Substructure studies of Fornax have found a number of significant overdensities, which provide further signs of its complex formation scenario and possible accretion of (enriched) gas or subhalos. Two young (<2 Gyr) overdensities were found by Coleman et al. (2004, 2005), one of which was later found to be due to background galaxies (Bate et al. 2015). A further three overdensities were found more toward the center of Fornax, all containing young stellar populations and metallicities (as derived from their color–magnitude diagrams [CMDs]) similar to the young Fornax field stars (de Boer et al. 2013; Bate et al. 2015). This may indicate a connection between the overdense components and the field stars (Olszewski et al. 2006), providing tentative clues to the formation of Fornax, possibly through reaccretion of previously expelled gas (e.g., Pasetto et al. 2011; Nichols et al. 2012; Yozin & Bekki 2012).

Despite the copious amount of work done in classifying the stellar content of Fornax, a detailed study of its extratidal structure has been lacking owing to the need for deep ($g \approx 24$), wide-field, and homogeneous data in multiple filters, necessary to constrain the distribution of very low surface brightness features in the galaxy outskirts. A study covering the outskirts of Fornax was conducted recently using data from the VST/ATLAS survey (Bate et al. 2015), which found no evidence of strong extratidal features. However, their data were only marginally able to reach the MSTO of Fornax, resulting in only a low number of tracer stars beyond the tidal radius. In contrast, the Fornax system and its environs have been covered recently by the Dark Energy Survey (DES) using deep and homogeneous imaging from the Dark Energy Camera (DECam; Flaugher et al. 2015) mounted on the 4 m Blanco Telescope at the Cerro Tololo Inter-American Observatory in Chile (Flaugher et al. 2015). The resulting data are deeper than any previous wide-field studies such as Bate et al. (2015) and therefore allow for a better selection of different stellar tracer features from the CMD to reduce contamination in its low-luminosity outskirts.

In this article we present an extensive wide-field study of Fornax’s stellar population using deep $g$- and $r$-band photometry reaching down to the MSTO stars from DES Y3 data that cover 25 deg$^2$ of the Fornax dSph. We perform a detailed and quantitative analysis of the structural properties of this galaxy.
In particular, we investigate the extratidal stellar content of Fornax to probe for the existence of tidal features. Given its brightness and stellar kinematics, Fornax might be one of the most likely galaxies in the LG to be affected by MW tidal forces and likely to contain distorted outer isophotes and possible tidal tails (Wang et al. 2017), although different conclusions are reached in different works (Battaglia et al. 2015).

The article is organized as follows. In Section 2 we will describe the specifics of the DES survey and the reduction of the data. Section 3 will discuss the structural parameters of the Fornax dwarf galaxy, followed by a discussion of the features in the CMD in Section 4. The surface density maps and search for tidal features are described in Section 5, and the spatial distribution of different stellar populations is studied in Section 6. Finally, Section 7 discusses the search for stellar substructures, followed by the summary and conclusions in Section 8.

## 2. Data

DES is a deep, wide-area imaging survey in the grizY bands performed with the DECam mounted on the 4 m Blanco Telescope at the Cerro Tololo Inter-American Observatory in Chile. Here we use wide-field imaging data from the first 3 yr of DES operation (DES Y3). DES Y3 is the first DES data set to cover the full DES footprint, which covers ~5000 deg² of the southern Galactic cap, including the Fornax dSph galaxy area.

The details of the DES Y3 data reduction can be found in Morganson et al. (2018), and the definition of the star–galaxy classification model can be found in Shipp et al. (2018). Here we briefly describe some of the procedures. The DES Y3 images were processed by the DES data management pipeline. Individual exposures were remapped to a consistent pixel grid and co-added to increase imaging depth (Morganson et al. 2018). Object detection was performed on a combination of the r+i+g co-added images using the SExtractor toolkit (Bertin & Arnouts 1996; Bertin et al. 2002) with an object detection threshold of ~5σ (Morganson et al. 2018). The ngmix method (Sheldon 2014), which constructed and modeled the point-spread function (PSF) of each object simultaneously to all available epochs and bands to derive the quoted magnitudes and errors, was applied to overcome the difficulties of performing precise photometric and morphological measurements in the co-adds.

To further select a high-quality stellar sample, the ngmix method was used to fit a composite galaxy model (bulge plus disk) to each source in all bands simultaneously (Drlica-Wagner et al. 2018). The best-fit size, CM_T, and associated uncertainty, CM_T_ERR, from this galaxy-model fit were then used to distinguish point-like objects from those that are spatially extended. Specifically, we defined an extended classification variable, NGMIX_CLASS, based on the sum of three selection criteria,

$$ \text{NGMIX\_CLASS} = \left( \begin{aligned} & (\text{CM\_T} + 5 \text{CM\_T\_ERR}) > 0.1 \\ & + ((\text{CM\_T} + \text{CM\_T\_ERR}) > 0.05) \\ & + ((\text{CM\_T} - \text{CM\_T\_ERR}) > 0.02). \end{aligned} \right) \quad (1) $$

Our final stellar sample used NGMIX_CLASS ≤ 1 criteria. This stellar selection was designed to yield a compromise between completeness and purity in the resulting stellar sample. It was compared with deeper imaging data from Hyper Suprime

Cam DR1 (Aihara et al. 2018; see Figure 1 in Shipp et al. 2018). It is found that, in uncrowded regions, this selection is >90% complete for g = 23.5 with a galaxy contamination rising from ≲5% at g ≲ 22.5 to ~30% by g ~ 23.5. We constrained our stellar sample to the range 19 < g, r < 23.5. The faint-end limit was imposed to avoid spurious density fluctuations resulting from inhomogeneous survey depth and galaxy contamination.

To account for interstellar dust extinction, we followed the procedure described in Abbott et al. (2018). We started with E (B − V) values from the reddening map of Schlegel et al. (1998). We computed fiducial interstellar extinction coefficients and integrated over the DES standard bandpasses considering a fixed source spectrum that is constant in spectral flux density per unit wavelength. The resulting multiplicative coefficients for the g and r band are R_g = 3.185 and R_r = 2.140. Throughout this paper, all magnitudes refer to extinction-corrected PSF magnitudes derived by ngmix.

To assess the crowding issues in Fornax dSph galaxy in the DES photometry, which may affect studies such as galaxy central profile fitting, we implement a pipeline to run DAOphot (Stetson 1987) photometry. The DAOphot approach is designed to better recover stellar sources in crowded fields than SExtractor, although the drawback is low confidence for star/galaxy separation. The DAOphot pipeline was implemented for the DES Y3 co-added images (Morganson et al. 2018) in g, r, and i bands, in a total area of 6.572 deg² on the sky. After the PSF photometry using DAOphot tasks was run, the sources were matched using data from three bands, with a maximum deviation of 1". We also add a zero-point with no color terms. Typical values for the standard deviation of the zero-point are 30 mmag. Similar zero-points (with a maximum deviation of 0.03 mag) were found compared to the ngmix PSF photometry from DES Y3 catalogs. Comparing the sources from the DAOphot catalog to the DES Y3 sources within magnitude range g, r, i = 17–21, the magnitude errors are less than 0.03. The final DAOphot catalog presents a continuous coverage of the sky and successfully recovers the five known globular clusters of the Fornax dSph galaxy, for which we find that some of them have severe crowding issues in the DES ngmix photometry. We find a mild difference, ~30%, in the Fornax central stellar density profile (r ≲ 10') between the DAOphot and the DES ngmix catalogs, and the comparison is shown in the Appendix. We also find that the DAOphot catalog is slightly more complete than the ngmix catalog in the inner 10' of Fornax, though both catalogs still suffer from crowding in this region. For example, we perform artificial star tests for the DAOphot pipeline on the central tiles of Fornax and assess the completeness level of our catalog. We find that within the central 10' a significant fraction of the input artificial stars are not recovered with the correct magnitude and location owing to blending, and that within 10' <50% completeness is reached at g < 22 if both criteria of matching spatially with 1" precision and magnitude within 0.5 mag are applied. The incompleteness level decreases for areas farther away from the galaxy center and appears to be symmetric on the east and west sides of Fornax. We therefore mask the inner 10' region when we derive structural parameters. A detailed analysis of the inner region of Fornax is reserved for future work.
<table><caption>**Table 1**

Fornax Structural Parameters Derived with the MCMC Method for Different Profile Models.</caption><thead><tr><th>Parameter</th><th>King</th><th>Exponential</th><th>Plummer</th><th>Sérsic</th><th>Bate15</th><th>B06</th></tr></thead><tbody><tr><td>α<sub>2000</sub></td><td>2<sup>h</sup>39<sup>m</sup>53<sup>s</sup></td><td>2<sup>h</sup>39<sup>m</sup>53<sup>s</sup></td><td>2<sup>h</sup>39<sup>m</sup>53<sup>s</sup></td><td>2<sup>h</sup>39<sup>m</sup>53<sup>s</sup></td><td>2<sup>h</sup>39<sup>m</sup>52<sup>s</sup></td><td>2<sup>h</sup>39<sup>m</sup>52<sup>s</sup></td></tr><tr><td>δ<sub>2000</sub></td><td>-34°30′32″</td><td>-34°30′21″</td><td>-34°30′25″</td><td>-34°30′32″</td><td>-34°30′39″</td><td>-34°30′49″</td></tr><tr><td>Ellipticity</td><td>0.31 ± 0.002</td><td>0.31 ± 0.002</td><td>0.29 ± 0.002</td><td>0.31 ± 0.002</td><td>0.31 ± 0.01</td><td>0.30 ± 0.01</td></tr><tr><td>Position angle</td><td>42°2 ± 0°2</td><td>43°2 ± 0°2</td><td>42°7 ± 0°3</td><td>42°3 ± 0°2</td><td>41°6 ± 0°2</td><td>46°8 ± 1°6</td></tr><tr><td>Sérsic r<sub>s</sub></td><td>...</td><td>...</td><td>...</td><td>16.4 ± 0.2</td><td>14.5 ± 0.1</td><td>17.3 ± 0.2</td></tr><tr><td>Sérsic index m</td><td>...</td><td>...</td><td>...</td><td>0.80 ± 0.006</td><td>0.78 ± 0.005</td><td>0.71 ± 0.01</td></tr><tr><td>Exponential r<sub>e</sub></td><td>...</td><td>11.8 ± 0.04</td><td>...</td><td>...</td><td>...</td><td>11.0 ± 0.1</td></tr><tr><td>Plummer r<sub>p</sub></td><td>...</td><td>...</td><td>19.9 ± 0.06</td><td>...</td><td>...</td><td>19.6 ± 0.1</td></tr><tr><td>King r<sub>t</sub></td><td>77.5 ± 0.4</td><td>...</td><td>...</td><td>...</td><td>69.7 ± 0.3</td><td>69.1 ± 0.4</td></tr><tr><td>King r<sub>c</sub></td><td>20.3 ± 0.1</td><td>...</td><td>...</td><td>...</td><td>14.6 ± 0.1</td><td>17.3 ± 0.2</td></tr><tr><td>2D Half-light r<sub>h</sub></td><td>20.8</td><td>19.8</td><td>19.4</td><td>20.1</td><td>17.0 (King)<br>17.2 (Sérsic)</td><td>18.1 (King)<br>18.6 (Sérsic)<br>18.5 (Exp.)</td></tr><tr><td>χ<sup>2</sup><sub>reduced</sub></td><td>1.89</td><td>2.04</td><td>2.54</td><td>1.41</td><td></td><td></td></tr></tbody></table>

**Note.**

The last two columns refers to the structural parameters derived by Bate et al. (2015) (Bate15) and Battaglia et al. (2006) (B06).

### 3. Structural Parameters

Several previous works have analyzed the radial distribution of stars in Fornax (Irwin & Hatzidimitriou 1995; Battaglia et al. 2006; Bate et al. 2015; del Pino et al. 2015). Profile models can be fitted to provide important constraints on the size and mass of the galaxy. We rederive the structural parameters of Fornax dSph using star objects from the uniform, wide-coverage DES Y3 data to compare with previous values. We use several different models, including an empirical King profile (King 1962), an exponential profile, a Sérsic profile (Sérsic 1968), and a Plummer profile (Plummer 1911). We fit a surface density profile to the data in a region of radius 2.5 around the center of the galaxy. Magnitude cuts at g, r < 23.5 are applied to ensure good-quality photometry, and stars with g - r > 1.2 are excluded. We define the ellipticity $\epsilon$ to be $\epsilon = 1 - b/a$, where b and a are the semiminor and semimajor axes of the galaxy, respectively. The position angle is taken from north toward east. We use the 2D unbinned maximum likelihood algorithm described in Martin et al. (2008). The likelihood is evaluated with an affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) of the form proposed by Goodman & Weare (2010), which is implemented as the MCMC Hammer emcee package (Foreman-Mackey et al. 2013) for Python. We allow the following parameters to vary freely: the center of Fornax (R.A. and decl.), the position angle, the ellipticity, the contamination stellar density, and the parameters for the radial density profile. We ran the MCMC with at least an ensemble of 50 walkers, each of which did approximate 1000 steps to ensure chain convergence. We note that the central stellar distribution, as discussed in the Appendix, has crowding issues. To reduce the potential effects from crowding on derived structural parameters, we exclude stars within the inner 10' radius at the galaxy center while we fit the observed data.

The structural parameters from fitting various profiles are shown in Table 1, where the uncertainties represent the marginalized errors of the 16th and 84th percentiles in the fitting procedure. The results from Battaglia et al. (2006) and Bate et al. (2015) are also shown here for comparison and are in general agreement with parameters derived here. However, we note a few specific differences when we compare with results

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Surface density profile for Fornax dSph galaxy with overlaid best-fitting exponential (blue dotted line), Plummer (purple dashed line), Sérsic (red dashed-dotted line), and King models (green solid line). The data are binned in elliptical annuli using the best-fit position angle and ellipticity $\epsilon = 0.31$, and the fits were performed to the full unbinned data allowing all parameters to vary freely, as described in Section 3. We note that our fits are derived excluding stars within a 10' radius (shown with the gray band). Parameter values of the fits are provided in Table 1.</figcaption></figure>

from those previous works. For example, our fits show systematically noticeably larger values in the effective radius, or the derived 2D half-light radius in all profile types. For example, the resulting Sérsic 2D half-light radius $r_h$ is higher by 1.5 (7.5%) than the derived Battaglia et al. (2006) value and 2.9 (14.4%) higher than the Bate et al. (2015) result. This is likely due to deeper and better-quality photometry in our data that have resolved a fainter and older main-sequence (MS) population that has more diffused stellar distribution. Our King tidal radius ($r_t = 77.5$) is also 7.8 (8.4) larger than the Bate et al. (2015) (Battaglia et al. 2006) values. Otherwise, our position angle, ellipticity, and other structural parameter values are consistent with differences between the literature values.

Figure 1 provides the best-fit models for each profile compared to the data binned in elliptical annuli with the best-fit position angle and ellipticity. These profiles were obtained using the parameter values from the unbinned fits in Table 1. We also compute the reduced chi-square $\chi^2_{reduced}$ as a goodness-of-fit indicator. It is calculated from the observed surface
<figure><img src="image_4.png" /><figcaption>**Figure 2.** CMDs in *g* and *r* bands for three regions around Fornax and one background region (rightmost panel). From left to right, we show all stars within an annulus with inner and outer elliptical radius *r*<sub>ell</sub> of *r*<sub>ell</sub> < 0.2, 0.2 < *r*<sub>ell</sub> < 0.8, 0.8 < *r*<sub>ell</sub> < 1.5, and the background region in a rectangle with -2°5 < *ξ* < -2°0 and -2°0 < *η* < 2°0, where *ξ* and *η* are coordinates centered on the Fornax centroid. Several different stellar populations are marked by red boxes: blue RGB (B-RGB), red RGB (R-RGB), red clump (RC), horizontal branch (HB), blue plume (BP), and blue loop (BL) stars. The black points with error bars show the average photometry errors at each magnitude bin.</figcaption></figure>

number density profile derived from elliptical, concentric
annuli ranging from 10' to 120' from the galaxy center and
the best-fitting profiles, by assuming that star counts follow
Poisson distributions with means equal to the star counts from
the best-fit model. Similar to Battaglia et al. (2006) and Bate
et al. (2015), we also found that a Sérsic profile provides a
much better fit around the tidal radius, whereas the King profile
generates slightly better fits at small radii. Overall a Sérsic
profile provides a best fit both at the tidal radius and in the inner
regions, as reflected by smaller $\chi^2_{\text{reduced}}$ compared to other
models. We note that from the $\chi^2_{\text{reduced}}$ values it is indicated that
most of the profile models we consider here are not very good
fits to the Fornax light profile. As we will show later in
Sections 5 and 6, near the galaxy center many substructures are
present, and the inner distribution is slightly off-centered for
some stellar populations, resulting in a complex profile shape
for Fornax that may not be captured by simple models.

4. The CMDs of Fornax

In Figure 2 we show *g*, *g* − *r* CMDs of several regions of Fornax and one background region. From left to right, we show all stars within an annulus with elliptical radius *r*<sub>ell</sub> (with ellipticity = 0.31) of *r*<sub>ell</sub> < 0.2, 0.2 < *r*<sub>ell</sub> < 0.8, and 0.8 < *r*<sub>ell</sub> < 1.5 and a contamination field (rightmost panel) consisting of data within −2.5° < *ξ* < −2.0° and −2.0° < *η* < 2.0°, where *ξ* and *η* are coordinates centered on the Fornax centroid. The CMD of Fornax shows a multitude of features, indicating that it has experienced a prolonged period of star formation. The well-sampled stellar features are clearly different from MW contamination features, as traced by the surrounding background region. The bright part of the CMD is dominated by the RGB feature, which shows a clear bifurcation (also see discussions in Bate et al. 2015; Battaglia et al. 2006). Following Battaglia et al. (2006) and Bate et al. (2015), we will refer to the two branches of RGBs as the blue RGB (B-RGB) and red RGB (R-RGB). The two features probe stars of different metallicities and ages and show different spatial distribution, pointing to a radial gradient of metallicity within Fornax (see also Battaglia et al. 2006; de Boer et al. 2012). Blue loop (BL) stars are also visible around *g* − *r* ≈ 0.5, indicating the presence of helium-core burning stars with ages ranging from a few million years to 1 Gyr. Lower down in the CMD, a strongly populated horizontal branch (HB) and red

clump (RC) are visible, both of which change distribution at
different radii. Ancient (>10 Gyr), metal-poor blue HB stars
are strongly present in the 0°.8 < *r*<sub>ell</sub> < 1°.5 bin, while the RC
stars dominate in the innermost bin. This indicates that the
radial age gradient is also linked to a metallicity gradient, as
also seen from previous spectroscopic measurements at
different radii (Battaglia et al. 2006). The rich, composite red
clump is composed of intermediate-age, metal-enriched stars
that have been used as a standard candle to determine the
distance to Fornax (Bersier 2000; Rizzi et al. 2007).

The faint part of the CMD is dominated by a composite
subgiant branch and MSTO feature, as a result of multiple
overlapping stellar populations covering a large range in age.
This is further supplemented by a strongly populated blue
plume (BP) population extending to blue color and bright
magnitudes. The CMD of BP is largely mixed with the
population of blue straggler stars, which can also be found on
the blueward side of old MSTOs, except at the very bright
magnitude end, where it is composed of genuine young BP
stars. Those bright BP stars, which are only reachable by young
(<1 Gyr) MS stars, are mostly visible in the innermost Fornax
sample. The composition of the CMD features also changes at
different radii, with younger populations gradually disappear-
ing as the radius increases until only the oldest populations are
left over in the outskirts. This once again points to younger
stars being preferentially found in the center of Fornax.

To provide nonparametric, model-independent evidence of
the spatial distribution difference among different subpopula-
tions, we calculate the elliptical distance of the stars and
compare their cumulative distributions using the Kolmogorov–
Smirnov (K-S) test. We find that they all return *p*-values <0.01
no matter which comparison combination we try, indicating
that they all have statistically significantly different spatial
distribution. The composition of features in the CMD makes it
clear that Fornax has formed stars across most of cosmic time
without significant interruptions and just recently (≈100
million years ago) stopped star formation (de Boer et al. 2013).

5. Surface Density Maps

The Fornax dwarf galaxy is among the systems most likely
to be affected by MW tidal forces, given its relatively cold
kinematics for its brightness (Wang et al. 2017). If present,
these tidal features are likely to have a very low surface
density, making them hard to detect against a screen of MW contamination. Therefore, efficient decontamination techniques, such as “matched-filtering” methods (e.g., Rockosi et al. 2002; Grillmair 2009), are needed to enhance their signal over the numerous contaminants.

To investigate the presence of tidal features in Fornax, we construct matched-filtered maps of different stellar populations. To ensure that the filtered data are still Poissonian distributed, we forgo using a classical matched-filter technique and instead adopt a Boolean mask procedure. In our case, the Hess diagram of the galaxy dense regions, where the ratio of contamination/ source densities is low, is used to build a “source” filter, defining the shape of the dSph stellar populations in the color–magnitude plane. The source filter region is derived within an ellipse with a size of 0.8 times the tidal radius. A “contamination filter” is obtained from a region far enough away to be free of galaxy members, i.e., the background region with $-2.5 < \xi < -2.0$ and $-2.0 < \eta < 2.0$ that is shown in the rightmost panel of Figure 2.

We describe the source filter, in the form of a probability density function $P_{str}(g, g-r)$ in color–magnitude space. The contamination is described by $P_{bg}(g, g-r)$. Both the source and background filter masks are then normalized, and the fraction $P_{str}(g, g-r)/P_{bg}(g, g-r)$ is used to select CMD pixels above a certain threshold. The filter threshold is chosen by finding the value that maximizes the signal-to-noise ratio. We then assign a weight of 0 to values of $P_{str}(g, g-r)/P_{bg}(g, g-r)$ lower than the threshold and a weight of 1 for values higher than the threshold.

The decontaminated Fornax stellar spatial distribution derived from the matched-filtering method is displayed in the top left panel of Figure 3. In the top right panel of Figure 3 it is shown in the form of isodensity contours with pixel size $3' \times 3'$. There are no significant signs of tidal features or distortion within the 25 deg$^2$ area that we have examined. However, there are mildly extended low surface density excesses in the distribution of stars on the east (left-hand side of the plot) and the west side of the galaxy, indicating a slightly skewed isophote shape in the outskirts.

We also derive a 2D map of residuals between the decontaminated surface density map of Fornax stars and the surface density predicted by the best-fitting Sérsic model. It is shown in the bottom two panels of Figure 3. Except for some over- and underdensities in the Fornax inner region that we will discuss later in Section 7, it shows reasonable agreement between the decontaminated map of Fornax stars and the models. Based on those maps, we analyze how frequent those residuals appear in the search area and adopt the $1\sigma$ confidence interval of the distribution as our detection limit; this corresponds to ~1.8 stars pixel$^{-1}$ (~32.1 mag arcsec$^{-2}$). We note that the inner $10'$ region, which is shown as a gray ellipse in the bottom right panel of Figure 3, is excluded in this analysis owing to effects from crowding. The lack of significant residuals in the periphery of the Fornax dSph galaxy in the residual map therefore indicates no statistically significant tidal debris detection down to a surface brightness limit of ~32.1 mag arcsec$^{-2}$.

# 6. Spatial Distribution of Different Stellar Populations

The analysis of the spatial distribution of stars in different evolutionary phases is a useful tool to study spatial variations of the stellar population mix as a function of age and/or

metallicity. There are multiple examples in the literature of this type of analysis, through which age gradients and star formation history in the Fornax dwarf galaxy have been quantified (Battaglia et al. 2006; Bate et al. 2015; del Pino et al. 2015).

In this section, we perform a full matched-filtering method (described in Section 5) and structural analysis of the spatial distribution of Fornax R-RGB, B-RGB, HB, and BP stars. We have chosen not to include an analysis of the RC feature, since the stellar population mix of the RC feature is similar to that of the R-RGB feature. Filters were generated from $(g-r, g)$ Hess diagrams, as displayed in Figure 2. The filters were built by truncating the Hess diagram at appropriate CMD boxes.

## 6.1. Matched-filtered Maps

The matched-filtered maps for four subpopulations are shown in Figure 4: B-RGB, R-RGB, HB, and BP. Each panel is normalized separately, based on their respective rms values of star counts in the contamination area with $-2.5 < \xi < -2.0$ and $-2.0 < \eta < 2.0$, where $\xi$ and $\eta$ are coordinates centered on the Fornax centroid. The rms values for each subpopulation are 0.14 star pixel$^{-1}$ for B-RGB, 0.14 star pixel$^{-1}$ for R-RGB, 0.20 star pixel$^{-1}$ for HB, and 0.08 star pixel$^{-1}$ for BP. Contour levels are shown at 10, 25, 50, 100, 200, 400, 600, 1000, and 1500 times the rms values calculated in the contamination region. The tidal radius (blue ellipses with solid lines) and the half-light radius (blue ellipses with dashed lines) of the overall Fornax distribution derived by this work are plotted in each panel. The red ellipses mark the half-light radius of each subpopulation. The center of the overall Fornax distribution is marked by the blue cross.

The B-RGB is more diffuse than the R-RGB, as pointed out in Bate et al. (2015). B-RGB, B-RGB, and RC distributions appear to align with the overall Fornax orientation but show slightly off-center positions and are slightly more extended toward the southwest side (lower right corners in Figure 4). The HB is more diffuse than other subpopulations and shows a mildly extended distribution toward the southwest direction as well. For BP, it is clear that the orientation and the ellipticity of the isodensity contours change to be more aligned with the R. A. direction and become more circularized toward the center. This is also shown in del Pino et al. (2015), where it is argued that the young MS stars do not follow the overall Fornax star distribution in terms of both orientation and ellipticity.

## 6.2. Structural Parameters

In order to statistically analyze the compatibility between the spatial distributions of stars in the different evolutionary phases, we obtained their structural parameters by fitting with Plummer models. For the BP population, we model it with two Plummer models superposed on each other to capture the change of spatial orientation. Since our aim here is not to determine which functional form best represents the surface relative density of the various populations, but to quantify relative differences in their spatial distribution, we fit only a Plummer profile, in order to restrict the number of free parameters. Even in the regime of low number statistics for some of the Fornax subpopulation stars (e.g., B-RGB), the affine-invariant ensemble sampler method provides relatively well constrained position angle, ellipticity, and half-light radius estimates. We note that unlike the overall structural parameter
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top left panel: spatial distribution of Fornax dSph matched-filtered stars. Top right panel: surface number density map of the matched-filtered stars overlaid with isodensity contours. The light-blue cross marks the galaxy centroid. Bottom left panel: map of residuals between the Fornax surface number density profile from the matched-filtered stars and the surface density from the best-fitting Sérsic model (see Table 1 for the parameter values). Pixel size is 2' × 2', smoothed by a Gaussian kernel with 4' dispersion. Bottom right panel: residual map, with blue (orange) and dark-blue (dark-orange) contours showing 2σ and 3σ detections above (below) the mean residuals. The red ellipses in each panel show the nominal King tidal radius derived from this work. The gray ellipse at the center with a radius of 10' marks the region that is excluded in the residual analysis.</figcaption></figure>

fits in Section 3, we do not mask the inner 10' for the subpopulation analysis here. Most of the stellar subpopulations we consider here are brighter than $g < 22$ and therefore are much less subject to the incompleteness due to crowding. The faintest subpopulation, BP, actually has the most compact distribution (see the discussion below). This is opposite of what we would expect from a completeness artifact. Still, this may impose greater uncertainties on the BP structural parameters than the values quoted here, and readers should proceed with caution with those results.

In Table 2 we summarize the Plummer model fitting results of these various subpopulations. The results listed in Table 2 reflect what are shown in the matched-filtered contour maps in Figure 4 and provide quantitative descriptions. The red ellipses in each panel in Figure 4 also show the orientation, ellipticity, and half-light radius from the Plummer model fits. The

centroids of all four subpopulations show very mild deviations from the Fornax center. It is clear in Table 2 that the HB has a larger Plummer radius $r_p$ ($r_p = 29.3$) than other subpopulations and the overall population ($r_p = 19.4$). The central compact BP distribution, which roughly consists of ~16% of the whole BP population, is the most compact (with the smallest $r_p$ of 12.0) and has a very different orientation ($\theta = 105.8$) than other populations and its own diffuse components. The diffuse BP component has orientation very similar to the overall distribution but with slightly rougher shape (ellipticity = 0.25). The R-RGB has profile parameters most similar to the overall Fornax population (see Table 1). For B-RGB and RC, they are also similar to the overall Fornax population, except for a slightly different orientation or ellipticity.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Matched-filtered stellar density contours for the blue RGB (B-RGB; top left panel), red RGB (R-RGB; top middle panel), red clump (RC; top right panel), horizontal branch (HB; bottom left panel), and blue plume (BP; bottom middle panel) stars. Contour levels are shown at 10, 25, 50, 100, 200, 400, 600, 1000, and 1500 times the rms values calculated in the contamination region. The blue ellipses show the nominal King tidal radius, and the small blue ellipses with dashed lines mark the 2D half-light radius of the overall Fornax population. The red ellipses display the half-light radius from the Plummer profile fits (see Table 2) of individual subpopulations. The blue crosses mark the galaxy centroid.</figcaption></figure>

## 7. Inner Overdensities

In Figure 5 we show a residual map by subtracting a pixelized number density map with the smooth best-fit Sérsic model zoomed in within the 4 deg² area around the Fornax center. We mark several previously known and newly discovered features by this work. Assuming Poisson statistics, the blue (orange) and dark-blue (dark-orange) contours, respectively, show the 3σ and 5σ significance regions above (below) the smooth Sérsic model. We note that the additional 5% uncertainty in smooth background counts is added to take into account the propagated uncertainties from the best-fit model parameters. The light-gray and brown contours show the 2σ significance regions (above and below the smooth background). The light-green diamonds show locations of the known globular clusters, cyan pentagons indicate the position of previously identified overdensities from the literature (Coleman et al. 2004; de Boer et al. 2013; Bate et al. 2015), and yellow triangles indicate the position of two overdensities with significance ≳3 in the young star distribution discussed in Section 7.1.

We also mark three regions of residual excesses with red lines on the southwest side and northeast side that resemble shell-like features wrapping around the galaxy. These low surface brightness features are ~20′–40′ away from the galaxy center, with a few regions along those features having ≳3σ

above the smooth background. The nature and true significance of these features are unclear, as we note that the significance of those residuals may be further reduced if a more complex density profile model is considered, such as a superposition of two Sérsic profiles with one of them slightly off-center. Still, there are hints in the literature that suggest that Fornax has experienced one or more merger events in the past. We speculate that these overdensities shown in Figure 5 might be the debris of those tidally disrupted stellar systems.

For example, this scenario is supported by chemo-kinematic analyses of Fornax, which reveal that Fornax may have complex dynamical patterns, including multirotation components with different metallicity distributions or directions of angular momentum (Amorisco & Evans 2012; del Pino et al. 2017). This is ascribed to a merger between two or more stellar systems that each carry a particular angular momentum when they merged with Fornax. Furthermore, several stellar substructures have been linked to possible merger origins (Coleman et al. 2004; Yozin & Bekki 2012; del Pino et al. 2015). If some of those globular clusters in Fornax are brought in by merger events in the past, this might explain Fornax’s unusually high number of globular clusters.
<table><caption>**Table 2**
Fornax Subpopulation Structural Parameters (Plummer Model)</caption><thead><tr><th>Parameter</th><th>B-RGB</th><th>R-RGB</th><th>RC</th><th>HB</th><th>BP</th></tr></thead><tbody><tr><td>Δα<sub>2000</sub></td><td>0°02</td><td>0°01</td><td>0°005</td><td>0°014</td><td>-0°002</td></tr><tr><td>Δδ<sub>2000</sub></td><td>-0°018</td><td>0°000</td><td>0°005</td><td>0°015</td><td>-0°005</td></tr><tr><td>Ellipticity</td><td>0.31 ± 0.01</td><td>0.31 ± 0.01</td><td>0.29 ± 0.01</td><td>0.32 ± 0.01</td><td>0.39 ± 0.1</td></tr><tr><td>θ</td><td>38°9 ± 1°4</td><td>42°1 ± 0°8</td><td>43°2 ± 0°4</td><td>37°8 ± 1°1</td><td>105°8 ± 7°0</td></tr><tr><td>r<sub>p</sub></td><td>22°6 ± 0°4</td><td>18°1 ± 0°2</td><td>18°0 ± 0°1</td><td>29°3 ± 0°3</td><td>12°0 ± 0°7</td></tr><tr><td>No. of stars</td><td>3268</td><td>9117</td><td>22751</td><td>4889</td><td>fraction = 0.16 ± 0.05</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Residual significance map within the 4 deg² area of Fornax. The blue (orange) and dark-blue (dark-orange) contours show the 3σ and 5σ significance regions above (below) the smooth Sérsic model. The light-gray and brown contours show the 2σ significance regions (above and below the smooth background). The black solid ellipses show the tidal radius, and the dashed ellipses are with radii of 0.7 and 0.5 times the tidal radius. The light-green diamonds show locations of the five known globular clusters, the cyan pentagons indicate the position of previously identified overdensities from the literature (Coleman et al. 2004; de Boer et al. 2013; Bate et al. 2015), and the yellow triangles indicate the position of two additional overdensities with significance >3σ from the young star distribution discussed in Section 7. The pink triangle marks the location of the shell discussed by Coleman et al. (2004). The red dashed lines on the southwest side and northeast side mark three regions that resemble possible shell-like features >3σ. North is up, and east is to the left.</figcaption></figure>

## 7.1. Young Stellar Overdensities

The inner regions of Fornax (i.e., $r_{ell} < 0°.5$) are dominated by stars of intermediate age and show a large fraction of stars younger than ≈2 Gyr on the young MS. These BP stars are readily visible in Figure 2 for blue colors with $g - r < 0$ extending up to $g ≈ 20$. Previous studies of the young stars have discovered several overdensities, some of which are composed of stars younger than what are found in the Fornax field (Coleman et al. 2004; de Boer et al. 2013; Bate et al. 2015).

To study the intermediate-age and young population in DES Y3 data, we select stars on the BP region with a color cut shown in Figure 2. We exclude stars fainter than $g, r = 23.5$ to avoid stars on the old MS that may scatter into our selection box. In Figure 6, we show the spatial density of BP stars, with colors representing the density in each bin. The first thing to notice in Figure 6 is that the young stars are oriented quite differently than the old stellar populations on the RGB, which

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Spatial density distribution of Fornax BP stars, selected from the $g$ and $r$ bands using cuts of $g - r < 0$ and $r < 23.5$. Black pentagons indicate the position of previously identified overdensities from the literature (Coleman et al. 2004; de Boer et al. 2013; Bate et al. 2015). Large dashed circles in the central region indicate the centers of the roughly bimodal distribution making up the core of the young Fornax population. Finally, solid numbered circles indicate the position of five overdensities or “regions of interests” detached from the central young star distribution. North is up, and east is to the left.</figcaption></figure>

is also shown in Figure 4 and the position angle from the Plummer model fits (see Table 2). The central distribution within 0°.2 is aligned very close to the R.A. direction, as pointed out in Section 6. However, there is a diffuse and extended component that is more aligned with overall orientation and ellipticity of Fornax, and therefore the derived position angle for BP stars in Table 2 only differs by 7°.8 compared to the overall Fornax position angle. The origin of this orientation mismatch is unclear but may be linked to the episodic formation history of Fornax and the possibility of gas accretion fueling the creation of the young stars (e.g., Yuan et al. 2016). The distribution of young stars can be roughly characterized as two large overdensities on either side of the north-south axis, with the eastern one being slightly more extended.

To further investigate the central populations, we extract the CMD of the central features within the large dashed circles shown in Figure 6. These two circles are positioned on the peaks of the two main density distributions, after doing a simple Gaussian fit. For comparison, we show the expected distribution of young, intermediate-age, and old stars by overlaying stellar isochrones from the MIST isochrone library (Choi et al. 2016; Dotter 2016). The CMDs show the presence
<figure><img src="image_9.png" /><figcaption>**Figure 7.** CMDs of the two main overdense regions within the young Fornax populations. These two overdensities show the presence of a wide range of stellar populations, ranging from ancient and metal-poor to intermediate-age (≈5 Gyr) stars and young (<1 Gyr) metal-rich stars. The blue ([Fe/H] = −2.50 dex, age = 13 Gyr) and red ([Fe/H] = −1.00 dex, age = 4 Gyr) lines show reference isochrones for the main Fornax field population, while the green line indicates a stellar population with solar metallicity at 0.5 Gyr. Comparison between the two panels shows that the eastern region has younger stars (with $g-r \lesssim 0.0$) extending all the way up to $r = 19.5$, while the young stars in the western region only extend up to $r = 20.5$ and are therefore slightly older.</figcaption></figure>

of a wide range of stellar populations, including ancient metal-poor stars on the blue RGB and the dominant intermediate-age (≈5 Gyr), metal-rich population that results in a red RGB and dense red clump. Furthermore, as expected, the CMDs show a clear presence of young stars on the BP population. The density and morphology of the young BPs are different between the two regions, with stars in the eastern region (top panel) extending to brighter magnitudes. Comparison to stellar isochrones, which is shown in Figure 7, indicates that the western region is well matched to a population with solar metallicity at an age of 0.5 Gyr, while the eastern region has brighter stars corresponding to a younger MS age.

Figure 6 shows several other overdense features of young stars surrounding the central features. The positions of known overdensities are labeled by pentagons and clearly correlate with regions with a density higher than the local background. The known overdensities are described in more detail in Coleman et al. (2004) (#1; bottom left pentagon), de Boer et al. (2013) (#2; top left pentagon), and Bate et al. (2015) (#3 and #4; top right and bottom right pentagons, respectively). We extract stars of the known features using regions with a radius of 2' and show the resulting CMDs in Figure 8. Compared to the VST/ATLAS data from Bate et al. (2015), the

DES Y3 photometry is deeper and allows us to characterize the stellar populations of these features with greater accuracy. Therefore, we determine the best fit to the overdense population using a modified version of the Talos code (de Boer et al. 2012), which fits a single isochrone to a cut-out of the CMD centered on the blue part of the CMD. These fits are made to generate an indication of the age of the overdensities, and no attempt is made here to fit the detailed shape of the density distribution. Metallicities are assumed to follow the Fornax young age–metallicity relation as determined in de Boer et al. (2013) for the field population. Given their brightness, no spectroscopic metallicities are known for the young stars of Fornax. We find that the overdensity discovered by Coleman et al. (2004) (#1) is the oldest, with an age of 1.5 Gyr, in good agreement with the literature. The two features discovered by Bate et al. (2015) (#3 and #4) have ages of 0.5 and 1 Gyr, respectively, while the overdensity found previously by de Boer et al. (2013) (#2) has a very young age of 0.2 Gyr. It is interesting that the youngest overdensities are found on the eastern side of Fornax, for which the field population is also younger than the western side. This might be an indication that this side of the Fornax stellar distribution might be younger as a whole, possibly as the result of an external event.

The distribution of young stars in Figure 6 shows several other high-density regions away from the main central core. These might simply be extensions of the main Fornax field distribution but could also be separate structures with a distinct population. However, due to the complex distribution of young stars, we do not attempt to identify individual substructures, which are very likely noise peaks. Instead, we select five regions of interest highlighted in Figure 6 and study their CMDs for acquiring additional information. Those regions are placed on local overdense peaks with $N > 10$ per pixel that are sufficiently far from previously known high overdensities. The location is defined using a simple Gaussian fit to the region surrounding the hot pixel. We select all stars around high-density features within a radius of 3' as shown in Figure 6. For region #6 we adopt a smaller radius of 2' to avoid overlapping with the main Fornax field population. By examining their CMDs, some regions are mostly populated by blue horizontal branch stars (e.g., #5), and some are likely extensions of the young Fornax field population (e.g., #7) and therefore are likely not genuine. Several regions exhibit similar age and metallicity fitting results and also show spatial proximity, indicating hints of possible links between them, for example, #6 and the shell of Coleman et al. (2004) (#1), region #4 from Bate et al. (2015) and our region #9, and #8 and the nearby overdensity of de Boer et al. (2013) (#2).

## 8. Conclusions and Discussions

In this work we utilize the DES Y3 data that have homogeneous wide-field coverage of the Fornax dSph galaxy down to a depth of $g, r \sim 23.5$ to study several properties of the galaxy.

We have fitted the Fornax structural parameters with several different models, including Plummer, Sérsic, exponential, and King profiles, and found that the Sérsic model gives the best fit to the data, while the other models seem to fail to capture the complex light profile shape of Fornax. We also found that our data indicate that Fornax is more extended than what was shown in previous studies. For example, our results show slightly larger half-light radius than the literature values, which
<figure><img src="image_10.png" /><figcaption>**Figure 8.** CMDs of previously identified overdensities in Fornax from the literature (#1 is from Coleman et al. 2004, #2 from de Boer et al. 2013, and #3 and #4 are from Bate et al. 2015). Overlaid on the CMDs are reference isochrones for the Fornax field population (blue: [Fe/H] = -2.50 dex, age = 13 Gyr; red: [Fe/H] = -1.00 dex, age = 4 Gyr), as well as the locus of the main overdense populations (in green with metallicity and age marked in each panel). We note that the magnitude distribution of the young stars (with $g-r \lesssim 0.0$) reflects on their inferred age.</figcaption></figure>

were derived from either shallower data or incomplete coverage. Our derived King tidal radius $r_t$, which has a value of 77.5, is also larger than previous results ($r_t$ = 69.7 from Bate et al. 2015 and $r_t$ = 69.1 from Battaglia et al. 2006). Nevertheless, in general the rest of the properties, such as ellipticity, position angle, and profile characteristic radii, agree with previous results.

To examine the possible low surface density tidal feature around Fornax, we also applied the matched-filtering method to enhance the signal over contamination. We examined a 25 deg² area around Fornax and found no significant tidal features or distortion down to ~32.1 mag arcsec⁻² with a 1σ confidence interval by analyzing the frequency distribution of the smooth model subtraction residuals. This could imply that either the tidal debris are too diffuse to be reached with current detection limits or the results align with predictions of some previous N-body simulation works that suggest a lack of tidal distortions (Barber et al. 2015; Battaglia et al. 2015).

We further examine the spatial distribution of several different stellar populations, including R-RGB, B-RGB, HB, and BP stars, by checking their matched-filtered maps and fitting with the Plummer profile model. We find that the HB stars are more spatially extended than the RGB stars, in accordance with previous findings. Our matched-filtered maps of R-RGB, B-RGB, and BP stars show mild departure from symmetric distribution and have distorted or extended features toward the southwest side. The orientation of BP stars in the center region deviates significantly from the overall Fornax orientation and gradually changes farther out from the center.

When we subtract a smooth Sérsic model from the observed stellar distribution and derive the significance of the residuals above the smooth background, we find roughly three low surface density features at ~3σ above the background. These shell-like residual features are located ~20'-40' away from the galaxy center and wrap around Fornax at the northeast side and southwest side. However, we note that the significance of those residuals may be further reduced if a more complex density profile model is considered, such as a superposition of two Sérsic profiles with one of them slightly off-center. Although the low significance of the detected signal, their proximity to the region masked in our analysis, and their suspicious symmetry avoid any concrete conclusion of their nature, those features might be debris due to past merger events.

We have studied the young populations in Fornax by focusing on the distribution and morphology of blue stars ($g-r < 0$) on the young MS. Stars that occupy this region of the CMD have ages typically below 2 Gyr and can be used to study the spatial distribution of Fornax in a regime where previous studies have seen hints of interactions or accretion (Coleman et al. 2004, 2005). The depth of DES Y3 data allows us to study the stellar population content in young stars with greater accuracy than was previously possible using other wide-field surveys such as VST/ATLAS.

The spatial density of young stars (see Figure 6) shows a very different morphology than the distribution of old stars in Fornax. While the dominant 4–7 Gyr old stars have a position angle of nearly 45°, the inner young stars are in bimodal distribution that is more aligned along the R.A. direction. This
indicates that the change in spatial distribution occurred in a time frame of at most several gigayears. This is a promising opportunity for deciphering the formation history of Fornax, especially given that the Gaia satellite (Gaia Collaboration et al. 2016) will produce accurate proper motions for the brightest stars in this sample.

The young star distribution in Figure 6 shows several high-density regions, some of which are new and some of which were previously discovered (Coleman et al. 2004; de Boer et al. 2013; Bate et al. 2015). We use the deep DES data to assign an age to the overdense population based on their CMD morphology (see Figure 8). We find a range of ages between 0.2 and 1.5 Gyr, in good agreement with literature values. Given the discrete sets of populations we find, the star formation happened in a stochastic way with inhomogeneous mixing across the galaxy.

It is striking that the youngest overdensities are all found on the eastern side of Fornax, where the field population itself is slightly younger than in the west (see Figure 7). Furthermore, we find tantalizing evidence that the previously discovered overdensities are linked to other overdense features farther in toward the bulk of the Fornax young field population. This gives the impression that these features are in fact shells or clumps of stars stripped out from the central young distribution. It is unclear what implications this has for the recent evolution of Fornax, especially in the minor merger scenario (del Pino et al. 2015). The morphology and ages of the features imply that Fornax may have undergone an encounter of some sort during the past 2 Gyr that radically changed the spatial distribution of newly formed stars but left the distribution of already-formed stars intact. An accurate kinematic study of these relatively faint stars (down to $r \approx 22.5$) would be invaluable in this context.

We would like to thank Vasily Belokurov for helpful discussions. We also would like to thank the anonymous referee for insightful and constructive comments that greatly improved this manuscript. M.Y.W. acknowledges the support of the McWilliams Postdoctoral Fellowship. T.d.B. acknowledges support from the European Research Council (ERC StG-335936). S.K. is supported by National Science Foundation grant AST-1813881. This paper has gone through internal review by the DES Collaboration.

This research has made use of NASA's Astrophysics Data System Bibliographic Services.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at The Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungs-gemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at the Cerro Tololo Inter-American Observatory and the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grant Nos. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDE funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013), including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project No. CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC, under contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States government retains and the publisher, by accepting the article for publication, acknowledges that the United States government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States government purposes.

Facility: Blanco (DECam).

Software: SExtractor (Bertin et al. 2002), DAOPhot (Stetson 1987), astropy (Astropy Collaboration et al. 2013), matplotlib (Hunter 2007), numpy (Van Der Walt et al. 2011), scipy, ngmix (Sheldon 2014), emcee (Foreman-Mackey et al. 2013).

## Appendix

### Surface Density Profile from DAOphot Data

To evaluate the impact of crowding on the Fornax DES Y3 data and how much the derived structural parameters could
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left panel: surface density profile for Fornax dSph galaxy with stars of *g* < 23.5 overlaid with best-fitting Sérsic models from DES ngmix photometry data (model: black dashed-dotted line; data: black points) and DAOphot data (model: red dashed-dotted line; data: red points). The data are binned in elliptical annuli using the best-fit position angle and ellipticity, and the fitting procedure is described in Section 3. Right panel: surface density profile from DES ngmix and DAOphot data with stars of *g* < 22.0.</figcaption></figure>

potentially be affected, we generate *DAOphot* (see Section 2) catalogs and compare the resulting surface density profiles between these two data sets. The binned 1D density profiles from DES Y3 data (black points) and the *DAOphot* data (red points) are displayed in Figure 9. In general, the DES Y3 profile agrees well with the *DAOphot* profile. However, it is noticeable that at the center region (*r* < 10′) the profile amplitude from the DES Y3 ngmix catalog is slightly lower than the *DAOphot* result. The density discrepancy at the center is more prominent for faint magnitude limits. This is shown in Figure 9 when we compare profiles in the left panel (stars with *g* < 23.5) to the ones in the right panel (stars with *g* < 22.0). This indicates that the DES Y3 data do present crowding issues at the center of Fornax. Also as mentioned in Section 2, we perform artificial start tests and find that the *DAOphot* catalog is <50% complete within the inner 10′ for stars with *g* < 22. In order to reduce the impacts of crowding on the derived structural parameters, we therefore mask the central 10′ region when we perform the fits.

## ORCID iDs

M. Y. Wang https://orcid.org/0000-0002-8226-6237

A. Drlica-Wagner https://orcid.org/0000-0001-8251-933X

S. E. Koposov https://orcid.org/0000-0003-2644-135X

A. K. Vivas https://orcid.org/0000-0003-4341-6172

A. B. Pace https://orcid.org/0000-0002-6021-8760

A. R. Walker https://orcid.org/0000-0002-7123-8943

L. Strigari https://orcid.org/0000-0001-5672-6079

J. L. Marshall https://orcid.org/0000-0003-0710-9474

B. Yanny https://orcid.org/0000-0002-9541-2678

A. Roodman https://orcid.org/0000-0001-5326-3486

T. M. C. Abbott https://orcid.org/0000-0003-1587-3931

J. Annis https://orcid.org/0000-0002-0609-3987

D. Brooks https://orcid.org/0000-0002-8458-5047

D. L. Burke https://orcid.org/0000-0003-1866-1950

M. Carrasco Kind https://orcid.org/0000-0002-4802-3194

C. B. D'Andrea https://orcid.org/0000-0002-8198-0332

J. García-Bellido https://orcid.org/0000-0002-9370-8360

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

D. Gruen https://orcid.org/0000-0003-3270-7644

R. A. Gruendl https://orcid.org/0000-0002-4588-6517

D. L. Hollowood https://orcid.org/0000-0002-9369-4157

D. J. James https://orcid.org/0000-0001-5160-4486

K. Kuehn https://orcid.org/0000-0003-0120-0808

R. Miquel https://orcid.org/0000-0002-6610-4836

E. Sanchez https://orcid.org/0000-0002-9646-8198

M. Smith https://orcid.org/0000-0002-3321-1432

F. Sobreira https://orcid.org/0000-0002-7822-0658

M. E. C. Swanson https://orcid.org/0000-0002-1488-8552

G. Tarle https://orcid.org/0000-0003-1704-0781

## References

Abbott, T. M. C., Abdalla, F. B., Allam, S., et al. 2018, ApJS, 239, 18

Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4

Amorisco, N. C., & Evans, N. W. 2012, ApJL, 756, L2

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Barber, C., Starkenburg, E., Navarro, J. F., & McConnachie, A. W. 2015, MNRAS, 447, 1112

Bate, N. F., McMonigal, B., Lewis, G. F., et al. 2015, MNRAS, 453, 690

Battaglia, G., Sollima, A., & Nipoti, C. 2015, MNRAS, 454, 2401

Battaglia, G., Tolstoy, E., Helmi, A., et al. 2006, A&A, 459, 423

Bersier, D. 2000, ApJL, 543, L23

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in ASP Conf. Ser. 281, Astronomical Data Analysis Software and Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley (San Francisco, CA: ASP), 228

Buonanno, R., Corsi, C. E., Zinn, R., et al. 1998, ApJL, 501, L33

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Coleman, M., Da Costa, G. S., Bland-Hawthorn, J., et al. 2004, AJ, 127, 832

Coleman, M. G., Da Costa, G. S., Bland-Hawthorn, J., & Freeman, K. C. 2005, AJ, 129, 1443

Coleman, M. G., & de Jong, J. T. A. 2008, ApJ, 685, 933

de Boer, T. J. L., Tolstoy, E., Hill, V., et al. 2012, A&A, 544, A73

de Boer, T. J. L., Tolstoy, E., Saha, A., & Olszewski, E. W. 2013, A&A, 551, A103

de Vaucouleurs, G., & Ables, H. D. 1968, ApJ, 151, 105

del Pino, A., Aparicio, A., & Hidalgo, S. L. 2015, MNRAS, 454, 3996

del Pino, A., Aparicio, A., Hidalgo, S. L., & Łokas, E. L. 2017, MNRAS, 465, 3708

del Pino, A., Hidalgo, S. L., Aparicio, A., et al. 2013, MNRAS, 433, 1505

Dotter, A. 2016, ApJS, 222, 8

Drlica-Wagner, A., Sevilla-Noarbe, I., Rykoff, E. S., et al. 2018, ApJS, 235, 33

Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
Gallart, C., Aparicio, A., Zinn, R., et al. 2005, Near-fields Cosmology with Dwarf Elliptical Galaxies, in IAU Coll. 198, ed. H. Jerjen & B. Binggeli (Cambridge: Cambridge Univ. Press), 25

Goodman, J., & Weare, J. 2010, CAMCS, 5, 65

Grillmair, C. J. 2009, ApJ, 693, 1118

Hodge, P. W. 1961a, AJ, 66, 83

Hodge, P. W. 1961b, AJ, 66, 249

Hodge, P. W., & Smith, D. W. 1974, ApJ, 188, 19

Hunter, J. D. 2007, CSE, 9, 90

Irwin, M., & Hatzidimitriou, D. 1995, MNRAS, 277, 1354

Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201

King, I. 1962, AJ, 67, 471

Larsen, S. S., Strader, J., & Brodie, J. P. 2012, A&A, 544, L14

Letarte, B., Hill, V., Jablonka, P., et al. 2006, A&A, 453, 547

Letarte, B., Hill, V., Tolstoy, E., et al. 2010, A&A, 523, A17

Łokas, E. L. 2009, MNRAS, 394, L102

Martin, N. F., de Jong, J. T. A., & Rix, H.-W. 2008, ApJ, 684, 1075

Mateo, M., Olszewski, E., Welch, D. L., Fischer, P., & Kunkel, W. 1991, AJ, 102, 914

McConnachie, A. W. 2012, AJ, 144, 4

Morganson, E., Gruendl, R. A., Menanteau, F., et al. 2018, PASP, 130, 074501

Nichols, M., Lin, D., & Bland-Hawthorn, J. 2012, ApJ, 748, 149

Olszewski, E. W., Mateo, M., Harris, J., et al. 2006, AJ, 131, 912

Pasetto, S., Grebel, E. K., Berczik, P., Chiosi, C., & Spurzem, R. 2011, A&A, 525, A99

Pietrzyński, G., Górski, M., Gieren, W., et al. 2009, AJ, 138, 459

Plummer, H. C. 1911, MNRAS, 71, 460

Pont, F., Zinn, R., Gallart, C., Hardy, E., & Winnick, R. 2004, AJ, 127, 840

Rizzi, L., Held, E. V., Saviane, I., Tully, R. B., & Gullieuszik, M. 2007, MNRAS, 380, 1255

Rockosi, C. M., Odenkirchen, M., Grebel, E. K., et al. 2002, AJ, 124, 349

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Sérsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico)

Shapley, H. 1938, Natur, 142, 715

Sheldon, E. S. 2014, MNRAS, 444, L25

Shetrone, M., Venn, K. A., Tolstoy, E., et al. 2003, AJ, 125, 684

Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862, 114

Stetson, P. B. 1987, PASP, 99, 191

Suda, T., Hidaka, J., Aoki, W., et al. 2017, PASJ, 69, 76

Tolstoy, E., Irwin, M. J., Cole, A. A., et al. 2001, MNRAS, 327, 918

Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, arXiv:1102.1523

Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2006, AJ, 131, 2114

Wang, M.-Y., Fattahi, A., Cooper, A. P., et al. 2017, MNRAS, 468, 4887

Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 147

Yozin, C., & Bekki, K. 2012, ApJL, 756, L18

Yuan, Z., Qian, Y.-Z., & Jing, Y. P. 2016, MNRAS, 456, 3253