© ESO 2018

1. Introduction

The finding by Meléndez et al. (2009) that the Sun departs from most solar twins in the solar neigbourhood in that the composition of its surface layers is comparatively rich in volatile and poor in refractory elements, has led to considerable discussion (for a recent review see Adibekyan et al. 2017). Although the finding, a correlation between the chemical abundances for the Sun relative to the twins with the condensation temperatures of the elements in a gas, has been verified by several independent studies, for example Ramírez et al. (2009), Gonzalez et al. (2010), Adibekyan et al. (2014), Nissen (2015), and Spina et al. (2016), its interpretation is still disputed. The early suggestion by Meléndez et al. (2009) that it could be due to the formation of terrestrial planets in the protoplanetary nebula before the latter was dumped onto the Sun, polluting its convection zone by dust-cleansed gas, requires either an unusually long-lived nebula or a very rapid retraction of the solar convection zone to the surface (cf. Gustafsson et al. 2010). As suggested by Önehag et al. (2011) this early retraction might result from the episodic accretion scenario suggested by Baraffe & Chabrier (2010), if applied to the Sun but not to the twins. The alternative – that the effect is due to Galactic chemical evolution and migration of the Sun from galactic regions other than those of the twins – was suggested and partly supported by tendencies found by Adibekyan et al. (2014), Maldonado et al. (2015), Maldonado & Villaver (2016), and Nissen (2015) for some correlation of the abundance effects to depend on stellar ages. Still another possibility was proposed by Gustafsson et al. (2010) that the presolar nebula was cleansed early on from dust by radiation pressure from hot stars in the neighbourhood. There are indications from the isotopic composition (primarily daughters of ${}^{26}$Al and ${}^{41}$Ca) that the Sun was formed in a relatively rich cluster with at least one nearby core-collapse supernova affecting the solar composition (cf. Adams 2010 for a summary, see, however, Fujimoto et al. 2018 for an alternative scenario), and thus with some O- or B-type stars that could have partly cleansed the solar cloud before the solar system was formed. In order to test this hypothesis Önehag et al. (2011 and 2014) carried out detailed analyses of solar twins and solar analogues in the rich old cluster M 67, which is known to have an overall chemical composition close to solar and an age of the order of 4 Gyr. It was found that the cluster stars indeed have an abundance profile more similar to that of the Sun than most solar twins in the neighbourhood, which then seemed to support the assumption that the protocluster and proto-sun were indeed cleansed by radiation from hot stars.

Since the pioneering work by Schalén (1939) and Spitzer (1941) on radiatively accelerated interstellar dust grains, various effects of such grains on interstellar gas have been explored, for example on the dynamics of clouds in the Galaxy (Franco et al. 1991) or on the internal structure of star-forming clouds with H II regions hollowing their centres (Mathews 1967, Cochran & Ostriker 1977). Several relatively detailed numerical models of H II regions with dust-driven flows have been published. In addition to the static models by Draine (2011a), where the dust drift in the region was explored, resulting in piled-up dust in the surrounding outer shell of the region, dynamic simulations have been made in studies such as Arthur et al. (2004), Martínez-González et al. (2014), and Kim et al. (2016). Several studies have been focussed on the structure and evolution of density gradients and shock fronts in these objects (Franco et al. 1990, Shu et al. 2002) and the structure of the photoionization fronts (Diaz-Miller et al. 1998). Skinner & Ostriker (2015) have explored the feedback via dust grains of radiation from massive star clusters embedded in giant molecular clouds.

The issue studied in the present paper is whether the cleansing mechanism proposed could be efficient enough to deplete the dust, not only in the protosolar cloud, but even in the cloud in which M 67, with its initially about 20 000 stars (Hurley et al. 2005), once formed. In this study we use a simple model of a giant molecular cloud with massive stars to set upper limits of the amount of dust-cleansed neutral gas that could be expected. Thus, our aim is not to present a realistic model of this complex physical situation but sooner to obtain constraints on the possibilities of cleansing large gas masses by the mechanism proposed. Our study is also a simple preparatory work for possible more detailed numerical simulations provided that the present results suggest that such studies may be warranted.

Our model represents a star-forming cloud by an initially homogeneous sphere of dense neutral atomic or molecular gas with dust and with a bright light source in its centre. Initially, we follow how the dust is affected by the radiative pressure from the star(s) (Sect. 2.1). The friction force of the dust moving relative to the gas has to be considered, and this friction also moves the gas. However, we see (in Sect. 2.2) that the complex effects of this gas dynamics are to a considerable extent compensated for when the total amount of cleansed gas is calculated. Finally, the motion of the ionization front around the star also has to be considered; this ionization reaches the gas cleansed, presses it outwards, and may stop the possibilities for stars to form easily from that gas (Sect. 2.3). Unlike the model studies referred to above, the present work focusses on the dust in the neutral gas, outside the H II region. In Sect. 3 we draw conclusions concerning how great masses of cleansed gas that could possibly result and be available for forming stars.

2. Model

2.1. Dust front

In interstellar clouds, the mean free path of the molecules is orders of magnitude longer than the size of the dust grains. In this case the equation of motion for a single spherical dust grain, illuminated by a star with luminosity L, can be written (see Draine 2011a)$$\begin{array}{c}\hfill m·\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}=\frac{L{e}^{-\tau }\mathrm{\pi }{a}^2\overline{q_{\text{pr}}}}{4\mathrm{\pi }{d}^{2}c}-2\mathrm{\pi }{a}^2\mathit{nkT} F(s),\end{array}$$(1)

with$$\begin{array}{c}\hfill F(s)\approx \frac{s}{\alpha }{[1+{(\alpha ·s)}^2]}^{1/2}+{\left(\frac{\mathit{\text{eU}}}{\mathit{\text{kT}}}\right)}^2\text{ln}\Lambda \frac{s}{3\sqrt{\mathrm{\pi }}/4+s^3},\end{array}$$(2) $$\begin{array}{c}\hfill s\equiv \frac{v_{\text{drift}}}{\sqrt{2kT/m_{\mathrm{H}}}},\hspace{1em}\alpha \equiv 3\sqrt{\mathrm{\pi }}/8,\hspace{1em}\Lambda =\frac{3\mathit{\text{kT}}}{2\mathit{ae}|\mathit{eU}|}{\left(\frac{\mathit{kT}}{\mathrm{\pi }{n}_{\mathrm{e}}}\right)}^{1/2}.\end{array}$$(3)

Here, the motion of the dust grain, with mass m and radius a, in the x direction is along the light rays from the star with luminosity L. The light is attenuated by extinction and the grain is at time t located at an optical depth $\tau$ and a distance d from the star. The quantity $q_{\mathrm{pr}}$ is a numerical factor, which is dependent on the grain composition, shape, and wavelength. For the relevant dust and stellar light discussed here a mean $\overline{q_{\mathrm{pr}}}$ of about 1.5 may be adopted (Draine 2011a). The speed of light is c and the mass of the hydrogen atoms or molecules, whatever species assumed to dominate the gas, is $m_{\mathrm{H}}$. The dust grain is retarded by a drag force, described by the second term of the right-hand side of Eq. (1). The variable $v_{\mathrm{drift}}$ is the drift velocity of the grain relative to the gas with the number density of atoms or of molecules n and the kinetic temperature T. In the expression for F(s) in Eq. (2) the last term represents the Coulomb forces, where U is the grain potential, e the electron charge, and $n_{\mathrm{e}}$ the electron density. Following Draine (2011a), we set $|\mathit{\text{eU}}|/\mathit{\text{kT}}=2.5$. For our neutral gas (with degrees of ionization from Table 30.1 in Draine 2011b) we find the Coulomb force term to be negligible for cool (10–100 K) gas and even rather small for warm $\sim 5000$ K gas as compared with the drag force due to the neutral gas collisions. We note that the dependence of the drag force on the drift velocity is linear for small speeds – in typical cases with $T=10$ K the shift to a mainly quadratic dependence occurs around $v_{\mathrm{drift}} \sim 400$ m s${}^{-1}$, which is higher than typical sound speeds in the cool neutral gas. Below, we call the assumption of a linear dependence of the drag on the drift velocity the small velocity approximation.

For the cases to be discussed the acceleration term, the first term in Eq. (1), may often be neglected and we have$$\begin{array}{c}\hfill F(s)=F_1(v_{\text{drift}})=\frac{L{e}^{-\tau } \overline{q_{\text{pr}}}}{8\mathrm{\pi }{d}^{2}c \mathit{\text{nkT}}},\end{array}$$(4)

which then gives $v_{\mathrm{drift}}$ from Eqs. (2) and (3). We note that $v_{\mathrm{drift}}$ is then independent of the grain size a, except for a possible weak dependence through $\overline{q_{\mathrm{pr}}}$.

In the idealized 1D case of a fully homogeneous gas and an initially equally homogeneous dust distribution in the gas, a stellar flux that is suddenly switched on may be expected to cause the dust to assemble in a narrow dust front. Grains that are first hit by the incoming radiation are rapidly accelerated until they reach an equilibrium with the gas friction at drift velocities of perhaps several hundred metres per second; this is dependent on the illuminating star, its distance, the extinction of its light and the gas density. When the optical depths increase due to shielding by the dust passed, and the available momentum flux from the light source is distributed on an increasing number of swept-up dust grains, the velocities are successively retarded to typically just a few metres per second. Dust grains lagging behind the dust front are exposed to the full illumination by the star; these are immediately accelerated and soon reach the front where they lose much of the driving momentum. An interesting aspect of this is that the bigger grains, which have greater inertia, may be expected to reach further into or across the dense front before they are halted by friction. Thus, some separation of dust of various sizes may occur. Similarly, dust grains ahead of the front, which have initially zero drift relative to the gas, are affected differently. The bigger grains should be passed by the front and be more slowly accelerated towards it, again because of their inertia. Thus, in this simplified model the big grains are more dispersed across the front than the small grains.

To explore this further, we modelled the dust distribution by solving Eq. (1) for a stack of plane-parallel dust sheets of equal masses, which are pushed along in a homogeneous and steady gas column by a light source that has a flux perpendicular to the plane of the sheets. Each dust sheet is free to move through the others and the relevant changes in the optical depths $\tau$ are successively taken into account. To be specific, we set the dust mass density ${\rho }_{\mathrm{d}}=1$ g cm${}^{-3}$ and the gas-to-dust mass ratio $\varphi =100$. We switched on an illuminating star at a distance d from the column. Two different sizes of grains were explicitly considered; these have radii a of ${10}^{-4}$ and ${10}^{-5}$ cm, respectively and have equal total mass contributions, and their combined effects were included in the calculation of $\tau$. The density of H${}_2$ molecules $n={10}^5$ cm${}^{-3}$ and for the stellar luminosity and distance we set $L={10}^{38}$ erg s${}^{-1}$, and $d=$ 1 pc. Typical results are shown in Figs. 1 and 2. In the first Figure we illustrate the evolution of a dust region just 50 au thick (limited in size for clarity) and embedded in an extensive homogeneous and static gas cloud. This figure shows how the smaller (black) grains gather in a thin decelerating dust front while the bigger (red) grains start more slowly, which is why they tend to lag behind the smaller grains during the first 50 kyr. Then the bigger grains pass to the front and lead, although the small grains gradually catch up. In Fig. 2, representing the evolution of the illuminated surface layers of an extended dusty gas, the distribution across the front for the two types of grains is shown for a number of times with intervals of 50 kyr. As is seen, the leading more massive grains actually show a two-peak distribution: with a first peak ahead of the small-grain dust front by at first typically 5 au and next successively thinner and the other and more massive dust front more aligned with the front of the small grains. Obviously, the front may be totally about 10 au wide for big grains, which have diameters of $\sim 1$ micron in typical cases, while the front of the smallest grains ($<0.02$ microns) is just a few au.

The picture presented above is, however, too simplified. The front may have a speed relative to the gas of only some metres per second, and the grains that happen to be in front of it are affected by more dynamical forces and motions, for instance related to turbulence. Thus, the dust front may be much less well defined. In addition to this, it is easy to see that the dust front may be affected by instabilities; a spurious reduction of the number of grains in a section of the front admits radiation to push the dust (which may also drag the gas) in that section further ahead of the front. Similarly, an increase of the number of grains in a region causes it to move more slowly. If these perturbations occur early on, in regions with moderate optical depths, they may have noticeable consequences for the shape of the front such as the formation of finger-like structures beyond it.

Fig. 1. Motion of a stack of 100 dust sheets, initially extending across 50 au, moving as a result of perpendicular radiative pressure (from below in the figure) through a homogeneous gas. The motions of small (${10}^{-5}$ cm) and big (${10}^{-4}$ cm) grains are demonstrated by different colours.

Fig. 2. Dust density in the dust front shown in Fig. 1, at 9 different consecutive times with intervals between them of 50 kyr for two different grain sizes.

2.2. Motion of the dust front

We now estimate the velocity of the dust front through the neutral gas and subsequently explore the resulting cleansing effects. We make these estimates based on the assumption of spherical symmetry and the application of the momentum equation and conservation of the number of dust grains. We also make some simple approximations concerning the acceleration phase of the dust during its motion through optically thin layers. We show that the character of the narrow dust front discussed above substantially simplifies the discussion. A general schematic outline of our schematic model is shown in Fig. 3.

We first assume that the dust front becomes optically thick and we explore departures from this later. For the moving dust front at a distance r from the star or central group of stars, the total momentum transferred by the stellar radiation to the dust per second is represented by the LHS in Eq. (5), below. The variable $m_{\mathrm{cl}}(t)$ denotes the total gas mass swept over by the dust front, and thus cleansed, at time t. The radiative momentum partially goes into accelerating the new dust grains caught by the front from the gas to the dust-front velocity as well as accelerating the grains already present there, resulting in the acceleration $\ddot{r}$. The total mass of this dust swept up by the front is $m_{\mathrm{cl}}(t)/\varphi$, where $\varphi$ is the gas to dust mass ratio. If the speed of the dust front is $\mathrm{d}{r}_{\mathrm{d}}/\mathrm{d}t=v_{\mathrm{d}}$ the added momentum per second to the new grains in the front, in addition to the momentum they have when swept along by the gas, is $4\mathrm{\pi }·r{(t)}^{2}n(t)·m_{\mathrm{H}}\mu ·(v_{\mathrm{d}}-v_{\mathrm{g}})/\varphi$, where $v_{\mathrm{g}}$ is the velocity of the gas and $\mu =1+Y+Z,$ where Y and Z are the mass fractions of He and heavier elements in the gas, respectively. The value n(t) is the local number density of gas molecules (or atoms if the molecules are dissociated) at the dust front. The rest of the momentum is lost to the gas by friction, according to the last term in the following equation:$$\begin{array}{c}\hfill \frac{L \overline{q_{\mathit{\text{pr}}}}}{c}=\frac{m_{\mathit{\text{cl}}}(t)}{\varphi }·\ddot{r}+4\mathrm{\pi }r{(t)}^2·\frac{m_{\mathrm{H}}\mu }{\varphi }n(t)·v_{\mathrm{d}}·(v_{\mathrm{d}}-v_{\mathrm{g}})\\ \hfill +\frac{G·m_{\mathrm{cl}}(t)}{r{(t)}^2\varphi }(m_{\mathrm{cl}}(t)+m_{\ast })\\ \hfill +2\mathrm{\pi }kTn(t)·F_1(v_d-v_g){\int }_{a_{\mathrm{l}}}^{a_{\mathrm{u}}}N(a,t)a^2 \mathrm{d}a.\end{array}$$(5)

The value of $\overline{q_{\mathrm{pr}}}$ is dependent on the detailed radiative transfer in the dust front. Assuming that the photons are finally absorbed and give up their total momentum in the optically thick front, we can set $\overline{q_{\mathrm{pr}}}=1$. The fourth term represents the gravity from the gas mass on the dust shell and the gravity from the stars with mass $m_{\ast }$ formed inside the dust shell at early stages before the cleansing process started. One should note that the local number density of H${}_2$ molecules or free hydrogen atoms n(t) may be very different from the initial density, $n_0$. The temperature-dependent choice of hydrogen species is specified later.

The quantity N(a)da is the total number of dust grains in the dust front with radii in the interval ($a,a +$da), and $a_l$ and $a_u$ are the lower and upper limits, respectively, of this dust-size distribution. We have adopted the classical distribution $N(a,t)=N_t{a}^{-3.5}$ of Mathis et al. (1977), see also Casuso & Beckman (2010), in the range from $a={10}^{-7}$–${10}^{-4}$ cm. It is worth noting that, with this size distribution, the most important grains contributing extinction (scaling as $a^2$) are those with smaller sizes, while the mass contribution ($\propto a^3$) is more evenly distributed among the various sizes.

We can now normalize the dust size distribution, assuming spherical grains, to the proper dust mass ratio at the time t $$\begin{array}{c}\hfill {\int }_{a_l}^{a_u}{N}_t{a}^{-3.5}\frac{4\mathrm{\pi }{a}^3}{3}{\rho }_{\mathrm{d}} \mathrm{d}a=\frac{8\mathrm{\pi }}{3}{N}_t{\rho }_{\mathrm{d}}({\mathit{\text{au}}}_{1/2}^{}-a_l^{1/2})=\frac{m_{\mathrm{cl}}(t)}{\varphi }·\end{array}$$(6)

Since $a_u^{1/2}>>a_l^{1/2}$ we may neglect the $a_l^{1/2}$ term in Eq. (6) and obtain$$\begin{array}{c}\hfill {\mathit{\text{Nt}}}_{}=\frac{3}{8\mathrm{\pi }}·\frac{m_{\mathrm{cl}}(t)}{{\rho }_{\mathrm{d}}\varphi \sqrt{{\mathit{\text{au}}}_{}}}·\end{array}$$(7)

Correspondingly approximating the integral in Eq. (5) and defining $\overline{a}\equiv \sqrt{a_l{a}_u}$ we may now write the equation of motion for the dust front, as$$\ddot{r}+\lambda \frac{r{(t)}^{2}n(t)}{m_{\text{cl}}(t)}{\upsilon }_{\text{d}}\cdot ({\upsilon }_{\text{d}}-{\upsilon }_{\text{g}})+\kappa Tn(t)F_1({\upsilon }_{\text{d}}-{\upsilon }_{\text{g}})=-\frac{G}{r{(t)}^2}(m_{\text{cl}}(t)+m_{*})+\frac{\varphi }{c}\frac{L}{m_{\text{cl}}(t)}.$$(8)

Here$$\begin{array}{c}\hfill \kappa =\frac{3}{2}\frac{k}{{\rho }_{\mathrm{d}}\overline{a}}\end{array}$$(9)

and$$\begin{array}{c}\hfill \lambda =4\mathrm{\pi }{m}_{\mathrm{H}}\mu \end{array}$$(10)

are two constants, depending on the dust properties and the chemical composition of the gas, respectively.

Fig. 3. Schematic picture of our model. The bright star(s) is located in the middle, surrounded by the expanding H II region with radius $R_{\mathrm{II}}$. Outside the H II region is dust-cleansed neutral gas, surrounded by the dust front. The proportions in the figure are not realistic: the true thickness of the cleansed region is only a few percent of the radius of the H II region, and the thickness of the dust front is in turn about 1% of that of the cleansed region.

2.3. Cleansed gas mass

It is clear from Eq. (8) that for a full solution we must also supply information about the gas velocity $v_{\mathrm{g}}$ and gas density n(t) as functions of r and t. Another circumstance of importance is the handling of the radiative force term in the optically thin layers of the gas where this equation must be modified. We return to these complications but first develop a simple theory that partially circumvents these difficulties, keeping in mind that our main ambition is to obtain an estimate of the cleansed gas mass, $m_{\mathrm{cl}}(t)$, sooner than studying the detailed dynamics of the situation.

By direct solution with typical parameters we find that the first two terms of Eq. (8) after relatively short times, $t_{\mathrm{d}}$, become unimportant in comparison with the third term, i.e. the dust grains relatively rapidly acquire the drift velocity. For a linear dependence of $F(s) \sim s/\alpha$ in Eq. (2) we find$$\begin{array}{c}\hfill t_{\mathrm{d}}\sim \frac{{\rho }_{\mathrm{d}}\overline{a}}{\sqrt{kT{m}_{\mathrm{H}}}·n(t)},\end{array}$$(11)

while for the case $F(s) \sim s^2$ $$\begin{array}{c}\hfill t_{\mathrm{d}}\sim \frac{{\rho }_{\mathrm{d}}\overline{a}}{n(t)(v_{\mathrm{d}}-v_g)·m_{\mathrm{H}}}.\end{array}$$(12)

Typically, this leads in both cases to $t_{\mathrm{d}} \sim {10}^4$ years or less and total masses of cleansed gas during this time of a few solar masses. For the optically thin case, similar results are obtained from Eqs. (5) and (19).

The neglect of the two first terms in the equation of motion for the dust front is subsequently called the drift approximation. In practice, with realistic values of $m^{\ast }$, we also find the gravity term in Eq. (8) to be unimportant. The gravity term would only be comparable to the dominating terms for relatively compact but massive regions in which $m_{\mathrm{cl},\odot }/(r_{\mathrm{pc}}\times L_{38}^{1/2})>3\times {10}^3$ , where $L_{38}$ is the luminosity of the central star(s) in ${10}^{38}$ erg/s and $m_{\mathrm{cl},\odot }$ is the cleansed mass in solar masses. Also, in such cases we have found by numerical experimentation that the gravity term does not increase the maximum value of $m_{\mathrm{cl}}$. There is, however, the possibility that gravity could, in the long run, turn the expansion of the gas cloud to a collapse with a second phase of star formation. This is discussed later in this work. Thus, neglecting also the gravity term we get$$\begin{array}{c}\hfill \kappa T n(t)F_1(v_{\mathrm{d}}-v_{\mathrm{g}})\approx \frac{\varphi }{c}\frac{L}{m_{\mathrm{cl}}(t)}.\end{array}$$(13)

Defining F and s according to Eqs. (2) and (3) we find the solution for s from the resulting quadratic equation. Next, we obtain the gas mass passed (and thus cleansed) by the dust front in the time interval dt as$$\begin{array}{c}\mathrm{d}{m}_{cl}=({\upsilon }_d-{\upsilon }_g)n(t)\mu m_H\cdot 4\pi r{(}^{t}d\\ =4\sqrt{2\pi }\sqrt{kT{m}_H\mu }\cdot s\cdot n(t)r{(}^{t}dt.\end{array}$$(14)

We now combine the solution for s and Eq. (14) to a differential equation from which $m_{\mathrm{cl}}(t)$ may be obtained, if the variation of the local density at the dust front, n(t), and the distance of the dust front from the star, r(t), are known as follows:$$\frac{d{\text{m}}_{cl}}{d(t)}=\frac{4\pi }{\alpha }\mu \sqrt{kT{m}_H}\times {\left[{\left\{1+{\left(\frac{2{\alpha }^2\varphi \cdot L}{c\kappa \cdot T\cdot m_{cl}(t)\cdot n(t)}\right)}^2\right\}}^{1/2}-1\right]}^{1/2}n(t)\cdot r{(t)}^2.$$(15)

In the early phase of the expansion of the dust front the second term in the RHS of Eq. (15) may dominate. In this case we find$$\begin{array}{c}\hfill \frac{\mathrm{d}{m}_{\mathrm{cl}}}{\mathrm{d}t}=\frac{8\sqrt{\mathrm{\pi }}}{3}\mu ·{[\frac{\varphi {\rho }_{\mathrm{d}}\overline{a}{m}_{\mathrm{H}}}{c}·\frac{Ln(t)}{m_{\mathrm{cl}}(t)}]}^{1/2}·r{(t)}^2.\end{array}$$(16)

However, as the dust front accumulates more dust, so that the stellar radiation momentum is distributed on a greater dust mass, the second term on the RHS of Eq. (15), gets $\ll 1$. Putting in realistic numbers, we find that in this case, implying a shift from a quadratic to a linear dependence of the friction force on the drift velocity, occurs when$$\begin{array}{c}\hfill \frac{T_{10}{m}_{\mathrm{c},\odot }n{(t)}_4}{L_{38}}>100,\end{array}$$(17)

where $T_{10}$ is the temperature in the dust shell in units of 10 K and $n{(t)}_4$ is the local gas number density in ${10}^4$ cm${}^{-3}$. This case thus applies if considerable masses ($m_{\mathrm{cl}}\gg 1$ solar mass) are to be cleansed. It obviously corresponds to the small velocity approximation as discussed above. We have then$$\begin{array}{c}\hfill \frac{\mathrm{d}{m}_{\mathrm{cl}}}{\mathrm{d}t}=\sqrt{2}{\mathrm{\pi }}^{3/2}\mu \sqrt{\frac{m_{\mathrm{H}}}{\mathit{kT}}}·\frac{\varphi {\rho }_{\mathrm{d}}\overline{a}}{c}·\frac{L}{m_{\mathrm{cl}}(t)}·r{(t)}^2.\end{array}$$(18)

We note that the problem is significantly simplified when this approximation can be used. The explicit dependence on density fluctuations through n(t) is then cancelled. Hydrodynamics beyond the simple application of momentum conservation only appears indirectly through the radial distance of the dust front from the star, r(t), and a possible variation of the gas temperature T. Also in the high-velocity case (Eq. (16)) we note that the dependence on density only goes as $n{(t)}^{1/2}$. Since r(t) is very much governed by the expansion of the H II region, the growth of $m_{\mathrm{cl}}$ is thus in general not very sensitive to the details of the structure of the neutral gas shell.

For the case of an optically thin gas Eq. (1) applies to every individual dust grain. In the extreme completely transparent case, with $\tau =0$, the grains initially closest to the star have been pushed out to a certain distance from the star by radiative pressure and there at a given time define the limits of the cleansed region. With the small velocity approximation and applying Eq. (14) we find$$\begin{array}{c}\hfill \frac{\mathrm{d}{m}_{\mathrm{cl}}}{\mathrm{d}t}=\frac{3\sqrt{2\mathrm{\pi }}}{16}\sqrt{\frac{m_{\mathrm{H}}}{\mathit{kT}}}\mu \frac{\overline{q_{\mathrm{pr}}}L}{c}.\end{array}$$(19)

In this case, not only $v_{\mathrm{d}}-v_{\mathrm{g}}$ and n(t) but also r are cancelled, and, with a constant T and L, the cleansed mass grows linearly with time.

It is easily seen by inserting typical numbers into Eqs. (18) and (19) that the growth of the cleansed mass may be orders of magnitude more rapid in the optically thin case than in the thick case. This is an obvious consequence of the shadowing in the thick case of radiation by the grains closest to the star, radiation that would otherwise directly strike further grains ahead. The optical depth $\tau$ of the dust collected from the cleansed mass is$$\begin{array}{c}\hfill \tau ={\int }_{r_0}^r\mathrm{d}r{\int }_{a_l}^{a_u}{N}_v(a,r)\sigma (a)\mathrm{d}a=\frac{3}{16\mathrm{\pi }}\frac{\omega }{{\rho }_{\mathrm{d}}\varphi \overline{a}}·\frac{m_{\mathrm{cl}}}{r{(t)}^2},\end{array}$$(20)

where $\sigma (a)$, which we here set $=\mathrm{\pi }{a}^2$, is the absorption cross section per grain and $\omega$, $1\le \omega \le 3$, depends on the r distribution of dust. The variable $N_v(a,r)$ is the number of dust grains with radii in the interval $(a,a+\mathrm{d}a)$ per volume unit. The lowest value of $\omega$, chosen subsequently, corresponds to a thin dust shell with radius r and the highest to the original homogeneous distribution.

It is important to consider how far the dust front proceeds through the gas before it gets optically thick. If we assume that the initial number density of H molecules in the gas was $n_0$ and that all dust that was initially within the radius r(t) has been collected in the dust front when it reaches this radius at time t, we can use the condition that $\tau <1$ and Eq. (20) to estimate an upper limit of the maximum mass cleansed during this first phase when the dust front is optically thin. We find$$\begin{array}{c}\hfill m_{\mathrm{cl}}<\frac{256\mathrm{\pi }}{3}{\left(\frac{{\rho }_{\mathrm{d}}\varphi \overline{a}}{\omega }\right)}^3·{\left(\frac{1}{m_{\mathrm{H}}\mu }\right)}^2·\frac{1}{n_0^2},\end{array}$$(21)

which, with characteristic numbers for the parameters leads to$$\begin{array}{c}\hfill m_{cl,\odot }<2\times {10}^7\frac{1}{n_0^2},\end{array}$$(22)

where $n_0$ is in cm${}^{-3}$. Since the original molecular number density in the cores of Giant Molecular Clouds is typically of the order of ${10}^4$ cm${}^{-3}$, we conclude that the mass cleansed in such a cloud during the optically thin phase is small.

As a check we modelled the dust motion in the optically thin layer of the dusty molecular cloud closest to the star by subdividing the region into ten different subshells, each with thickness corresponding to 0.2 in optical depth. For these subshells, we applied Eq. (8), properly modified with a factor of $\exp (-\tau )$ multiplied to L. In solving the equation for each subshell individually, we found that these subshells pile up and merge in a dust front close to a common geometrical depth point P, and from that point and onwards we integrated Eq. (8) for the full front. In calculating the swept gas mass we added the different sub-layers for depths outside P. As expected, the effects of this more detailed treatment of the radiative transfer is minor for realistic shells.

It is easy to see that the optical depth of the dust shell increases linearly with its radius. The effect of the increase of the number of swept-up grains, varying as $r^3$, combines with that of the dispersion of the dust which varies as $r^{-2}$. This ends when the border of the cloud is reached. After that the optical depth decreases in proportion to $r{(t)}^{-2}$. However, except for very small clouds, this occurs at such late stages in the evolution of the dust shell, when the speed of the shell is so slow and the star(s) in the centre have declined so much in brightness, that it does not affect the resulting cleansed gas mass much.

Regarding the motion of the dust front, we conclude that the front may be assumed to be optically thick and the drift approximation may be applied (i.e. Eq. (15) may be used), while the small velocity approximation with a linear dependence of the friction force on the drift velocity is not valid in early phases (small r) and strong luminosities (cf. Eq. (17)). It remains to allow for the effects of the functions r(t) and n(t) on the motion of the dust front. This requires considerations concerning the dynamics of the H II region and some discussion of the pressure equilibrium of the neutral gas.

2.4. Ionization front

The evolution of the region close to the bright illuminating star (or stars) is of key significance for the dust cleansing of the surrounding cool gas. The expanding H II region may well rapidly penetrate into the molecular cloud, heat it, and thus erode the possibilities for star formation. Also, this expansion, owing to the great pressure in the hot ionized gas ($T \sim 10 000$ K) and radiation pressure (Krumholz & Matzner 2009, Akimkin et al. 2015, Akimkin et al. 2017), pushes the cool gas with its dust front to a rapidly increasing distance from the star. The gas just outside the H II region is compressed by the outgoing shock wave to high densities (typically $n \sim {10}^{4-5}$ cm${}^{-3}$ just preceding the ionization front (Draine 2011a, Kim et al. 2016). We apply a simple approach to modelling the H II region to obtain a lower limit of r(t) in Eq. (15).

A bright and hot star that forms in a hydrogen cloud with initial mass density ${\rho }_{\mathrm{g}}$ very rapidly ionizes a region around it, generating an initial Strömgren sphere (Strömgren 1939) with a radius$$\begin{array}{c}\hfill R_{\mathrm{Str}}={\left(\frac{3 L f m_{\mathrm{H} \mathrm{I}}^2}{4\mathrm{\pi } {\alpha }_r X{\rho }_{\mathrm{g}}^2}\right)}^{1/3}.\end{array}$$(23)

Here, ${\alpha }_r=2.6\times {10}^{-13}$ cm${}^3$s${}^{-1}$ denotes the Case B recombination coefficient for hydrogen in the on-the-spot approximation (cf. e.g. Walch et al. 2012). The value X is the mass fraction of hydrogen in the gas, while f is the fraction of the stellar flux that is directly able to photoionize hydrogen from the ground state (i.e. with a wavelength $<91.2$ nm). The value of f is typically 10–20% for O-type stars. Next, the great overpressure in this sphere, owing to its high temperature of $\sim {10}^4$ K with respect to its surrounding, causes the gas to expand (D-type expansion). Bisbas et al. (2015) have found that the time development of the radius of the region, $R_{\mathrm{H} \mathrm{II}}(t)$ may be approximated by a combination of two functions, ${\mathcal{R}}_{\mathrm{I}}(t)$ and ${\mathcal{R}}_{\mathrm{II}}(t)$, according to$$R_{\mathrm{H}\hspace{0.17em}\text{II}}(t)={ℛ}_{\text{II}}(t)+(1-0.733\times {\text{e}}^{-t/\text{Myr}})\times ({ℛ}_{\text{I}}(t)-{ℛ}_{\text{II}}(t)),$$(24)

where ${\mathcal{R}}_{\mathrm{I}}$ and ${\mathcal{R}}_{\mathrm{II}}$ are solutions to two ordinary differential equations (see Bisbas et al. 2015, their Eqs. (8) and (11), respectively). The authors compared this approximation with more detailed numerical simulations and found it to agree excellently with the detailed results obtained for $t\ge 0.05$ Myr. We solved those ordinary differential equations adapted to our model parameters and consider that the cleansing process stops if or when the HII front reaches the dust front.

The effects on the dust cleansing in the molecular cloud are dependent on the early evolution of the bright star and its H II region, i.e. on the phase when the radiation, but not the ionization front, has penetrated into the neutral gas. The early evolution of massive stars and their surrounding ionized gas is much dependent on their accretion history. This has been investigated recently in a number of studies; for references see Haemmerlé & Peters (2016). These authors followed accretion histories of up to 10${}^5$ years, where in some cases the model did not approach the zero-age main sequence until late, and found that convective pre-main-sequence (PMS) objects may delay the growth of the H II region by up to 10${}^4$ years. We schematically explored the consequences of this by assuming a linear growth of the Strömgren radius from an initial value corresponding to L / 10 to L within the time $t_{\mathrm{br}}$, which we varied from 0 to ${10}^5$ yr.

The assumed overall arrangement of the star, its expanding H II region into an H I region, and the H${}_2$ cloud where cleansing is expected to occur are illustrated in Fig. 3. In our calculations, we alternatively assumed all neutral gas to be molecular, or all to be atomic, with corresponding temperature corrections; see below. In our calculation of the cleansed mass available for new stars we subtracted this ionized gas. From Bisbas et al. (2015), their Eq. (7), we obtain the following mass estimate for the ionized gas $m_{\mathrm{i}}(t)$:$$\begin{array}{c}\hfill m_{\mathrm{i}}(t)=\frac{4\mathrm{\pi }}{3}{n}_0{m}_{\mathrm{H}}\mu R_{\mathrm{Str}}^{3/2}{R}_{\mathrm{H} \mathrm{II}}^{3/2},\end{array}$$(25)

and thus get the reduced cleansed neutral mass $m_{\mathrm{c}}(t)$ available for star formation as follows:$$\begin{array}{c}\hfill m_{\mathrm{c}}(t)=m_{\mathrm{cl}}(t)-m_{\mathrm{i}}(t).\end{array}$$(26)

The composition and fate of the dust in the H II region is an important issue. How much of the dust is destroyed in the hot region, how much is piled up at the border and then gradually pushed into the neutral gas, and how much dust is shot far outside the H II region and maybe even outside the molecular cloud? There have been works carried out on related situations (see e.g. Krumholz & Matzner 2009, Draine 2011a, Martínez-González et al. 2014, Kim et al. 2016), but not explicitly on the effects on the dust density in the gas outside the H II region.

2.5. Gas density and radius of the dust front

Obviously the gas plays an important role in regulating the cleansing effects, primarily through its velocity in the friction term in the equation of motion for the dust (Eq. (8)), but also more indirectly in affecting the density distribution and transporting mass and momentum through gas flows, sound waves, and shocks. Since our interest is focussed on determining an upper limit on the cleansed mass and the drift approximation may be used, we do not need to determine the gas velocity directly; what matters is primarily the local gas density at the dust front, n(t). In order to estimate that we applied the quasi-hydrostatic-equilibrium condition in the gas shell, following (Rahner et al. 2017) as suggested by the observations by Pellegrini et al. (2007). Thus, the gradient of the internal pressure in the gas, P, is assumed to be in equilibrium with the gradient of the deposited radiative momentum per surface area and time unit. For any given time t, we may then obtain the pressure structure by solving$$\begin{array}{c}\hfill \frac{\mathrm{d}P(r)}{\mathrm{d}r}=\frac{L}{4\mathrm{\pi }c{r}^2}·{\mathrm{e}}^{-\tau }\frac{\mathrm{d}\tau }{\mathrm{d}r},\end{array}$$(27)

with the gradient of $\tau$ from Eq. (20),$$\begin{array}{c}\hfill \frac{\mathrm{d}{\tau }_d}{\mathrm{d}r}=\frac{3}{8\mathrm{\pi }}\frac{\omega }{{\rho }_{\mathrm{d}}\varphi \overline{a}}·\frac{m_{\mathrm{cl}}(t)}{r^3}.\end{array}$$(28)

We solve Eqs. (13) and (28) together with the proper differential equation (Eqs. (15) or (19)) and with the boundary conditions $\tau (r_0)=0$, $P(r_0)=P_0$ and the initial condition $m_{\mathrm{cl}}(0)=0$. We combine the solution P(r) with the equation of state to get $n(t)=P(t)/(kT)$ (for values of the temperature T, see Sect. 2.6 below). We find an almost linear increase of log n(t) with log$\tau$, up to a constant value at $\tau \sim 1$. This is found in practice to be independent on the values of $r_0$, which are set to $R_{\mathrm{Str}}$, and of $P_0$. This value, $n(\tau =1),$ is taken as representative for n(t), the density in the dust shell.

We also need to estimate the radius of the dust shell. This is done by noting that the gas pressure should dominate in the warm ionized gas (Lopez et al. 2014) and assuming that the gas pressures on each side of the ionization frontier are close to identical. Thus, $n_{\mathrm{i}}{T}_{\mathrm{i}}=n_{\mathrm{n}}{T}_{\mathrm{n}}$, where $n_{\mathrm{i}}$ and $T_{\mathrm{i}}$ are the number density of particles (nuclei and electrons) and temperature, respectively, in the ionized zone and $n_{\mathrm{n}}$ and $T_{\mathrm{n}}$ are the corresponding quantities in the neutral zone, extending from the ionization front to the inner border of the dust front. Kim et al. (2016) have shown that the piling up of ionized gas at the ionization front from radiative pressure may increase the thermal pressure at the border by as much as a factor of two or more for high values of L and $n_0$, while the gas gets highly depleted inside the H II region. We find, however that such an increase of $n_{\mathrm{i}}{T}_{\mathrm{i}}$ in our boundary condition only leads to small effects on the calculated cleansed mass. We obtain$$\begin{array}{c}\hfill \frac{m_{\mathrm{i}}(t)·T_i}{{\mu }_{\mathrm{i}}·\frac{4\mathrm{\pi }}{3}{R}_{\mathrm{H} \mathrm{II}}{(t)}^3}=\frac{m_{\mathrm{c}}(t)·T_n}{{\mu }_{\mathrm{n}}·\frac{4\mathrm{\pi }}{3}(r{(t)}^3-R_{\mathrm{H} \mathrm{II}}{(t)}^3)},\end{array}$$(29)

where ${\mu }_{\mathrm{i}}$ and ${\mu }_{\mathrm{n}}$ are mean molecular weights in the ionized and neutral regions, respectively. From this we find$$\begin{array}{c}\hfill r(t)=R_{\mathrm{H} \mathrm{II}}(t)·{[1+0.65·\frac{m_{\mathrm{c}}(t)}{m_{\mathrm{i}}(t)}\frac{T_n}{T_i}]}^{1/3},\end{array}$$(30)

and may thus estimate r(t) from $R_{\mathrm{H} \mathrm{II}}(t)$ and $m_{\mathrm{i}}(t)$ as obtained from Eqs. (24) and (25).

2.6. Model

In the model calculations the equation of motion for the dust front as transformed to the equation for the growth of the cleansed mass, alternatively Eqs. (15) and (19) depending on the optical depth, was solved together with the equation for the ionization front, Eq. (24), and its two underlying differential equations; these equations were supplied with the equations for estimating the density, optical depth, and radius of the dust front (Eqs. (27), (28), and (30)). The system of equations is efficiently solved by the Matlab ODE solver ODE23.

In combining the elements discussed above into a relevant model, a number of parameters are to be considered and set:

L: The luminosity of the illuminating bright star (or sum of luminosities for an ensemble of such stars). Values of L ranging from ${10}^{37}$ erg s${}^{-1}$ to ${10}^{40}$ erg s${}^{-1}$ have been chosen. A flux of ${10}^{37}$ erg s${}^{-1}$ corresponds to a main-sequence star of 6–9 $M$${}_{\odot }$, ${10}^{38}$ erg s${}^{-1}$ to 14–16 $M$${}_{\odot }$, and ${10}^{39}$ to 35–80 $M$${}_{\odot }$, according to the mass-luminosity relations of Salaris & Cassisi (2005) and Eker et al. (2015), see also the evolution for a 60 $M$${}_{\odot }$ model by Groh et al. (2014). If we adopt a Salpeter initial mass function (IMF; Salpeter 1955) with $\mathcal{N}(m)\propto m^{-2.35}$, we find that the summed-up luminosity distribution from a young stellar generation peaks at stellar masses of $\sim 20$$M$${}_{\odot }$, and that the light from stars around that value dominates over that from the much fewer stars with $m >$ 40 $M$${}_{\odot }$. It also follows that the integrated luminosity from 30 000 stars (i.e. of the magnitude, including binaries, suggested by Hurley et al. (2005) to constitute the starting population of the cluster M 67), distributed according to the Salpeter IMF, produces a luminosity of $1.8\times {10}^{39}$ erg s${}^{-1}$. The strong luminosity increase of the stars in connection with their supernova phases turns out to be of less importance because of the short duration of those as compared with the main-sequence phases; we neglect the effects of the supernova winds, which are not thought to be very selective in their effects on grains relative to gas. The luminosity of an individual star as well as its effective temperature may however vary significantly during the accreting PMS phase with timescales as short as 100–1000 years (cf. Haemmerlé & Peters 2016); this leads to variations in the number of photoionizing photons by about a factor of two, which affects our f parameter.

f: The fraction of the stellar flux that is emitted at wavelengths shorter than 91.2 nm, i.e. that is able to photoionize hydrogen atoms from the ground state. Values chosen here are 0.1 and 0.2, as estimated from the model atmospheres of hot stars by Castelli & Kurucz (2004). We note however that according to the tables supplied by Rahner et al. (2017), values as high as 0.4 may be valid for early stellar populations.

$t_{\mathrm{br}}$: The characteristic time for the H II region around the still forming bright star to grow to the initial Strömgren radius. As a standard, we set the parameter to 0, but also explored the effects of setting it to values up to 100 000 years. This, following Haemmerlé & Peters (2016), probably overestimates the effects of the delayed growth of the region.

$t_{\mathrm{end}}$: The ending time for the evolution of the model system is set to 6 Myr. During this time the central star/s, after an initial brightening, is assumed to have a constant total luminosity. A value of 6 Myr is significantly longer than the lifetimes of the most massive stars, but since the stars with masses around 20 solar masses collectively contribute most ionizing radiation, and the formation time of such stars in a cluster may be somewhat different, we have chosen this value of $t_{\mathrm{end}}$. Also, we note that the cleansed mass $m_{\mathrm{c}}(t)$ tends to level out beyond 1 Myr (see Fig. 5), which makes the choice of $t_{\mathrm{end}}$ less significant.

$n_0$: The initial density in the gas determines the expansion speeds of the H II region and the dust front and the rate of cleansing, $m_{\mathrm{c}}(t)$. The density of H${}_2$ molecules in the cloud, or alternatively of H I if a temperature of 100 K is chosen (see below), is usually estimated to be in the range 10${}^2$–10${}^5$ cm${}^{-3}$ (Bergin & Tafalla 2007) this range reflects the highly clumpy character of the medium. The various values in this range have been chosen alternatively.

$T_{\mathrm{n}}$: The temperature of the neutral gas, alternatively chosen to be 10 and 100 K, for an adopted composition of H${}_2$ molecules and H atoms, respectively. The value $m_{\mathrm{H}}$ is changed from the mass of one hydrogen molecule to the mass of one hydrogen atom, accordingly. As a standard, we chose the molecular alternative, i.e. $T_0=10$ K and $m_{\mathrm{H}}=2\times 1.67\times {10}^{-24}$ g, and after that studied the consequences of the alternative choice.

$T_{\mathrm{i}}$: The temperature of the gas with fully ionized hydrogen is set to 10 000 K.

$\varphi$: The gas/dust ratio by mass in star-forming regions is usually set to 100 (see e.g. Liseau et al. 2015) and this value has been chosen here.

$\mu$: The mass density of the gas, as normalized on the mass density of hydrogen, $\mu =1+(Y+Z)/X$, is set equal to 1.35.

${\rho }_{\mathrm{d}}$: The mass density of the dust grains is set to 1.0 g cm${}^{-3}$.

The maximum cleansed masses in 6 Myr, as a function of the stellar luminosity, are illustrated for different $n_0$ in Fig. 4. The time developments of the cleansed mass and radius of the dust front are shown for the various choices of the parameters in Fig. 5.

Fig. 4. Logarithmic cleansed mass, $m_{\mathrm{cl}}(t)$, after 6 Gyr as a function of the luminosity of the central light source for three different initial molecular number densities $n_0$.

Fig. 5. Logarithmic cleansed mass, $m_{\mathrm{c}}(t)$ (black and blue lines), and the radius of the dust front, r(t) (red lines), as functions of time for two different luminosities of the central source, L, and two different initial molecular number densities $n_0$.

3. Results and discussion

3.1. Results

Figure 5 shows how the rapid growth of $m_{\mathrm{c}}$, while the optical depth is small, decelerates when the dust front gets optically thick and then follows Eq. (18). We can easily understand the main behaviour in this figure by approximating r(t) by $R_{\mathrm{H} \mathrm{II}}$ according to the Hosokawa & Inutsuka approximation (see Bisbas et al. 2015, their Eq. (13)) which for $R_{\mathrm{St}}\ll {\mathcal{R}}_{\mathrm{I}}\ll {(T_i/T_0)}^{3/4}$ leads to$$\begin{array}{c}\hfill R_{\mathrm{HII}}\sim R_{\mathit{St}}^{3/7}{t}^{4/7}.\end{array}$$(31)

Then, we find from the Eqs. (23) and (18) by integration$$\begin{array}{c}\hfill m_{\mathrm{cl}}(t)\sim L^{1/2}{n}_0^{-1/4}·t,\end{array}$$(32)

and$$\begin{array}{c}\hfill r(t)\sim n_0^{-2/7}·t^{4/7},\end{array}$$(33)

which agrees well with the tendencies in Figs. 4 and 5. It is interesting to note that the cleansed mass $m_{\mathrm{c}}$ is only substantially affected by the density in the cloud via the effects of the initial density $n_0$ onto the radius of the ionized region and that an increasing density thus leads to a reduction of $m_{\mathrm{c}}$.

A measure of some, but limited, value of the ability for stars to form by gravitational instability out of the cleansed gas would be the thickness of the cleansed shell D relative to the Jeans length, $l_{\mathrm{J}}$,$$\begin{array}{c}\hfill l_{\mathrm{J}}=\sqrt{\frac{15kT}{4\mathrm{\pi }G\mu \rho }}.\end{array}$$(34)

We have used the last term of Eq. (30) to obtain D and estimated the gas density $\rho (t)$ as$$\begin{array}{c}\hfill \rho (t)=\frac{m_{\mathrm{c}}(t)}{4\mathrm{\pi }{(R_{\mathrm{II}}+D)}^2·D}.\end{array}$$(35)

The results are illustrated in Fig. 6. It is seen that the thickness of the shell of cleansed gas is relatively small, as compared with the Jeans length, which seems to indicate that effective star formation in this shell may be problematic.

3.2. Discussion

It is clear from the results above that an efficient production of cleansed gas is favoured by a high stellar luminosity L and a relatively low initial gas density $n_0$. In order to obtain the observed abundance effect of Meléndez et al. (2009) complete dust cleansing is not necessary; a reduction limited to about 20$\%$ of the refractories relative to volatiles would suffice. Thus, a cleansed mass of just 0.2 $M$${}_{\odot }$ would in principle be sufficient to produce the effect observed for the Sun, relative to solar twins in our neighbourhood in the Galaxy. However, this relaxing circumstance is counterbalanced by the fact that the star-fomation efficiency would probably not be much higher than 20% although local efficiencies of 50% or more cannot be excluded (see e.g. Kim et al. 2016 and references given there).

It is important to consider whether the rough approximations behind the present estimates may have led to severe underestimates of the cleansed masses. One such approximation may be the early evolution of the bright star with a successive brightening instead of, as supposed here, a direct switch on to a given constant L. We tested the effects of a gradual brightening of the star, from 1/10 L to L in ${10}^5$ yr, but found only a small effect on the cleansed masses of less than a few percent. The effect of our schematic assumption of all the neutral gas to be molecular and at a low temperature was studied by adopting the alternative choice of atomic gas at 100 K. It was found that the reduction of $m_{\mathrm{H}}$ by a factor of two only caused small changes in $m_{\mathrm{c}}(t),$ while the increase of the temperature typically reduced $m_{\mathrm{c}}(t)$ by a factor of 2.2. This may also be seen directly from Eq. (18). We also considered the possibility that the gas at the dust front and beyond that is molecular, while the cleansed gas is hotter and atomic again with small consequences. More important are the effects of inhomogeneities on the expansion rates of the H II regions, found by Bisbas et al. (2015) to slow the mean expansion by up to 40%. This would lead to a systematic increase of the possible cleansed masses by up to a factor of about 2.

A severe problem for the cleansing mechanism studied here is that mixing may even out the differences of gas to dust ratios between cleansed and uncleansed regions. Turbulence with typical velocities of 1–2 km s${}^{-1}$ is observed in molecular clouds (see Hennebelle & Falgarone 2012 for a review). The extent to which turbulence is generated by self-gravity sooner than the result of, for example radiative feedback, from star-forming activity (for reviews, see Dobbs et al. 2014, Padoan et al. 2014, and also Raskutti et al. 2016) is still not clear. However, there are observational indications that turbulence may be strong in regions of a cloud in the early phases of star formation and even be the result of gravitational collapse (see Palau et al. 2015). It may also be driven by Rayleigh–Taylor and Kelvin–Helmholtz instabilities resulting from large-scale flows when the bubble of H II gas collides with the surrounding dense molecular gas; see Ochsendorf et al. (2014). Such turbulence could be assumed to erase the effects of radiative cleansing and fractionation discussed in this work; the velocities caused by the latter effects are about one to two orders of magnitude less than the observed turbulent velocities. As shown in Fig. 7, the drift velocities are only higher than 0.1 km s${}^{-1}$ in the very first phases, that is when $m_{\mathrm{c}}(t)$ is less than $10\%$ of its final value, the mass of the dust in the dust front is small and the front is close to the radiatively accelerating star/s. It is clear from the Figure that the highest drift velocities, which could perhaps compete with the turbulent mixing, would be in the high-luminosity low-density case, where $n_0 \sim {10}^2 {\mathrm{cm}}^{-3}$. If the dusty gas were still optically thin high velocities could possibly be expected even after a considerable gas mass were cleansed. Applying Eqs. (2), (3) and (4) it is easy to see that the key factor is $L/(d^{2}n),$ where d is the distance from the star/s to the dust grain. Setting $\tau \approx 0$ we find from Fig. 8 that$$\begin{array}{c}\hfill \psi =\frac{(L/{10}^{38}\mathrm{erg}{\mathrm{s}}^{-1})}{{(d/1 \mathrm{pc})}^2(n/{10}^4 {\mathrm{cm}}^{-3})}>10,\end{array}$$(36)

is the condition for $v_{\mathrm{drift}}$ to be equal or greater than 1 km s${}^{-1}$ in the optically thin case, for both the case of cool molecular and a somewhat hotter atomic hydrogen gas. Obviously, this condition on $\psi$ may be fulfilled, for example by the early M 67 cloud in which L may have been about $2\times {10}^{39}$, and in particular if n had been small.

However, the requirement of a small optical depth is indeed not met for this scenario. In order to obtain a satisfactory yield of cleansed gas to explain the M 67 abundances, we find from Eq. (20) that $\tau \sim 400 m_{c,4}/d_{\mathrm{pc}}^2$, where $m_{\mathrm{c},4}$ is the cleansed gas mass in units of ${10}^4$ solar masses. This estimate brings us back to the estimates made above for optically thick dust fronts. Obviously, it is impossible to meet both the requirements of a high drift velocity of the order of 1 km s${}^{-1}$ and a high cleansed mass.

The question is then whether the small drift velocities, and the observed high turbulence velocities, make the radiative cleansing hypothesis untenable. If the minimum spatial scale of the turbulence is significantly smaller than our typical values of the thickness of the cleansed gas shell D, of the order of 10 000 au in typical cases, the cleansing is inhibited or at least reduced considerably. Although the spacial resolution in the observations are limited, there are no signs in present data suggesting that this turbulence would be confined to large scales. Instead, contemporary simulations suggest smaller scales (Seifried et al. 2017).

It remains to explore whether efficient small-scale mixing close to the dust front, and further successive small-scale mixing in the cleansed gas would be able to totally inhibit the effects of the cleansing.

The Giant Molecular Clouds are known to have complex structures that have dense clumps of H${}_2$ gas interfoliated with regions with cool molecular gas with much lower densities (see further Hennebelle & Falgarone 2012). If such a cloud is exposed to hot stellar radiation, the less dense regions may be rapidly cleansed. It can be argued that the inhomogeneities, including those in the dust distribution (cf. e.g. Hopkins & Lee 2016) could lead to instabilities in the neutral gas that could possibly amplify the radiative cleansing. Instabilities of various types are known to occur during the evolution of H II regions, not the least at their outer boundaries (Garcia-Segura & Franco 1996, Whalen & Norman 2008a,b). Of a certain interest are the “resonant drag instabilities” (RDIs) of dust-gas mixtures (Squire & Hopkins 2018), which are unstable at a wide range of spatial scales and rapidly growing, in particular when the phase velocity of a wave in the gas agrees with the drift velocity component along that wave propagation vector (Hopkins & Squire 2018a). As pointed out by the referee of the present paper, these instabilities might affect the degrees of radiative cleansing. It should be noted, however, that the cleansed gas borders the dust front which in itself is very over-rich in dust. It remains to study how great variations in the gas to dust ratio that may occur in these regions, and in particular what the possibilities are for dust-cleansed gas to remain so.

An argument against the effects of turbulent mixing of a cleansed and a more dust-rich gas could be the existence of magnetic fields in the gas. If the fields were structured perpendicularly relative to the direction towards the illuminating stars, the fields could hamper turbulence provided that the ionization of the gas is non-negligible. However, such fields would also reduce the drift of the photoelectrically charged dust grains and thus diminish radiative cleansing considerably; see Draine (2011a). As recently shown by Hopkins & Squire (2018b), the wealth of the newly discovered RDIs in magnetized dusty gases may, however, lead to modifications of this conclusion due to more or less vastly growing magneto-sonic RDIs. At any rate, the inhomogeneities and local instabilities discussed in this work would probably not contribute a systematic global dust cleansing of as much mass as needed for a cluster such as M 67.

An interesting issue is whether the small possibilities for global radiative cleansing could be increased by a second gravitational collapse of the star-forming cloud, such as that which may have affected the core of 30 Doradus where two different star generations have been found (Sabbi et al. 2012). Rahner et al. (2017) have studied the circumstances under which re-collapse may occur and found that the star-formation efficiency must be below a few percent for low-mass clouds and below about 25% for clouds with masses of the order of ${10}^6$$M$${}_{\odot }$ to prevent stellar feedback due to the first generation of stars from totally dispersing the cloud. The characteristic time from a star-forming episode to the next would be of the order of some free-fall times (i.e. 3–6 Myr). If we could assume that the large-scale spherical symmetry could survive these episodes, such that the dust shell is not fully mixed into the collapsing cleansed gas within it, we might envisage a gradual dust-cleansing in several steps. For this to work, however, the effects of turbulence must again be limited. It is highly questionable whether this is possible. Certainly, during a collapse the drift speed increases with the decreasing distance d to the cloud centre in proportion to $d^{-2}$, but it also decreases with an increasing gas density $n \sim d^3$; at a uniform contraction this gas density dependence combines with the $d^{-2}$ dependence to a linear dependence of the drift speed on d (cf. Eq. (1)). Thus, the drift speed in a cloud contracting towards newly formed stars in fact decreases to values far below turbulence velocities of the order of 1 km s${}^{-1}$.

Fig. 6. Ratio of the thickness of the cleansed shell relative to the Jeans length, $D/l_{\mathrm{J}}$, for models with various central luminosities L and initial densities $n_0$.

Fig. 7. Logarithmic drift velocity $v_{\mathrm{drift}}$ in m s${}^{-1}$ a function of the cleansed mass, $m_{\mathrm{c}}(t)$ for a number of different stellar luminosities and initial number densities of H${}_2$ molecules. The horizontal line indicates a typical value of the sound speed.

Fig. 8. Logarithmic drift velocity in m s${}^{-1}$ vs. the quantity $\psi =L_{38}/d_{\mathit{pc}}^2·n_4^{-1}$ (for explanation, see text). Blue curve: atomic gas and $T_0=100$ K, red curve: molecular gas, $T_0=10$ K.

4. Conclusions

Using simple arguments, essentially based on momentum conservation in a dusty neutral and homogeneous gas affected by radiative pressure from nearby hot stars, we have demonstrated that radiative cleansing of star-forming clouds might in principle explain composition differences between some single stars, such as the Sun as compared with solar twins. However, as shown in Fig. 6 the cleansed gas shell is so geometrically thin that considerable inhomogeneities are needed to induce gravitational instabilities of cleansed gas. In order to produce such abundance effects for a large percentage of stars in a rich cluster, such as those observed for M 67 (Önehag et al. 2011, 2014) an unrealistic radiative flux is required. As discussed above, the luminosity of all the cluster stars in its youth is expected to be about $2\times {10}^{39}$ erg/s. This may produce, following the results presented in Fig. 5, well-cleansed gas for maximally a few thousand stars. To this we should also apply the modifications due to the limits of 20% cleansing observed for the refractory elements and the counteracting star-forming efficiency discussed above. Another difficulty is posed by the relative chemical homogeneity found for the 14 cluster members measured by Önehag et al. (2014). It is hard to understand why the progression of the dust front in a molecular cloud, surrounding a massive association of hot and bright stars, would leave behind a gas that is systematically cleansed from dust and corresponding heavy elements to the observed degree. One order of magnitude more extensive cleansing than found in the present study is needed, or even more than that, since our estimates are believed to give upper boundaries to the real effects. We must conclude that the mechanism discussed in this work does not explain the composition of M 67. Also, shocks with density fluctuations, magnetic fields, and not the least turbulence should reduce the cleansing very considerably. The drift speeds of the dust relative to the gas in most of our models are an order of magnitude smaller than the observed turbulence velocities., which also casts considerable doubt on whether this is the real explanation for the odd solar composition.

More detailed 3D hydrodynamical simulations of dust cleansing at the rim of expanding H II regions with due consideration of shocks and instabilities would be interesting. We do not, however, consider it probable that such simulations will lead to results that enhance the cleansing effects enough to explain our observations of M 67. The present result urges continued detailed studies of accurate stellar abundances in M 67 and other rich clusters to verify our present observational finding of global enhancement of volatile elements (or missing refractories) in the cluster relative to the majority of field stars.

Looking back at the development of this line of research, our finding in itself may be reflected on with some irony. The discovery of the solar-composition peculiarity by Meléndez et al. (2009) stimulated various proposed explanations. One of these was the dust-cleansing hypothesis, which, in fact, was the reason for our decision to study M 67 in detail. We then thought that a demonstration of solar-like abundance anomalies for the cluster would be a clear indication that the solar composition is the result of dust cleansing. The anomaly was indeed found for the cluster by Önehag et al. (2011 and 2014) and later verified by Liu et al. (2016). But this fact, we now find, does not verify the hypothesis but just elicits further questions. It turns out to be too difficult to cleanse enough of gas to get a sufficiently cleansed protocluster of this magnitude. So, we cannot explain the composition of M 67, as the result of self-cleansing of the gas by the first hot stars in the forming cluster. Other possibilities, like those of the chemo-dynamical evolution of the Galaxy, have to be considered further.

Since the discovery by Meléndez et al. (2009), a trend contrary to the more or less uniform tendency for dust depletion of the M 67 stars has been discovered: several physical binary stars show significant abundance differences between their ingoing components and differences for the various elements are correlated with the condensation temperature of the element in several cases. Notable such cases are the binaries 16 Cyg (Tucci Maia et al. 2014, Nissen et al. 2017), XO-2N/XO-2S (Ramírez et al. 2015, Biazzo et al. 2015), WASP 94 A/B (Teske et al. 2016a), HD 133131 A/B (Teske et al. 2016b), HAT-P-4 (Saffe et al. 2017), and the co-moving pair HD 240430/HD 1240429 (Oh et al. 2018). These findings certainly suggest more local effects than those explored here for M 67. They are not in total disagreement with the hypothesis of dust cleansing by bright stars; the results of the present study indicate that if such cleansing occurs, it would be more local as a result of the thinness of the cleansed shell that is probably deformed by instabilities. Turbulent dust or gas separation, such as that discussed by Liu et al. (2016), may possibly also play a role here. Some further efforts should thus be spent on modelling the possible cleansing of individual smaller gas clouds, like the protosolar cloud, by neighbouring hot stars. Similar simulations, in 3D with proper treatment of instabilities in the neighbourhood of the HII front and its preceding shock front expanding into a clumpy molecular cloud, are already carried out by several groups, but might also include the radiative fields and turbulence and their effects on the dust drift to illuminate cleansing possibilities.

Other suggested reasons for the odd solar atmospheric composition, such as those presented in Sect. 1 including the early suggested planetary connection and the possibility of radiative self-cleansing from dust by protostars of their late accretion flows, should be further scrutinized.

Acknowledgments

Bengt Edvardsson, Kjell Eriksson, Bernd Freytag, Anders Johansen, Lars Mattsson, Jorge Meléndez, Dhruba Mitra, Poul Erik Nissen, and an unknown referee are thanked for many valuable suggestions and Amy Mednick for comments on the language of the manuscript.

References

All Figures

	Fig. 1. Motion of a stack of 100 dust sheets, initially extending across 50 au, moving as a result of perpendicular radiative pressure (from below in the figure) through a homogeneous gas. The motions of small (${10}^{-5}$ cm) and big (${10}^{-4}$ cm) grains are demonstrated by different colours.
In the text

	Fig. 2. Dust density in the dust front shown in Fig. 1, at 9 different consecutive times with intervals between them of 50 kyr for two different grain sizes.
In the text

Fig. 3. Schematic picture of our model. The bright star(s) is located in the middle, surrounded by the expanding H II region with radius $R_{\mathrm{II}}$. Outside the H II region is dust-cleansed neutral gas, surrounded by the dust front. The proportions in the figure are not realistic: the true thickness of the cleansed region is only a few percent of the radius of the H II region, and the thickness of the dust front is in turn about 1% of that of the cleansed region.

In the text

	Fig. 4. Logarithmic cleansed mass, $m_{\mathrm{cl}}(t)$, after 6 Gyr as a function of the luminosity of the central light source for three different initial molecular number densities $n_0$.
In the text

	Fig. 5. Logarithmic cleansed mass, $m_{\mathrm{c}}(t)$ (black and blue lines), and the radius of the dust front, r(t) (red lines), as functions of time for two different luminosities of the central source, L, and two different initial molecular number densities $n_0$.
In the text

	Fig. 6. Ratio of the thickness of the cleansed shell relative to the Jeans length, $D/l_{\mathrm{J}}$, for models with various central luminosities L and initial densities $n_0$.
In the text

	Fig. 7. Logarithmic drift velocity $v_{\mathrm{drift}}$ in m s${}^{-1}$ a function of the cleansed mass, $m_{\mathrm{c}}(t)$ for a number of different stellar luminosities and initial number densities of H${}_2$ molecules. The horizontal line indicates a typical value of the sound speed.
In the text

	Fig. 8. Logarithmic drift velocity in m s${}^{-1}$ vs. the quantity $\psi =L_{38}/d_{\mathit{pc}}^2·n_4^{-1}$ (for explanation, see text). Blue curve: atomic gas and $T_0=100$ K, red curve: molecular gas, $T_0=10$ K.
In the text