© The Authors 2025

1. Introduction

The chemical abundance of star-forming regions and stellar populations store relevant information about the history of galaxy formation (Tinsley 1980; Maiolino & Mannucci 2019; Jara-Ferreira et al. 2024). Most of the heavy chemical elements are produced by stellar nucleosynthesis and injected into the interstellar medium (ISM) by supernovae explosions and stellar winds (Matteucci 2012; Naab & Ostriker 2017). In the current cosmological paradigm, the Lambda cold dark matter (ΛCDM) scenario, the structure is assembled hierarchically, with the gas cooling and collapsing within dark matter halos, which grow by accretion and mergers. Within halos, the cool gas settles, and eventually, star formation activity takes place (White & Rees 1978).

The formation of disk galaxies can be explained by theso-called inside-out formation model with overall angular momentum conservation (Fall & Efstathiou 1980), in which the central regions formed first. The global conservation of angular momentum results from a balance between the infalling and the expelled material (e.g., Sales et al. 2012; Pedrosa & Tissera 2015; Somerville & Davé 2015). Within this scenario, the central regions of galaxies are enriched before galaxy outskirts and for a longer time, resulting in a negative metallicity gradient. Observations of HII regions in disk galaxies have been detected since the 1970s (Searle 1971; Peimbert 1979). Recent integral field unit (IFU) observations, such as MaNGA (Bundy et al. 2015), CALIFA (García-Benito et al. 2015) and MUSE-Wide (Urrutia et al. 2019), have provided high quality data useful for the study of the spatial distribution of metallicity in galaxies in star-forming regions (e.g., Sánchez et al. 2014; Sánchez-Menguiano et al. 2016).

Observational results have shown the existence of a mass (luminosity)-metallicity gradient relation (MZGR) in the Local Universe (Vila-Costas & Edmunds 1992; Zaritsky et al. 1994; Garnett et al. 1997). However, with the advent of large databases, the MZGR has been reshaped. Belfiore et al. (2017) found a dependence of the gas-phase metallicity gradient on stellar mass in non-interacting nearby galaxies from the MaNGA survey, suggesting that low-mass galaxies (log(M_⋆/M_⊙)∼9.5) and high-mass galaxies (log(M_⋆/M_⊙)∼11) exhibit shallower gradients than intermediate-mass galaxies with log(M_⋆/M_⊙)∼10.3 (see also Khoram & Belfiore 2025). Flat gradients are expected in high-mass galaxies as the result of their more advanced evolutionary stage (Mollá et al. 2019). Galactic winds have also been proposed to contribute to that behavior at the low mass end since their effects are greater in galaxies with shallow potential wells (Dekel & Silk 1986), where they can more easily remove metal-rich material from the inner regions of the galaxy. Belfiore et al. (2017) also found a systematic flattening in the center of the most massive galaxies and in the outskirts of some galaxies (even an inversion of the gradient in the latter). According to the baryon cycle (i.e., the interplay between gas inflows, outflows and star formation, Péroux & Howk 2020), this could be consistent with an scenario where metal-rich gas is ejected from the inner regions and feeds a metal-rich gas halo, and simultaneously the outer regions of the galaxy re-accrete this gas (e.g. Diaz et al. 1989; Perez et al. 2011). If galaxies are in interactions, the transfer of material between them and the opening of the spiral arms in response to the tidal torque generated during the interaction could also contribute to changing the metallicity gradients in the outskirts (Sillero et al. 2017).

Observations have also reported significant deviations from a steep negative gradient in interacting galaxies and mergers compared to isolated ones (Kewley et al. 2010). Inverted or positive gradients (i.e., central regions are less enriched than the outer ones) have also been detected with a increasing frequency with increasing redshift (Troncoso et al. 2014). However, new observations show a large diversity of metallicity gradients, even up to z ∼ 8 (Venturi et al. 2024; Li et al. 2025). Several numerical works have reproduced shallower and inverted metallicity gradients as a result of mergers, interactions, or even strong energetic feedback (e.g., Rupke et al. 2010a; Tissera et al. 2019).

In terms of the evolution of the gas metallicity gradients, there are several studies that have used different simulations and subgrid models that predict different levels of evolution mainly depending on the strength of feedback (Pilkington et al. 2012; Gibson et al. 2013; Stinson et al. 2013; Tissera et al. 2016; Ma et al. 2017; Ibrahim & Kobayashi 2025; Garcia et al. 2025). Recently Tissera et al. (2022, hereafter T22), analyzing simulated star-forming galaxies with M_⋆ ≥ 10⁹ M_⊙ from the EAGLE project in a redshift range z = [0.0, 2.5], found that the median oxygen abundance gradient is close to zero at all z but exhibits an increase in the scatter at higher redshift. Such a mild evolution is predicted in simulations with enhanced supernova feedback (such as the MaGICC simulation, Gibson et al. 2013). T22 also explored the stellar mass dependence of the metallicity gradients, finding that these gradients exhibit a weak positive dependence on stellar mass at z ∼ 0, a trend that becomes stronger with increasing z. Observations of high redshift galaxies using the Very Large Telescope (VLT) (Wuyts et al. 2016; Carton et al. 2018), Hubble Space Telescope (HST), (Jones et al. 2015; Wang et al. 2020; Simons et al. 2021; Li et al. 2022) and James Webb Space Telescope (JWST) (Wang et al. 2022; Arribas et al. 2024; Venturi et al. 2024; Ju et al. 2025) often report shallow metallicity gradients with significant scatter and an increasing fraction of positive metallicity gradients with redshift. Globally, simulations are able to reproduce these observations, however, at a very high redshift, it is difficult to model the physical processes with enough resolution.

Most observational and numerical studies of metallicity gradients assume a single linear relation in order to describe the radial metallicity profiles, in general, within [0.5, 1.5] r₅₀, where r₅₀ is the half-mass radius. However, the presence of inner and outer breaks of the metallicity profile in certain galaxies have been reported for decades (e.g., Martin & Roy 1995; Diaz et al. 1989, respectively) and confirmed by the larger statistical results obtained with IFU observations (Sánchez et al. 2014). Sánchez-Menguiano et al. (2018) (hereafter SM18) looked for deviations from single gradients by automatically fitting the metallicity profiles in nearby galaxies using IFU observations provided by MUSE. They found that apart from the negative trend, inner drops and outer flattenings are common in disk galaxies (see also Easeman et al. 2022; Cardoso et al. 2025).

Numerical simulations of interacting galaxies have shown that gas inflows can be triggered by tidal torques generated during these events, which retrieve angular momentum and induce bar formation (e.g., Barnes & Hernquist 1996; Mihos & Hernquist 1996; Tissera et al. 2016). These gas inflows dilute the central abundances flattening the metallicity gradients (Perez et al. 2011; Torrey et al. 2012). However, they can also feed new star formation and consequently increase the level of oxygen, thus making the metallicity profiles more negative again (Sillero et al. 2017; Acharyya et al. 2025; Pan et al. 2025), possibly associated with the formation of bars (Chen et al. 2023). The newborn stars will produce supernova feedback which could trigger metal-loaded outflows that transport material out of the galaxy (Pilkington et al. 2012). Depending on the strength of the star formation activity and the potential well of the galaxies, the ejection of material might invert the metallicity gradients, as shown by numerical simulations (Stinson et al. 2013; Tissera et al. 2019). The metallicity gradients could again recover a weak negative metallicity gradient within around [1.4 − 2] Gyr (T22). Even if galactic outflows are not powerful enough to imprint strong changes on the slope of the metallicity profiles, they can weaken them and induce galactic fountains that transport the enriched material to the disk and outer parts, thus changing the slope of the metallicity profiles (e.g., Rupke et al. 2010b; Perez et al. 2011; Torrey et al. 2012). Within a cosmological context, Garcia et al. (2023) carried out a statistical study of the metallicity gradients in the Illustris-TNG50 simulation and found a break that separates an inner highly enriched region and a mixing dominated outer disk region.

In this work, we aim to characterize the shape of the metallicity gradients of simulated galaxies in a cosmological context and to analyze the physical process involved in shaping them. For this purpose, we present a new fitting algorithm, DB-A, to automatically classify metallicity profiles. In this paper, we use the suite of cosmological simulations of the Chemo-dynamIc propErties of gaLaxies and the cOsmic web project, (CIELO; Tissera et al. 2025). We characterized the metallicity profiles of the star-forming regions of 45 central galaxies with stellar masses within M_⋆ = [10^8.5, 10^10.5] M_⊙ at z = 0. We estimated the simulated break radius and compared them with available observations. To study the different physical processes acting to shape the distribution of chemical elements in the ISM, we analyzed the metallicity profiles of the star-forming gas in the two main galaxies of a Local Group (LG)-type environment as case studies. The simulation has an adequate cadence (∼ 0.16 Gyr) to follow the impact of different process on the shape of the metallicity distribution. We selected this LG analog as a case study because it was also analyzed in two previous papers. Rodríguez et al. (2022) studied the impact of the environment on infall disk satellites, which also interacted with the central galaxies. Cataldi et al. (2023) focused on the impact of galaxy assembly on the shape of the dark matter distribution. Both works followed the infalling material and their impact. Hence, these works complement each other in relation to the study of galaxy formation and the impact of environment on baryons and dark matter.

This paper is organized as follows. In Section 2 we introduce the CIELO simulations and the DB-A algorithm used to automatically fit the metallicity profiles. In Section 3 we report the occurrence of breaks in metallicity profiles at z = 0 as well as the relation between the metallicity gradients and stellar mass, star formation, and galaxy size. We also explore the break radius dependence on stellar mass and star formation. In Section 4 we present the evolution of the metallicity gradients of our two case studies, including the identification of star formation rates (SFRs), gas accretion rates, satellite interactions, and mergers. In Section 5 we discuss in detail the physical origin of the inner and outer breaks along the evolutionary history of the case studies in the context of the previous results. Finally, in Section 6 we summarize our main findings.

2. CIELO simulations

The CIELO project aims to study galaxy formation and the chemical evolution of baryons in different field and filamentary environments. The initial conditions (ICs) have been chosen to be centered in halos with virial masses in the range M₂₀₀ = [10¹⁰, 10¹²] M_⊙ (Pehuen halos), and includes two LG analogs.

The CIELO ICs are taken from a dark matter-only run of a cosmological periodic cubic box of side length L = 100 Mpc h⁻¹ consistent with a ΛCDM universe model with Ω₀ = 0.3170, Ω_Λ = 0.6825, Ω_B = 0.049 and h = 0.6711. Zoom-in initial conditions were generated using the MUSIC code (Hahn & Abel 2011). The L12 level corresponds to high resolution with m_dm = 1.36 × 10⁵ h⁻¹ M_⊙, while the L11 level represents intermediate resolution with m_dm = 1.28 × 10⁶ h⁻¹ M_⊙. The initial gas masses are m_gas = 2.1 × 10⁴ h⁻¹ M_⊙ for L12 and m_gas = 2.0 × 10⁵ h⁻¹ M_⊙ for L11. At L11 resolution, the gravitational softening values are ϵ_g = 400 pc for gas and stellar particles and ϵ_dm = 800 pc for dark matter particles. In contrast, the L12 runs have ϵ_g = 250 pc for gas and stellar particles and ϵ_dm = 500 pc for dark matter particles. The detailed description of the initial conditions and the comparison of the main properties of the CIELO galaxies with observations can be found in Tissera et al. (2025).

Briefly, the CIELO simulations were run using a version of the P-GADGET-3 code, which includes a multiphase model for the gas component, metal-dependent cooling, star formation and supernova feedback, as described in detail in Mosconi et al. (2001), Scannapieco et al. (2005) and Scannapieco et al. (2006). Here, we only provide a summary, highlighting the main aspects which are relevant for our work. A more detailed discussion can be found in Tissera et al. (2025).

The star formation algorithm is implemented in a stochastic fashion depending on the gas density and considering gas clouds with a density, ρ_gas > 10⁻²⁶g cm⁻³, temperature T < 15 000 K, and in a convergent flow, ∇ · v < 0. A detailed description is given by Scannapieco et al. (2005). Initially, baryons are in the form of gas with primordial abundance, i.e., X_H = 0.76, Y_He = 0.24 and Z = 0. A Chabrier (2003) Initial Mass Function is assumed in the mass range [0.1, 40] M_⊙ (Chabrier 2003). The chemical evolution model includes the enrichment by Type Ia (SNIa) and Type II (SNII) supernovae¹.

The SNIa events are assumed to originate from CO white dwarf (CO WD) binary systems, in which the explosion is triggered when the primary star, due to mass transference from its companion, exceeds the Chandrasekhar limit. For simplicity, the lifetime of the progenitor systems (delay times) are assumed to be randomly distributed over the range [0.7, 1.1] Gyr. To calculate the number of SNIa, an observationally motivated relative SNIa rate is adopted (see also Jiménez et al. 2015 for a detailed comparison with the single degenerated delay time distribution (DTD) model, Mosconi et al. 2001). The CIELO project assumes the SNIa yields by Iwamoto et al. (1999).

The SNII are assumed to originate from stars more massive than 8 M_⊙. The nucleosynthesis products are taken from the metal-dependent yields of Woosley & Weaver (1995). The metallicity-dependent lifetimes of SNII are estimated according to Raiteri et al. (1996). The code traces the following 12 different chemical elements: H, ⁴He, ¹²C, ¹⁶O, ²⁴Mg, ²⁸Si, ⁵⁶Fe, ¹⁴N, ²⁰Ne, ³²S, ⁴⁰Ca and ⁶²Zn.

The CIELO galaxies have been used to study the impact of LG environment onto infall satellites (Rodríguez et al. 2022), the effects of galaxy formation on the shape of dark matter halos (Cataldi et al. 2023), the possible contribution of primordial black holes to the dark matter component (Casanueva-Villarreal et al. 2024) and the origin of the stellar halos (Gonzalez-Jara et al. 2025).

2.1. The simulated galaxies

All CIELO simulations are analyzed using a consistent methodology as explained in detailed in Tissera et al. (2025). Briefly, halos are identified at their virial radius, r_vir, using the friends-of-friends (FoF) algorithm (Davis et al. 1985). Substructures are detected with a SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009), and merger trees are constructed using the AMIGA algorithm (Knollmann & Knebe 2009). We recalculated the center of mass of each substructure by applying the shirking sphere technique. For the gas component, we looked for cold and dense star-forming regions by selecting only gas particles with a temperature of T < 15 000 K and a density of ρ > 7 × 10⁻²⁶ g cm⁻³, hereafter referred to as SF particles. We note that these are the thresholds used by the subgrid physics model to convert gas into star particles together with the requirement for the gas particles to be in a convergent flow. Hence, we are taking into account all the cold and dense gas component because we aim at describing the chemical distribution in the gas that is broadly available for star formation. The larger advantage is that we can better trace the outskirts of discs. Additionally, these results will serve as a baseline to compare with a future work where we mock observations of HII more closely, particulary at high redshift (Tapia-Contreras et al., in prep.).

Galaxies are rotated considering the direction of the total angular momentum of the cold gas. Hence, all metallicityprofiles are estimated on the face-on projections. We consider as part of a galaxy all mass enclosed within 2 × r₈₃, where r₈₃ is the galactocentric radius that enclosed ∼83 percent of the total stellar mass of the system (Tissera et al. 2019). We also estimate the stellar half-mass radius, r₅₀, the SFR, considering young stars with ages younger than 0.5 Gyr, and the specific SFR (sSFR), defined as sSFR=SFR/M_⋆.

2.2. The oxygen profiles and the DB-A algorithm

We estimate the metallicity profiles of the gas component as follows. For each galaxy, we compute the oxygen to hydrogen abundance ratio² 12 + log(O/H) for all selected star-forming regions. Galaxies with fewer than 500 SF particles were discarded. Abundance radial profiles are constructed by computing a running median over radially ordered SF particles, using a window of N = 40 particles. This significantly reduces thenumerical noise associated with the stochasticity of the subgrid methods. However, we are applying azimuthal averages, and hence, azimuthal variations related to the 2D metallicity distribution are not considered (Ho et al. 2017; Kreckel et al. 2019; Solar et al. 2020; Orr et al. 2023).

Because our aim is to study the shape of the profiles and the physical mechanisms behind them, the detection of breaks in the radial profiles needs to be automatically done. This is not a trivial task considering that there is a significant dispersion in the chemical abundances at a given radius as shown in Fig. 1. A similar situation occurs in observations, as discussed by SM18 and Chen et al. (2023), for example.

Fig. 1. Examples of the variety of metallicity profiles in the CIELO galaxies. From left to right: (i) Face-on and edge-on (top and bottom subpanels) gas density projections in a large cubic volume of 200 kpc-side centered at a central galaxy. (ii) Face-on and edge-on gas density projections in a cubic volume of 50 kpc-side located at the center of mass of the galaxy. (iii) Face-on and edge-on projected metallicity maps in the 50 kpc-side volume. Arbitrary density contours are included in (ii) and (iii) to illustrate the morphology. (iv) Gas-phase metallicity profile of the central galaxy (background heat map). For illustrative purposes, median values within 0.5 kpc bins are shown as filled squares, with error bars indicating the standard deviation in each bin. The DB-A fit adopted in this work (see Sect. 2.2) is shown as a black solid line. For comparison, a linear fit within [0.5–1.5] r50 is also displayed (red dashed line). Vertical lines are included to indicate the gravitational softening ϵ (gray, dashed line), inner break radii (teal solid line) and/or outer break radii (orange solid line) when appropriate.

For this purpose, we developed an algorithm, hereafter referred to as DB-A, that allows for the possible existence of an arbitrary number of breaks. However, based on observational constraints, we set a maximum of two breaks in a given metallicity profile, i.e., an inner break and/or an outer break. No fixed bounds were placed for the break radius, although conflictive cases (i.e., very small or very large break radius or inner and outer break radius too close to each other) were manually revised, decreasing the number of breaks used to fit the profiles. The algorithm applies a Gaussian Kernel Density Estimation³ to estimate the probability density of every particle in the profile and to eliminate outliers if such a value is lower than 0.05. DB-A is based on a piecewise linear function and applies a least-squares method to minimize the error. The algorithm works as follows:

1. The maximum number of breaks of the piecewise function n_max and a tolerance parameter τ are defined. We selected n_max = 2, assuming that the limit case is a doubly broken linear profile; and τ = 0.99.

2. The code iterates over the number of breaks of the piecewise function, n_br, such that 0 ≤ n_br ≤ n_max. The respective piecewise function is fit in each case, considering scenarios of a single linear trend, a profile broken once, and a profile broken twice. The resultant error⁴, ε, is stored.

3. In every iteration we compare the relative errors (ε_n/ε_n − 1) with our τ parameter. If ε_n/ε_n − 1 < τ, the current n_br is chosen to be fit. Otherwise, the lower value of n_br is preferred.

We tested the robustness of our method against both the Bayesian Information Criterion and the Akaike InformationCriterion. The results obtained with each criterion were consistent, reinforcing the reliability of our approach.

We define inner, intermediate and outer regions according to the detection of breaks. The metallicity gradient, ∇_(O/H), is then defined as the slope of the linear function that fits the distribution of chemical abundances in each region, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$. We note that when the profile is doubly broken, both inner and outer gradients can therefore be directly defined. However, when the profile has a single break, r_br, this can be either an inner or an outer break. Based on the characteristic break radius reported by B17 and SM18, our criteria to classify them are the following:

1. If r_br < 1.5 × r₅₀, the break is considered an inner break, and only inner and middle ∇_(O/H) are defined.

2. If r_br > 1.5 × r₅₀, the break is considered an outer break, defining only middle and outer ∇_(O/H).

In the case that a profile is not broken, only the middle gradient is defined, which corresponds to a linear fit to the SF particles within [ϵ, 2 × r₈₃ ].

In Fig. 1 we display several examples to illustrate both the diversity of metallicity profiles and the performance of our fitting method. In the right panels of Fig. 1 there are examples of linear, inner break, outer break and doubly broken profiles. It is important to note that both inner and outer breaks can exhibit either negative or positive gradients. Gas density projections of the simulated galaxies and their environment suggest that the shape of the metallicity profile is affected by the interaction with nearby galaxies as expected. We also show the case of a linear inverted metallicity profile (positive metallicity gradient), which is associated with a clearly disturbed gas distribution (see Rupke et al. 2010a; Fragkoudi et al. 2016; Sillero et al. 2017; Tissera et al. 2019). We will come back to this point later.

For comparison, we also fit each profile a linear regression (TS)⁵ within a fixed radial range of r ∈ [0.5, 1.5]× r₅₀. Figure 2 shows a comparison between the slope estimated with the two methods, ${\mathrm{\nabla }}_{\text{(O/H)}}$ and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$. We note that for every galaxy we have one gradient estimated with the TS method and up to three gradients estimated with the DB-A method, corresponding to the slope in each region. Figure 2 displays ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ versus ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ for the four different types of profiles identified. In the case of broken profiles, we also added the corresponding inner or outer slopes. For every type of profile, there tends to be a good agreement between the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, favoring a direct comparison between both approaches. Small differences are likely to be related to the adopted radial ranges. When the profiles exhibit an inner break, they show the largest discrepancies between the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ and the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$.

Fig. 2. Comparison between the metallicity gradients (filled black circles) estimated by the DB-A fit, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, and by the Theil-Sen method, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$, for galaxies exhibiting, from top to bottom, (i) linear profiles, (ii) inner break profiles, (iii) outer break profiles, and (iv) double break profiles. When corresponds, inner (teal symbols) and outer (orange symbols) metallicity gradients are also shown. Vertical lines are included as a visual aid to identify gradients belonging to the same galaxy.

3. The gas-phase metallicity gradients in the Local Universe

From the analysis described in the previous section, we found 17 (38%) linear, 6 (13%) inner-broken, 18 (40%) outer-broken, and 4 (9%) doubly broken profiles at z = 0. Table A.1 summarizes the best fitting parameters to each profile, including the respective abundance gradients in each region and the break radii if corresponds. In the left, upper panel of Fig. 3 the MZGR is displayed using the metallicity gradients of the intermediate regions, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$. The different symbols and color code indicate the global shape of the profiles (i.e., inner, outer or no breaks). For comparison the metallicity gradients fit by SM18 for MUSE galaxies and by Cardoso et al.(2025) for CALIFA galaxies are also displayed in a similar fashion. Additionally, the metallicity gradients reported by Chen et al. (2023), obtained by fitting a broken profile to bar disk galaxies, are also included. The ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ (dex r${}_{50}^{-1}$) of the CIELO galaxies generally aligns with the observations. However, a larger variety is detected for simulated galaxies with an outer broken profile, a feature observed across the entire range of stellar masses. SM18 found 53% of their sample to show linear negative gradients while 21% have inner breaks, 10% only outer breaks and 17% are doubly broken. These frequencies differ from those obtained for the CIELO galaxies, where the majority exhibit breaks. We estimated the binomial probability of our results given the SM18 rates, and found a low probability (p < 0.01), which indicates that the frequency of breaks in CIELO is larger than expected from SM18 observations. We note that the stellar mass range covered by the CIELO galaxies is different from that of SM18 and this could be an important factor (Belfiore et al. 2017). Additionally, this difference could also reflect the fact that in the simulated galaxies we used the whole cold and dense gas distribution while observations are restricted to detections limits, which can prevent the measurement of the outer breaks, for example. Extending both the observed and the simulated samples to cover a similar mass range of galaxies could allow for a more robust assessment of the differences, which could also provide clues regarding galaxy formation models and the adopted subgrid physics.

Fig. 3. Stellar mass-, sSFR-, and size-metallicity gradient relation for the 45 star-forming central galaxies (left, middle, and right panels, respectively). In the upper panels, the relations were constructed by using only ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$. The colors of solid symbols encode the type of profile that was best fit to the whole metallicity distribution, inner break (IB; teal triangles), outer break (OB; orange triangles), and doubly broken (DB; gray hexagons). In the lower panels are the same relations but constructed using ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ (red circles). Observational data from Martin & Roy (1995), Ho et al. (2015), and SM18 (MUSE); Franchetto et al. (2021) (MaNGA+MUSE); Chen et al. (2023) (TYPHOON); Cardoso et al. (2025) (CALIFA); and MaNGA DR17 are also displayed in a similar fashion. The black dashed line indicates ∇(O/H) = 0 for reference purposes only. The two main galaxies LG1-4469 and LG1-4337 of the LG analog are highlighted, as they are later analyzed in detail as case studies.

In the lower-left panel of Fig. 3, the MZGR is displayed using ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ for the CIELO sample. For comparison, observations reported by SM18, i.e., those which follow a linear profile, and by Franchetto et al. (2021), Ho et al. (2015) and Belfiore et al. (2017), where a linear regression is fit. As can be seen, the simulated values are within the observed range. There are some galaxies which show steeper negative gradients but, nevertheless, within the observational range. In the next section we further discuss the origin of the inverted/positive slopes in the metallicity distributions.

In the middle and right panels of Fig. 3, we display the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ (upper panels) and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ (lower panels) as a function of sSFR and r₅₀, respectively. As can be seen the relation between ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ and both parameters are within the observed range. We note that CIELO galaxies are smaller compared to the SM18 as expected because of the differences in stellar mass as can be seen from the left panels of this figure. At z = 0 we find only three galaxies with positive ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ (see Table A.1). These galaxies show signals of recent interactions.

Regarding the inner and outer slopes, we performed similar relations and include available data to confront them. In Fig. 4, the corresponding relations for the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$ (upper panels) and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ (lower panels) are displayed. As can be seen the simulated values are within the observed range. In the case of the inner drops, they match the correlation signal reported by Cardoso et al. (2025). For low stellar mass galaxies, inner rises or negative slopes appear to be preferred. The simulated slopes are consistent with observations reported in Chen et al. (2023).

Fig. 4. Similar to Fig. 3 but featuring the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$ and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ in broken metallicity profiles. Observational data from Martin & Roy (1995), SM18 (MUSE), Chen et al. (2023) (TYPHOON) and Cardoso et al. (2025) (CALIFA) are included for comparison.

The lower panels of Fig. 4 displays the same relation but for the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$. There is a significant dispersion that decrease with stellar mass. Observations and simulations report similar slopes for galaxies with M_⋆ > 10^9.5 M_⊙. Simulated galaxies with positive ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ tend to have more active star formation. Positive slopes in the outer region could be linked to accretion ofpre-processed material or satellite interactions (Perez et al. 2011). We will discuss these different physical processes by using case studies in the following section.

Finally, in Fig. 5, we display the characteristic inner and outer radii as a function of stellar mass and sSFR. We also included observations from SM18. The simulated inner breaks have a median (25th – 75th percentiles) of 1.08 (0.96–1.26) r₅₀ and outer breaks of 2.62 (1.77–2.912) r₅₀. We do not detect any statistically significant trend with stellar mass or sSFR for the simulated or the observational data as quantified by the Spearman Correlation Coefficient (SCC), except for the observed inner breaks, which seems to define an anti-correlation with sSFR (see also Cardoso et al. 2025).

Fig. 5. Inner (teal symbols) and outer (orange symbols) break radii as a function of (i) galaxy stellar mass (left panel) and (ii) sSFR (right panel). The CIELO galaxies are displayed as filled symbols, while empty symbols represent observations from Martin & Roy (1995), SM18, Chen et al. (2023) and Cardoso et al. (2025) as indicated in the legend. The SCC was calculated to evaluate possible correlations between galaxy properties and inner and outer break radii. The upper part of each panel shows the SCC computed considering: CIELO galaxies (top rows) and observations (bottom rows). The p-value is shown in parenthesis in each case.

In Fig. 3 and Fig. 4 we have highlighted with larger symbols two galaxies, LG1-4469 and LG1-4337. These target galaxies will be used as examples to analyze the impact of different physical mechanisms on the shape of the metallicity profiles. At z = 0, LG1-4469 and LG1-4337 have a double break and an outer break, respectively.

4. Case studies on the evolution of gas-phase metallicity gradients

In this section, we trace LG1-4469 and LG1-4337 from the LG1 analog, back in time along their merger trees, using them as case studies. We analyze the details of their evolution to gain insights into the physical processes that can generate these breaks and the timescales over which these profiles change shape. We expect the evolution of the metallicity gradients to reflect their distinct evolutionary histories (e.g., Ma et al. 2017; Tissera et al. 2022), but to share common features. At each available snapshot, the same procedure applied to galaxies at z = 0 (see Sect. 2.2) is performed on their progenitors up to z ∼ 6. Hence, between z = 0 and z ∼ 6 this analysis was applied to 87 available snapshots⁶. We remind the reader that only systems with more than 500 SF particles are considered.

At all available redshifts, we rotate the gaseous discs along the direction of total angular momentum, estimate the projected metallicity gradients and applied the DB-A method to all of them. Hence, at each z we have the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ with the corresponding breaking radii. Additionally, we estimated r₈₃ and r₅₀ for the progenitors of both galaxies.

We define the regions associated with breaks or equivalent ones to evaluate the impact of the acting physical processes. For this purpose, we defined the inner, middle, and outer regions as follows:

1. Inner region: r ∈ [0, r_inner) kpc

2. Middle region: r ∈ [r_inner, r_outer) kpc

3. Outer region: r ∈ [r_outer, 2 × r₈₃ kpc

where

$$\begin{array}{c}\hfill {\mathrm{r}}_{\mathrm{inner}}=\{\begin{array}{cc}{\mathrm{r}}_{\mathrm{br}},\hfill & \hfill \text{if the profile has an inner break at}{\mathrm{r}}_{\mathrm{br}}\\ {\mathrm{r}}_{50},\hfill & \hfill \text{otherwise}\end{array}\end{array}$$(1)

$$\begin{array}{c}\hfill {\mathrm{r}}_{\mathrm{outer}}=\{\begin{array}{cc}{\mathrm{r}}_{\mathrm{br}},\hfill & \hfill \text{if the profile has an outer break at}{\mathrm{r}}_{\mathrm{br}}\\ {\mathrm{r}}_{83},\hfill & \hfill \text{otherwise}\end{array}.\end{array}$$(2)

In each of these regions, we computed the SFR in a similar way as explained in the previous section. Since inner, middle and outer regions can differ in size during the evolution of the galaxies, we also calculate the SFR surface density, Σ_SFR, within each of the three defined regions (see Appendix C).

4.1. Environmental and secular processes

To understand the physical processes modulating the shape of the metallicity profiles, we explored the relationship between the central galaxy and its environment by following the accretion of gas from different sources and the interaction and merger histories. Satellite interactions and merger events are mechanisms that affect the galaxy dynamics and distribution of both gas and stars, in addition to being able to contribute gas for star formation as discussed in the Introduction. We also search for gas particles accreted from the circumgalactic medium (CGM) (Péroux & Howk 2020) and from cold flows (Ceverino et al. 2016).

4.1.1. Surviving satellites

In order to study the effects on the abundance distributions associated with surviving satellites, we identified all the substructures surrounding the central galaxy at z = 0, and select as satellite galaxies those whose center of mass are within r_vir of the central galaxy, and have a stellar mass log(M_⋆/M_⊙) > 7. An interaction is defined from the time a satellite enters the virial radius of a central galaxy. We followed their evolution back time and obtained their stellar masses, gas fractions and distance to the center of mass (d_CM) of the progenitors of our two target galaxies up to z ∼ 6. These properties are displayed in detail in Appendix B as a function of time. The most significant periods of interaction, in terms of stellar mass fraction and gas content, are highlighted in the upper panels of Fig. B.1 (dashed, green regions).

4.1.2. Mergers

Furthermore, we followed the merger tree of each central galaxy to identify the most significant merger events that occur during each galaxy assembly process. We selected mergers with a stellar mass ratio ${\mathrm{M}}_{\star ,\mathrm{merger}}/{\mathrm{M}}_{\star ,\mathrm{central}}>0.01$ at the infall time (i.e., when the center of mass of the merging galaxy crosses r_vir of the central galaxy progenitor). We define a merger event when the infalling satellite is not longer identified as a separate substructure by the SUBFIND algorithm. In terms of the gas content of the merging galaxies, we considered both wet and dry mergers. A high amount of gas can be accreted into the central galaxy during wet mergers, while the torque exerted by dry mergers can also have an impact on the redistribution of gas within the central galaxy. The properties of such mergers are displayed in detail in Appendix B. Analogous to satellite interactions, the selected merger events are highlighted in the lower panels of Fig. B.1 (dashed, blue regions).

4.1.3. Accretion rates

In order to understand the gas motion and infall into the target galaxies, we estimated three gas accretion rates, ̇ℳ, one for each region defined at the beginning of Sect. 4. The accretion rates are calculated as the gaseous mass that enters the inner, middle or outer regions per unit Gyr (rates per unit area are provided in Appendix C). We also considered five different accretion sources. To do this, at a given time, the gas particles that currently belong to the central galaxy, but did not belong to its progenitor in the previous available snapshot, are followed back in time for ∼ 1 Gyr to check if: they belonged to another galaxy (hereafter, satellite gas); they were not gravitationally bound to any structure (hereafter, unbound gas); or they were part of the gas surrounding the progenitor or CGM, defined as the gas located further than 2 × r₈₃) but within the r_vir.

In addition, we also estimated the redistribution of gas within a galaxy (ISM) as the gaseous mass that moves from one region to another, in any direction. Radial and vertical displacements of gas particles between two consecutive snapshots would be noted as δr and δz, respectively. We discarded radial displacements smaller than the gravitational softening (ϵ_g ∼ 0.4 physical kpc for LG1 at z = 0). As a convention, we will use δr < 0 to refer to particles moving outward (hereafter, gas radial outflows), and in the opposite case, δr > 0 to refer to particles moving inward (hereafter, gas radial inflows).

4.2. Evolution of the star formation rate and gas accretion

Figure 6 shows the temporal evolution of both the SFR and $\dot{\mathcal{M}}$ in units of M_⊙ yr⁻¹. Fig. C.1 shows the evolution of the Σ_SFR and Σ_̇ℳ per region (i.e., SFR and $\dot{\mathcal{M}}$ normalized by the area of each region). As can be seen, both quantities are strongly correlated in the inner region, particularly when there is active star formation. This suggests that the accreted gas is efficiently condensed and converted into stars. In the middle and outer regions, both quantities evolve similarly but the amount of accreted gas tend to be ∼1 − 2 orders of magnitude higher than the SFR. This behavior can be explained by the lower gas density with respect to the central regions. This is consistent with the low star formation efficiency typically observed in the outskirts of galaxies compared to their central regions (Leroy et al. 2008; Bigiel et al. 2008; Krumholz et al. 2009). It suggests that gas accumulates in the disk and is steadily converted into stars over time. At low redshift, the outer region of LG1-4337 provides an example of active accretion with very low star formation activity, emphasizing the potential influence of additional physical processes beyond accretion in regulating star formation activity in galactic outskirts, thereby adding to the complexity (e.g., turbulence).

Fig. 6. Time evolution of the SFR (shaded areas) and the total accretion rates (̇ℳ; solid lines) per region (inner, left column; mid, central column; outer, right column) for each galaxy (LG1-4469, top row; LG1-4337, bottom row). The epochs where inner or outer breaks are identified are denoted by using teal or orange symbols, respectively (see Section 4 for a detailed explanation on the definition of the different regions). Figure C.1 shows the same quantities normalized by the area of each region.

In order to understand the origin of the gas accretion, in Fig. 7 we display contribution of each of the defined accretion sources. In the inner region, the dominant accretion source transitions from unbound gas (blue areas) to inflowing ISM (pink area), for decreasing redshift. this indicates that at later times, most of the gas that reaches the central regions is primarily accreted from the ISM reservoir in the discs. Peaks in accretion from satellites and the CGM can results in up to f_acc ∼ 0.4 − 0.5, often driven by interactions that either directly transfer gas from companions or trigger gas inflows through angular momentum loss from the disk ISM (pink areas) and CGM (yellow areas).

Fig. 7. Fractions of the total accreted gas (̇ℳ; see Fig. 6) in the inner (left column), middle (central column) and outer (right column) regions for each galaxy (LG1-4469: top row; LG1-4337; bottom row). The gas accreted fractions correspond to inflowing ISM (pink areas), outflowing ISM (red areas), the CGM (yellow areas), interacting/merging satellite (green areas) or unbound gas (blue areas).

The middle regions show a more diverse origin for accreted gas, although the infall of ISM from the outermost disk remains the dominant source (pink areas), typically contributing between 0.3 and 0.8 to the fraction of gas accretion. Surges of accreted gas from satellites and the CGM, comparable to those observed in the inner region, occur within the same time intervals. Additionally, outflowing material of the ISM (red areas), potentially redistributed by supernova feedback or the disturbance of the discs produced by tidal torques during close encounters, contributes up to f_acc ∼ 0.3 (see also Sillero et al. 2017).

The outer regions show the most diverse sources of accretion, with distinct patterns observed in the two case studies. The high-mass galaxy, LG1-4337, displays a smooth transition where unbound gas is gradually replaced by CGM gas as the dominant accretion source. This CGM gas has been likely processed in central region of the galaxy, later ejected by supernova feedback to feed the CGM, and subsequently re-accreted due by the deep gravitational potential of the galaxy (Tumlinson et al. 2017; Péroux & Howk 2020; Péroux et al. 2020). Peaks in accretion of radially outflowing gas correlate with mergers (see Fig. B.1) that can disrupt part of the external spiral arms and trigger star formation bursts (see Fig. 6) that subsequently drives galactic outflows (Miranda et al. in prep.).

For the low-mass galaxy, LG1-4469, the transition begins with unbound gas being replaced by CGM gas (lookback time of 10–7 Gyr) and later by outflowing ISM gas (lookback time: of 7–5 Gyr). During the first period, the galaxy is disrupted by a merger event (see Fig. B.1), with CGM accretion resulting from the relaxation of the galaxy and rapid infall of surrounding gas remnants.

In the second period, the contribution of the outflowing ISM becomes significant. This corresponds to radially heated particles and the opening of spiral arms caused by interactions with a surviving satellite (see Fig. B.1). Both galaxies also receive significant contributions from the accretion of satellite gas during close encounters, as previously observed in the inner and mid regions.

Globally, at earlier times, unbound gas, representing cold flows from cosmic filaments, contributes significantly across all regions. This source dominates (f_acc ≥ 0.7 ) until approximately z ∼ 1 in the outer regions and z ∼ 1.5 in the other two inner regions. Later on, f_acc decreases down to less than 0.05.

While each galaxy exhibits a unique accretion pattern in each region, reflecting the interplay between their internal properties and the characteristics of their environment over time, the physical processes contributing gas remain the same. Hence, what varies with redshift is their relative importance and the properties of the accreted gas. In the following section, we explore the impact of these physical processes on the shape of the metallicity gradients.

4.3. Evolution of the metallicity gradients

Figure 8 illustrates the evolution of the gradients for both targeted galaxies across time. We show, when corresponds, the gradients in the inner (teal symbols) and outer (orange symbols) regions defined by the DB-A method, as explained in Section 2.2. The gradients are normalized by the corresponding r₅₀ at every time. The general and consistent behavior observed in these galaxies is that the metallicity gradient undergoes substantial variability at high redshift, progressively transitioning to a more gradual evolution at lower redshifts, with this trend being particularly evident in the middle gradient. This behavior is consistent with previous works (e.g., Ma et al. 2017; Tissera et al. 2022), indicating a higher variability of the metallicity gradients associated with a higher frequency of merger and gas-richness of galaxies at higher redshift.

Fig. 8. Evolution of the metallicity gradients as a function of lookback time (and redshift) for the targeted galaxies LG1-4469 (upper panel) and LG1-4337 (lower panel). The metallicity gradients in the inner (teal symbols), middle (black solid lines), and outer (orange symbols) regions are displayed. The shaded regions depict the mergers and interaction events (blue and green, respectively).

In Fig. 8, the shaded regions indicate periods when the target galaxy is involved in an interaction with surviving satellites (green regions) or undergoing a merger event (blue regions). These events were identified and selected as described in Sects. 4.1.1 and 4.1.2, respectively. The full merger history is shown and discussed in the Appendix (see Fig. B.1), including the information about the gas content of the interacting galaxies and its mass ratio with respect to the central galaxy. It can be seen that the periods of greater variability in the metallicity gradient are linked to these events, with a stronger impact at high redshift due to the higher gas fraction and relative mass ratio of the galaxies involved.

Inner and outer gradients are also steeper at a higher redshift and tend to be more negative. It is clear from Fig. 8 that, at least for the target galaxies, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ tends to be negative at higher redshift, but as the galaxies evolved they became flatter or even positive. In the case of ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$, most of them are negative, and they became positive only in few occasions.

In the following sections, we focus the discussion on the physical processes behind the occurrence of inner and/or outer breaks. As shown in Fig. 2 and previously discussed in Sect. 3, the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ is similar to the gradient estimated with more conservative methods assuming a single linear function. Hence, we further investigate the global evolution of the gradients in this region, which has been extensively discussed in previous works (Gibson et al. 2013; Ma et al. 2017; Yates et al. 2021; Tissera et al. 2022; Hemler et al. 2021; Vallini et al. 2024; Venturi et al. 2024). Instead, we focus on the processes that could originate outer and inner breaks.

As shown in Fig. 8, inner breaks occur at different epochs throughout the entire analyzed redshift range. The two target galaxies exhibit a total of 46 inner breaks (21 for LG1-4469 and 25 for LG1-4337), representing 26% of the total analyzed metallicity profiles. These breaks can appear either as a single break or as part of doubly broken profiles.

1. Inner rise: The inner gradient is negative and steeper than the mid-region gradient (${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$ < ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ ).

2. Inner drop: The inner gradient is flatter or inverted compared to the mid-region gradient (${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$ > ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ ).

Using this classification, we found that 72% of inner breaks correspond to inner rise profiles, which are characterized by a high central metallicity and a steep negative slope in the inner region. The remaining 28% are inner drop profiles, which display diluted central metallicities, resulting in a flattened or inverted metallicity gradient in the inner region.

We find outer breaks to be the more common feature across time, being present in a total of 75 metallicity profiles (35 for LG1-4469 and 40 for LG1-4337), 43% of the complete sample. We found two scenarios for the outer region:

1. Outer rise: The outer gradient is flatter or inverted compared to the mid-region gradient (${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ > ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ ).

2. Outer drop: The outer gradient is negative and steeper than the mid-region gradient (${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ < ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ ).

Outer rises, characterized by high metallicities in the outskirts, are found in 68% of the cases. Outer drops are present in the remaining 32%, and they are found almost exclusively at high redshift (88% at z > 1; see Fig. 8).

5. Discussion on the origin of inner and outer breaks

In this section we discuss the origin of the inner and outer breaks using the analysis we have performed so far. We apply our analysis to our two case studies to illustrate the action of the mentioned physical processes. There might be variations in the response of the metallicity distribution to these effects depending on the specific conditions of the progenitors and their environment. We consider this detailed analysis as the foundation for a statistical study that will be conducted in a separate paper and will consider galaxies in different environments.

5.1. Inner rises

From Fig. 7 we can see that the inner region receives gas from two main sources, the ISM at low redshift and unbound gas at high redshift. Additionally, from Fig. C.1, we can appreciate that the inner region displays the highest accretion rates per unit area compared to the middle and outer regions.

Most of this material originates from ISM gas moving into the center of the galaxy, driven by tidal torques produced in galaxy interactions (see Fig. B.1). While the unbound gas is found to fall along filamentary structures.

In Fig. 9 we display the motions, the relative level of enrichment and temperature of the accreted gas according to their origin during a key period of evolution when the profiles exhibit an inner rise. The gas motion is quantified by δr and δz (see definition in Section 4.1.3). The relative level of enrichment, Δ(O/H), is defined as the difference between the metallicity of the accreted gas and that of the respective region. Negative values indicate that the accreted gas is less enriched than the preexisting gas in the region. Finally, the temperature distribution of the accreted material.

Fig. 9. Diagram of δr − δz illustrating the motion of the gas accreted into the inner region of LG1-4337 between lookback time ∼4.8 Gyr and ∼3 Gyr. During this period, the galaxy exhibits an inner break of its metallicity profile. Left panel: KDE distribution of the accreted particles. The colored numbers displayed at the top represent the fractional contribution of each gas accretion source, consistent with the color code used in Fig. 7 and the legend. The numbers below indicate the fraction of particles with δr < 0 (outflowing) and δr > 0 (inflowing). Central panel: Distribution of accreted particles, color coded by Δ(O/H), which represents the difference between the metallicity of the accreted gas and that of the respective region. Negative values indicate that the accreted gas is less enriched than the preexisting gas in the region. The median (med) value is computed over the entire time interval, while the maximum (max) value corresponds to the largest difference observed within that period. Right panel: Distribution of accreted gas particles colored by their temperature.

Overall, the dominant process that change the inner slopes is the infalling of cold gas from the disk ISM, with lower metallicity than the central regions (Kewley et al. 2006; Ellison et al. 2013; Bustamante et al. 2020). We find that this gas inflows first dilute the central metallicity but soon after they also trigger star formation which injects new elements producing a negative increase of the metallicity gradient as shown by previous works (Perez et al. 2011; Torrey et al. 2012). In fact, Sillero et al. (2017) observed a similar pattern in simulations of interacting galaxies, where interactions induced gas inflows, boosting SFR and diluting central metallicity, leading to flatter gradients.

The impact of mergers is detected at high and low redshift with the difference that there is larger frequency of wet mergers at high redshift (Appendix B). Wet mergers can trigger stronger starbursts and hence, increase more significantly the level of enrichment. At the same time, they induce stronger outflows which can flatten the metallicity gradients again. An example of the sequence of evolution in the shape of the metallicity gradients can be seen in Fig. 10 for LG1-4469 (upper two panels) and LG1-4337 (lower two panels). We find a difference in the lifetime of these features, which seems to depend on redshift. Galaxies at high redshit can recover their overall negative gradients faster because of higher probability of gas accretion. For these examples, the inner rise lasts about 150 Myr. While at lower redshift, mergers are less gas-rich and accretion from cold flows is very small. This can be clearly seen in Fig. 10, where we display examples of this temporal evolution of inner broken metallicity profiles. Both galaxies, with initially negative gradients, get their metallicity profiles flatten near the first apocenter due to the infall and accumulation of metal-poor gas, consistent with previous findings (e.g., Rupke et al. 2010a; Perez et al. 2011; Sillero et al. 2017). This is followed by a starburst (∼0.16 Gyr) that increases the central oxygen abundance by ∼0.2 − 0.3 dex, inducing a break in the profile. Chen et al. (2023) report similar metallicity profiles in observed barred galaxies with high central SFRs (see also Martin & Roy 1995) and Zahid & Bresolin (2011).

Fig. 10. Temporal sequence showing the evolution of the inner rise metallicity profiles for LG1-4469 (first and second rows) and LG1-4337 (third and fourth rows). The labeled images depict gas density projections at different snapshots, while the corresponding metallicity profiles are displayed below, following the same format as in Figure 1. We note that the timescales differ between the two galaxies.

The lower mass LG1-4469 retains broken profiles for ∼0.3 Gyr before supernova feedback ejects significant gas from the central region, flattening the profile, as can be followed in Fig. 8. In contrast, the high-mass galaxy retains metal-rich gas against galactic winds (Peeples & Shankar 2011; Tumlinson et al. 2017), with inner breaks surviving about ∼1.5 Gyr (Fig. 8). Over time, breaks migrate outward (in units of r₅₀ ) and the level of enrichment at the break radius becomes higher (∼0.1 dex), reflecting inside-out chemical evolution tied to the star formation history. After ∼1.5 Gyr, the break shifts to the outer regions, leaving the outskirts with weak or flat slopes, likely influenced by additional processes discussed in Sections 5.3 and 5.4.

5.2. Inner drops and inverted gradients: The effect of energetic feedback

Inner drops are identified in the metallicity profiles of the two analyzed CIELOgalaxies although they are less frequent than inner rises. Indeed, they are extensively reported in observational studies of more massive local galaxies (e.g., Sánchez-Menguiano et al. 2018; Cardoso et al. 2025). Observational challenges have been recognized in explaining the presence of inner drops in local galaxies, highlighting that the shape of the metallicity profile may depend on the chosen diagnostic method (Easeman et al. 2022). Additionally, contamination from active galactic nucleus (AGN) and diffuse ionized gas (DIG) can affect the identification of HII regions, impacting metallicity estimations (Cardoso et al. 2025). Hence, the comparison with observations should be always taken with caution (see also Acharyya et al. 2020; Cornejo-Cárdenas et al. 2025).

Despite these challenges, correlations between global galaxy properties and the presence of inner drops suggest a physical origin for this feature. SM18 find that inner drops are more common in high-mass galaxies, a result later confirmed in MaNGA galaxies by Easeman et al. (2022). Additionally, these galaxies exhibit systematically lower sSFR, suggesting that inner drops may be linked to quenching processes in the central regions. A large series of works reported lower metallicity and enhanced star formation in interacting galaxies, suggesting that tidal torques generate low metallicity inflows that trigger star formation (Kewley et al. 2006; Ellison et al. 2020). Recent results from Garay-Solis et al. (2024) found both increasing (decreasing) central abundances related with decreasing (increasing) sSFR. This is consistent with our results where we show that the metallicity evolves along the interaction and merger, modulated by the characteristics of the event and the gas available for star formation.

The lower occurrence of inner drops in the two analyzed CIELO galaxies compared to observations could arise from the differences in the simulated stellar mass range with respect to observed ones and the subgrid physics for the injection of α-elements by SNII, which is simultaneous with the injection of energy or the absence of AGN feedback, although our CIELO sample does not include very massive galaxies, which are expected to be the most affected by this feedback (e.g., Somerville & Davé 2015). We highlight that numerical resolution can also affect the detection of inner features at small radii.

Interestingly, Figs. 6 and 8 reveal that inner drops occur temporally close to strong SF bursts and inverted ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, suggesting a common physical origin to both features. During these periods, strong supernova feedback is expected to generate metal-loaded outflows. A visual inspection of the profiles supports this interpretation, as galaxies with inner drops are morphologically disrupted, further demonstrating the impact of energetic feedback on gas distribution. This disruption is also common in inverted ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ galaxies.

In Fig. 11 we present a set of examples of galaxies showing different signals of interaction at different redshifts (decreasing redshift from left to right), showcasing diverse patterns in abundance distribution (bottom row), including a double break with an inner drop, two inverted ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ and even a negative ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$. All of them display ring-like distributions in the gas density maps as can be seen from this figure (top row). We also present metallicity maps from an edge-on perspective (central row), including streamlines to depict the gas motions. These maps reveal the presence of bipolar, metal-loaded outflows at all redshifts. If the outflows are strong enough while star formation is concentrated in outer regions of the galaxy, a clear inverted gradients can be produced as can be seen in the third column of Fig. 11.

Fig. 11. Examples of morphologically disrupted galaxies showing evidence of central gas ejection driven by strong energetic feedback. The top row displays face-on gas density projections; the middle row shows edge-on metallicity maps with streamlines indicating gas velocity directions; and the bottom row presents metallicity profiles, formatted as in Figs.1 and 10. The corresponding galaxy and lookback time/redshift are labeled in the upper panels.

Finally, the fourth column showcases a galaxy with a negative ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, but with evidence of ejected gas, as seen from the lack of gas in the central region in the density projection and the metallicity profile, and proven by the strong bipolar metal-loaded outflows seen in the metallicity map. This galaxy actually corresponds to the later stage of an inner rise, previously discussed in Sect. 5.1 and shown in the right-most column of Fig. 10. In this case, the effect of feedback does not generate an inner drop or inverted gradient, but instead acts to erase the previously developed inner rise.

The sample of 13 inner drops identified in our case studies originate from the discussed processes at low and high redshift. In a future analyzes we expect to draw statistical results regarding an evolution of the efficiency and dependence on the characteristics of the involved physical mechanisms.

5.3. Outer rises

Observations of outer rises (also referred to as outer flattenings in the literature) have been extensively reported for galaxies in the local universe (Bresolin et al. 2012; Sánchez et al. 2014; Sánchez-Menguiano et al. 2016, 2018) and for the far side of the Milky Way (Minniti et al. 2020).

In Sect. 5.1 we discussed how an inner break, due to the continuous and extended SF activity, is able to move outward until it is classified as an outer break. A similar scenario was proposed by Wang & Lilly (2022), assuming a significant coplanar inflow of gas into the disc. In their model, the gradient inside the break radius is the result of the radial inflow of gas being enriched by the in-situ star formation, and the outer flattening is interpreted as a metallicity floor set by the level of enrichment of the inflow.

Bresolin et al. (2012), analysing a small sample of local galaxies, discarded the possibility that a significant fraction of the metals in the outer region are produced by in situ star formation. Alternative possibilities, such as radial gas mixing (Bresolin et al. 2009; Werk et al. 2011; Bresolin et al. 2012) and the re-accretion of enriched material in a galactic fountain scenario driven by supernova feedback (Davé et al. 2011; Bresolin et al. 2012; Belfiore et al. 2017) has been proposed to explain these enriched outskirts. Our previous discussion on the properties of the accreted gas supports this hypothesis (Sect. 4.1.3). Figure 12 (lower panels) showcases the motion, Δ(O/H) and temperature of the gas component in a galaxy with an outer break in the indicated time period. This figure shows that most of the accreted material that comes from the CGM is slightly more enriched than the preexisting one, exhibiting Δ(O/H) ∼ 0, while there is also a significant part of outflowing material from the ISM, which is hot and systematically more metal-rich. Both processes coexist and simultaneously affect the properties of the ISM, resulting in a slightly flatter slope in the outer region.

Fig. 12. Diagrams of δr − δz showing the motion of the gas accreted into the outer region of LG1-4337 from different sources (similar to Fig. 9) between lookback time 11.60 Gyr and 10.1 Gyr (upper panels) and between lookback time 6.5 Gyr and 4.9 Gyr (lower panels). The diagrams in the upper panels illustrate the accretion of a high-z galaxy that exhibits an outer drop, with a dominant contribution from metal-poor unbound gas (62% of the accreted gas). The diagrams in the lower panels illustrate a low-z galaxy with an outer rise, with the clear indication recycled metal-rich gas from the CGM being the main source of accretion (56% of the accreted gas).

Figure 13 shows examples of outer breaks for both case studies (upper panels: LG1-4469; lower panels: LG1-4337). Columns one to four correspond to temporal sequences where the outer break is persistent. It can be noted that during these periods, in both cases, the abundances in the outer region are similar and the position of the break mildly evolves to larger radii. We estimated a surviving timescale for this feature that ranges within the range 1 − 1.5 Gyr (see Fig. 8). The shape of the metallicity profile is also affected by galaxy interactions. As discussed in Sect. 4.1.3, close interactions with gas-rich satellites cause the direct accretion of gas from the companion. We found that, for LG1-4469, the sequence of outer rises shown in Fig. 13 is disturbed due to the interaction with a gas rich satellite (lookback time ∼4 Gyr; see Fig. 8 and B.1), whose accreted material (Fig. 7) dilutes the metallicity in the outskirts and the profile recovers the monotonic decreasing behavior. On the other hand, during the sequence shown in Fig. 13, LG1-4337 has a merger event and interactions with small satellites, most of them with small or null gas fractions. The accreted material from these companions is not enough to dilute the metallicity in the outer region, and the outer rise persists. It is also likely that the torques exerted by these satellites open up the spiral arms and favor gas mixing (Sillero et al. 2017). The cold counterpart in the outflowing accreted gas seen in lower panels of Fig. 12 supports this scenario.

Fig. 13. Similar to Fig. 10 but illustrating the evolution of two examples of outer rise metallicity profiles in LG1-4469 (upper panels) and LG1-4337 (lower panels).

5.4. Outer drops

In the target galaxies most outer drops occur at high redshift, where the variety and complexity of physical processes involved in disk assembly make their study particularly challenging. These profiles are typically short-lived. In Fig. 14, we present a sequence illustrating an outer drop that persists during a merger. For this galaxy, initially, the metallicity profile does not exhibit an outer drop, characterized instead by a flat slope in the outer region, with very low metallicities. Intense star formation raises the metallicity in the mid-region while the outer region accretes low-metallicity material, generating an outer drop, a process sustained for approximately 0.3 Gyr.

Fig. 14. Similar to Fig. 10 but illustrating the evolution of outer drop metallicity profiles identified in LG1-4337.

The merger induces a strong, centrally concentrated starburst (see Fig. C.1), ejecting most of the gas in the mid-region (fourth column). Despite the dispersion of ejected gas and the inversion of the gradient in the mid-region, the continuous accretion of metal-poor material maintains the outer drop (fourth and fifth columns). The metallicity map in the leftmost column of Fig. 11 displays the actual shape and pristine gas content of the filamentary structure involved in this sequence.

Form the upper panels of Fig. 12 we can see the characteristics of accreted material during a period when an outer drop is detected. In this case, most of the gas comes from cold flows with low metallicity as expected. We acknowledge that these are two examples which have given us a clear visualization of different processes acting together to modulate the distribution of chemical elements in the ISM. As mentioned before, these results will be used to define a strategy to perform a statical study in a forthcoming paper.

6. Conclusions

We studied the azimuthally averaged metallicity profiles of the SF gas in 45 central galaxies of the CIELO suite of simulations (Tissera et al. 2025). The sample of galaxies covers a stellar masses mass range within the M = [10^8.5, 10^10.5] at z = 0, representing sub-Milky Way galaxies in the field, filamentary environments and LGs.

A new algorithm, DB-A, was developed to self-consistently fit multiple power laws to the metallicity profiles. DB-A selects the best combination of power laws by minimizing chi-squared values, which allows for a flexible assessment of metallicity gradients within various regions of galaxies. We applied this algorithm to CIELO galaxies at z = 0 to study the global dependence of the gradients on galaxy properties. Additionally, the evolution of the metallicity gradients of two galaxies from a simulated LG analog were analyzed as case studies to individualize the physical processes responsible of breaking the metallicity profiles. We stress that these simulations do not include AGN feedback. Hence, we restricted the analysis to sub-Milky Way mass-size galaxies.

At z = 0, the analysis revealed a diversity in shapes for the metallicity profiles, whose properties can be summarized as follows:

• Approximately 38% of the CIELO galaxies exhibit metallicity profiles that can be well represented by a linear negative gradient within ≈[0.5, 2.5] r₅₀. The rest displayed broken profiles, with 40% having outer breaks, 13% having inner breaks, and 9% displaying double breaks in the metallicity profiles.

• Linear negative metallicity profiles are consistent with the inside-out formation model, where central regions are enriched earlier and for a longer time. However, galaxies are affected by the environment and internal processes that can mix gas and stars, shaping the slopes and redistributing chemical elements in the gas and stellar populations, modulating the metallicity profiles, and in some cases inducing breaks in these profiles.

• The mid-region gradients (${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$ ) obtained by the DB-A method are comparable to the gradients obtained with a traditional single linear fit within the range [0.5, 1.5] r₅₀.

• The ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$, and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ agree with observational values, and their relation with the stellar mass and the SFR are also consistent with observations. Smaller galaxies with lower masses tend to have flatter metallicity gradients, likely due to the influence of supernova feedback.

• At z = 0, the median galactocentric position of the inner breaks is at ∼1.2 r₅₀, and that of the outer breaks is at ∼2.5 r₅₀ for the CIELO galaxies. The outer break radii agree well with available observations, while the inner break positions tend to be larger.

The analysis of the two case study galaxies revealed that the metallicity gradient undergoes substantial variability at high redshift, transitioning to a more gradual evolution at lower redshifts. Periods of greater variability in the metallicity gradient can be linked to satellite interactions and merger events. Our study identifies inner breaks in the metallicity profiles, which may be associated with gas inflows, star formation, and energetic feedback, while outer breaks are linked to gas accretion patterns.

Inner breaks are identified in 26% of the metallicity profiles. Inner rises, which have a steep negative slope in the inner region, are observed in 72% of the inner break profiles. Inner drops, where the central metallicity is diluted, are found in 28% of the inner break profiles. The existence and characteristics of inner rises and inner drops in the two case studies can be explained as follows:

• Inner rises are often caused by the efficient triggering of star formation after gas inflows toward the center from the ISM at lower redshifts or from unbound gas at higher redshifts take place. As a consequence, the central metallicity increases.

• Inner drops can be linked to strong feedback, which causes metal-rich gas to be ejected from the center of the galaxy, evidenced by the lack of gas in the central region and disrupted gas distributions. The same physical processes can produce inverted metallicity gradients under similar conditions. At least, for the two histories of the evolutionary histories analysed in this work, inner drops seem to have shorter lifetimes than inner rises.

Outer breaks are found in 43% of the metallicity profiles, characterized by a change in slope in the outer regions of the galaxy. Outer rises, where the metallicity in the outskirts is higher, are found in 68% of outer break profiles. Outer drops, where the metallicity decreases more steeply in the outer region, occur in the remaining 32%, and are more common at high redshift. The outer breaks can be explained as follows:

• Outer rises can be formed by the re-accretion of enriched material through galactic fountains or through continuous and extended star formation in the outer regions, after a break has moved from inner regions. Enhanced gas mixing due to the tidal torques exerted in galaxy-galaxy interactions can also contribute to outer rises.

• Outer drops are often caused by the accretion of metal-poor gas from cold flows, especially at high redshift. Mergers and interactions can also cause outer drops by enhancing star formation in the mid-region and direct accretion of low-metallicity gas in the outer region.

Double breaks exhibit both inner and outer breaks. They represent the interplay and combination of the processes that cause the individual breaks in each case.

Our findings underscore the complex interplay between internal and external processes shaping the metallicity profiles of galaxies. Our detailed analysis highlights the role of gas inflows, feedback, and accretion in sculpting the distinctive features observed in metallicity gradients over cosmic time. Our results are consistent with previous simulations studying mergers, interactions, the CGM, and cold flows. However, our work advances the field by isolating the specific contributions of these processes during different stages of galaxy formation and evolution in a cosmological context, and, hence, underscore their potential as diagnostic tools for tracing the effects of various physical processes. Future research will extend the analysis to a larger sample to quantify and further understand the relative impact of each mechanism in shaping galaxy metallicity profiles.

Acknowledgments

We thank the anonymous referee for his/her constructive report which helped to improve this paper. BTC gratefully acknowledges funding by ANID (Beca Doctorado Nacional, Folio 21232155). PBT acknowledges partial funding by Fondecyt-ANID 1240465/2024 and Núcleo Milenio ERIS NCN2021_017. This project has received funding from the European Union Horizon 2020 Research and Innovation Programme under Marie Skłodowska-Curie Actions (MSCA) grant agreement No. 101086388-LACEGAL. We acknowledge partial support by ANID BASAL project FB210003. This project used the Ladgerda Cluster (Fondecyt 1200703/2020 hosted at the Institute for Astrophysics, Chile), the NLHPC (Centro de Modelamiento Matemático, Chile) and the Barcelona Supercomputer Center (Spain). This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. ES acknowledges funding by Fondecyt-ANID Postdoctoral 2024 Project N°3240644 and thanks the N°cleo Milenio ERIS. NP acknowledges support from Préstamo BID PICT Raices 2023 N°0002. RDT thanks the Ministerio de Ciencia e Innovación (Spain) for financial support under project grant PID2021-122603NB-C21. This project has made use of: Astropy (Astropy Collaboration 2013, 2018), Matplotlib (Hunter 2007) and Py-SPHViewer (Benitez-Llambay 2015). We acknowledge the developers for their work.

References

Appendix A: Physical properties of CIELO galaxies at z = 0.

Appendix B: Merger histories

In Sections 4.1.1 and 4.1.2, we discussed how interactions between the central galaxies and their satellites or mergers were traced back in time using their merger trees to identify key events contributing to their assembly histories. For illustration purposes we selected specific periods of either mergers or interactions that clearly highlight the effect of these physical processes on the metallicity distribution. Figure B.1 provide a detailed view of these interactions, with the upper panels depicting surviving satellites and the lower panels showing mergers. These interactions are tracked through the physical separation of their centers of mass, along with the evolution of their stellar mass ratios and gas fractions at various stages. The most significant events, based on proximity, stellar mass ratios, gas fractions, and coupled interactions, are highlighted in green and blue shaded regions.

The assembly histories of the two galaxies differ markedly. Wet mergers are frequent at high redshift for both, but LG1-4337, the more massive galaxy, experiences a greater number of massive mergers and completes its assembly earlier. At lower redshift, merger galaxies are typically 25% the mass of the central galaxy, with higher gas fractions observed in LG1-4469. This trend also applies to their surviving satellites. We note that LG1-4469 is undergoing a major interaction with a gas-rich companion at z = 0. Coupled interactions, in both merging galaxies and surviving satellites, are common in the evolution of both systems (see, e.g., Gámez-Marín et al. 2024).

Fig. B.1. Time evolution of the physical separation between the center of mass of the central galaxy (left panel: LG1-4469; right panel: LG1-4337) and its surviving satellites (upper panel) and merged galaxies (lower panel). Each continuous line represents a different satellite or merged galaxy, with line thickness proportional to their stellar mass ratio. The line color indicates the gas fraction of the satellite/merged galaxy at different stages of the interaction. Key events are highlighted with green and blue shaded regions, consistent with those shown in Fig. 8.

Appendix C: Normalized $\dot{\mathcal{M}}$ and star formation rate

In Fig. 6 we show the mass of accreted gas ($\dot{\mathcal{M}}$) and stars formed (SFR) per unit time. However, the inner, mid and outer regions will naturally have different sizes, and they will evolve in time as well. In Fig. C.1 we display the ($\dot{\mathcal{M}}$) and SFR normalized by the projected area of each region. The normalization reveals that the amount of accreted gas and stars formed per unit area is highest in the inner region, even though it tends to be the smallest in size. As discussed in Section 4.1.3, the normalized values also show a strong correlation between SFR and accreted gas in the inner region, while the mid and outer regions show an offset between these quantities.

Fig. C.1. Same as Fig. 6 but normalized by the projected area of each region.

All Tables

All Figures

Fig. 1. Examples of the variety of metallicity profiles in the CIELO galaxies. From left to right: (i) Face-on and edge-on (top and bottom subpanels) gas density projections in a large cubic volume of 200 kpc-side centered at a central galaxy. (ii) Face-on and edge-on gas density projections in a cubic volume of 50 kpc-side located at the center of mass of the galaxy. (iii) Face-on and edge-on projected metallicity maps in the 50 kpc-side volume. Arbitrary density contours are included in (ii) and (iii) to illustrate the morphology. (iv) Gas-phase metallicity profile of the central galaxy (background heat map). For illustrative purposes, median values within 0.5 kpc bins are shown as filled squares, with error bars indicating the standard deviation in each bin. The DB-A fit adopted in this work (see Sect. 2.2) is shown as a black solid line. For comparison, a linear fit within [0.5–1.5] r50 is also displayed (red dashed line). Vertical lines are included to indicate the gravitational softening ϵ (gray, dashed line), inner break radii (teal solid line) and/or outer break radii (orange solid line) when appropriate.

In the text

Fig. 2. Comparison between the metallicity gradients (filled black circles) estimated by the DB-A fit, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$, and by the Theil-Sen method, ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$, for galaxies exhibiting, from top to bottom, (i) linear profiles, (ii) inner break profiles, (iii) outer break profiles, and (iv) double break profiles. When corresponds, inner (teal symbols) and outer (orange symbols) metallicity gradients are also shown. Vertical lines are included as a visual aid to identify gradients belonging to the same galaxy.

In the text

Fig. 3. Stellar mass-, sSFR-, and size-metallicity gradient relation for the 45 star-forming central galaxies (left, middle, and right panels, respectively). In the upper panels, the relations were constructed by using only ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Mid}}$. The colors of solid symbols encode the type of profile that was best fit to the whole metallicity distribution, inner break (IB; teal triangles), outer break (OB; orange triangles), and doubly broken (DB; gray hexagons). In the lower panels are the same relations but constructed using ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{TS}}$ (red circles). Observational data from Martin & Roy (1995), Ho et al. (2015), and SM18 (MUSE); Franchetto et al. (2021) (MaNGA+MUSE); Chen et al. (2023) (TYPHOON); Cardoso et al. (2025) (CALIFA); and MaNGA DR17 are also displayed in a similar fashion. The black dashed line indicates ∇(O/H) = 0 for reference purposes only. The two main galaxies LG1-4469 and LG1-4337 of the LG analog are highlighted, as they are later analyzed in detail as case studies.

In the text

	Fig. 4. Similar to Fig. 3 but featuring the ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Inner}}$ and ${\mathrm{\nabla }}_{\text{(O/H)}}^{\mathrm{Outer}}$ in broken metallicity profiles. Observational data from Martin & Roy (1995), SM18 (MUSE), Chen et al. (2023) (TYPHOON) and Cardoso et al. (2025) (CALIFA) are included for comparison.
In the text

Fig. 5. Inner (teal symbols) and outer (orange symbols) break radii as a function of (i) galaxy stellar mass (left panel) and (ii) sSFR (right panel). The CIELO galaxies are displayed as filled symbols, while empty symbols represent observations from Martin & Roy (1995), SM18, Chen et al. (2023) and Cardoso et al. (2025) as indicated in the legend. The SCC was calculated to evaluate possible correlations between galaxy properties and inner and outer break radii. The upper part of each panel shows the SCC computed considering: CIELO galaxies (top rows) and observations (bottom rows). The p-value is shown in parenthesis in each case.

In the text

Fig. 6. Time evolution of the SFR (shaded areas) and the total accretion rates (̇ℳ; solid lines) per region (inner, left column; mid, central column; outer, right column) for each galaxy (LG1-4469, top row; LG1-4337, bottom row). The epochs where inner or outer breaks are identified are denoted by using teal or orange symbols, respectively (see Section 4 for a detailed explanation on the definition of the different regions). Figure C.1 shows the same quantities normalized by the area of each region.

In the text

	Fig. 7. Fractions of the total accreted gas (̇ℳ; see Fig. 6) in the inner (left column), middle (central column) and outer (right column) regions for each galaxy (LG1-4469: top row; LG1-4337; bottom row). The gas accreted fractions correspond to inflowing ISM (pink areas), outflowing ISM (red areas), the CGM (yellow areas), interacting/merging satellite (green areas) or unbound gas (blue areas).
In the text

	Fig. 8. Evolution of the metallicity gradients as a function of lookback time (and redshift) for the targeted galaxies LG1-4469 (upper panel) and LG1-4337 (lower panel). The metallicity gradients in the inner (teal symbols), middle (black solid lines), and outer (orange symbols) regions are displayed. The shaded regions depict the mergers and interaction events (blue and green, respectively).
In the text

Fig. 9. Diagram of δr − δz illustrating the motion of the gas accreted into the inner region of LG1-4337 between lookback time ∼4.8 Gyr and ∼3 Gyr. During this period, the galaxy exhibits an inner break of its metallicity profile. Left panel: KDE distribution of the accreted particles. The colored numbers displayed at the top represent the fractional contribution of each gas accretion source, consistent with the color code used in Fig. 7 and the legend. The numbers below indicate the fraction of particles with δr < 0 (outflowing) and δr > 0 (inflowing). Central panel: Distribution of accreted particles, color coded by Δ(O/H), which represents the difference between the metallicity of the accreted gas and that of the respective region. Negative values indicate that the accreted gas is less enriched than the preexisting gas in the region. The median (med) value is computed over the entire time interval, while the maximum (max) value corresponds to the largest difference observed within that period. Right panel: Distribution of accreted gas particles colored by their temperature.

In the text

Fig. 10. Temporal sequence showing the evolution of the inner rise metallicity profiles for LG1-4469 (first and second rows) and LG1-4337 (third and fourth rows). The labeled images depict gas density projections at different snapshots, while the corresponding metallicity profiles are displayed below, following the same format as in Figure 1. We note that the timescales differ between the two galaxies.

In the text

Fig. 11. Examples of morphologically disrupted galaxies showing evidence of central gas ejection driven by strong energetic feedback. The top row displays face-on gas density projections; the middle row shows edge-on metallicity maps with streamlines indicating gas velocity directions; and the bottom row presents metallicity profiles, formatted as in Figs.1 and 10. The corresponding galaxy and lookback time/redshift are labeled in the upper panels.

In the text

Fig. 12. Diagrams of δr − δz showing the motion of the gas accreted into the outer region of LG1-4337 from different sources (similar to Fig. 9) between lookback time 11.60 Gyr and 10.1 Gyr (upper panels) and between lookback time 6.5 Gyr and 4.9 Gyr (lower panels). The diagrams in the upper panels illustrate the accretion of a high-z galaxy that exhibits an outer drop, with a dominant contribution from metal-poor unbound gas (62% of the accreted gas). The diagrams in the lower panels illustrate a low-z galaxy with an outer rise, with the clear indication recycled metal-rich gas from the CGM being the main source of accretion (56% of the accreted gas).

In the text

	Fig. 13. Similar to Fig. 10 but illustrating the evolution of two examples of outer rise metallicity profiles in LG1-4469 (upper panels) and LG1-4337 (lower panels).
In the text

	Fig. 14. Similar to Fig. 10 but illustrating the evolution of outer drop metallicity profiles identified in LG1-4337.
In the text

Fig. B.1. Time evolution of the physical separation between the center of mass of the central galaxy (left panel: LG1-4469; right panel: LG1-4337) and its surviving satellites (upper panel) and merged galaxies (lower panel). Each continuous line represents a different satellite or merged galaxy, with line thickness proportional to their stellar mass ratio. The line color indicates the gas fraction of the satellite/merged galaxy at different stages of the interaction. Key events are highlighted with green and blue shaded regions, consistent with those shown in Fig. 8.

In the text

	Fig. C.1. Same as Fig. 6 but normalized by the projected area of each region.
In the text