Boundary conditions in hydrodynamic simulations of isolated galaxies and their impact on the gas-loss processes

Three-dimensional hydrodynamic simulations are commonly used to study the evolution of the gaseous content in isolated galaxies, besides its connection with galactic star formation histories. Stellar winds, supernova blasts, and black hole feedback are mechanisms usually invoked to drive galactic outflows and decrease the initial galactic gas reservoir. However, any simulation imposes the need of choosing the limits of the simulated volume, which depends, for instance, on the size of the galaxy and the required numerical resolution, besides the available computational capability to perform it. In this work, we discuss the effects of boundary conditions on the evolution of the gas fraction in a small-sized galaxy (tidal radius of about 1 kpc), like classical spheroidal galaxies in the Local Group. We found that open boundaries with sizes smaller than approximately 10 times the characteristic radius of the galactic dark-matter halo become unappropriated for this kind of simulation after about 0.6 Gyr of evolution, since they act as an infinity reservoir of gas due to dark-matter gravity. We also tested two different boundary conditions that avoid gas accretion from numerical frontiers: closed and selective boundary conditions. Our results indicate that the later condition (that uses a velocity threshold criterion to open or close frontiers) is preferable since minimizes the number of reversed shocks due to closed boundaries. Although the strategy of putting computational frontiers as far as possible from the galaxy itself is always desirable, simulations with selective boundary condition can lead to similar results at lower computational costs.


INTRODUCTION
Differential equations are widely used in different astrophysical contexts, from physical phenomena involving our solar system to the large-scale structures in the universe, as clusters and superclusters of galaxies. In particular, fluid-dynamic problems are formulated through differential equations that represent the conservation of mass, momentum, and energy of a fluid with a certain equation of state (e.g., Landau & Lifshitz 1987).
Corresponding author: Anderson Caproni anderson.caproni@cruzeirodosul.edu.br To find a particular solution of any differential equation, it is necessary to provide some initial condition and/or some boundary conditions (BCs; e.g., Arfken & Weber 2005). However, because of high complexity involving hydrodynamic (HD) problems in astrophysical systems, analytical solutions for the temporal behavior of a fluid are rare, leading to applications of numerical methods to solve the HD equations (e.g., Toro 2009). There are several numerical codes dedicated to dealing with astrophysical gas/particle dynamics (e.g., Stone & Norman 1992;Fryxell et al. 2000;Raga et al. 2000;Teyssier 2002;Gammie et al. 2003;Anninos et al. 2005;Springel 2005;Mignone et al. 2007;Bryan et al. 2014), adopting different strategies to numerically evolve arXiv:2302.04825v1 [astro-ph.GA] 9 Feb 2023 gas/particle flows. Distinct BCs are usually available in those codes, which are chosen according to the specific physical situation to be studied.
In this work, we focus on HD simulations planned to study the time evolution of the gas content inside an isolated galaxy under the influence of a dark-matter distribution and supernova feedback (e.g., Silich & Tenorio-Tagle 1998;Mac Low & Ferrara 1999;Fragile et al. 2003;Wada & Venkatesan 2003;Marcolini et al. 2006;Stinson et al. 2007;Revaz et al. 2009;Ruiz et al. 2013;Recchi 2014;Caproni et al. 2015;Emerick et al. 2016;Caproni et al. 2017;Emerick et al. 2020), aiming to verify the influence of the BCs on the gas removal efficiency.
This paper is structured as follows. In Section 2, we describe the three BCs analyzed in this work. Initial setup and the general results from the HD simulations performed in this work are presented in Section 3 and discussed in Section 4. The main conclusions obtained in this work are listed in Section 5.

THE SELECTIVE BOUNDARY CONDITION
Before introducing our selective boundary condition (SBC), it is useful to present the main characteristics of the additional boundary conditions (BCs) used in this work.
Let ρ, P , and v be the mass density, the thermal pressure, and the velocity of a fluid element at a position r measured in a given reference frame. The former three quantities are usually referred to as primitive variables in HD problems. For grid numerical simulations, the region of interest is discretized on computational cells, where the HD equations are evolved in space and time. The region of interest, or simply the computational domain, is enclosed by numerical boundaries. Boundary conditions are implemented numerically by the usage of guard or ghost cells adjacent to the boundaries of the computational domain.
Let us also definen as the unit vector orthogonal to the boundaries of a computational domain, always pointing outwards by convention. In the case of the open boundary condition (OBC), also known as outflow BC, the gradient of any primitive variable across the boundary alongn is set equal to zero (e.g., Mignone et al. 2007). Caproni et al. (2015) and Caproni et al. (2017) adopted closed boundary conditions (CBC) in their HD simulations. It differs from open boundaries only in terms of the values of v at the boundaries: all velocity components were set to zero in Caproni et al. (2015), while the null value was set only for the velocity component parallel ton, v n = v ·n, in Caproni et al. (2017). Those authors adopted such boundary condi-tions to avoid that the frontiers of the computational domains in their simulations behaved as an infinity reservoir of matter due to the dark-matter gravitational potential (see Section 3 for further discussion). A similar BC was also adopted by Fragile et al. (2003), where density and temperature at static boundaries (v = 0) is kept fixed at their initial values.
The SBC (also known as diode BC; e.g., Fryxell et al. 2000;Zingale et al. 2002) is a variant of the CBC adopted in Caproni et al. (2017), in the sense that if the fluid element that reaches the boundary is moving outwards, as well as it having a speed higher than a predefined threshold value, v th , the CBC is switched to OBC at that location. In other words, the selective boundaries allow those fluids that are moving fast enough to leave the computational domain; otherwise SBC blocks their passage, keeping them inside the domain. Thus, the SBC can be defined as follows where SBC is the boundary condition to be used for a given position at the boundaries in a given time step, and |v ·n| (= |v n |) is the the absolute value of v n in a given cell adjacent to the boundary.

Initial setup
Aiming to test the impact of the boundaries on the evolution of the gas content inside an isolated (dwarf spheroidal) galaxy, we decided to use in our simulations a similar initial gas configuration found in Caproni et al. (2017). In a few words, an isothermal gas is put in hydrostatic equilibrium with a cored, static dark-matter gravitational potential (e.g., Equation 6 in Caproni et al. 2017), so that its density distribution is peaked at the center of the gravitational potential well, decreasing radially as the galactocentric distance increases.
Adopting a dwarf galaxy as a proxy for an isolated galaxy in our simulations avoids working with large computational domains, since dwarf galaxies are relatively small in size, with a tidal radius roughly below of some thousands of parsecs (e.g., Mateo 1998). Consequently, it helps to conduct high numerical resolution experiments without a large number of computational cells, decreasing substantially the involved execution times.
We show in Table 1 the main physical characteristics of the isolated galaxy used in our numerical simulations. These values are compatible with those inferred for the classical dwarf spheroidal galaxy Ursa Minor. c At the center of the galaxy.

Perturbing the galactic gas: types Ia and II supernovae feedback
The initial gas distribution is perturbed by supernova (SN) blasts in our simulations. We followed basically Caproni et al. (2017) for the SN feedback recipe, even though the new version of our code used in this work follows independently types Ia and II supernovae 1 . In a few words, the rates of types Ia and II SNe in our simulations were constrained by the chemical evolution model for Ursa Minor galaxy (Lanfranchi & Matteucci 2004, a typical classical dwarf spheroidal galaxy in the Local Group. The imposed types Ia and II SNe rates are strictly respected during the whole of the simulations, telling to the code when an SN event must occur. On the other hand, where a SN event must take place depends on its type: denser regions are more prone to be selected for harboring a type II SN blast, while type Ia SNe are distributed randomly inside the galaxy. Independent of the type of supernova, an internal energy of 10 51 erg is added into the computational cell elected as an SN site. The SN feedback injects momentum into the interstellar medium, producing a net motion of the gas that is directed outward the galaxy. These galactic winds drive the gas losses that the simulated galaxy will experience as the time evolves. A portion of this galactic wind can reach the boundaries after a given interval, so that we 1 Further details of this new approach will be provided in a future paper in preparation (Lanfranchi et al. 2023).
must be concerned about the influence of the chosen boundaries on it.

Boundary conditions and instantaneous gas-loss rates
All the numerical HD simulations performed in this work made use of the PLUTO code 2 (Mignone et al. 2007) in its version 4.2. The classical hydrodynamic differential equations are evolved in time by a third-order Runge-Kutta algorithm, while the primitive variable reconstruction is done by a piecewise parabolic method (Colella & Woodward 1984). The flux computation among numerical cells was done by the advection upstream splitting method (AUSM+; Liou 1996). We also assumed that gas respects the ideal equation of state, and it is under influence of a cored, dark-matter gravitational potential and cooling processes (see Caproni et al. 2017 for additional details).
We performed 16 HD numerical simulations to study the impact of the boundaries on the gas-loss rates. The main characteristics of these simulations are listed in Table 2. They include two simulations adopting OBC (OBL60N170 and OBL3N100), one simulation with CBC (CBL3N100), and the remaining 13 simulations used to test the behavior of the SBC.

Comparing open, closed and selective boundary condition simulations
Following (Caproni et al. 2015), we estimated the instantaneous gas mass inside a galactocentric radius R gal = 950 pc (compatible with the tidal radius of the Ursa Minor dSph galaxy; e.g., Irwin & Hatzidimitriou 1995), after integrating numerically the mass density distribution obtained in the simulations where M gas is the total gas mass at a time t inside a spherical volume V gal with a radius of R gal . The instantaneous mass fraction of the gas inside R gal , f gas , is calculated from where M gas,0 is the initial gas mass inside R gal (see Table  1). We show in Figure 1 the behavior of f gas for three distinct simulations: OBL3N100, CBL3N100, and V2L3N100 (see Table 2 for further details). They show similar decreasing rates in the gas mass fraction due to SN feedback considering the first 600 Myr of evolution. After this interval, the situation changes dramatically: the OBC induces the breaking of the previous monotonic trend due to the rising of strong inflows of matter, which leads to extremely high (and nonphysical!) masses in comparison with what would be expected for a dwarf galaxy. This issue was already found by Caproni et al. (2015) in their simulations: the OBC acts as an infinite reservoir of matter, which provides gas whenever the pressure equilibrium within the computational domain is broken due to the domain discretization (e.g., Zingale et al. 2002).
On the other hand, the decrease of the amount of gas still remains after 600 Myr for both CBL3N100 and V2L3N100 runs, even though the loss rates lower gradually until they become almost null after ∼1.6 Gyr. There is also a systematic offset between the gas mass fractions from CBL3N100 and V2L3N100 runs, the former presenting a substantially higher value after an elapsed time of 3 Gyr (∼0.40 5 against ∼0.18 for V2L3N100). The reason is that gas flows reaching boundaries with speeds higher than v th = 2 km s −1 are allowed to leave the computational domain, while the CBC retains them, independent of how fast the gas flow is.
The quasi-saturation in the gas removal is detected in both CBL3N100 and V2L3N100 runs. The same result was also found by Caproni et al. (2017) (see their Figure  4). Caproni et al. (2017) also pointed out that reverse shocks generated when the SN-driven gas reaches the computational boundaries could decrease the inferred gas losses due to the gas retention. From a simple analytical calculation using the escape velocity associated with the dark-matter (DM) halo, these authors estimated that the CBC decreased the gas losses by a factor of ∼2.5 after 3 Gyr of evolution.
An alternative way to estimate the influence of the computational frontiers on the results is to put them farther and compare the amount of gas left inside the galaxy after a certain elapsed time. Thus, we run two additional simulations, V64L6N200 and V64L12N200, in which the computational domain was extended, respectively, by a factor of 2 and 4 in relation to the original box size. These two simulations were compared with the V64L3N100 run, in which v th was set to 64 km s −1 (the DM escape velocity used by Caproni et al. (2017) in their analytical calculations). The time evolution of the gas mass fraction inside R gal for these three simulations is shown in Figure 2.
Until ∼500 Myr, these three simulations produced the same instantaneous gas mass fractions, as it can be seen in Figure 2. After that, the discrepancy between V64L3N100 and the other two runs increases monotonically until approximately 1.5 Gyr, when such differences stabilize approximately around a factor of about 2.7. Note this factor is compatible with the value estimated previously by Caproni et al. (2017) from arguments based on the comparison between the escape velocity of the dark-matter halo and the velocities in the adjacent cells at the computational boundaries.
Simulations V64L6N200 and V64L12N200 agree with each other, indicating that a computational domain with a size of 6 kpc (∼6 times the tidal radius of our fiducial galaxy) is enough to minimize boundary effects on the gas losses in similar simulations performed in this Figure 1. Instantaneous gas mass fraction inside a galactocentric radius of 950 pc (tidal radius of the galaxy) for the simulations using OBC (OBL3N100, blue circles), CBC (CBL3N100, green triangles), and SBC with v th = 2 km s −1 (V2L3N100, red squares). These three simulations were made considering a cubic domain of 3 3 kpc 3 . Figure 2. Instantaneous gas mass fraction inside a galactocentric radius of 950 pc (tidal radius of the galaxy) for SBC simulations with different sizes of the computational domain: L = 12 kpc (V64L12N200, blue circles), L = 6 kpc (V64L6N200, green triangles), and L = 3 kpc (V64L3N100, red squares). All these three runs adopt v th = 64 km s −1 . Dashed black line represents the results from an OBC simulation with L = 60 kpc (OBL60N170).
work. Besides, a low-amplitude oscillatory behavior in the curves concerning these two simulations is seen in Figure 2. We have interpreted this as the result of the competition between outward galactic winds driven by SNe and the DM halo's gravity, which tries to push gas back into the galaxy. In other words, these oscillations could be the numerical realization of the "keeping the gas spread and heated" suggested previously in the literature (e.g., Read et al. 2006;Caproni et al. 2017).
Even though Figure 2 indicates some convergence in the gas mass fraction derived from the simulations V64L6N200 and V64L12N200, their domain dimensions are small in comparison to R 200 (∼30 kpc; see Table  1) of the dark-matter halo of our fiducial galaxy. To verify whether these results are indeed representative of the gas losses driven by supernovae, we run an additional OBC simulation, OBL60N170, with a length of 60 kpc (∼2R 200 ) in each Cartesian direction. Its derived instantaneous mass fraction is shown by a dashed black line in Figure 2, indicating a quite similar behavior inferred from the simulations V64L6N200 and V64L12N200. Thus, a domain size of about 6 kpc seems to be enough for HD simulations of isolated galaxies with similar properties listed in Table 1.

The impact of the threshold velocity on the selective boundary condition simulations
As it was discussed in the previous section, the usage of v th = 64 km s −1 in SBC simulation V64L3N100 Figure 3. Instantaneous gas mass fraction inside a galactocentric radius of 950 pc (tidal radius of the galaxy) for SBC simulations with different threshold speeds (from 64 to 1 km s −1 ). Simulation V64L12N200 is also plotted (purple circles). . Instantaneous gas mass fraction inside a galactocentric radius of 950 pc for the SBC simulations V2L3N200 (blue circles), V2L3N100 (green triangles), and V2L3N50 (red squares). All these simulations use the same value for v th (= 2 km s −1 ), but differing in terms of numerical resolution. The solid black line shows the results from the simulation V64L6N250.
increased the amount of gas left after 3 Gyr of evolution by a factor of 2.7 in comparison to the simulations V64L6N200 and V64L12N200, which made use of larger computational domains. A question that arises is whether it is possible to recover the results found in larger box simulations just varying the value of v th in the SBC simulations. Thus, we run seven additional simulations with the same initial setup and resolution of V64L3N100, but decreasing the value of v th mostly by multiples of 2. A comparison among these simulations, as well with the larger box simulation V64L12N200 can be seen in Figure 3.
Again, no apparent differences among all simulations in Figure 3 are noted until approximately 500 Myr of evolution. After this interval, the instantaneous amount of gas left inside R gal decreases systematically as v th is lowered from 64 to 1 km s −1 . The largest differences occur when v th 4 km s −1 , indicating that only a small portion of the gas that is pushed away by the SNe reaches the boundaries with speeds higher than ∼4 km s −1 .
It can be seen in Figure 3 that simulations V2L3N100 and V1.5L3N100 led to gas mass fractions closer to that obtained in the simulation V64L12N200, indicating that the appropriated value of v th must be roughly between 1.5 km s −1 and 2.0 km s −1 for simulations with a cubic domain of 3×3×3 kpc 3 .

Effects of the numerical resolution
Adopting simulation V2L3N100 as a reference, we multiplied (divided) by a factor of 2 the number of computational cells, but keeping all additional parameters fixed, generating simulation V2L3N200 (V2L3N50). It means a change in the numerical resolution from 30 pc cell −1 to 15 pc cell −1 in the case of V2L3N200, while a resolution of 60 pc cell −1 is attained for the simulation V2L3N50. We show in Figure 4 the influence of the numerical resolution on the instantaneous amount of gas left inside the galaxy. No difference in the mass fraction among these three simulations is seen during the first 200 Myr, when V2L3N50 begins to show a higher massloss rate in comparison to the other ones. This trend is inverted after about 1 Gyr and remains so until the end of the simulations.
Concerning simulations V2L3N100 and V2L3N200, there is no significant difference between them until ∼1 Gyr, but after 1.5 Gyr, f gas decreases slowly in V2L3N200, in contrast with simulation V2L3N100 that presents a small-amplitude oscillations around f gas ∼ 0.185. The increment in the numerical resolution from 30 to 15 pc cell −1 is allowed to solve the snowplow transition radius for number densities as low as 1 cm −3 (e.g., Cioffi et al. 1988;Ostriker & McKee 1988), avoiding over cooling issues that weaken the kinetic feedback from supernovae (e.g., Creasey et al. 2011;Simpson et al. 2015;Caproni et al. 2017). Thus, a larger fraction of gas reached speeds higher than the threshold speed of 2 km s −1 in simulation V2L3N200, leaving definitively the computational domain.
At this point, it is interesting to verify whether the monotonic decrease of f gas in V2L3N200 is due to a rather low value of v th in SBC. For this aim, we run an additional simulation, V64L6N250, where we doubled the size of the computational domain but keeping the numerical resolution of 15 pc per cell inside a cubic subdomain of 3 kpc in size (see Table 2 for further details). The behavior of f gas as a function of time is shown in Figure 4 by the black solid line. No monotonic decrease of f gas after 1.5 Gyr is seen but there are smallamplitude variations around f gas ∼ 0.185 instead, as in simulation V2L3N100. It suggests that v th = 2 km s −1 adopted in V2L3N200 is subestimated somehow. Based on the results shown in Figure 3, an increment of about 0.5-1.0 km s −1 in v th might be enough to reconcile simulations V2L3N200 and V64L6N250.
Finally, we can also note that the increment of numerical resolution led to a low amount of gas left inside the galaxy after 3 Gyr of evolution, a factor of ∼2.5 between the lowest and highest numerical simulations in Figure 4 (V2L3N50 and V2L3N200, respectively). This difference is not too big if we consider the usual uncertainties regarding the estimates of the mass in stars, gas and dark matter in galaxies, as well as a relatively poor knowledge concerning the individual efficiencies of the feedback mechanisms to remove gas in those systems.
As it was already mentioned, a small fine tuning in the value of v th can diminish or even eliminate those discrepancies.

DISCUSSION
Our results showed that no influence of the BCs on the instantaneous gas-loss rates is observed until ∼600 Myr in the simulations discussed in Section 3.3. It suggests that noncosmological grid-based HD simulations involving isolated galaxies will not be substantially influenced by the choice of a particular BC if the simulated time is less or of the order of some hundreds of Myr, as it is the case of several previous works involving different types of galaxies (e.g., Mac Low & Ferrara 1999;Fragile et al. 2003;Wada & Venkatesan 2003;Fragile et al. 2004;Melioli & de Gouveia Dal Pino 2015;Emerick et al. 2019Emerick et al. , 2020. For analogous HD simulations but involving longer timescales of evolution, the usage of SBC may be a useful alternative without sacrificing the numerical resolution and/or increasing the computational costs in the case of putting the numerical frontiers far from the galaxy (e.g., Marcolini et al. 2006;Mori & Burkert 2000;Emerick et al. 2016). To provide a sense of the gain in CPU time, we run four additional simulations evolved during 200 Myr in a workstation equipped with 128 2.2 GHz processors. Two of these simulations adopt SBC with v th = 2 km s −1 , while the two complementary ones use OBC. Besides, domain volumes of 3×3×3 and 60×60×60 kpc 3 were built for both SBC and OBC. In the case of the small volume domain (3×3×3 kpc 3 ), 60 computational cells per Cartesian axis were generated, implying a numerical resolution of 50 pc per cell. For the larger volume simulations, we kept the same numerical resolution of 50 pc per cell between -1.5 and 1.5 kpc, but decreasing nonmonotonically this resolution until it reaches the numerical boundaries at -30 and 30 kpc, leading to 102 cells per Cartesian direction. The results can be summarized as follows: • The size of the computational domain is fixes, and the elapsed time to complete a simulation does not depend strongly on the assumed BC: in the case of a domain size of 3 kpc, ∼12.8 and 13.1 hours for SBC and OBC, respectively; for a 60-kpc box, the elapsed times were ∼48.7 and 48.6 hr for SBC and OBC, respectively; • However, a larger domain implies in a substantial longer time for the completion of the simulation. For instance, a larger computational domain with OBC led to a longer execution time by a factor of ∼3.7. Even though this factor is smaller than the ratio between the total number of the cells used in the simulations, (102/60) 3 ∼ 4.9, it shows that larger domains imply higher computational costs that could become prohibitive if high-performance computational resources are not accessible in practice.
Besides the avoidance of the frontiers of the computational domain behaving as an infinity reservoir of matter in simulations with gravity, the SBC can decrease (or even eliminate) the occurrence of reversing flows of mat-ter due to a pure CBC. As any flow colliding with a CBC will have its normal velocity reversed, it induces spurious backflows of matter that could modify the previous gas motions at interacting zones, as well as the physical conditions (density and temperature) of the gas (mainly if the created backflows become strong shocks). Note also that these reversing flows are expected to occur even in the absence of gravity forces.
The choice of a particular BC may also influence on the stability of the initial gas configuration under hydrostatic equilibrium with a gravitational potential. To quantify this effect on grid-based simulations of isolated galaxies, we rerun simulations OBL3N100 and V2L3N100 for 500 Myr turning off the SN feedback during the whole simulation. The results are shown in Figure 5. We note in the case of SBC (upper panels in Figure 5) that the initial gas distribution is well preserved during the whole simulation, with spurious speeds being lower than 0.25 km s −1 (∼50 percent of the cells have speeds lower than 20 m s −1 ). On the other hand, the usage of OBC in a relatively small computational domain (lower panels in Figure 5) induces catastrophic inflows of gas that destroyed the initial spherically symmetric cored distribution of the gas, producing spurious speeds as high as some tens of kilometers per second. To reduce such spurious motions using OBC, larger computational domains are needed.
The magnitude of the spurious accretion also depends on the numerical resolution, as pointed out previously by Zingale et al. (2002). We analyzed the impact of the numerical resolution on the time stability of the initial gas configuration rerunning simulations V2L3N50, V2L3N100, V2L3N200 without SN feedback, including also an extra simulation with a lower numerical resolution in comparison with the previous ones (l = 150 pc cell −1 ). We show in Figure 6 the behavior of the spurious speeds in terms of numerical resolution after 500 Myr of evolution considering SBC. Trends of the increase of the mean and maximum spurious speed with the decrease of the numerical resolution are clearly seen in Figure 6, even though their values are always very small ( 1 km s −1 ) in comparison to the OBC simulation shown in Figure 5. These results suggest that SBC may be also useful in numerical problems somehow involving the hydrostatic equilibrium condition.

CONCLUSIONS
In this work, we studied the influence of the computational frontiers on the gas removal process in (small) galaxies. The option for using an initial configuration compatible with a typical dwarf galaxy (tidal radius of about 1 kpc) is justified by keeping the computational domain as small as possible without sacrificing substantially the numerical resolution, and keeping the computational costs relatively low as well. Three different boundary conditions were employed in this work: open (or outflow), closed, and selective boundary conditions. The 16 hydrodynamic simulations with types Ia and II supernovae feedback performed in this work adopted a cubic domain where the galactic center coincides with the center of the computational box. The majority of these simulations have frontiers put at a galactocentric distance corresponding to ∼ 1.6R gal . Our main results are summarized as follows.
• No difference in the gas mass fraction left inside the galaxy is noted until about 600 Myr of evolution, independent of the three boundary conditions analyzed in this work. It suggests that similar simulations involving short periods of time can adopt open boundary conditions without any loss of integrity of the results; • After 600 Myr of evolution, open boundary conditions for a relatively small computational box (sizes smaller than ∼ 3R gal or about 10 times the characteristic radius of the galactic dark-matter halo) act as an infinity reservoir of gas due to darkmatter gravity whenever the pressure equilibrium within the computational domain is broken due to the domain discretization (e.g., Zingale et al. 2002). In this case, closed or selective boundary conditions are preferable if the increase of the computational edges are somehow unfeasible; • As it was already expected (e.g, Caproni et al. 2015), closed frontiers tend to retain more gas in comparison to the selective boundary condition, impacting on the amount of mass left inside the galaxy: a factor of 2 approximately (see Figure 1); • Concerning the influence of the value of v th used in the selective boundary condition simulations, no difference in f gas is seen until approximately 500 Myr of evolution. It remains true until 3 Gyr for the simulations using v th 8 km s −1 , coinciding with the results from the closed boundary simulation. For the simulations with v th 4 km s −1 , the instantaneous amount of gas left inside the galaxy decreases systematically as v th is lowered; • For v th 1.5 km s −1 , f gas decreases with time, in contrast with SBC simulations with higher v th that present a plateau-like behavior after ∼1.5 Gyr of evolution. Numerical simulations with larger computational domains ( 6R gal ) show similar plateau-like behavior, but showing also a smallamplitude oscillation around f gas ∼ 0.185 possibly produced by the competition between the pull from the dark-matter gravitational potential and the push due to the supernova feedback; • In terms of numerical resolution, our results show no difference in the mass fraction during the first 200 Myr when l is varied from 60 to 15 pc cell −1 . This interval is extended to about 1 Gyr considering simulations with 30 and 15 pc cell −1 only. The monotonic decreasing of f gas seen in V2L3N200 is not present in V64L6N250 with a larger computational domain, indicating that v th = 2 km s −1 adopted in V2L3N200 is subestimated somehow. Based on the results shown in Figure 3, a small increment of about 0.5-1.0 km s −1 in v th might be enough to reconcile simulations V2L3N200 and V64L6N250.
• Although the strategy of putting computational frontiers as far as possible from the galaxy is al-ways desirable, our simulations with a selective boundary condition can lead to similar results but at less expensive demands regarding computational resources.
As a final remark, even though we have analyzed the influence of the boundary conditions over the gas-loss rates using a dwarf spheroidal galaxy, the SBC strategy can be adopted for any type of galaxy or astrophysical system that demands closed numerical frontiers (e.g., see Lanfranchi et al. 2021 for an application involving SBC in the context of intermediate-mass black hole feedback in dwarf spheroidal galaxies).