Numerical models for the dust in RCW 120

Context. The interstellar bubble RCW 120 seen around a type O runaway star is driven by the stellar wind and the ionising radiation emitted by the star. The boundary between the stellar wind and interstellar medium (ISM) is associated with the arc-shaped mid-infrared dust emission around the star within the HII region. Aims. We aim to investigate the arc-shaped bow shock in RCW 120 by means of numerical simulations, including the radiation, dust, HII region, and wind bubble. Methods. We performed 3D radiation-hydrodynamic simulations including dust using the GUACHO code. Our model includes a detailed treatment of dust grains in the ISM and takes into account the drag forces between dust and gas and the effect of radiation pres- sure on the gas and dust. The dust is treated as a pressureless gas component. The simulation uses typical properties of RCW 120. We analyse ﬁve simulations to deduce the effect of the ionising radiation and dust on both the emission intensity and the shape of the shock. Results. The interaction of the wind and the ionising radiation from a runaway star with the ISM forms an arc-shaped bow shock where the dust from the ISM accumulates in front of the moving star. Moreover, the dust forms a second small arc-shaped structure within the rareﬁed region at the back of the star inside the bubble. In order to obtain the decoupling between the gas and the dust, it is necessary to include the radiation-hydrodynamic equations together with the dust and the stellar motion. In this work all these elements are considered together, and we show that the decoupling between gas and dust obtained in the simulation is in agreement with the morphology of the infrared observations of RCW120.


Introduction
The Infrared Astronomical Satellite (IRAS, Neugebauer et al. 1984) all-sky survey detected many extended arc-like structures associated with OB-runaway stars and paved the way for further studies to reveal the bow shocks that are present around these powerful stellar winds (van Buren et al. 1995;Peri et al. 2012). Furthermore, the mid-infrared (MIR) observations from Spitzer-MIPSGAL (Carey et al. 2009) of interstellar bubbles along the Galactic plane show that the interiors of HII regions contain dust, and some of them present 24 µm dust emission near to the central ionising source Kendrew et al. 2012;Simpson et al. 2012). The dust grains are expected to evacuate the interiors of the hot bubbles as a result of the mechanical energy and radiation injected by the massive stars. However, recently, infrared observations show a significant amount of dust in the inner part of the bubble, particularly sitting near to the stellar feedback source. Indeed, high-resolution infrared observations show that HII bubbles contain a significant amount of dust in their interiors Martins et al. 2010;Anderson et al. 2012). Everett & Churchwell (2010) explored the possibility that the evaporation of small dense cloudlets could resupply the hot gas with a new generation of dust grains. While this mechanism can explain the presence of dust inside HII bubbles, it is not capable of producing the morphology of the MIR radiation: an arc-shaped and a peaking close to the ionising source. (Ochsendorf et al. 2014, see also van Buren & McCray 1988) proposed that dust grains are expelled from ionised stars by radiation pressure and Pavlyuchenkov et al. (2013) performed 1D numerical calculations on the dust in the boundary between the stellar wind and the interstellar medium (ISM), and they concluded that dust deviation by itself is not a viable mechanism to adequately explain the arcs observed in 24 µm IR emission in HII regions. One of the more studied interstellar bubbles with a 24 µm arc is RCW 120. The RCW 120 bubble is located at 1.3 kpc (Zavagno et al. 2007) at l = 348. • 3, b = +0. • 5 (α J200 = 17 h 09 m 0 s , δ J200 = −38 • 24 (Rodgers et al. 1960) with a radius of 1.9 pc (Anderson et al. 2015). Its ionising star is CD-38 • 11636 or LSS 3959, an O8 type star (Georgelin & Georgelin 1970;Avedisova & Kondratenko 1984;Russeil 2003;Zavagno et al. 2007) located at α J200 = 17 h 12 m 20.6 s , δ J200 = −38 • 29 26 . This star is found in a Galactic HII region embedded in a molecular cloud with numerical density of n a ∼ 3000 cm −3 (Zavagno et al. 2007). The age of the massive star is <3 Myr, and the dynamical age of the shell is ∼0.2 Myr (Arthur et al. 2011). More recently, Torii et al. (2015) found that the dynamical age is ∼0.4 Myr. Ochsendorf et al. (2014) proposed a pressure and density gradient in the ISM around the ionising source of RCW 120, in order to explain the asymmetry of the HII region and the MIR arcs. However the diameter of the RCW 120 is around 3.8 pc (Anderson et al. 2015) and a pressure or density gradient must be related to a stellar motion (Weaver et al. 1977;Raga 1986;Mac Low et al. 1991;Arthur & Hoare 2006;Mackey et al. 2013;Steggles et al. 2017). Most probably, the stars, particularly the massive stars, are born in motion with respect to their surroundings because they keep the momentum gained during their star formation process, where the turbulence is needed, with velocities of 2−5 km s −1 (e.g., Peters et al. 2010;Dale & Bonnell 2011).
Many authors have modelled bow shocks from runaway O and B stars (e.g. Mac Low et al. 1991;Comeron & Kaper 1998;Meyer et al. 2014) and there have been a few studies of HII regions for v ≥ 20 km s −1 , (Tenorio-Tagle 1979;Raga et al. 1997;Mackey et al. 2013). In both cases a partial shocked shell forms, upstream for the bow shock and in the lateral direction for the HII region, and can be unstable in certain circumstances. In the case of RCW 120, Mackey et al. (2016) performed 2D radiation hydrodynamical numerical simulations of stellar wind bubbles and post-processing calculations to obtain synthetic intensity maps at infrared wavelengths. They found good agreement between their model and the RCW 120 observations. In particular, they obtained an infrared arc with the same shape and size as the arc observed in RCW 120, with a brightness comparable to the observed one. They concluded that their results could be showing the outer edges of an asymmetric stellar wind bubble. On the other hand, the respective effects of the smaller and larger dust grains has been studied by van Marle et al. (2011). These latter authors modelled the gas and dust dynamics using high-resolution numerical simulations of a fast-moving red supergiant and concluded that the smaller sized grains follow the motion of the gas and the larger grains are weakly coupled to the gas. Therefore, they proposed that the small dust grains dominate the emission in the infrared observations. Moreover, De Colle & Raga (2006) and Meyer et al. (2017) showed that a small magnetic field can affect the predicted emission from the post-shock recombination regions by one or two orders of magnitude. However, the work of these latter authors was focused on the optical emission prescriptions and they did not study the infrared emission mainly produced by the dust grains. Indeed, Meyer et al. (2017) made postprocessing calculations of the dust continuum infrared emission and showed that the presence of the magnetic field decreases the infrared emission. Finally, using 3D magneto-hydrodynamical simulations Katushkina et al. (2018) showed that the drag force cannot be neglected in models with ISM densities larger than 10 cm −3 .
More recently, Henney & Arthur (2019a,b,c), studied the formation of IR-emission arcs around luminous stars that move supersonically relative to their surrounding medium. These latter authors found four different regimes depending on the type of coupling between dust grains and gas, and taking into consideration the relative values between the stellar wind and the ultraviolet optical depth of the bow shell. They presented cases in which a radiation-supported dust wave decouples from the gas to form an infrared emission arc outside of any hydrodynamic bow shock, considering the magnetic fields in their numerical models. They considered the drag force because of the dust and the gas interactions assuming that the fraction of the dust is zero in the wind.
In the present study, we modelled the RCW 120 physical properties, performing a series of numerical simulations that solve the radiation-hydrodynamics and dust equations in a fully 3D fixed Cartesian mesh. This paper is organized as follows: the numerical simulation setup is described in Sect. 2, the simulation results are presented in Sect. 3, and finally we conclude in Sect. 4.

Governing equations
We performed 3D simulations using the Guacho 1 code (Esquivel et al. 2009;Esquivel & Raga 2013). The Guacho code uses a Cartesian grid to solve the near-conservation laws governing the gas-dynamics and the neutral hydrogen rate. It includes the radiative transfer at the Lyman limit of ionising photons emitted from an isotropic point source photoionising a neutral hydrogen distribution. Also, it uses the ray-tracing scheme where the diffuse radiation is considered only in the B recombination case. The code is described in detail in Esquivel & Raga (2013) and in Schneiter et al. (2016). In this project we incorporate the dustgrain dynamics coupled to gas dynamics and the momentum feedback of radiation on the dust.
The radiative cooling is included as described by Raga & Reipurth (2004) for the atomic gas, and for lower temperatures we have included the parametric molecular cooling function presented by Kosiński & Hanasz (2007) and Rivera-Ortiz et al. (2019). The complete set of equations is as follows: where, ρ, u, p, and E are the density, the velocity, the thermal pressure and the total energy density of the gas, respectively. G rad and L rad are the energy gains and losses due to radiation. We assume a standard ideal gas law for closure, and therefore the gas total energy density and the gas thermal pressure are related by E = ρu 2 /2 + P/(γ − 1), where the ratio between specific capacities is γ = 5/3. The factors ρ d and u d are the density and the velocity of the dust. The drag force f d is given by a combination of Epstein's drag law for the subsonic regime and Stokes' drag law for the supersonic regime (see, van Marle et al. 2011), where c s = γ p ρ is the gas sound speed. The test of the hydrodynamics equation coupled with dust using the drag forces in Eq. (6) shows the same results as those obtained by van Marle et al. (2011). Finally, the radiation pressure gradient force due to photoionisation is and the radiation pressure gradient force due to grain dust absorption is where φ H is the Lyman α photoioinzation rate, χ H is the hydrogen ionization potential, σ d is the grain dust cross-section (the geometrical cross section), F * is the photons flux, and c is the light speed. More details about the radiative force are given in Rodríguez-Ramírez et al. (2016).
In the ray tracing scheme, the photon package emitted by an isotropic point source decays by a factor of e −(∆τ H +∆τ d ) when travelling through neutral gas and dust. Along a path-length ∆l, the opacities induced by neutral gas and dust are ∆τ H = a 0 n HI ∆l and ∆τ d = a d n d ∆l, respectively. Here, a 0 is the photoionisation cross-section constant, a d is the geometrical dust cross-section, and n HI and n d are the HI and dust number density. The photoionisation rate φ HI is computed by equating the ionising photon rate, S, to the ionization per unit time in each cell, that is, S = n HI φ HI dV.
Simulation parameters and setup. The coupled gas-dust and radiation pressure equations are advanced using the HLL scheme from Harten et al. (1983) with the second-order spatial re-construction minmod slope limiter coupled with the secondorder time step.
The ISM has a uniform number density of n a = 3000 cm −3 , temperature T 0 = 100 K, and is moving in the z direction with a velocity of 5 km s −1 (to emulate the stellar motion of RCW 120). Our numerical models considered appropriate values of the stellar wind velocity (2300 km s −1 ) and a mass injection rate (2.7 × 10 −7 M yr −1 ) for one single O8V star (Sternberg et al. 2003). The steady spherical stellar wind was imposed in a sphere of radius 0.09 pc at the simulation grid centre.
We used a photon ionisation rate of 1.6 × 10 48 phot s −1 (see Zavagno et al. 2007) and the photons are injected in the surface of a sphere with a radius of 1.5 times the wind source radius.
We assume the dust to be composed of silicates, for which the internal particle density is ρ grain = 3.3 g cm −3 (Draine & Lee 1984). Concerning the dust shape, we assume spherical grains with a radius of a dust . The simulations are performed with two dust radii, a dust = 0.1 and a dust = 1 µm. The dust to gas fraction is χ = ρ d ρ 0 · ρ grain * V grain , where V grain is the volume of the grains. The dust is only included into the uniform ISM.
We computed simulations of RCW 120, varying the dust fraction and size, the momentum of the wind, the number of ionising photons, and the density in the ISM (see Table 1

The gas and dust decoupling
To understand the stellar wind and the photon pressure effect in the gas and dust dynamics and the dust-gas interaction, we  performed five numerical models (Table 1). A cut in the Z−X plane of the model M3 is shown in Fig. 1. As we can see, two shock regions are formed: one strong shock in front of the runaway star, the nearest one; and the second, a weaker shock behind the runaway star. In the tail of the HII region, the cometary HII region, the shock wave forms a shell of less-dense, shocked ISM. In order to analyse the behaviour of these two regions, we show in Fig. 2 a 1D cut along two directions: a cut perpendicular to the moving star direction, the X-axis, and a cut in the moving direction, the Z-axis; at time 0.4 Myr (left an right panels of Fig. 2, respectively). It is clear that the two regions mentioned above are observed. For the M1 model, where only the stellar wind is considered without including the photons flux (the diamonds line in Fig. 2), the gas and dust are well coupled and are stacked at the external shell. Inside the bubble within the shell of the shocked ISM, there is a rarefied region with low gas and dust density. However, the gas is hot in this region. In Fig. 3 we show the fraction of ionised gas for cuts on the X and Z axes (left and right panels, respectively). In the cut on the X axis, one can see that the gas is completely ionised within a narrow region where the ram pressure and the over-pressure of the hot bubble reach equilibrium. This bubble of ionised (hot) gas is much larger in the Z direction than in the X-axis (1.2 pc.), namely about 3 pc extending in an opposite direction to the movement of the star, giving the cometary shape in this region. For this model, the ionised region has an ellipsoidal (egg-like) shape. In the other models, M2-M5, the ionising photons are included (Table 1). In those models, the ISM is heated by the ionization front which decreases the shock Mach number resulting in a lower compression in the external shell. At the tail of the cometary HII region only a compression wave forms and the shock remains only at the front and the side of the runaway star (Fig. 2). As a result, at the back of the moving star inside the bubble, the density increases and the temperature decreases. Models M3-M4, with dust in ISM and photoionisation, show a greater difference behind the runaway star inside the bubble, where more ISM gas and dust remains. Indeed, the dust decreases the compression wave strength. Moreover, the structure of this rarefied zone changes under the effect of the dust (Fig. 2). There is an increment in the gas density to 0.25 pc from the star and a minimum (local) at 1 pc for M3 and M4 models in comparison with 1.1 pc for the M2 model. Concerning the dust distribution in the M1, M3, M4, and M5 models, in front of the runaway star (Fig. 2) the compression shock is strong and the dust is tightly coupled with the gas. However, inside the bubble at the back of the runaway star in the rarefied region the dust is decoupled from the gas. This decoupling between the dust and gas is more important in model M3 because the dust grains are bigger. In models M3 and M4 there are zones where the dust accumulates. The denser zone is at the forward shock and the second is inside the bubble in the rarefied region (Fig. 2). As a result, there are two regions with a significant dust grain emission: (a) an outer region of swept up ISM and (b) a region more internal and close to the ionising star, at approximately 0.2 pc, which is formed with ISM gas filling the decreased pressure region due to the star motion.

Dust distribution in RCW 120
A visual comparison of our models is described in this section. Figure 4 shows the emission of the Galactic RCW 120 bubble at 24 µm from Spitzer- MIPSGAL Carey et al. (2009). MIPSGAL is a survey of the same region as Spitzer-GLIMPSE surveys (see Churchwell et al. 2009 for more details) using the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) instrument at 24 µm. The MIPSGAL resolution is 6 at 24 µm. This image shows the extended emission of small dust grains at 24 µm in the region of the ionised gas and the emission in the envelope associated with young stellar objects (YSOs). As one can see, the emission at 24 µm of dust is mainly localised in two places: (a) in an external and thin shell in the front part of the cometary HII region, and (b) near the ionising source in a central region of the bubble that has a significant dust fraction. Figure 1 shows a cut on the X-Y axis for the M3 model; the grey map is the gas density and the dust density contours. The blue contours correspond to the normalized density of dust larger or equal to 0.8 (the outer shell) and the contour corresponds to densities lower than 0.8 (the inner dust structure). As we can see in the Fig. 1, the HII region has an ellipsoidal (egglike) structure and does not look like a spherical shape as shown in the observation in Fig. 4. In addition, although the distribution of the dust is in the external shell, because it was dragged by the leading shock, it is also seen that the internal structure of dust is behind the position of the source (marked with an asterisk in Fig. 1). Marsh & Whitworth (2019) determined the shape of the external shell of RCW 120 and concluded that it is a circular, and not elongated, ring with an internal radius of 1.66 pc and a depth of 0.33 pc. The elongation presented in the current numerical models with a star motion cannot explain the shape of the RCW 120. Torii et al. (2015) studied the proper motion in RCW 120 and found a problem with the determination of the molecular shell velocity with respect to the observer. It is therefore likely that the object has a projection angle that does not allow us to see an elongated shape in this HII region. Finally, in order to compare our model with observations, we built column density maps, N(i, j) columnar = n(i, j) * dR, where a rotation has been included. Figure 5 shows the column density for the dust for a rotating angle of 72 degrees of the projection plane X-Z. If the line of vision is perpendicular to the axis of movement, the external shell has the shape of a ring, with the ionising star at the centre of it. However, the rotation of 72 degrees that we have used results in an external shell with a spherical section in the direction of the movement of the star. Also, the internal region with dust now coincides with the position of the star.

Conclusions
We present several numerical models for the evolution of the HII region, RCW 120. Our models consider the photoionisation effect, the wind injected by the central star, and the effect of dust on the dynamical evolution of this object and the stellar motion. With respect to the M1 model, (without photoionisation), the gas and dust do not seem to uncouple, since they are dragged by the leading shock and both of them are stacked in an external shell. We performed a numerical simulation for a weak wind model (where the stellar wind transfers a momentum reduced by a factor of 100 compared with the other models). This model, M5, produces a bubble with gas and dust and a significant decoupling of the dust and gas is not seen.
Our strong stellar winds models show an internal region with a high density of gas and dust very close to the central star. On the other hand, the distribution of dust and gas inside the bubble is different, with a partial decoupling between them. The model in which the dust grains are smaller maintains the coupling or partially coupling of gas and dust for a longer time despite the interactions with the ionisation and shock fronts. The decoupling is greater in models M3 and M4, with larger dust grains 1 µm. In these models, the dust is mainly concentrated in the external shell and in the central structure near the position of the star. That is, the gas and dust density profiles are found with a very similar distribution as inside the RCW 120 bubble.
We show that when making a rotation of 72 degrees, the morphology results in an external shell with a spherical ring shape. Surprisingly, the internal region with dust now coincides with the position of the star. This coincides with the observations which suggest that there is an angle with respect to the line of sight that presents us from observing an elongated shape in RCW 120.
As far as we are aware, this is the first time that 3D gasdynamic equations, coupled with the photoionisation and the dust together with the stellar motion, have been solved for RCW 120. Therefore, according to our models, all of these are of great importance to understand the distribution of dust in the RCW 120 bubble.
Our results show good agreement with those obtained by Mackey et al. (2016), and are in accordance with the general analysis of Henney & Arthur (2019b) for the case of a dust wave where the dust is decoupled from the gas. Indeed, the distribution of dust is mainly located in a region found immediately after the bow shock, as in the two works mentioned previously.
We highlight the fact that we have not included the magnetic field in our simulations coupled with radiation hydrodynamics and dust, and that this is important in order to compute synthetic maps with optical and infrared emission. The inclusion of the magnetic field in simulations and the effect of this on the emission is left for future work.
Finally, we have shown the importance of strong stellar wind in the decoupling of the gas and dust during the lifetime of the bubble.