The effect of radiation pressure on spatial distribution of dust inside HII regions
MMNRAS , 1–10 (2015) Preprint 16 October 2018 Compiled using MNRAS L A TEX style file v3.0
The effect of radiation pressure on spatial distribution ofdust inside H II regions Shohei Ishiki, (cid:63) Takashi Okamoto, and Akio K. Inoue Department of Cosmoscience, Hokkaido University, N10 W8, Kitaku, Sapporo, 060-0810, Japan Department of Environmental Science and Technology, Faculty of Design Technology,Osaka Sangyo University, 3-1-1 Nakagaito, Daito, Osaka 574-8530, Japan
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We investigate the impact of radiation pressure on spatial dust distribution insideH ii regions using one-dimensional radiation hydrodynamic simulations, which includeabsorption and re-emission of photons by dust. In order to investigate grain size effectsas well, we introduce two additional fluid components describing large and small dustgrains in the simulations. Relative velocity between dust and gas strongly depends onthe drag force. We include collisional drag force and coulomb drag force. We find that,in a compact H ii region, a dust cavity region is formed by radiation pressure. Resultingdust cavity sizes ( ∼ . ii region is strongly charged, relative velocity between dust and gas ismainly determined by the coulomb drag force. Strength of the coulomb drag force isabout 2-order of magnitude larger than that of the collisional drag force. In addition,in a cloud of mass 10 M (cid:12) , we find that the radiation pressure changes the grain sizedistribution inside H ii regions. Since large (0.1 µ m) dust grains are accelerated moreefficiently than small (0.01 µ m) grains, the large to small grain mass ratio becomessmaller by an order of magnitude compared with the initial one. Resulting dust sizedistributions depend on the luminosity of the radiation source. The large and smallgrain segregation becomes weaker when we assume stronger radiation source, sincedust grain charges become larger under stronger radiation and hence coulomb dragforce becomes stronger. Key words: radiative transfer – methods: numerical – ISM: clouds – H ii regions Radiation from young massive stars plays a crucial role instar forming regions, and its effect on spatial dust distribu-tion inside H ii regions is also non-negligible. O’dell & Hub-bard (1965) firstly observed dust inside the H ii region andmany other observations found dust in H ii regions (O’dellet al. 1966; Ishida & Kawajiri 1968; Harper & Low 1971).O’dell & Hubbard (1965) observationally estimated the dis-tribution of dust inside H ii regions, concluding that gas-to-dust mass ratio decreases as a function of distance from thecentre of the nebulae. Nakano et al. (1983) and Chini et al.(1987) observationally suggested the existence of dust cavityregions. There have been some theoretical attempts to revealdust distribution inside H ii regions (Mathews 1967; Gail &Sedlmayr 1979a,b). Gail & Sedlmayr (1979b) suggested thata dust cavity can be created by radiation pressure. (cid:63) E-mail: [email protected]
Radiation pressure may also produce spatial variationsin the grain size distribution inside H ii regions as sug-gested by recent observational data of IR bubbles. From theGalactic Legacy Ingrared Mid-Plane Survey Extraordinaire(GLIMPSE; Benjamin et al. 2003), Churchwell et al. (2006)found that about 25% of IR bubbles are associated withknown H ii regions and they claimed that the IR bubblesare primarily formed around hot young stars. Deharvenget al. (2010) then pointed out that 86% of IR bubbles areassociated with ionzed gas. Since Churchwell et al. (2006)missed the large ( >
10 arcmin) and small ( < µ m continuum emission appears further from radi-ation source than that of 8 µ m continuum emission. Sincethey assumed that 250 µ m continuum emission traces thebig grains (BGs) and 8 µ m continuum emission traces thepolycyclic aromatic hydrocarbons (PAHs), they argued that c (cid:13) a r X i v : . [ a s t r o - ph . GA ] N ov S. Ishiki, T. Okamoto & A. K. Inoue the dust size distribution depends on the distance from aradiation source.Inoue (2002) argued the presence of the central dust de-pleted region — dust cavity — in compact/ultra-compactH ii regions in the Galaxy by comparing the observedinfrared-to-radio flux ratios with a simple spherical radia-tion transfer model. The dust cavity radius is estimated tobe 30% of the Stromgren radius on average, which is toolarge to be explained by dust sublimation. The formationmechanism of the cavity is still an open question, while theradiation pressure and/or the stellar wind from the excita-tion stars have been suggested as responsible mechanisms.We will examine whether the radiation pressure can producethe cavity in this paper. By considering the effect of radi-ation pressure on dust and assuming steady H ii regions,Draine (2011) theoretically explained the dust cavity sizethat Inoue (2002) estimated from observational data.Akimkin et al. (2015, 2017) estimated dust size dis-tribution by solving motion of dust and gas respectively,and they concluded that radiation pressure preferentially re-moves large dust from H ii regions. Their simulations have,however, assumed a single OB star as a radiation source. Asmentioned by Akimkin et al. (2015), grain electric potentialis the main factor that affects the dust size distribution. Ifwe assume a stronger radiation source, such as a star cluster,dust would been more strongly charged and their conclusionsmight change.In this paper, we investigate the effect of radiation pres-sure on spatial dust distribution inside compact H ii re-gions and compare it with the observational estimates (In-oue 2002). In addition, we perform multi-dust-size simula-tions and study the effect of the luminosity of the radiationsource on dust size distribution inside H ii regions.The structure of this paper is as follows: In Section 2, wedescribe our simulations. In Section 3, we describe our sim-ulation setup. In Section 4, we present simulation results. InSection 5, we discuss the results and present our conclusions. We place a radiation source at the centre of a sphericallysymmetric gas distribution. The species we include in oursimulations are H i , H ii , He i , He ii , He iii , electrons, anddust. We assume the dust-to-gas mass ratio to be 6 . × − corresponding to a half of the abundance of elements heavierthan He (so-called ’metal’) in the Sun (Asplund et al. 2009).We neglect gas-phase metal elements in this paper. We solvethe radiation hydrodynamic equations at each timestep asfollows:step 1 Hydrodynamic equationsstep 2 Radiative transfer and other related processessubstep 2.1 Static radiative transfer equationssubstep 2.2 Chemical reactionssubstep 2.3 Radiative heating and coolingsubstep 2.4 Grain electric potentialThe methods we use for radiation transfer, chemical reac-tions, radiative heating, cooling and time stepping are thesame as Ishiki & Okamoto (2017, hereafter paper I). We include absorption and thermal emission of photons bydust grains in our simulations. To convert the dust mass den-sity to the grain number density, we assume a graphite grainwhose material density is 2.26 g cm − (Draine & Salpeter1979). We employ the cross-sections of dust in Draine & Lee(1984) and Laor & Draine (1993) . Dust sizes we assume are0.1 µ m or 0.01 µ m. Dust temperature is determined by theradiative equilibrium, and thus the dust temperature is in-dependent from gas temperature. We assume that the dustsublimation temperature is 1500 K; however, dust never beheated to this temperature in our simulations. We do notinclude photon scattering by dust grains for simplicity. In our simulations, we solve hydrodynamics including thecoulomb drag force which depends on grain electric po-tential. In order to determine the grain electric potential,we consider following processes: primary photoelectric emis-sion, auger electron emission, secondary electron emission,and electron and ion collisions (Weingartner & Draine 2001;Weingartner et al. 2006). The effect of auger electron emis-sion and secondary electron emission is, however, almostnegligible in our simulations, because high energy photons( > eV) responsible for the two processes are negligiblein the radiation sources considered in this paper. Since thetime scale of dust charging processes is so small ( (cid:46) In our simulations, we calculate the effect of drag force F drag on a dust of charge Z d and radius a d (Draine & Salpeter1979) as follows: F drag = 2 πa kT g (cid:34)(cid:88) i n i (cid:0) G ( s i ) + z φ ln(Λ /z i ) G ( s i ) (cid:1)(cid:35) , where s i ≡ (cid:112) m i v / (2 kT g ) ,G ( s i ) ≈ s i / (3 √ π ) (cid:113) πs / ,G ( s i ) ≈ s i / (3 √ π/ s ) ,φ ≡ Z d e / ( a d kT g ) , Λ ≡ / (2 a d e | φ | ) (cid:112) kT g /πn e ,k is the Boltzmann constant, T g is the temperature of gas, n i is the number density of i th gas species, n e is the numberdensity of electron, z i is the charge of i th gas species ( i =H i , H ii , He i , He ii , He iii ), and m i is the mass of i th species.000
10 arcmin) and small ( < µ m continuum emission appears further from radi-ation source than that of 8 µ m continuum emission. Sincethey assumed that 250 µ m continuum emission traces thebig grains (BGs) and 8 µ m continuum emission traces thepolycyclic aromatic hydrocarbons (PAHs), they argued that c (cid:13) a r X i v : . [ a s t r o - ph . GA ] N ov S. Ishiki, T. Okamoto & A. K. Inoue the dust size distribution depends on the distance from aradiation source.Inoue (2002) argued the presence of the central dust de-pleted region — dust cavity — in compact/ultra-compactH ii regions in the Galaxy by comparing the observedinfrared-to-radio flux ratios with a simple spherical radia-tion transfer model. The dust cavity radius is estimated tobe 30% of the Stromgren radius on average, which is toolarge to be explained by dust sublimation. The formationmechanism of the cavity is still an open question, while theradiation pressure and/or the stellar wind from the excita-tion stars have been suggested as responsible mechanisms.We will examine whether the radiation pressure can producethe cavity in this paper. By considering the effect of radi-ation pressure on dust and assuming steady H ii regions,Draine (2011) theoretically explained the dust cavity sizethat Inoue (2002) estimated from observational data.Akimkin et al. (2015, 2017) estimated dust size dis-tribution by solving motion of dust and gas respectively,and they concluded that radiation pressure preferentially re-moves large dust from H ii regions. Their simulations have,however, assumed a single OB star as a radiation source. Asmentioned by Akimkin et al. (2015), grain electric potentialis the main factor that affects the dust size distribution. Ifwe assume a stronger radiation source, such as a star cluster,dust would been more strongly charged and their conclusionsmight change.In this paper, we investigate the effect of radiation pres-sure on spatial dust distribution inside compact H ii re-gions and compare it with the observational estimates (In-oue 2002). In addition, we perform multi-dust-size simula-tions and study the effect of the luminosity of the radiationsource on dust size distribution inside H ii regions.The structure of this paper is as follows: In Section 2, wedescribe our simulations. In Section 3, we describe our sim-ulation setup. In Section 4, we present simulation results. InSection 5, we discuss the results and present our conclusions. We place a radiation source at the centre of a sphericallysymmetric gas distribution. The species we include in oursimulations are H i , H ii , He i , He ii , He iii , electrons, anddust. We assume the dust-to-gas mass ratio to be 6 . × − corresponding to a half of the abundance of elements heavierthan He (so-called ’metal’) in the Sun (Asplund et al. 2009).We neglect gas-phase metal elements in this paper. We solvethe radiation hydrodynamic equations at each timestep asfollows:step 1 Hydrodynamic equationsstep 2 Radiative transfer and other related processessubstep 2.1 Static radiative transfer equationssubstep 2.2 Chemical reactionssubstep 2.3 Radiative heating and coolingsubstep 2.4 Grain electric potentialThe methods we use for radiation transfer, chemical reac-tions, radiative heating, cooling and time stepping are thesame as Ishiki & Okamoto (2017, hereafter paper I). We include absorption and thermal emission of photons bydust grains in our simulations. To convert the dust mass den-sity to the grain number density, we assume a graphite grainwhose material density is 2.26 g cm − (Draine & Salpeter1979). We employ the cross-sections of dust in Draine & Lee(1984) and Laor & Draine (1993) . Dust sizes we assume are0.1 µ m or 0.01 µ m. Dust temperature is determined by theradiative equilibrium, and thus the dust temperature is in-dependent from gas temperature. We assume that the dustsublimation temperature is 1500 K; however, dust never beheated to this temperature in our simulations. We do notinclude photon scattering by dust grains for simplicity. In our simulations, we solve hydrodynamics including thecoulomb drag force which depends on grain electric po-tential. In order to determine the grain electric potential,we consider following processes: primary photoelectric emis-sion, auger electron emission, secondary electron emission,and electron and ion collisions (Weingartner & Draine 2001;Weingartner et al. 2006). The effect of auger electron emis-sion and secondary electron emission is, however, almostnegligible in our simulations, because high energy photons( > eV) responsible for the two processes are negligiblein the radiation sources considered in this paper. Since thetime scale of dust charging processes is so small ( (cid:46) In our simulations, we calculate the effect of drag force F drag on a dust of charge Z d and radius a d (Draine & Salpeter1979) as follows: F drag = 2 πa kT g (cid:34)(cid:88) i n i (cid:0) G ( s i ) + z φ ln(Λ /z i ) G ( s i ) (cid:1)(cid:35) , where s i ≡ (cid:112) m i v / (2 kT g ) ,G ( s i ) ≈ s i / (3 √ π ) (cid:113) πs / ,G ( s i ) ≈ s i / (3 √ π/ s ) ,φ ≡ Z d e / ( a d kT g ) , Λ ≡ / (2 a d e | φ | ) (cid:112) kT g /πn e ,k is the Boltzmann constant, T g is the temperature of gas, n i is the number density of i th gas species, n e is the numberdensity of electron, z i is the charge of i th gas species ( i =H i , H ii , He i , He ii , He iii ), and m i is the mass of i th species.000 , 1–10 (2015) ffect of radiation pressure In this section we describe the procedure to solve the set ofhydrodynamic equations: ∂∂t ρ g + ∂∂x ρ g v g = 0 ∂∂t ρ d + ∂∂x ρ d v d = 0 ∂∂t ρ g v g + ∂∂x ρ g v = ρ g a gra + f rad , g − ∂∂x P g + K d ( v d − v g ) ∂∂t ρ d v d + ∂∂x ρ d v = ρ d a gra + f rad , d + K d ( v g − v d ) ∂∂t (cid:18) ρ g v + 12 ρ d v + e g (cid:19) + ∂∂x (cid:20)(cid:18) ρ g v + h g (cid:19) v g + 12 ρ d v (cid:21) = ( ρ g v g + ρ d v d ) a gra + f rad , g v g + f rad , d v d where ρ g is the mass density of gas, ρ d is the mass densityof dust, v g is the velocity of gas, v d is the velocity of dust, a gra is the gravitational acceleration, f rad , g is the radiationpressure gradient force on gas, f rad , d is the radiation pres-sure gradient force on dust, P g is the gas pressure, e g is theinternal energy of gas, h g is the enthalpy of gas, and K d isthe drag coefficient between gas and dust defined as follows: K d ≡ n d F drag | v d − v g | , where n d is the number density of dust grains.In order to solve the dust drag force stably, we use fol-lowing algorithm for the momentum equations: (cid:20) p ∗ d (cid:0) = ρ t+∆td v ∗ d (cid:1) p ∗ g (cid:0) = ρ t+∆tg v ∗ g (cid:1)(cid:21) = (cid:20) p t d p tg (cid:21) + (cid:20) F p , d ( ρ t d , v t d ) F p , g ( ρ tg , v t g , e tg ) (cid:21) ∆ t, (1) (cid:20) p t+∆td p t+∆tg (cid:21) = (cid:20) ρ t+∆td ρ t+∆tg (cid:21) p ∗ d + p ∗ g ρ t+∆td + ρ t+∆t g + (cid:20) ρ t+∆td ρ t+∆tg (cid:21) (cid:20) a gra + f d + f g ρ t+∆td + ρ t+∆t g (cid:21) ∆ t + (cid:20) − (cid:21) ρ t+∆t g ρ t+∆t d ρ t+∆t g + ρ t+∆t d (cid:0) v ∗ g − v ∗ d (cid:1) e − ∆ ttd + (cid:20) − (cid:21) t d ρ t+∆t d f g − ρ t+∆t g f d ρ g + ρ d (1 − e − ∆ ttd ) , (2)where ∆ t is the time step, ρ ti is the mass density of i th speciesat time t , p ti is the momentum of i th species at time t , e tg is the internal energy of gas at time t , F X , i is the advectionof the physical quantity X of the i th species, f d is the forceon dust ( f d = f rad , d ), f g is the force on gas ( f g = f rad , g − ∂P g /∂x ), and the inverse of the drag stopping time, t d , is t − = K d ρ t+∆td + ρ t+∆tg ρ t+∆td ρ t+∆tg . (3)Equation (2) that determines the relative velocity betweendust and gas is the exact solution of the following equations: ρ g ddt v g = f rad , g − ∂∂x P g + ρ g a gra + K d ( v d − v g ) ,ρ d ddt v d = f rad , d + ρ d a gra + K d ( v g − v d ) . (4)Momentum advection and other hydrodynamic equa-tions are solved by using AUSM + (Liou 1996). We solve thehydrodynamics in the second order accuracy in space andtime. In order to prevent cell density from becoming zeroor a negative value, we set the minimum number density, n H (cid:39) − cm − . We have confirmed that our results arenot sensitive to the choice of the threshold density as longas the threshold density is sufficiently low.In order to investigate whether our method is reliable,we perform shock tube tests in Appendix A. ln AppendixB, we describe how we deal with the dust grains with twosizes. In the first simulation, in order to investigate whether oursimulation derives a consistent result with the observa-tional estimate for compact/ultra-compact H ii regions (In-oue 2002), we model a constant density cloud of hydrogennumber density 4 × cm − and radius 1.2 pc. As a ra-diation source, we place a single star (i.e. black body) atthe centre of the sphere. Since we are interested in the for-mation of a dust cavity, we neglect the gravity which doesnot affect the relative velocity between dust grains and gas(see equation (2)). We assume a single dust grain size in thissimulation.In the second set of simulations, in order to investigatethe effect of radiation pressure on the dust grain size dis-tribution inside a large gas cloud, we model a cloud as aBonner-Ebert sphere of mass 10 M (cid:12) and radius 17 pc. Asthe radiation source, we consider a single star (black body,BB) or a star cluster (a simple stellar population, SSP) andwe change the luminosity of the radiation source to inves-tigate the dependence of the dust size distribution on theluminosity of the radiation source. We compute its luminos-ity and spectral-energy distribution as a function of timeby using a population synthesis code, P´EGASE.2 (Fioc &Rocca-Volmerange 1997, 1999), assuming the Salpeter ini-tial mass function (Salpeter 1955) and the solar metallicity.We set the mass range of the initial mass function to be 0.1to 120 M (cid:12) .Materials at radius, r , feel the radial gravitational ac-celeration, a gra ( r ) = − G M star r ( r + r ) / − G M ( < r ) r where M( < r) represents the total mass of gas inside r and M star is the mass of the central radiation source, which is50 M (cid:12) for the single star case and 2 × or 2 × M (cid:12) for the two star cluster cases. Since the gravity from theradiation source has a non-negligible effect on simulationresults and causes numerical instability in the case of SSP,we need to introduce softening length, r soft . We set it to0.5 pc for the SSP. Since the gravity from a single star isnegligible effect on simulation results, we set 0 pc for thesingle star. MNRAS , 1–10 (2015)
S. Ishiki, T. Okamoto & A. K. Inoue
Table 1.
Initial conditions and numerical setup for the simulations. The radius of each cloud are shown as r cloud . The number densities, n H , n He , and n d indicate the initial number densities of hydrogen, helium, and dust in the innermost cell, respectively. Spatial distributionof dust and gas are indicated by ‘C’ (constant density) and ‘BE’ (Bonnor-Ebert sphere) in the 7th column. The spectrum of a radiationsource is shown in the 8th column, where ‘SSP’ and ‘BB’ respectively indicate a simple stellar population with solar metallicity and theblack body spectrum of given temperature. ˙ N ion represents the number of ionized photon emitted from a radiation source per unit time.The initial temperature of gas and dust are represented by T g and T d , respectively. The mass of a central radiation source is indicatedby M star . The number of dust sizes is shown in the 13th column. When ‘Gravity’ is on/off, we include/ignore gravitation force in thesimulations. Cloud r cloud n H n He n d , Large n d , Small
Distribution Source ˙ N ion T g T d M star Dust Gravity(pc) (cm − ) (cm − ) (10 − cm − ) (10 − cm − ) (10 s − ) (K) (K) (10 M (cid:12) )Cloud 1 1.2 4 × × × Following the dust size distribution of Mathis et al.(1977), so-called MRN distribution, we assume two dust sizein these simulations. We assume the initial number ratio oflarge to small dust as n d , Large : n d , Small = 1 : 10 . , where n d , Large and n d , Small are the number density of dustgrains of 0.1 µ m and 0.01 µ m in size, respectively.The details of initial conditions are listed in Table 1.We use linearly spaced 128 meshes in radial direction, 128meshes in angular direction, and 256 meshes in frequencydirection in all simulations to solve radiation hydrodynam-ics. We present density, gas temperature, dust-to-gas mass ratio,grain electric potential ( V d ≡ eZ d /a d ), and relative velocitybetween dust and gas as functions of radius in Fig 1. In thetop panels, the hydrogen number density is indicated by thered solid line. The number density of H ii is indicated by theblue dash-dotted line. The initial state of the simulation isshown by the black dotted line.The average electron number density within an H ii re-gion, n e , the H ii region radius, r H ii , the dust cavity ra-dius, r d , and the ratio between the radius of the dust cav-ity to the Str¨omgren radius, y d , obtained by our simulation( t = 0 .
42 Myr) and the observational estimate are shownin Table 2. We find that our simulation results are in broadagreement with the observational estimate. The dust cavi-ties, hence, could be created by radiation pressure. The pa-rameter y d obtained by the simulation is somewhat smallerthan the observational estimate. However, we could find abetter agreement if we tuned the initial condition such asthe gas density. In addition, the agreement would be betterif we included the effect of stellar winds, which was neglectedin this paper.Since dust inside the H ii region is strongly charged,relative velocity between dust and gas is determined bycoulomb drag force. Magnitude of the coulomb drag forceis about 2-order of magnitude larger than that of the colli-sional drag force. The relative velocity, thus, becomes largestwhen the dust charge is neutral.Grain electric potential gradually decreases with radius and then suddenly drops to negative value. Near the ion-ization front, the number of ionized photons decreases andhence collisional charging becomes important. This is thereason behind the sudden decrease of the grain electric po-tential. In the neutral region, there is no photon which is ableto ionize the gas and hence there is few electrons that collidewith dust grains. On the other hand, there are photons thatphotoelectrically charge dust grains. Therefore, the grainelectric potential becomes positive again at just outside ofthe H ii region. Then, the UV photons are consumed and theelectron collisional charging becomes dominant again in theneutral region. We present densities, gas temperature, dust-to-gas mass ra-tios for large and small grains, large-dust-to-small-dust massratios ( ρ d , Large /ρ d , Small ), the grain electric potential, and rel-ative velocity between dust and gas as functions of radius inFig 2. In order to compare the simulation results on the dustsize distribution with each other, we present the results atthe time when the shock front reaches to ∼
15 pc. In order tostudy the dependence of the dust size distribution on timeand the luminosity of the radiation source, we also presentthe simulation result of Cloud 2 at t =1.1 My: the sameirradiation time as Cloud 4. In the top panels, the hydro-gen number density is indicated by the red solid lines andthat of H ii is indicated by the blue dot-dashed lines. Initialstates of the simulations are shown by black dotted lines.In the fifth row, the charges of dust grains with size 0.1 µ mand 0.01 µ m are indicated by the red solid and blue dot-dashed lines, respectively. The black dotted lines show theinitial profiles (i.e. 0 V). In the bottom panels, the relativevelocity between dust grains with size 0.1 µ m and gas andbetween dust grains with size 0.01 µ m and gas are indicatedby the red solid and blue dot-dashed lines, respectively. Notethat the radiation source becomes stronger from Cloud 2 toCloud 4.We find that radiation pressure affects the dust distri-bution within an H ii region depending on the grain size. InFig 2, we divide them into the following four regions:(a) From the central part, radiation pressure removesboth large and small dust grains and creates a dust cavity(the yellow shaded region). MNRAS000
15 pc. In order tostudy the dependence of the dust size distribution on timeand the luminosity of the radiation source, we also presentthe simulation result of Cloud 2 at t =1.1 My: the sameirradiation time as Cloud 4. In the top panels, the hydro-gen number density is indicated by the red solid lines andthat of H ii is indicated by the blue dot-dashed lines. Initialstates of the simulations are shown by black dotted lines.In the fifth row, the charges of dust grains with size 0.1 µ mand 0.01 µ m are indicated by the red solid and blue dot-dashed lines, respectively. The black dotted lines show theinitial profiles (i.e. 0 V). In the bottom panels, the relativevelocity between dust grains with size 0.1 µ m and gas andbetween dust grains with size 0.01 µ m and gas are indicatedby the red solid and blue dot-dashed lines, respectively. Notethat the radiation source becomes stronger from Cloud 2 toCloud 4.We find that radiation pressure affects the dust distri-bution within an H ii region depending on the grain size. InFig 2, we divide them into the following four regions:(a) From the central part, radiation pressure removesboth large and small dust grains and creates a dust cavity(the yellow shaded region). MNRAS000 , 1–10 (2015) ffect of radiation pressure n H [ c m − ] Cloud 1 t =0.42My HI+HIIHIIinitial state T g [ K ] -3 -2 ρ d ρ g -1 -1 -10 V d [ V ] r [pc] -2 -1 v d − v g [ k m s − ] Figure 1.
Density (top-row), gas temperature (second-row), dust-to-gas mass ratio (third-row), grain electric potential (fourth-row),and relative velocity between dust and gas (bottom-row) profile at t = 0 .
42 Myr. The black dotted lines show the initial profiles. Thered solid lines represent the simulation results. The blue dashed lines in the top panel shows the ionized hydrogen density profile.
Table 2.
Comparison between the simulation results ( t = 0 .
42 Myr) and the observational estimates. The number densities of ionizedelectrons inside H ii regions are represented by n e . The number of ionized photons emitted from radiation sources per unit time isrepresented by ˙ N ion . The radius of H ii region is represented by r H ii . The radius of dust cavity radius is represented by r d . The parameter y d is defined by y d ≡ r d /R St , where R St is the Str¨omgren radius. Since, observationally, the number density is driven from the columndensity, the electron number density in the simulation is defined as by n e ≡ (cid:82) r H ii n e dx/r H ii . In addition, in order to match the definitionof r d with Inoue (2002), the dust cavity radius is defined by r d ≡ r H ii (cid:0) − (cid:82) r H ii ρ d dx/ (cid:82) r H ii ρ g ( ρ d /ρ g ) initial dx (cid:1) , where ( ρ d /ρ g ) initial represents the initial dust-to-gas mass ratio. n e ˙ N ion r H ii r d y d (cm − ) (10 s − ) (pc) (pc) ≡ r d /R St this work ( t = 0 .
42 Myr) 1247 6.2 0.73 0.21 0.20Inoue (2002) 1200 ±
400 6.8 ± ± ± ± (b) Within an H ii region, ρ d , Large /ρ d , Small has a peak.Between the region ‘a’ and this peak, there is a region where ρ d , Large /ρ d , Small takes the local minimum value (the cyanshaded region), for example, at r ∼ t = 2 . ρ d , Large /ρ d , Small is also reduced just behind theionization front (the gray shaded region), for example, at r ∼ t = 2 . MNRAS , 1–10 (2015)
S. Ishiki, T. Okamoto & A. K. Inoue n H [ c m − ] IF SF
Cloud 2 t =1.1My HI+HIIHIIinitial state T g [ K ] -4 -3 -2 ρ d , L a r g e ρ g -4 -3 -2 ρ d , S m a ll ρ g -1 ρ d , L a r g e ρ d , S m a ll a b cd -10 -1 -1 -10 V d [ V ] µ m µ m r [pc] -10 -1 -10 -2 -1 -2 -10 v d − v g [ k m s − ] µ m µ m − IF SF
Cloud 2 t =2.9My a b cd r [pc] IF SF Cloud 3 t =1.9My a b c d r [pc] IF SF Cloud 4 t =1.1My a b c d r [pc] Figure 2.
Density (top-row), gas temperature (second-row), large-dust-to-gas-mass ratio (third-row), small-dust-to-gas-mass ratio(fourth-row), large-dust-to-small-dust-mass ratio (fifth-row), grain electric potential (sixth-row), and relative velocity between dust andgas (bottom-row) profiles. From left to right, we show the results for Clouds 2 at t = 1 . t = 2 . t = 1 . t = 1 .000
Density (top-row), gas temperature (second-row), large-dust-to-gas-mass ratio (third-row), small-dust-to-gas-mass ratio(fourth-row), large-dust-to-small-dust-mass ratio (fifth-row), grain electric potential (sixth-row), and relative velocity between dust andgas (bottom-row) profiles. From left to right, we show the results for Clouds 2 at t = 1 . t = 2 . t = 1 . t = 1 .000 , 1–10 (2015) ffect of radiation pressure the radiation source becomes brighter (the regions ‘a’). Thereasons are as follows. grain electric potential of the dustgrains with the same size within r = 2 pc is almost the sameamong all simulations and the number density of the gasbecomes smaller for stronger radiation source. Since the dustdrag force strongly depends on the grain electric potential,the number density of gas, and radiation pressure on dust,relative velocity between dust and gas becomes larger forthe brighter source.In the region ‘b’ and ‘d’, the ratio ρ d , Large /ρ d , Small is de-creased from the initial condition when the radiation sourceis a single OB star (Cloud 2). Except in the ionization front(vertical brown dashed lines) which is contained in the re-gion ‘d’, radiation pressure preferentially removes large dustgrains from these regions. The photoelectric yield of the largedust grains is smaller than that of the small dust grains, andhence grain electric potential of the large dust grains be-comes smaller than that of the small dust grains. Coulombdrag between large dust grains and gas therefore becomesweaker than that between small dust grains and gas. On theother hand, since Cloud 4 has the strongest radiation sourceand hence it makes grain electric potentials largest amongthe simulations, the dust size segregation in the regions ‘b’and ‘d’ is less prominent. Even when we compare Cloud 2and 4 at the same irradiation time, t =1.1 Myr, the dust sizedistributions inside H ii regions are different. Luminosity ofthe radiation source must be the main cause of the dust sizesegregation.The ratio ρ d , Large /ρ d , Small in all simulations has a peakin the regions ‘c’. Since dust grains have large negativecharge in the regions ‘c’ and ‘d’, the coulomb drag forcebetween dust and gas is strong and hence dust and gas aretightly coupled each other. Large dust grains are, therefore,removed from the regions ‘a’ and ‘b’ and gathered in theregions ‘c’.At the ionization front and the shock front (verticalgreen dot-dot-dashed lines), the relative velocity, v d − v g hasdownward peaks. In theses fronts, gas pressure force exceedsradiation pressure force. Since the dust drag time dependson the dust grain size, the dust-gas relative velocity alsodepends on the grain size. As a result, ρ d , Large /ρ d , Small isslightly reduced in these fronts.
We have investigated radiation feedback in dusty clouds byone-dimensional multi-fluid hydrodynamic simulations. Inorder to study spatial dust distribution inside H ii regions,we solve gas and dust motion self-consistently. We also in-vestigate dust size distribution within H ii regions by con-sidering dust grains with two different sizes.We find that radiation pressure creates dust cavity re-gions. We confirm that the size of the dust cavity region isbroadly agree with the observational estimate (Inoue 2002).We also find that radiation pressure preferentially re-moves large dust from H ii regions in the case of a singleOB star. This result is almost the same as in Akimkin et al.(2015). The dust size distribution is, however, less affectedwhen the radiation source is a star cluster, in other word, amore luminous case. Resulting dust size distributions largelydepend on the luminosity of the radiation source. We assume dust is graphite. There are, however, otherforms of dust such as silicate. Since the photoelectric yieldand the absorption coefficient depend on a dust model, spa-tial dust distribution of dust grains may become differentwhen we use a different dust model. For example, since sili-cate has a larger work function and a smaller absorption co-efficient than graphite, the cavity size in the silicate case maybecome larger than that in the graphite case (see Akimkinet al. 2015, 2017 for details).In our simulations, we neglect the effect of sputteringthat changes the dust grain size. We estimate this accordingto Nozawa et al. (2006), and confirm that sputtering effectis negligible in our simulations. However, if we consider thesmaller dust grains, we may have to include the sputtering. ACKNOWLEDGEMENTS
We are grateful to Takashi Kozasa, Takashi Hosokawa, andShu-ichiro Inutsuka for helpful discussion. SI acknowledgesGrant-in-Aid for JSPS Research Fellow (17J04872) and TOacknowledge the financial support of MEXT KAKENHIGrant (16H01085). Numerical simulations were partly car-ried out with Cray XC30 in CfCA at NAOJ.
REFERENCES
Akimkin V. V., Kirsanova M. S., Pavlyuchenkov Y. N., WiebeD. S., 2015, MNRAS, 449, 440Akimkin V. V., Kirsanova M. S., Pavlyuchenkov Y. N., WiebeD. S., 2017, MNRAS, 469, 630Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A,47, 481Benjamin R. A., et al., 2003, PASP, 115, 953Chini R., Kruegel E., Wargau W., 1987, A&A, 181, 378Churchwell E., et al., 2006, ApJ, 649, 759Deharveng L., et al., 2010, A&A, 523, A6Draine B. T., 2011, ApJ, 732, 100Draine B. T., Lee H. M., 1984, ApJ, 285, 89Draine B. T., Salpeter E. E., 1979, ApJ, 231, 77Fioc M., Rocca-Volmerange B., 1997, A&A, 326, 950Fioc M., Rocca-Volmerange B., 1999, astro-ph, 9912179Gail H.-P., Sedlmayr E., 1979a, A&A, 76, 158Gail H. P., Sedlmayr E., 1979b, A&A, 77, 165Harper D. A., Low F. J., 1971, ApJ, 165, L9Inoue A. K., 2002, ApJ, 570, 688Ishida K., Kawajiri K., 1968, PASJ, 20, 95Ishiki S., Okamoto T., 2017, MNRAS, 466, L123Laor A., Draine B. T., 1993, ApJ, 402, 441Liou M.-S., 1996, Journal of Computational Physics, 129, 364Mathews W. G., 1967, ApJ, 147, 965Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425Nakano M., Kogure T., Sasaki T., Mizuno S., Sakka K., Wirami-hardja S. D., 1983, Ap&SS, 89, 407Nozawa T., Kozasa T., Habe A., 2006, ApJ, 648, 435O’dell C. R., Hubbard W. B., 1965, ApJ, 142, 591O’dell C. R., Hubbard W. B., Peimbert M., 1966, ApJ, 143, 743Paladini R., et al., 2012, ApJ, 760, 149Salpeter E. E., 1955, ApJ, 121, 161Simpson R. J., et al., 2012, MNRAS, 424, 2442Weingartner J. C., Draine B. T., 2001, ApJS, 134, 263Weingartner J. C., Draine B. T., Barr D. K., 2006, ApJ, 645, 1188MNRAS , 1–10 (2015)
S. Ishiki, T. Okamoto & A. K. Inoue
APPENDIX A: SHOCK TUBE TESTS
In order to investigate whether our method is reliable, weperform shock tube tests. Since the effect of the dust be-comes almost negligible in shock tube tests if we assumedust-to-gas mass ratio as 6.7 × − (the value we assumein our simulations) and hence we will not be able to in-vestigate whether the numerical code is reliable or not, weassume dust-to-gas mass ratio as 1 in the shock tube tests.The initial condition of the shock tube problem is as follows: ρ g = ρ d = (cid:40) , ( x < . , . , ( x > . ,P g = (cid:40) , ( x < . , . , ( x > . ,γ = 1 . , where γ is heat capacity ratio. Since the analytic solutionsare known for K d = 0 and ∞ , we perform test calculationsfor K d = 0 and K d = 10 (∆ t sim (cid:29) ( ρ g ρ d ) / ( ρ d + ρ g ) K − ≡ t d ), where ∆ t sim is the time scale of the shock tube problemand t d is the drag stopping time. We use linearly spaced400 meshes between x = 0 and 1. Time steps we use forthese simulations are ∆ t = 2 . × − for K d = 0 and ∆ t =4 . × − for K d = 10 . The results are shown in FigA1. We confirm that the numerical results agree with theanalytic solutions. APPENDIX B: DUST GRAINS WITH TWOSIZES AND GAS DYNAMICS
In order to investigate the spatial variation of the grain sizedistribution inside H ii regions, we solve following hydrody-namics equations, where we consider dust grains with twosizes (dust-1 and dust-2): ∂∂t ρ g + ∂∂x ρ g v g = 0 ∂∂t ρ d1 + ∂∂x ρ d1 v d1 = 0 ∂∂t ρ d2 + ∂∂x ρ d2 v d2 = 0 ∂∂t ρ g v g + ∂∂x ρ g v = ρ g a gra + f rad , g − ∂∂x P g + K d1 ( v d1 − v g ) + K d2 ( v d1 − v g ) ∂∂t ρ d1 v d1 + ∂∂x ρ d1 v = ρ d1 a gra + f rad , d1 + K d1 ( v g − v d1 ) ∂∂t ρ d2 v d2 + ∂∂x ρ d2 v = ρ d2 a gra + f rad , d2 + K d2 ( v g − v d2 ) ∂∂t (cid:18) ρ g v + e g (cid:19) + ∂∂t (cid:18) ρ d1 v + 12 ρ d2 v (cid:19) + ∂∂x (cid:18) ρ g v + h g (cid:19) v g + ∂∂x (cid:18) ρ d1 v + 12 ρ d2 v (cid:19) = ( ρ g v g + ρ d1 v d1 + ρ d2 v d2 ) a gra + f rad , g v g + f rad , d1 v d1 + f rad , d2 v d2 where ρ d1 is the mass density of dust-1, ρ d2 is the massdensity of dust-2, v d1 is the velocity of dust-1, v d2 is thevelocity of dust-2, f rad , d1 is the radiation pressure gradient force on dust-1, f rad , d2 is the radiation pressure gradientforce on dust-2, K d1 is the drag coefficient between gas anddust-1, and K d2 is the drag coefficient between gas and dust-2. In order to solve dust drag force stably, we use followingalgorithm for the equation of momentum: p ∗ d1 (cid:0) = ρ t+∆td1 v ∗ d1 (cid:1) p ∗ g (cid:0) = ρ t+∆tg v ∗ g (cid:1) p ∗ d2 (cid:0) = ρ t+∆td2 v ∗ d2 (cid:1) = p t d1 p tg p t d2 + F p , d1 ( ρ t d1 , v t d1 ) F p , g ( ρ tg , v t g , e tg ) F p , d2 ( ρ t d2 , v t d2 ) ∆ t, (B1) p t+∆td1 p t+∆t g p t+∆td2 = ρ t +∆ t d1 ρ t +∆ t g ρ t +∆ t d2 p ∗ d1 + p ∗ g + p ∗ d2 ρ t +∆ t d1 + ρ t +∆ tg + ρ t +∆ t d2 + ρ t+∆td1 ρ t+∆t g ρ t+∆td2 (cid:20) a gra + f d1 + f g + f d2 ρ t+∆td1 + ρ t+∆t g + ρ t+∆td2 (cid:21) ∆ t + ba + x cd + x ρ g e x ∆ t x ( x − y ) (cid:2) b ( d + x ) v ∗ d1 + ( a + x )( d + x ) v ∗ g + c ( a + x ) v ∗ d2 (cid:3) + ba + y cd + y ρ g e y ∆ t y ( y − x ) (cid:2) b ( d + y ) v ∗ d1 + ( a + y )( d + y ) v ∗ g + c ( a + y ) v ∗ d2 (cid:3) + ba + x cd + x e x ∆ t − x ( x − y ) [ a ( d + x ) f d1 + ( a + x )( d + x ) f g + d ( a + x ) f d2 ]+ ba + y cd + y e y ∆ t − y ( y − x ) [ a ( d + y ) f d1 + ( a + y )( d + y ) f g + d ( a + y ) f d2 ] , (B2)where f d1 is the force on dust-1 ( f d1 = f rad , d1 ), f d2 is theforce on dust-2 ( f d2 = f rad , d2 ), a = K d1 ρ t+∆td1 ,b = K d1 ρ t+∆tg ,c = K d2 ρ t+∆tg ,d = K d2 ρ t+∆td2 ,x = − (cid:104) ( a + b + c + d ) + (cid:112) ( a + b + c + d ) − ad + ac + bd ) (cid:105) , and y = ( ad + ac + bd ) x . As in section 2.4.1, in order to determine the relative velocitybetween gas and dust, we use equation (B2) which is the
MNRAS000
MNRAS000 , 1–10 (2015) ffect of radiation pressure ρ g K d = 0 exactcode ρ d v g v d P g K d = 10 Figure A1.
Gas mass density (top), dust mass density (second from the top), gas velocity (middle), dust velocity (second from thebottom), and pressure (bottom) profiles. In the left and right panels, we show the results for K d = 0 and K d = 10 , respectively. Thered dots indicate the numerical results. The black solid lines at the left panels represent the exact solutions for K d = 0. The black solidlines at the right panels represent the exact solutions for K d = ∞ . For both simulations, we show the results with the adiabatic index γ = 1.67 at t = 0.2.MNRAS , 1–10 (2015) S. Ishiki, T. Okamoto & A. K. Inoue exact solution of the following equations: ρ d2 ddt v d2 = f rad , d2 + ρ d2 a gra + K d2 ( v g − v d2 ) ,ρ g ddt v g = f rad , g − ∂∂x P g + ρ g a gra + K d1 ( v d1 − v g ) + K d2 ( v d2 − v g ) ,ρ d1 ddt v d1 = f rad , d1 + ρ d1 a gra + K d1 ( v g − v d1 ) . (B3)In order to solve momentum equations, we therefore firstsolve the momentum advection (B1), and then we solve theexact solution of the equation (B3) by equation (B2).In the case for | x | ∆ t (cid:28) | y | ∆ t (cid:28)
1, we use Taylarexpansion, e x ∆ t ≈ x ∆ t or e y ∆ t ≈ y ∆ t , and prevent thenumerical error in calculating (e x ∆ t − /x from becomingtoo large. APPENDIX C: THE TERMINAL VELOCITYAPPROXIMATION
We here show that the terminal velocity approximation maygive an unphysical result when the simulation time step ∆ t is shorter than the drag stopping time t d . In order to derivedust velocity and gas velocity, we have used the equation (2).On the other hand, Akimkin et al. (2017) used the termi-nal velocity approximation. When we employ the terminalvelocity approximation, the equation (2) is transforms intothe following form: (cid:20) p t+∆td p t+∆tg (cid:21) = (cid:20) ρ t+∆td ρ t+∆tg (cid:21) p ∗ d + p ∗ g ρ t+∆td + ρ t+∆t g + (cid:20) ρ t+∆td ρ t+∆tg (cid:21) a gra ∆ t + (cid:20) ( ρ t+∆td ∆ t + ρ t+∆tg t d ) f d ( ρ t+∆tg ∆ t + ρ t+∆td t d ) f g (cid:21) ρ t+∆tg + ρ t+∆td + (cid:20) ρ t+∆td f g ρ t+∆tg f d (cid:21) (∆ t − t d ) ρ t+∆tg + ρ t+∆td . (C1)The advantage of the equation (2) is that it is accurate evenfor ∆ t < t d . In contrast, the equation (C1) becomes inaccu-rate for ∆ t (cid:28) t d , since the relation of ∆ t and t d should be∆ t (cid:29) t d in order the terminal velocity approximation to bevalid. For example, the direction of f d on gas and that of f g on dust in the equation (C1) becomes opposit for ∆ t < t d .We perform simulations by using equation (C1) in steadof equation (2) and compare the simulation results. Simula-tion results do not largely change for Cloud 2, 3, and 4. Thenumerical simulation of Cloud 1, however, is crashed, sincethe timestep becomes ∆ t (cid:28) t d at some steps.The relation between ∆ t and t d becomes ∆ t < t d when the drag stopping time t d is larger than the chemicaltimestep ∆ t chem or the timestep ∆ t CFL defined by Clourant-Friedrichs-Lewy condition. The chemical timestep is definedin equation (7) in paper I.In Fig C1, we present the condition for t CFL ( ≡ α ∆ x/v ) > t d in the case of the intergalactic medium (IGM),the H ii region, and the H i region, where α is constant (weassume α = 0 . x is the mesh size, and v is velocity. Thedetails of numerical setup for the IGM, the H ii region, andthe H i region are listed in Tab C1. The green, red, and bluehatched regions represent the condition for t CFL > t d forthe IGM, the H ii region, and the H i region, respectively. If Table C1.
Numerical setup for the IGM, the H ii region, andthe H i region. The number densities of hydrogen and ionizedhydrogen are represented by n H and n H ii . The temperature of gasis represented by T g . The radius of a grain is represented by a dust .The grain electric potential of dust grains is represented by V d .The relative velocity between a dust grain and gas is representedby ∆ v . n H n H ii T g a dust V d ∆ v (cm − ) (cm − ) (K) ( µ m) (V) (km s − )IGM 10 − − ii region 10 10 10 i region 10 -2 -1 v [km s − ] -2 -1 ∆ x [ p c ] IGMHII regionHI region
Figure C1.
The green, red, and blue hatched regions representthe condition of t CFL > t d for the IGM, the H ii region, and theH i region, respectively. the relation between t CFL and t d becomes t CFL (cid:28) t d , thesimulation may become unstable. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000