# Dust-Evacuated Zones Near Massive Stars: Consequences of Dust Dynamics on
Star-forming Regions

Nadine H. Soliman TAPIR, Mailcode 350-17, California Institute of Technology,
Pasadena, CA 91125, USA  Philip F. Hopkins  TAPIR, Mailcode 350-17, California
Institute of Technology, Pasadena, CA 91125, USA  Michael Y. Grudić  Carnegie
Observatories, 813 Santa Barbara St, Pasadena, CA 91101, USA

###### Abstract

Stars form within dense cores composed of both gas and dust within molecular
clouds. However, despite the crucial role that dust plays in the star
formation process, its dynamics is frequently overlooked, with the common
assumption being a constant, spatially uniform dust-to-gas ratio and grain
size spectrum. In this study, we introduce a set of radiation-dust-
magnetohydrodynamic simulations of star forming molecular clouds from the
STARFORGE project. These simulations expand upon the earlier radiation MHD
models, which included cooling, individual star formation, and feedback.
Notably, they explicitly address the dynamics of dust grains, considering
radiation, drag, and Lorentz forces acting on a diverse size spectrum of live
dust grains. We find that interactions between radiation and dust
significantly influence the properties of gas surrounding and accreting onto
massive stars. Specifically, we find that once stars exceed a certain mass
threshold ($\sim 2M_{\odot}$), their emitted radiation can evacuate dust
grains from their vicinity, giving rise to a dust-suppressed zone of size
$\sim 100$ AU. Commencing during the early accretion stages and preceding the
Main-sequence phase, this process results in a mass-dependent depletion in the
accreted dust-to-gas (ADG) mass ratio within both the circumstellar disc and
the star. We predict massive stars ($\gtrsim 10M_{\odot}$) would exhibit ADG
ratios that are approximately one order of magnitude lower than that of their
parent clouds. Consequently, stars, their discs, and circumstellar
environments would display notable deviations in the abundances of elements
commonly associated with dust grains, such as carbon and oxygen.

Star formation(1569) — Stellar abundances(1577) — Massive stars(732) —
Interstellar dynamics(839) — Interstellar dust(836) — Interstellar dust
extinction(837)

STARFORGE: Dust-Evacuated ZonesSTARFORGE: Dust-Evacuated Zones

##  1 Introduction

Star formation and stellar evolution are complex processes influenced by a
multitude of physical mechanisms and environmental factors within Giant
Molecular Clouds (GMCs) (McKee & Ostriker, 2007; Girichidis et al., 2020). The
presence of magnetized, supersonic, and turbulent flows within these clouds
drives strong density fluctuations, giving rise to regions with varying
densities and spatial dimensions (Larson, 1981; Mac Low & Klessen, 2004). As
these density fluctuations develop, some reach a critical point where their
gravitational force begins to dominate, initiating protostellar collapse. The
collapse of these gas and dust over-densities gives birth to protostars, which
accrete nearby material and evolve into mature stars. During this process,
they dynamically interact with the surrounding cloud through feedback
mechanisms such as radiation, jets, radiatively-driven stellar winds, and
supernova explosions (Krause et al., 2020).

Dust, an essential component of GMCs, heavily influences all of these
processes. These particles, known as “grains”, are created as a byproduct of
stellar evolution, within the atmospheres of evolved stars and in supernova
remnants. Dust grains reprocess stellar radiation, absorbing far-ultraviolet
(FUV), optical and infrared (IR) photons, re-emitting them in the IR and
submillimeter wavelengths (Draine & Lee, 1984; Mathis, 1990; Li & Draine,
2001; Tielens, 2005). Furthermore, they scatter background radiation and emit
a thermal continuum in the IR. In addition to its radiative interactions, dust
governs the thermodynamics and chemistry of gas within the ISM, exerting a
substantial influence on the intricate processes triggering star formation
(Whittet et al., 1993; Salpeter, 1977; Weingartner & Draine, 2001, 2001;
Draine, 2003; Dorschner, 2003; Spitzer Jr, 2008; Watanabe & Kouchi, 2008;
Minissale et al., 2016). Moreover, dust acts as a reservoir, confining heavy
elements within a solid phase, which can subsequently become integrated into
stars and planets, thereby significantly influencing their overall
compositions. Once luminous sources form within the GMC, the presence of dust
can also initiate radiatively driven outflows in the surrounding regions
(Murray et al., 2005; King & Pounds, 2015; Höfner & Olofsson, 2018).

In addition, the dynamic interplay between dust and gas is pivotal in shaping
the evolution of star-forming regions. Turbulence within the cloud can
generate variations in the dust density on scales comparable to the turbulent
eddy turnover scale, leading to substantial fluctuations in the dust-to-gas
(DTG) ratio within prestellar cores (Thoraval et al., 1997, 1999; Abergel et
al., 2002; Flagey et al., 2009; Boogert et al., 2013). As demonstrated by
Hopkins (2014), this phenomenon has significant implications within large
GMCs, where dust dynamics alone could potentially lead to stellar populations
that exhibit notable variations in abundances of elements commonly found in
dust grains, including CNO, Si, Mg and Fe (Hopkins & Conroy, 2017).

Despite its well-recognized importance, most simulations that model star-
formation often assume a perfect coupling between the dust and the surrounding
gas. This treatment involves considering both components as moving together,
with dust essentially acting as an “additional opacity" to the gas. However,
the dynamics of dust is far more complex. Dust particles, ranging from a few
angstroms to several micrometres in size, are inherently aerodynamic and often
charged. When dust grains are accelerated by forces like radiation pressure,
they drift through the gas, and the imparted momentum couples to the
surrounding gas through electromagnetic and drag forces, redistributing
momentum in the environment. The efficiency of this coupling relies on both
the characteristics of the grains and the surrounding environment. For typical
densities of GMCs, where $n_{\rm H}\sim 10^{2}\,\rm cm^{-3}$, the collisional
interactions between dust and gas are relatively infrequent, resulting in a
relatively weak coupling of the dust dynamics with that of the gas. This
decoupling between grain and gas populations has notably been observed,
particularly in the case of larger dust grains (Krüger et al., 2001; Frisch &
Slavin, 2003; Meisel et al., 2002; Altobelli et al., 2006, 2007; Poppe et al.,
2010). To improve our understanding of dust dynamics within this environment,
this study serves as an initial investigation that relaxes some assumptions
regarding the perfect coupling of gas and dust dynamics.

In this work, we present the first radiation-dust-magnetohydrodynamic (RDMHD)
simulations of star-forming molecular clouds that explicitly account for dust
dynamics. We utilize simulations from the STAR FORmation in Gaseous
Environments (STARFORGE) project, which provide a comprehensive representation
of individual star formation. This encompasses the stages of proto-stellar
collapse, subsequent accretion and stellar feedback, main-sequence evolution,
and stellar dynamics, with a thorough consideration of the relevant physical
processes (as described in Grudić et al. (2022)). Within the context of this
paper, our focus is directed towards investigating the impacts of explicit
dust-gas radiation dynamics on the properties of the stellar populations that
emerge within the cloud.

The paper is structured as follows: In Section 2, a concise outline of the
code is presented (for a detailed description of the numerical methods and
implementations, refer to Grudić et al. (2021)) along with a description of
the initial conditions (ICs) of the simulations. We discuss the results
obtained from our fiducial simulation runs in Section 3 and compare them to
runs with simplified dust physics. Finally, we discuss the implications of our
findings in Section 4 and present our conclusions in Section 5.

##  2 Simulations

###  2.1 STARFORGE Physics

We run a set of RDMHD STARFORGE simulations to study star formation in GMCs,
using the GIZMO code (Hopkins, 2015). These simulations adopt the complete
physics setup from the “full” STARFORGE model, as described in Grudić et al.
(2022). For solving the equations of ideal magnetohydrodynamics (MHD), we rely
on the GIZMO Meshless Finite Mass MHD solver, detailed in Hopkins & Raives
(2016) and Hopkins (2016). In this study, we extend the standard STARFORGE
physics by explicitly modeling the dynamics of the dust particles, where the
details of the implementation are outlined in Section 2.2.

We utilize the GIZMO meshless frequency-integrated M1 solver (Lupi et al.,
2017, 2018; Hopkins & Grudić, 2019a; Hopkins et al., 2020; Grudić et al.,
2021; Hopkins & Grudić, 2019a; Hopkins et al., 2020) to evolve the time-
dependent, frequency-integrated radiative transfer (RT) equation adopting a
reduced speed of light $c=30\,\rm kms^{-1}$. The radiation is discretised into
five distinct frequency bands, which cover a range extending from the Lyman
continuum to the far infrared (FIR). As radiation interacts with matter, it
triggers processes such as photoionisation, photodissociation, photoelectric
heating, and dust absorption. Additionally, the matter undergoes radiative
cooling and heating, accounting for metal lines, molecular lines, fine
structure lines, as well as continuum and dust collisional processes, as
outlined in Hopkins et al. (2023). An especially noteworthy feature of our
radiative transfer approach is the direct coupling of each radiation band to
the dust particles. A more comprehensive description of this integration is
presented in Section 2.2. The radiation field is initialised with an external
heating source at the boundaries, representing the interstellar radiation
field (ISRF). As local sources (stars) emerge, they also contribute to the
radiation field.

Individual stars in the simulations are represented by sink particles. These
sink particles arise from gas cells that meet the criteria for runaway
gravitational collapse and follow the protostellar evolution model introduced
by McKee & Offner (2010). Sink particles can then accrete bound gas and dust
elements, where dust accretion adheres to the same critera as gas accretion.
Each sink particle separately tracks the quantities of both dust and gas it
has accumulated since its initial sink formation. To ensure realistic
accretion irrespective of resolution, the sink particles first accrete
particles onto an unresolved disk reservoir, the material in which is then
smoothly accreted onto the sink particle.

The sink particles undergo growth through accretion while progressing along
the main sequence. Their luminosity and radius follow the relationships
outlined in Tout et al. (1996), and they emit a black-body spectrum with an
effective temperature of $T_{\rm eff}=5780,\rm
K\left(L_{\star}/R_{\star}^{2}\right)^{1/4}$. Beyond radiation, sinks interact
with their surrounding medium via protostellar jets, stellar winds, and the
potential occurrence of supernovae. To calculate gravitational forces, we
employ an adaptive gravitational softening approach, which spans a range
extending down to approximately $\sim 2\times 10^{-5}$ AU for the gas cells.
Additionally, a fixed softening of $6$ AU is utilized for the dust and sink
particles. For a more comprehensive description of the modelled physics and
numerical techniques, readers are encouraged to refer to Grudić et al. (2021,
2022).

###  2.2 Dust Physics

A detailed description of the numerical methods used for dust modelling in our
simulations, albeit focussing on more idealised scenarios, can be found in
Hopkins & Lee (2016); Lee et al. (2017); Moseley et al. (2019); Soliman &
Hopkins (2023). Following a Monte Carlo sampling approach, we depict dust
grains in our simulation as "super-particles" (Carballido et al., 2008;
Johansen et al., 2009; Bai & Stone, 2010; Pan et al., 2011; McKinnon et al.,
2018). Each simulated “dust particle” encapsulates an ensemble of dust grains
characterized by similar attributes such as grain size
$(\epsilon_{\text{grain}})$, mass $(m_{\text{grain}})$, and charge
$(q_{\text{grain}})$ determined through live collisional, photoelectric, and
cosmic ray charging (Draine & Sutin, 1987; Tielens, 2005).

The motion of each dust grain is governed by the following equation:

$\displaystyle\frac{\mathrm{d}{\bf v}_{d}}{\mathrm{dt}}$ $\displaystyle={\bf
a}_{\rm gas,\,dust}+{\bf a}_{\rm grav}+{\bf a}_{\rm rad}$ (1)
$\displaystyle=-\frac{{\bf w}_{s}}{t_{s}}-\frac{{\bf w}_{s}\times\hat{\bf
B}}{t_{L}}+{\bf g}+\frac{\pi\,\epsilon_{\rm grain}^{2}}{m_{\rm
grain}\,c}\,\langle Q\rangle_{\rm ext}\,{\bf G}_{\rm rad},$

where ${\bf v}_{d}$ represents the velocity of the grain; ${\bf a}_{\rm
gas,\,dust}=-{{\bf w}_{s}}/{t_{s}}-{{\bf w}_{s}\times\hat{\bf B}}/{t_{L}}$
takes into account the forces exerted by the gas on the dust, including drag
(quantified by the "stopping time" $t_{s}$) and Lorentz forces (characterised
by the gyro/Larmor time $t_{L}\equiv m_{\rm grain}c/|q_{\rm grain},{\bf B}|$);
${\bf w}_{s}\equiv{\bf v}_{d}-{\bf u}_{g}$ corresponds to the drift velocity
of a dust grain relative to the gas velocity ${\bf u}_{g}$ at the same
position ${\bf x}$; ${\bf B}$ denotes the local magnetic field; ${\bf a}_{\rm
grav}={\bf g}$ is the gravitational force due to a local gravitational field;
and ${\bf a}_{\rm rad}$ is the force due to an incident radiation field ${\bf
G}_{\rm rad}\equiv{\bf F}_{\rm rad}-{\bf v}_{d}\cdot(e_{\rm
rad}\,\mathbb{I}+\mathbb{P}_{\rm rad})$ in terms of the radiation flux/energy
density/pressure density ${\bf F}_{\rm rad}$, $e_{\rm rad}$, $\mathbb{P}_{\rm
rad}$ for a grain of size $\epsilon_{\rm grain}$ and mass $m_{\rm
grain}\equiv(4\pi/3)\,\bar{\rho}_{\rm grain}^{\,i}\,\epsilon_{\rm grain}^{3}$,
where $\bar{\rho}_{\rm grain}^{\,i}$ is the internal grain mass density, and
dimensionless absorption+scattering efficiency $\langle Q\rangle_{\rm ext}$;
$c$ is the speed of light.

Our radiation treatment closely follows the methodology described in Hopkins
et al. (2022); Soliman & Hopkins (2023), which employed a gray-band opacity
assumption. In our current study, we build upon this approach by introducing a
five-band opacity treatment for each radiation band evolved by the RT solver.
We determine an effective opacity, expressed as $\langle Q\rangle_{\rm
ext}(\epsilon_{\rm grain},\lambda_{\rm eff})=\min\left(2\pi\epsilon_{\rm
grain}/\lambda_{\rm eff},1\right)$, where $\lambda_{\rm
eff}\equiv\sqrt{\lambda_{\min}\lambda_{\max}}$, representing the geometric
mean of the wavelength boundaries for the respective radiation bands.

The radiation pressure force on the dust is determined through the M1
radiative transfer method described above, encompassing terms up to
$\mathcal{O}(v^{2}/c^{2})$: $\partial_{t}e_{\rm rad}+\nabla\cdot{\bf F}_{\rm
rad}=-R_{\rm dust}\,{\bf v}_{d}\cdot{\bf G}_{\rm rad}/c^{2}$,
$\partial_{t}{\bf F}_{\rm rad}+c^{2}\,\nabla\cdot\mathbb{P}_{\rm rad}=-R_{\rm
dust}\,{\bf G}_{\rm rad}$. Specifically, it involves equations that govern the
temporal evolution of quantities related to the radiation field, where the
absorption and scattering coefficients, denoted as $R_{\rm dust}$, are
calculated locally from the dust grain distribution.

However, resolving the photon mean free path in global simulations is not
always feasible (Hopkins & Grudić, 2019b; Krumholz, 2018). To address this
limitation, we apply a correction factor as outlined by Grudić et al. (2021)
to estimate radiation flux. Additionally, in this work we account for multiple
scattering, where photons absorbed by dust are presumed to be re-emitted
within the infrared band.

To examine the effects of this implementation, we conduct a series of
numerical experiments, systematically introducing and removing the correction
factors for each radiation band. Our analysis reveals that while these
adjustments can influence the results quantitatively, their effects generally
remain within the interquartile range observed in this study.

For all forces originating from the interaction between gas and dust, denoted
$a_{\rm gas,dust}$, an equal and opposite force acts on the gas (referred to
as “back-reaction”). For the physical conditions under consideration, the drag
experienced by the dust is modeled using the Epstein drag formulation,
expressed as:

$\displaystyle t_{s}$
$\displaystyle\equiv\sqrt{\frac{\pi\gamma}{8}}\frac{\bar{\rho}_{\rm
grain}^{\,i}\,\epsilon_{\rm
grain}}{\rho_{g}\,c_{s}}\,\bigg{(}1+\frac{9\pi\gamma}{128}\frac{|{\bf
w}_{s}|^{2}}{c_{s}^{2}}\bigg{)}^{-1/2},$ (2)

where $\gamma$ is the adiabatic index, $\rho_{g}$ is the gas density, and
$c_{s}$ is the local sound speed. Additionally, Coulomb drag is considered,
though it typically constitutes a minor correction.

Furthermore, the thermochemistry of the medium is influenced by the presence
of dust through both indirect effects, arising from extinction by grains, and
direct effects, through dust heating and cooling terms (refer to Grudić et al.
(2022) for comprehensive details). These effects depend on the DTG ratio
and/or mean grain size. Instead of assuming fixed values, we estimate these
parameters by interpolating the local distribution of dust grains to the gas
cells, ensuring a self-consistent consideration of the thermodynamics
involved. We examine the implications of our explicit dust treatment on the
thermodynamics of star-forming regions and delve into the specific influence
of varying dust properties, such as grain size, in more detail in Soliman &
Hopkins (in prep.).

The simulations at hand do not resolve the detailed physics of the generation
and collimation of protostellar jets. Instead, jets are introduced by spawning
high velocity gas cells within a narrow cone (see Grudić et al. (2022) for
implementation details). The cone’s radius encompasses a few gas cells within
the sink radius, and is thus poorly resolved. To mitigate potential spurious
effects, such as the ejection of dust grains by the spawned jet particles, we
temporarily disable the interaction between dust particles and newly spawned
cells.This involves interpolating the gas properties from non-jet cells only
to dust for a specific duration, allowing the jets to escape the poorly
resolved sink radius. We experimented with the duration of the decoupling, and
find that results converge as long as the decoupling persists until the jets
can escape the sink radius or become resolved. We emphasize that we adopt this
approach as a precautionary measure to avoid overestimating the extent of dust
evacuation. Allowing the interaction to proceed would only lead to further
dust ejection near sink particles. Consequently, this approach does not
qualitatively affect our findings; instead, it strengthens our conclusions.

###  2.3 Initial conditions

For our fiducial simulations, we consider an initially uniform-density
turbulent molecular cloud with a mass of $M_{\rm cloud}=2\times 10^{3}\rm
M_{\odot}$ and a radius of $R=3\rm pc$. The cloud is enveloped within a $30\rm
pc$ periodic box with an ambient medium with density $\sim 10^{3}$ times lower
than that of the cloud. The initial velocity distribution follows Gaussian
random field with with an initial virial parameter $\alpha_{\rm
turb}=5\sigma^{2}R/(3GM_{\rm cloud})=2$. The initial magnetic field
configuration is uniform, and set to establish a mass-to-flux ratio equivalent
to 4.2 times the critical value within the cloud. For typical gas cells, the
mass resolution in our simulations is $\Delta m\sim 10^{-3}\rm M_{\odot}$,
while for the grain super-particles, we refine it to $2.5\times 10^{-6}\rm
M_{\odot}$ ($4\times$ higher resolution for the dust). Furthermore, cells
associated with protostellar jets and stellar winds have a higher mass
resolution of $10^{-4}\rm M_{\odot}$. In addition, we extend our sample to
include a larger cloud configuration with $M=2\times 10^{4}\rm M_{\odot}$ and
radius $R=10\rm pc$ also enveloped within a 10$\times$ larger box filled with
diffuse material. For the larger cloud, we use a coarser mass resolution,
scaled down by a factor of 10.

The dust component is initialized with a net zero drift velocity relative to
the surrounding gas. It is spatially distributed according to an initially
uniform DTG ratio $\rho_{d}^{0}=\mu^{\rm dg}\rho_{g}^{0}$, where $\mu^{\rm
dg}=0.01$ reflects the galactic value. The grain sizes follow an empirical
Mathis, Rumpl, & Nordsieck (MRN) power-law distribution with a differential
number density $dN_{\rm d}/d\epsilon_{\rm grain}\propto\epsilon_{\rm
grain}^{-3.5}$ (Mathis et al., 1977). The distribution spans a dynamic range
of $\epsilon_{\rm grain}^{\rm max}=100\epsilon_{\rm grain}^{\rm min}$, with
$\epsilon_{\rm grain}^{\rm max}=0.1\rm\mu m$ for our fiducial simulations. We
adopt the classic MRN mixture of carbonaceous and silicate composition,
assuming a uniform internal density and composition across different grain
sizes. This corresponds to grains with a sublimation temperature of
approximately 1500 K and an internal density of $\bar{\rho}_{\rm
grain}^{\,i}\sim 2.25\,\rm g/cm^{3}$. We do not model grain growth or
destruction; consequently, we adhere to a fixed size distribution, maintaining
constant sizes for particles throughout the simulation.

We provide a summary of the initial conditions for all simulations discussed
in this work in Table LABEL:table.

Name |  $M_{\rm cloud}$ [$\rm M_{\odot}$] |  $R_{\rm cloud}$ [pc] |  $\epsilon_{\rm grain}^{\rm max}$ [$\rm\mu m$] |  $\Delta m$ [$\rm M_{\odot}$] | Notes
---|---|---|---|---|---
m2e3_0.1 | $2\times 10^{3}$ | 3 | 0.1 | $10^{-2}$ | fiducial run
m2e3_0.1_hires | $2\times 10^{3}$ | 3 | 0.1 | $10^{-3}$ |  $\times$ 10 finer resolution
m2e3_0.1_nomag | $2\times 10^{3}$ | 3 | 0.1 | $10^{-2}$ | no Lorentz Forces on grains
m2e3_0.1_norad | $2\times 10^{3}$ | 3 | 0.1 | $10^{-2}$ | no radiation pressure forces on grains
m2e3_0.1_pass | $2\times 10^{3}$ | 3 | 0.1 | $10^{-2}$ | grains only feel drag & do not exert a back-reaction force
m2e3_0.1_nofb | $2\times 10^{3}$ | 3 | 0.1 | $10^{-2}$ | no stellar winds or protostellar jets
m2e3_1 | $2\times 10^{3}$ | 3 | 1.0 | $10^{-2}$ |  $\times$ 10 larger grains
m2e3_10 | $2\times 10^{3}$ | 3 | 10 | $10^{-2}$ |  $\times 100$ larger grains
m2e3_1_hires | $2\times 10^{3}$ | 3 | 1.0 | $10^{-3}$ |  $\times$ 10 larger grains & $\times$10 finer resolution
m2e3_10_hires | $2\times 10^{3}$ | 3 | 10 | $10^{-3}$ |  $\times 100$ larger grains & $\times$10 coarser resolution
m2e4_0.1 | $2\times 10^{4}$ | 10 | 0.1 | $10^{-2}$ | fiducial run
m2e4_0.1_hires | $2\times 10^{4}$ | 10 | 0.1 | $10^{-3}$ |  $\times 10$ finer resolution
m2e4_0.1_nomag | $2\times 10^{4}$ | 10 | 0.1 | $10^{-2}$ | no Lorentz forces on grains
m2e4_0.1_norad | $2\times 10^{4}$ | 10 | 0.1 | $10^{-2}$ | no radiation pressure forces on grains
m2e4_0.1_pass | $2\times 10^{4}$ | 10 | 0.1 | $10^{-2}$ | grains only feel drag & do not exert a back-reaction force
m2e4_1 | $2\times 10^{4}$ | 10 | 1.0 | $10^{-2}$ |  $\times$10 larger grains
m2e4_10 | $2\times 10^{4}$ | 10 | 10 | $10^{-2}$ |  $\times$100 larger grains
Table 1: The initial conditions for the simulations used in this study. The
columns include: (1) Simulation name. (2) Cloud mass $M_{\rm cloud}$. (3)
Cloud Radius $R_{\rm cloud}$. (4) Maximum grain size $\epsilon_{\rm
grain}^{\rm max}$. (5) Mass resolution $\Delta m$. (6) Notes indicating the
main variations from the fiducial run.

##  3 Results

###  3.1 Cloud Morphology

In Figure 1, we present the 2D integrated surface density of gas, denoted as
$\Sigma_{\rm gas}$ displayed on the left and dust, denoted as $\Sigma_{\rm
dust}$ shown on the right. These visualizations are from our m2e3_0.1_hires
simulation, which employs our comprehensive physics model to simulate a cloud
with an initial mass of $M_{\rm cloud}\sim 2\times 10^{3}\rm M_{\odot}$
evolved for a duration of $t\sim 3.5$ Myrs. We note that the dust surface
density has been scaled by a factor of $1/\mu^{\rm dg}$ to facilitate a
straightforward comparison between the distributions of gas and dust. Sink
particles, representing stars in the system, are portrayed as circles with
their sizes proportional to their respective masses. The initially spherically
uniform cloud undergoes a process of gravitational collapse and fragmentation.
This leads to the emergence of a stellar cluster near the central region with
a small number of sinks dispersed at greater distances from the cluster’s
center. When we compare the dust and gas distributions on parsec scales, we
observe that the distribution of the dust aligns with that of the gas. This
indicates that, on large scales, the dynamics of dust and gas are effectively
well-coupled.

Figure 1: The 2D integrated surface density of gas $\Sigma_{\rm gas}$ (left)
and dust $\Sigma-{\rm dust}$ (right) at $t\sim 3.5$ Myrs for the
m2e3_0.1_hires simulation, corresponding to a cloud with $M_{\rm cloud}\sim
2\times 10^{3}\rm M_{\rm\odot}$, $\Delta m\sim 10^{-3}\rm M_{\rm\odot}$
resolution, and $\epsilon_{\rm grain}^{\rm max}=0.1\rm\mu m$. This run employs
our full physics package and accounts for explicit dust dynamics, including
physical phenomena such as dust drag, dust back-reaction, Lorentz forces
acting on dust grains, and the explicit computation of dust opacity from grain
particles. Note that the dust surface density is scaled by $1/\mu^{\rm dg}$
for ease of comparison. Sink particles, representing stars, are depicted as
circles, with their radius reflecting their mass. We note that both the gas
and the dust exhibit similar large-scale structural features. For structural
differences at smaller scales, refer to Figure 5.  Figure 2: The rolling
median of the accreted dust-to-gas (ADG) mass ratio $\mu^{\rm adg}$ calculated
for sink particles formed within a cloud of initial mass $M_{\rm
cloud}=2\times 10^{3}\rm M_{\odot}$. The results are computed within
logarithmic mass bins and normalised to the initial mean dust-to-gas (DTG)
ratio of the cloud $\langle\mu_{0}^{\rm dg}\rangle$. The darkly shaded and
lightly shaded regions correspond to the interquartile and interdecile ranges
respectively. The plot contrasts the full-physics runs with a run with limited
dust physics, where m2e4_0.1 represents the full-physics fiducial run,
m2e4_0.1_nomag corresponds to a run with neutral dust grains (no Lorentz
forces), m2e4_0.1_norad corresponds to a run without radiation pressure forces
on grains, and m2e4_0.1_pass represents a ’passive’ grain run with only drag
forces and no back-reaction forces on the gas. Subsolar mass stars exhibit a
scattered ADG ratio centered around the average value of $\sim\langle\mu^{\rm
dg}_{0}\rangle$. The fiducial run reveals a clear trend of reduced ADG ratios
at higher sink masses influenced by radiation forces on dust, which diminishes
when radiative pressure forces are disabled.

While dust and gas exhibit coherence on cloud-size scales, stars primarily
grow by accreting the local mixture of dust and gas, a process that determines
various final stellar properties. Consequently, variations in the local DTG
ratio can influence the properties of the final stellar population. Due to the
simulation’s limited spatial resolution for sink accretion ($\sim 10$ AU), the
fate of dust, whether it gets accreted by a star or is ejected, remains
ambiguous. To address this uncertainty, we calculate a more robust measure:
the accreted dust-to-gas (ADG) ratio, which represents the mass ratio of dust
and gas elements satisfying accretion at the current resolution. We
acknowledge that this measure has the potential for overestimation,
considering the possibility of dust expulsion if accretion were tracked at
smaller spatial scales.

###  3.2 Reduced Dust Accretion in Massive Stars

In Figure 2, we present the moving median ADG ratio for individual sink
particles, calculated within logarithmic stellar mass bins, for a stellar
population. We present the sinks formed with our larger cloud with a mass of
$M_{\rm cloud}=2\times 10^{4}\rm M_{\odot}$ to ensure a robust statistical
dataset. To discern the contributions of various physical processes to the
results, we compare the ADG ratios in our fiducial full-physics run with runs
involving limited physics. Specifically, we experiment with runs, including a
“passive” grain run where grains experience only drag forces, without Lorentz
or radiation pressure forces and without exerting any back-reaction forces. We
also examine a run in which grains experience both drag and Lorentz forces,
and in turn induce back-reaction forces on the gas. However, in this specific
run, the dust particles neither explicitly experience radiation forces nor is
the dust opacity explicitly computed directly from the dust distribution.
Instead, the opacity of the gas is calculated by assuming a constant DTG
ratio. We also consider a scenario where we include all dust physics except
Lorentz forces, i.e. the case of neutral dust grains. Finally, we contrast
these scenarios with a comprehensive full-physics fiducial run, which
encompasses all the aforementioned physical effects, along with explicit
coupling of radiation to the dust grains.

As shown in the figure, we note a larger dispersions in the ADG ratio among
sub-solar mass sinks for all simulation runs, with an average value converging
to $\sim 1$. This dispersion is primarily ascribed to the inherent limitations
in resolution, set at around $\sim 10^{-2}\,\rm M_{\odot}$, which results in
low-mass sinks accreting a relatively limited number of gas and dust elements,
and hence exhibiting greater stochastic noise.

For sinks of higher masses ($\geq 2-8\rm M_{\odot}$), a trend emerges: the ADG
ratio decreases with increasing stellar mass. This phenomenon has implications
for stellar properties and the abundance levels of elements commonly found in
dust, including CNO, Si, Mg, and Fe, which we anticipate to be significantly
reduced. These levels can plummet to nearly an order of magnitude lower than
the original solid-phase abundances observed within the cloud. Strikingly,
this trend becomes negligible when the dust opacity is not directly computed
from the dust distribution, and the grains are not subject to radiation
pressure forces. Furthermore, this pattern is virtually indiscernible when the
grains are treated as passive entities within the system. Hence, we infer that
this phenomenon is likely driven by radiation pressure which expels dust
grains away from the vicinity of the sink particle. Massive stars, emitting
intense radiation, exert substantial radiation pressure on the surrounding
dust and gas, thereby significantly affecting their dynamics. This effect is
expected to dominate the dynamics when radiation pressure surpasses
gravitational forces, a condition expected to be met at specific stellar
luminosities and/or masses.

In addition to the aforementioned tests, we conducted simulations to assess
the impact of various stellar feedback mechanisms, including protostellar
jets, stellar winds, and supernova feedback, by systematically enabling and
disabling each in turn. We also explored the influence of magnetic fields
through radiative hydrodynamical simulations, distinct from our tests
involving neutral grains. We found that the identified trend remains robust
across different choices of these model parameters.

###  3.3 Dust Evacuation

####  3.3.1 Toy Model

We can gain a simple understanding into whether a star would accrete a dust
grain by examining the local dust dynamics through a set of simple
approximations. Consider a scenario involving an inflowing shell of gas and
dust influenced by a combination of drag and radiation from the star pushing
the shell outward. The grains should reach a terminal velocity, where they
drift relative to the gas at a velocity $w_{s}$. For simplicity, let us assume
sub-sonic drift with $v_{\rm inflow}\sim\epsilon c_{s}$ where $\epsilon$
characterises the efficiency of gravitational in-fall with respect to the
sound speed, and that the system is in the Rayleigh limit, characterised by
$\langle Q\rangle_{\rm ext}\propto\epsilon_{\rm grain}$, allowing for a
closed-form analytic solution. The amount of dust accreted by a star would
then depend on the ratio of the terminal drift velocity of dust grains away
from the star to their velocity through the in-falling gas $v_{\rm inflow}$.
This relationship can be approximated as follows:

$\displaystyle\frac{w_{\rm s}}{v_{\rm inflow}}$ $\displaystyle\sim\frac{a_{\rm
rad}t_{s}}{\epsilon c_{s}}\sim\frac{\langle Q\rangle_{\rm ext}F_{\star}}{4\pi
c\epsilon\rho_{g}c_{s}^{2}}\sim\frac{\langle Q\rangle_{\rm ext}\,u_{\rm
rad}}{\epsilon u_{\rm thermal}}$ (3) $\displaystyle\sim\frac{\langle
Q\rangle_{\rm ext}L_{\star}}{4\pi cr^{2}\epsilon\rho_{g}c_{s}^{2}}$
$\displaystyle\sim 0.1\left(\frac{\langle Q\rangle_{\rm
ext}}{0.2}\right)\left(\frac{L_{\star}}{\rm
L_{\odot}}\right)\left(\frac{0.01\,\rm
pc}{r}\right)^{2}\left(\frac{1}{\epsilon}\right)\left(\frac{10^{3}\,\rm
cm^{-3}}{n_{\rm H}}\right)\left(\frac{20\rm K}{T}\right),$

where $F_{\star}=L_{\star}/r^{2}$ is the incident flux for a star of
luminosity $L_{\star}$ on a grain at some radial distance $r$ situated within
a region of number density $n_{\rm H}$. Assuming that $L_{\star}\propto
M_{\star}^{3.5}$, we determine that this ratio reaches unity when
$M_{\star}\sim 2\rm M_{\odot}$ given the parameters we consider above.

We emphasize to the readers that this model is intentionally simplistic,
designed to offer qualitative physical insight into the relevant parameters
for this process. The model does not incorporate detailed models of accretion,
density profiles of dust and/or gas, an accurate description of the object’s
luminosity, or turbulence in the accretion flow, among other factors, all of
which are expected to vary spatially and temporally as the star evolves. An
additional point worth highlighting is that the ratio does not need to exceed
1 for the ADG ratio to be lower than the DTG ratio of the cloud; even modest
values of outward dust drift relative to gas inflow would lead to reduced dust
accretion.

We present a test of this hypothesis in Figure 3, where we use the full
expression for $t_{s}$ from Equation 2 to account for the subsonic and
supersonic limits. On the left, we show the mean dust drift velocity,
calculated as the mean radial velocity of dust relative to a nearby sink
subtracted from the mean radial velocity of the gas computed for a discretized
3D grid. This is denoted as $\langle v_{\rm dust}\rangle-\langle v_{\rm
gas}\rangle\equiv w_{s}$, and the value is then normalized to the local sound
speed at the position of the shell. Additionally, each data point is color-
coded based on the mean stellar accretion rate, averaged over approximately 1
million years. We analyze the relationship between these parameters and the
value of $L_{\star}/n_{\rm H}r^{2}$ within the grid element.

As expected from the simple relationship described in Equation 3 plotted in
the dashed black line, the dust drift velocity is proportional to
$L_{\star}/n_{\rm H}r^{2}$, maintaining this proportionality until
$w_{s}/c_{s}\sim 1$. Beyond this point, as the flow transitions to the
supersonic regime, the correlation transitions to a square root relationship
with $L_{\star}/n_{\rm H}r^{2}$. We extend this assessment by examining the
relationship between mean dust drift velocity and the mass of the central sink
particle. Similarly, $w_{s}$ positively correlates with the star’s mass, or in
other words, with the time-integrated accretion rate, which, in turn, is
associated with the integrated luminosity of the sink. Additionally, we
observe a correlation with the time-averaged sink accretion rate, where
periods of high accretion are associated with elevated dust drift velocities
and a higher value of $L_{\star}/n_{\rm H}r^{2}$. This correlation is expected
since it corresponds to increased luminosity and a higher gas density in the
vicinity of the star.

Figure 3: Left: The mean dust drift velocity $\langle w_{s}\rangle=\langle
v_{\rm dust}\rangle-\langle v_{\rm gas}\rangle$ normalized to the sound speed
$c_{s}$ around sink particles with final $M_{\star}\gtrsim 2M_{\odot}$ across
time. Each data point, color-coded to represent the time-averaged mean
accretion rate over a 1 Myr period, corresponds to mean values measured within
a discretized 3D grid centered around the sinks. We plot $\langle v_{\rm
dust}\rangle-\langle v_{\rm gas}\rangle$, against the luminosity of the sink
particle in its proximity, divided by its gas number density and the square of
its radial distance from the sink, $L_{\star}/n_{\rm H}r^{2}$. The dashed
black line represents the theoretical prediction from Equation 3,
incorporating both subsonic $\left(\langle w_{s}\rangle\propto
L_{\star}/n_{\rm H}r^{2}\right)$ and supersonic regimes $\left(\langle
w_{s}\rangle\propto\sqrt{L_{\star}/n_{\rm H}r^{2}}\right)$. The dust drift
exhibits a positive correlation with $L_{\star}/n_{\rm H}r^{2}$, as
theoretically predicted. There is scatter away from the theoretical prediction
due to the model’s simplicity and the influence of other forces and turbulence
on the dynamics. Right: The mean drift velocity against the mass of the sink
particle. Note that the discreteness in stellar mass corresponds to the
simulation resolution. The dust drift positively correlates with the sink’s
mass and mean accretion rate, in agreement with predictions.

####  3.3.2 A Case Study

In Figure 4, we present a case study of this phenomenon, focusing on a region
around a selected sink particle. We chose this particular candidate after
examining several others, as it demonstrates the median behaviour with minimal
noise. We also carefully chose specific timesteps to minimise the impact of
noise on the data and to clarify the dust evacuation process. Each line in the
plot corresponds to a snapshot in time and is colour-coded accordingly. All
values are computed as the mean values within narrow radial shells centered
around the star. The top four plots represent the state of the system during
the initial creation of the dust-evacuated region, while the lower section
corresponds to the subsequent dispersal of the dust-evacuated region.

At early times, there is only a minimal reduction in the DTG near the sink
particle. For instance, at $10^{-4}$ pc from the sink, the DTG ratio is only
half of the average initial value within the cloud. As we move farther from
the sink, we observe an accumulation of dust, leading to a DTG ratio of
approximately $\mu^{\rm dg}/\mu^{\rm dg}_{0}\sim 2$, which gradually
approaches the mean value as we reach a distance of about 0.1 parsecs from the
sink.

As the sink continues to accrete matter and increase in mass, we observe
corresponding changes in the gas environment. The gas number density increases
from its peak value of $10^{7}\rm cm^{-3}$ to $10^{9}\rm cm^{-3}$, and the gas
temperature transitions from an initial value of $T\sim 20$K to $T\sim 100$K.
As material accretes onto the sink, the dust is entrained alongside the gas,
and the peak in the DTG ratio shifts closer to the sink. Additionally, we
observe that the peak DTG ratio increases. This effect is primarily attributed
to the gas density dropping more rapidly than the dust density in that
specific region. We propose that this effect is driven by radiation pressure
forces that induce a net outward motion of the dust, in contrast to the net
inflow of gas. This process gives rise to the formation of a dust-evacuated
region surrounding the star. As a consequence, a dust shell, approximately
$10^{-2}$ pc in thickness, emerges, where $\langle\mu^{\rm
dg}_{0}\rangle/\langle\mu^{\rm dg}_{0}\rangle\sim 5$. In the immediate
vicinity of the star, this effect is accompanied by values dropping as low as
$\langle\mu^{\rm dg}\rangle/\mu^{\rm dg}_{0}\sim 10^{-2}$ in the central
region.

In the four bottom plots, we present the same parameters during a phase
characterised by a decline in the accretion rate, as indicated by the sink
mass-time plot. We observe a reduction in the gas number density and
temperature over time, consistent with expectations during a reduced accretion
phase. Additionally, we note a gradual increase in the DTG ratio within the
innermost regions, attributable to the diminishing gas density near the
central region. The reduced gas number density also leads to weaker drag
forces experienced by the dust grains, resulting in the outward shift of the
peak DTG ratio and the reduction of its maximum value as the dust shell
disperses.

Figure 4: The temporal evolution of environmental properties detailing the
creation and subsequent dispersion of a dust-evacuated region around a
reference sink particle. We only show specific timesteps that most clearly
highlight the evolution of this process. Each line in the plot corresponds to
a given timestep and is color-coded accordingly. The upper four plots, in the
hot colormap, depict the formation of a dust-evacuated region during a high
accretion state, while the lower four, in the cold colormap, show its
dispersal during reduced accretion. The properties are evaluated within narrow
radial shells. Top Left: The mean dust-to-gas (DTG) $\langle\mu^{\rm
dg}(r)\rangle$ ratio normalized to the initial mean value of the cloud
$\langle\mu^{\rm dg}_{0}\rangle$. Top Right: The mean number density of the
gas $n_{\rm H}(r)$. Bottom Left: The mean gas temperature $T(r)$. Bottom
Right: The sink particle’s mass evolution. The plots illustrate the in-fall of
material towards the sink, raising the central gas density. Ongoing accretion
and heightened luminosity increases the gas temperature due to radiation. As
gas flows inward, it carries dust closer, but the enhanced luminosity
counteracts the motion, propelling dust outward to form a $\sim 10^{-2}$ pc
thick dustshell, with a peak DTG ratio of $\langle\mu^{\rm dg}\rangle\sim
5\langle\mu^{\rm dg}_{0}\rangle$. As accretion slows down, the gas density and
temperature decrease, while the DTG ratio rises in the inner regions, and the
peak shifts outward. The plots illustrate variations in $\mu^{\rm dg}$ around
a sink particle induced by radiation pressure, and highlight the development
of a dust-evacuated zone encircled by a dust-rich shell.

While we demonstrated that both dust and gas exhibit similar features on large
scales, we will now discuss the differences that arise on smaller scales. In
Figure 5, we transition from the overview of the gas surface density of the
entire cloud in our m2e3_0.1_hires run ($6\text{pc}\times 6\text{pc}\times
6\text{pc}$) to a reduced volume. Specifically, the middle panel focuses on a
smaller volume with dimensions of ($0.01\text{pc}\times 0.01\text{pc}\times
0.005\text{pc}$), presenting the integrated 2D surface density of the gas
through a slice in that region. Individual dust particles, color-coded to
denote their respective grain sizes, are overlaid on the gas distribution. The
right panel zooms in on an even more compact space with dimensions of
($0.005\text{pc}\times 0.005\text{pc}\,\times 0.0025\text{pc}$), showcasing
the mean DTG ratio within the volume. Both zoomed-in plots are centered around
a specifically chosen sink particle. It is worth noting that the selection of
this particular sink particle and snapshot is intentional, as they effectively
illustrate the dust evacuation and dust pile-up phenomena.

The plot highlights a key observation: on sub-parsec scales, particularly in
the proximity of stars, the spatial distribution of dust particles deviates
from that of the gas. Further as indicated by the plots in Figure 4, we
observe a dust-suppressed zone near the sink particle followed by dust-rich
region, indicative of the presence of a dust pile-up or a dust “shell”
enveloping these stellar objects.

Figure 5: The 2D integrated gas surface density of the m2e3_0.1_hires
simulation. Left: The gas surface density of the cloud, with white circles
representing sink particles, sized according to their mass. Middle: The gas
surface density centered around a specific sink particle, with individual dust
particles color-coded by grain size. Right: The dust-to-gas mass ratio (DTG)
normalized to the mean DTG of the cloud within a 0.01 pc region centered
around the sink particle. This figure illustrates the presence of a dust-
evacuated zone around a sink particle, featuring a dust pile-up or “shell”
surrounding the sink. Figure 6: The rolling median of the accreted dust-to-gas
(ADG) ratio $\mu^{\rm adg}$ normalized to the initial mean DTG ratio of the
cloud $\langle\mu_{0}^{\rm dg}\rangle$ for clouds with $M_{\rm cloud}\sim
2\times 10^{3}\rm M_{\odot}$ and $M_{\rm cloud}\sim 2\times 10^{4}\rm
M_{\odot}$. We show both clouds simulated at a resolution of $\Delta m\sim
10^{-2}\rm M_{\odot}$. Darkly shaded and lightly shaded regions indicate the
interquartile and interdecile percentile ranges, respectively. The relatively
consistent trend suggests that, within the mass range investigated here, the
cloud mass does not affect the trend. Figure 7: The normalized rolling median
of the accreted dust-to-gas ratio ($\mu^{\rm adg}$) in simulations with
different resolutions. Shading indicates interquartile and interdecile ranges.
Dark blue line shows a $2\times 10^{3}\rm M_{\odot}$ cloud at resolutions of
$\Delta m\sim 10^{-3}$ while the pink line shows the same cloud at $\Delta
m\sim 10^{-2}\rm M_{\odot}$. Higher resolution reduces scatter, and aligns
$\mu^{\rm adg}$ closer to $\langle\mu_{0}^{\rm dg}\rangle$ before declining at
higher sink masses.

###  3.4 Effects of Cloud Mass

We consider how this phenomenon is affected by cloud properties. In Figure 6,
we compare the ADG ratio in two distinct clouds, both simulated with our full
physics setup. Both clouds are simulated with a common resolution of $\Delta
m\sim 10^{-2}\rm M_{\odot}$, with one having a mass of $2\times 10^{3}\rm
M_{\odot}$ and the other $2\times 10^{4}\rm M_{\odot}$ The results obtained
reveal a consistent trend in the ADG ratio, which remains unaffected by
variations in cloud mass. Across this mass range, the influence of radiation
pressure, inherently a localised phenomenon, exerts comparable effects across
diverse cloud masses. Nevertheless, although this phenomenon primarily impacts
the ADG ratio of a specific star, its influence could extend to larger radii.
As more massive stars form, their high luminosities could also lead to a
reduction in dust content for neighbouring stars.

While the primary emphasis of this work centers on reduced ADG ratios for
high-mass stars, our simulations also unveil a subset of low-mass stars that
are rich in dust content. Specifically, within the subset of solar and
subsolar mass stars, we observe instances of heightened ADG ratios. However,
we refrain from attributing this phenomenon to a single cause. Instead, we
postulate that it likely stems from the intricate interplay between turbulent
motions and radiation pressure emanating from neighboring stars, leading to
the redistribution of dust. Hence, stars emerging from regions with dust over-
densities are prone to experiencing increased levels of dust accretion.

###  3.5 Effects of Simulation Resolution

In Figure 7, we present a comparative analysis of the ADG ratios within
stellar populations formed in a cloud with a mass of $M_{\rm cloud}\sim
2\times 10^{3}\rm M_{\odot}$, utilising simulations at different resolutions
($\Delta m\sim 10^{-3}$ and $\Delta m\sim 10^{-2}\rm M_{\odot}$). As
anticipated, our observations indicate that increasing the simulation
resolution decreases the scatter in the distribution. In addition, for the
higher resolution cloud, the ADG ratio tends to converge to $\mu^{\rm dg}_{0}$
for relatively low-mass stars (those with masses less than a few solar
masses). This aligns with the expectation that stars with lower luminosities
lack the necessary conditions to evacuate dust, and therefore would have ADG
ratios that mirror the average DTG ratio within the cloud. These improvements
can be attributed to the enhanced sampling capabilities, allowing for a more
precise resolution of accretion and sink formation.

However, as shown in the plot, our findings are sensitive to resolution. This
observation is unsurprising given the inherent challenges in accurately
modelling accretion around luminous stars (Krumholz et al., 2009; Krumholz &
Thompson, 2012; Rosen et al., 2016). Achieving an accurate depiction of
relative dust-to-gas accretion onto a star, accounting for associated
radiation effects, requires significantly higher resolutions than currently
feasible with our simulations.

To illustrate, considering a 1 solar mass star that would have accreted
approximately $\sim 0.01M_{\odot}$ of dust mass, assuming $\mu^{\rm dg}\sim
10^{-2}$, at the $10^{-2}$ resolution (bearing in mind that the dust is up-
sampled by a ratio of 4 times the gas resolution), this roughly corresponds
accretion of $\sim$ 4 dust particles. Consequently, we do not necessarily
anticipate convergent results at our low resolutions, particularly at low sink
masses. However, this approach still provides valuable insights when comparing
simulations conducted at the same resolution. We hope that our results
motivate detailed zoom-in simulations of dusty accretion around protostellar
objects and onto massive stars, to provide a more detailed study of the dust
evacuation phenomenon.

To evaluate the effects of resolution on our findings, we conducted a series
of idealized tests of singular Shu (1977) collapse scenarios at resolutions of
$10^{-2}$, $10^{-3}$, and $10^{-4}\rm M_{\odot}$. These tests track the
problem hydrodynamically, allowing sink particle formation with subsequent
accretion and radiation. To simplify our analysis, we exclude the effects of
magnetic fields and protostellar jets and stellar winds. We find that the
observed phenomena qualitatively persist across different resolutions;
however, convergence is not achieved. In particular, neither the spatial
extent of the dust evacuation zone nor the precise magnitude of the stellar
ADG ratio converge with increasing resolution.

Considering the complexities of radiation-limited accretion, it remains
uncertain whether the system, as modeled in our simulations, should converge
within our resolution range. Higher resolution simulations capture finer-scale
structures, unveiling more intricate dynamics around the sink which thereby
altering the magnitude and scale of dust evacuation. Additionally, the
mechanisms of dust and gas accretion, whether via spherical collapse or disk
accretion, and their dynamics are resolution-dependent and influenced by the
specific physics in our simulations.

###  3.6 Effects of Altered Grain Size Distributions

Figure 8: The rolling median of the accreted dust-to-gas (ADG) ratio
($\mu^{\rm adg}$) for sink particles formed within a cloud of initial mass
$M_{\rm cloud}=2\times 10^{3}\rm M_{\odot}$ (top), $M_{\rm cloud}=2\times
10^{4}\rm M_{\odot}$ (bottom) with different initial grain-size distributions
($\epsilon_{\rm grain}^{\rm max}=0.1\rm\mu m$, $\epsilon_{\rm grain}^{\rm
max}=1\rm\mu m$, and $\epsilon_{\rm grain}^{\rm max}=10\rm\mu m$). The darkly
shaded and lightly shaded regions represent the interquartile and interdecile
percentile ranges, respectively. The m2e4_10 run exhibits notable scatter and
a low count of formed sink particles. To account for this, we scatter the ADG
ratio for each sink particle. Refer to Section 3.5 for a discussion on
potential drivers of this scatter. Changes in the grain size do not drive
significant deviations from the reduced $\mu^{\rm adg}$ trend for high mass
stars, but larger grain sizes lead to fewer high-mass stars due to reduced
star formation efficiency.

An additional variable to consider is the size of the dust grains. At
extremely small grain sizes, the grains should closely trace the dynamics of
the gas, whereas at larger grain sizes, the grains will be entirely decoupled
from gas dynamics. In Figure 8, we explore the impact of various plausible
grain sizes within the cloud environment, incorporating grains with maximum
sizes of $\epsilon_{\rm grain}^{\rm max}\sim 10\rm\mu m$ and $\epsilon_{\rm
grain}^{\rm max}\sim 1\rm\mu m$, alongside our baseline assumption of
$\epsilon_{\rm grain}^{\rm max}\sim 0.1\rm\mu m$. This is conducted while
ensuring a fixed total dust mass within the simulation and maintaining a
dynamic range of $\sim 100$ in grain sizes. We show the results for our higher
resolution $M_{\rm cloud}\sim 2\times 10^{3}\rm M_{\odot}$ simulations, in
addition to a $M_{\rm cloud}\sim 2\times 10^{4}\rm M_{\odot}$ cloud at lower
resolution. While in principle, grain size could influence our results, our
findings for these more reasonable sizes are consistent with Equation 3, which
demonstrates no grain-size dependence. It is worth noting that the grain size
does have an impact on the total number of sink particles that form, with
larger grains leading to a reduced overall sink count. This effect and related
changes will be studied in detail in Soliman & Hopkins (in prep.).

However, we begin to see some deviations from the reported trend for our
largest grain simulations. Note that we present individual $\mu^{\rm adg}$
values for each sink particle in the simulation with $M_{\rm cloud}\sim
2\times 10^{4}\rm M_{\odot}$ and $\epsilon_{\rm grain}^{\rm max}=10\rm\mu m$,
opting against plotting the median due to a limited number of sink particles
and a substantial dispersion in their $\mu^{\rm adg}$. The dispersion is
likely influenced by several factors, including the relatively lower
resolution. However, we posit that the increased grain size inherently
contributes to a higher level of scatter. Recall that these simulation assume
a fixed a grain size spectrum, where individual grain particles maintain a
fixed size throughout the simulation. Consequently, as grain size increases,
the individual grain count decreases, while each grain becomes more massive.
This suggests that whether a particle undergoes accretion or expulsion has a
more pronounced impact on the overall dust mass accreted by a sink, especially
when compared to its smaller size/mass grain counterpart. Further, Equation 3
assumes that $\langle Q\rangle_{\rm ext}$ is proportional to
$2\pi\epsilon_{\rm grain}/\lambda$, where $\lambda$ denotes the radiation
wavelength. However, for grains larger than $2\pi/\lambda$, $\langle
Q\rangle_{\rm ext}$ approaches 1. As a result, $w_{s}\propto\epsilon_{\rm
grain}^{-1}$, introducing an inverse dependence on grain size into Equation 3.
Considering the case of clouds with $\epsilon_{\rm grain}^{\rm max}=10\rm\mu
m$, all grains fall within this regime for radiation with $\lambda\lesssim
600\rm\mu m$, and larger grains fall into this regime for even shorter
wavelengths. This wavelength range corresponds to the optical/UV range, where
stars—especially more massive and hotter ones—typically emit peak radiation.
Therefore, we would expect to see increases in $\mu^{\rm adg}$ for such
massive grains, although such large grains are not expected to be abundant in
most GMCs.

In reality, dust grains undergo processes such as coagulation, accretion,
sputtering, photodestruction, and shattering which result in both changes in
the grain size distribution and the overall dust mass. These mechanisms,
currently not captured in our simulations, could impact the interpretation of
our results. For instance, if dust sublimation proceeds efficiently and the
elements comprising the grains become well-mixed in the gas phase near the
stars, they may be incorporated into the stars, potentially erasing the
effects of dust evacuation on stellar abundances. However, the dust-evacuated
zones around the stars would still persist or undergo further evacuation. We
acknowledge these limitations and plan to address them in future work.

In Figure 9, we present the bivariate distribution showing the correlation
between the local 3D gas density ($\rho_{g}$) and the dust density
($\rho_{d}$). This distribution is weighted by the dust mass and computed at a
spatial resolution of approximately $10^{-2}$ pc for the different grain size
runs ($\epsilon^{\rm max}_{\rm grain}=0.1,1,10,\mu\text{m}$) at a time of
approximately $2\,t_{\rm dyn}$. For the smallest grains, the distribution
exhibits a relatively narrow distribution, primarily centered around the
theoretically predicted perfect coupling line $\rho_{d}=\mu^{\rm dg}\rho_{g}$.
However, as the grain size increases, the distribution broadens, signifying a
decreased level of coupling between the gas and dust components. Nevertheless,
on average, the fluid remains sufficiently well-coupled to prevent large
fluctuations in the DTG ratio.

Figure 9: Bivariate distribution of the local 3D gas density ($\rho_{g}$) and
dust density ($\rho_{d}$), weighted by dust mass, illustrating the probability
distribution around a grain at the resolution scale of approximately $10^{-2}$
pc. From top to bottom, we show the distribution for a run with $\epsilon_{\rm
grain}^{\rm max}=0.1,1,10\,\rm\mu m$ respectively. We show the $1-\sigma$
(green), $2-\sigma$ (navy), $3-\sigma$ (orange), $4-\sigma$ (plum) contours.
The diagonal dotted lines represent perfect dust-gas coupling
($\rho_{d}=\mu^{\rm dg}\rho_{g}$), while the horizontal line denotes uniform
dust density. The distribution becomes progressively broader indicative of
weaker coupling for larger grain sizes.

##  4 Discussion

###  4.1 Implications

Our findings, which reveal a diminished accretion of dust onto massive stars
under the influence of radiation pressure, bear notable implications for the
stellar abundances of CNO, Mg, Fe, and Si, elements known to preferentially
deplete onto dust grains. This phenomenon was previously suggested in the
literature as a potential explanation for the anomalous abundance ratios
observed in nearby stars (Melendez et al., 2009; Gustafsson, 2018a, b;
Mathews, 1967; Cochran & Ostriker, 1977). A distinctive feature of this model
is its localized effect, with the evacuation zone extending up to $\sim
10^{-3}$ pc or a few hundred AUs. Consequently, this preferential accretion
mechanism may contribute to elucidating the anomalous abundance ratios
observed in specific nearby stars and the variations in abundance ratios among
stellar pairs (Maia et al., 2014; Nissen et al., 2017; Ramírez et al., 2015;
Biazzo et al., 2015; Teske et al., 2016a, b; Saffe et al., 2017; Oh et al.,
2018, e.g.,).

Further, the observed decrease in dust content has significant implications
for the gas dynamics and thermochemistry near the star. Lower dust content,
associated with diminished opacities and cooling rates, would accelerate the
expansion rate of HII regions (Ali, 2021). It would likely result in a larger
ionized volume, owing to the reduced UV absorption by dust.

Moreover, protoplanetary disks would also potentially affected by this
phenomenon. The DTG ratio plays a crucial role in determining the efficiency
of dust radial drift, influencing the structural evolution of the disk (Toci
et al., 2021). Moreover, as dust constitutes the fundamental building block of
planets, diminished dust content is likely to influence the types of planets
that can form and their compositions.

###  4.2 Caveats

The current study has inherent limitations that warrant careful consideration.
A major constraint is the absence of a model for dust evolution; our approach
relies on assuming a constant grain size distribution and a constant total
dust mass throughout the simulation. Consequently, should the grain size
distribution undergo significant changes, specific conclusions within our
study may require reevaluation. In particular, if dust grains experience
substantial destruction, whether through sublimation, shattering and/or
sputtering, evacuated regions would likely persist or even intensify. However,
if the elements released from the grains transition to the gas phase and
continue to accrete onto stars, spatial variations in dust distribution might
not necessarily indicate corresponding variations in stellar abundances.

Furthermore, our fiducial grain size distribution incorporates nanometer-sized
grains. The existence of such small grains in these conditions is disputed,
given the potential for coagulation or destruction due to the conditions in
molecular clouds. Additionally, it is uncertain whether such grains can be
modeled aerodynamically and whether they would adhere to the charge and mass
scalings we assume. However, it is important to note that the results
presented here predominantly depend on the grains holding the most mass,
which, under an MRN spectrum, are the largest grains. Nevertheless, these
considerations prompt us to explore different ranges of grain size in Section
3.6 to assess how varying grain properties might influence our results.

An additional constraint in our study is the resolution limit, which is
discussed in more detail in Section 3.5. To ensure a statistically robust
sample size for identifying the phenomenon of dust evacuation, we conducted
global star formation simulations. However, the most interesting behavior
occurs at the resolution limit of our study. Achieving detailed predictions
for the implications of dust dynamics on stellar metallicities, protostellar
envelope properties, and planet formation requires more intricate small-scale
physics. Addressing these complexities realistically demands a substantial
increase in resolution, which we plan to explore in future work.

##  5 Conclusions

This study introduces a set of RDMHD simulations of star forming GMCs as part
of the STARFORGE project. These simulations encompass a detailed
representation of individual star formation, accretion, and feedback
mechanisms, while also explicitly considering the influence of the dynamics of
dust grains. Our investigation focuses on the implications of these dynamics,
specifically radiation-dust interactions, on the emergent properties of
stellar populations.

Through our analysis, we find that when stars surpass a critical mass
threshold ($\sim 2M_{\odot}$), their luminosity exerts sufficient radiation
pressure on neighbouring dust grains, ultimately resulting in their expulsion
from the star’s accretion radius. This drives the formation of a dust-
evacuated region of size $\sim 100$ AU. Consequently, this process results in
a mass-dependent adjustment in the ADG mass ratio incorporated into these
stars via the accretion process. Furthermore, we investigate the potential
implications of varying cloud mass, grain sizes, and other physical parameters
such as decoupling the grains from magnetic fields and deactivating stellar
winds and jets. However, our findings indicated that these variations had
negligible effects within the parameter space of our study.

In summary, our findings shed light on the interplay between radiation, dust
dynamics, and star formation, offering valuable insights into the complex
processes that shape stars, their environments, and compositions within
molecular clouds.

Support for for NS and PFH was provided by NSF Research Grants 1911233,
20009234, 2108318, NSF CAREER grant 1455342, NASA grants 80NSSC18K0562, HST-
AR-15800. Support for MYG was provided by NASA through the NASA Hubble
Fellowship grant #HST-HF2-51479 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. Numerical
calculations were run on the TACC compute cluster “Frontera,” allocations
AST21010, AST20016, and AST21002 supported by the NSF and TACC, and NASA HEC
SMD-16-7592. This research is part of the Frontera computing project at the
Texas Advanced Computing Center. Frontera is made possible by National Science
Foundation award OAC-1818253.

## Data Availability Statement

The data supporting this article are available on reasonable request to the
corresponding author.

## References

  * Abergel et al. (2002) Abergel, A., Bernard, J. P., Boulanger, F., et al. 2002, A&A, 389, 239, doi: 10.1051/0004-6361:20020324
  * Ali (2021) Ali, A. A. 2021, Monthly Notices of the Royal Astronomical Society, 501, 4136
  * Altobelli et al. (2007) Altobelli, N., Dikarev, V., Kempf, S., et al. 2007, Journal of Geophysical Research (Space Physics), 112, 7105, doi: 10.1029/2006JA011978
  * Altobelli et al. (2006) Altobelli, N., Grün, E., & Landgraf, M. 2006, A&A, 448, 243, doi: 10.1051/0004-6361:20053909
  * Bai & Stone (2010) Bai, X.-N., & Stone, J. M. 2010, ApJ, 722, L220, doi: 10.1088/2041-8205/722/2/L220
  * Biazzo et al. (2015) Biazzo, K., Gratton, R., Desidera, S., et al. 2015, Astronomy & Astrophysics, 583, A135
  * Boogert et al. (2013) Boogert, A., Chiar, J., Knez, C., et al. 2013, The Astrophysical Journal, 777, 73
  * Carballido et al. (2008) Carballido, A., Stone, J. M., & Turner, N. J. 2008, MNRAS, 386, 145, doi: 10.1111/j.1365-2966.2008.13014.x
  * Cochran & Ostriker (1977) Cochran, W. D., & Ostriker, J. P. 1977, Astrophysical Journal, vol. 211, Jan. 15, 1977, pt. 1, p. 392-399. Research supported by the Aspen Center for Physics, 211, 392
  * Dorschner (2003) Dorschner, J. 2003, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 609, Astromineralogy; University Observatory Schillergässchen 3, D-07745 Jena, Germany, ed. T. K. Henning, 1–54
  * Draine & Lee (1984) Draine, B., & Lee, H. M. 1984, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 285, Oct. 1, 1984, p. 89-108. Research supported by the Alfred P. Sloan Foundation., 285, 89
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Draine & Sutin (1987) Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803, doi: 10.1086/165596
  * Flagey et al. (2009) Flagey, N., Noriega-Crespo, A., Boulanger, F., et al. 2009, The Astrophysical Journal, 701, 1450
  * Frisch & Slavin (2003) Frisch, P. C., & Slavin, J. D. 2003, ApJ, 594, 844, doi: 10.1086/376689
  * Girichidis et al. (2020) Girichidis, P., Offner, S. S., Kritsuk, A. G., et al. 2020, Space Science Reviews, 216, 1
  * Grudić et al. (2022) Grudić, M. Y., Guszejnov, D., Offner, S. S. R., et al. 2022, MNRAS, 512, 216, doi: 10.1093/mnras/stac526
  * Grudić et al. (2021) Grudić, M. Y., Guszejnov, D., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2021, Monthly Notices of the Royal Astronomical Society, 506, 2199, doi: 10.1093/mnras/stab1347
  * Gustafsson (2018a) Gustafsson, B. 2018a, Astronomy & Astrophysics, 616, A91
  * Gustafsson (2018b) —. 2018b, Astronomy & Astrophysics, 620, A53
  * Höfner & Olofsson (2018) Höfner, S., & Olofsson, H. 2018, A&A Rev., 26, 1, doi: 10.1007/s00159-017-0106-5
  * Hopkins (2014) Hopkins, P. F. 2014, ApJ, 797, 59, doi: 10.1088/0004-637X/797/1/59
  * Hopkins (2015) —. 2015, MNRAS, 450, 53, doi: 10.1093/mnras/stv195
  * Hopkins (2016) —. 2016, MNRAS, 462, 576, doi: 10.1093/mnras/stw1578
  * Hopkins & Conroy (2017) Hopkins, P. F., & Conroy, C. 2017, ApJ, 835, 154, doi: 10.3847/1538-4357/835/2/154
  * Hopkins & Grudić (2019a) Hopkins, P. F., & Grudić, M. Y. 2019a, MNRAS, 483, 4187, doi: 10.1093/mnras/sty3089
  * Hopkins & Grudić (2019b) —. 2019b, MNRAS, 483, 4187, doi: 10.1093/mnras/sty3089
  * Hopkins et al. (2020) Hopkins, P. F., Grudić, M. Y., Wetzel, A., et al. 2020, MNRAS, 491, 3702, doi: 10.1093/mnras/stz3129
  * Hopkins & Lee (2016) Hopkins, P. F., & Lee, H. 2016, MNRAS, 456, 4174, doi: 10.1093/mnras/stv2745
  * Hopkins & Raives (2016) Hopkins, P. F., & Raives, M. J. 2016, MNRAS, 455, 51, doi: 10.1093/mnras/stv2180
  * Hopkins et al. (2022) Hopkins, P. F., Rosen, A. L., Squire, J., et al. 2022, Monthly Notices of the Royal Astronomical Society, 517, 1491
  * Hopkins et al. (2023) Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154, doi: 10.1093/mnras/stac3489
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75, doi: 10.1088/0004-637X/704/2/L75
  * King & Pounds (2015) King, A., & Pounds, K. 2015, Annual Review of Astronomy and Astrophysics, 53, 115
  * Krause et al. (2020) Krause, M. G. H., Offner, S. S. R., Charbonnel, C., et al. 2020, Space Sci. Rev., 216, 64, doi: 10.1007/s11214-020-00689-4
  * Krüger et al. (2001) Krüger, H., et al. 2001, Planetary and Space Science, 49, 1303, doi: 10.1016/S0032-0633(01)00054-X
  * Krumholz (2018) Krumholz, M. R. 2018, MNRAS, 480, 3468, doi: 10.1093/mnras/sty2105
  * Krumholz et al. (2009) Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S., & Cunningham, A. J. 2009, Science, 323, 754
  * Krumholz & Thompson (2012) Krumholz, M. R., & Thompson, T. A. 2012, ApJ, 760, 155, doi: 10.1088/0004-637X/760/2/155
  * Larson (1981) Larson, R. B. 1981, Monthly Notices of the Royal Astronomical Society, 194, 809
  * Lee et al. (2017) Lee, H., Hopkins, P. F., & Squire, J. 2017, MNRAS, 469, 3532, doi: 10.1093/mnras/stx1097
  * Li & Draine (2001) Li, A., & Draine, B. 2001, The Astrophysical Journal, 554, 778
  * Lupi et al. (2018) Lupi, A., Bovino, S., Capelo, P. R., Volonteri, M., & Silk, J. 2018, MNRAS, 474, 2884, doi: 10.1093/mnras/stx2874
  * Lupi et al. (2017) Lupi, A., Volonteri, M., & Silk, J. 2017, MNRAS, 470, 1673, doi: 10.1093/mnras/stx1313
  * Mac Low & Klessen (2004) Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of modern physics, 76, 125
  * Maia et al. (2014) Maia, M. T., Meléndez, J., & Ramírez, I. 2014, The Astrophysical Journal Letters, 790, L25
  * Mathews (1967) Mathews, W. G. 1967, The Astrophysical Journal, 147, 965
  * Mathis (1990) Mathis, J. S. 1990, Annual Review of Astronomy and Astrophysics, 28, 37
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * McKee & Offner (2010) McKee, C. F., & Offner, S. S. R. 2010, ApJ, 716, 167, doi: 10.1088/0004-637X/716/1/167
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, Annu. Rev. Astron. Astrophys., 45, 565
  * McKinnon et al. (2018) McKinnon, R., Vogelsberger, M., Torrey, P., Marinacci, F., & Kannan, R. 2018, MNRAS, 478, 2851, doi: 10.1093/mnras/sty1248
  * Meisel et al. (2002) Meisel, D. D., Janches, D., & Mathews, J. D. 2002, ApJ, 579, 895, doi: 10.1086/342919
  * Melendez et al. (2009) Melendez, J., Asplund, M., Gustafsson, B., & Yong, D. 2009, The Astrophysical Journal, 704, L66
  * Minissale et al. (2016) Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, Astronomy & Astrophysics, 585, A24
  * Moseley et al. (2019) Moseley, E. R., Squire, J., & Hopkins, P. F. 2019, MNRAS, 489, 325, doi: 10.1093/mnras/stz2128
  * Murray et al. (2005) Murray, N., Quataert, E., & Thompson, T. A. 2005, The Astrophysical Journal, 618, 569
  * Nissen et al. (2017) Nissen, P., Aguirre, V. S., Christensen-Dalsgaard, J., et al. 2017, Astronomy & Astrophysics, 608, A112
  * Oh et al. (2018) Oh, S., Price-Whelan, A. M., Brewer, J. M., et al. 2018, The Astrophysical Journal, 854, 138
  * Pan et al. (2011) Pan, L., Padoan, P., Scalo, J., Kritsuk, A. G., & Norman, M. L. 2011, ApJ, 740, 6, doi: 10.1088/0004-637X/740/1/6
  * Poppe et al. (2010) Poppe, A., James, D., Jacobsmeyer, B., & Horányi, M. 2010, Geophysical Research Letters, 37, 11101, doi: 10.1029/2010GL043300
  * Ramírez et al. (2015) Ramírez, I., Khanal, S., Aleo, P., et al. 2015, The Astrophysical Journal, 808, 13
  * Rosen et al. (2016) Rosen, A. L., Krumholz, M. R., McKee, C. F., & Klein, R. I. 2016, Monthly Notices of the Royal Astronomical Society, 463, 2553
  * Saffe et al. (2017) Saffe, C., Jofré, E., Martioli, E., et al. 2017, Astronomy & Astrophysics, 604, L4
  * Salpeter (1977) Salpeter, E. 1977, Annual Review of Astronomy and Astrophysics, 15, 267
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488, doi: 10.1086/155274
  * Soliman & Hopkins (2023) Soliman, N. H., & Hopkins, P. F. 2023, Monthly Notices of the Royal Astronomical Society, 525, 2668
  * Spitzer Jr (2008) Spitzer Jr, L. 2008, Physical processes in the interstellar medium (John Wiley & Sons)
  * Teske et al. (2016a) Teske, J. K., Khanal, S., & Ramírez, I. 2016a, The Astrophysical Journal, 819, 19
  * Teske et al. (2016b) Teske, J. K., Shectman, S. A., Vogt, S. S., et al. 2016b, The Astronomical Journal, 152, 167
  * Thoraval et al. (1997) Thoraval, S., Boisse, P., & Duvert, G. 1997, A&A, 319, 948
  * Thoraval et al. (1999) Thoraval, S., Boissé, P., & Duvert, G. 1999, A&A, 351, 1051
  * Tielens (2005) Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge, UK: Cambridge University Press)
  * Toci et al. (2021) Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, Monthly Notices of the Royal Astronomical Society, 507, 818
  * Tout et al. (1996) Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257, doi: 10.1093/mnras/281.1.257
  * Watanabe & Kouchi (2008) Watanabe, N., & Kouchi, A. 2008, Progress In Surface Science, 83, 439
  * Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. 2001, The Astrophysical Journal, 563, 842
  * Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJS, 134, 263, doi: 10.1086/320852
  * Whittet et al. (1993) Whittet, D., Millar, T., & Williams, D. 1993, Dust and Chemistry in Astronomy, Institute of Physics Publishing, Bristol and Philadelphia