Molecular clumps photoevaporation in ionized regions
Davide Decataldo Andrea Ferrara Andrea Pallottini Simona Gallerani Livia Vallini
astro-ph.GA

 
 
 
firstpage–lastpage 2017
 
The [Y/Mg] clock works for evolved solar metallicity stars
Based on spectroscopic observations made with two telescopes: the Nordic Optical Telescope operated by NOTSA at the Observatorio del Roque de los Muchachos (La Palma, Spain) of the Instituto de Astrofísica de Canarias and the Keck I Telescope at the W.M. Keck Observatory (Mauna Kea, Hawaii, USA) operated by the California Institute of Technology, the University of California and the National Aeronautics and Space Administration.

    D. Slumstrup 1 F. Grundahl 1 K. Brogaard 1,2 A. O. Thygesen 3 P. E. Nissen 1 J. Jessen-Hansen 1 V. Van Eylen 4 
 M. G. Pedersen 5
    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




    
 
 
We study the photoevaporation of molecular clumps exposed to a UV radiation field including hydrogen-ionizing photons (hν > 13.6 eV) produced by massive stars or quasars. We follow the propagation and collision of shock waves inside clumps and take into account self-shielding effects, determining the evolution of clump size and density with time. The structure of the ionization-photodissociation region (iPDR) is obtained for different initial clump masses (M=0.01 - 10^4) and impinging fluxes (G_0=10^2 - 10^5 in units of the Habing flux). The cases of molecular clumps engulfed in the HII region of an OB star and clumps carried within quasar outflows are treated separately. We find that the clump undergoes in both cases an initial shock-contraction phase and a following expansion phase, which lets the radiation penetrate in until the clump is completely evaporated. Typical evaporation time-scales are ≃ 0.01 Myr in the stellar case and 0.1 Myr in the quasar case, where the clump mass is 0.1  and 10^3 respectively. We find that clump lifetimes in quasar outflows are compatible with their observed extension, suggesting that photoevaporation is the main mechanism regulating the size of molecular outflows.



ISM: clouds, evolution, photodissociation region - quasars: general
 






§ INTRODUCTION



The diffuse interstellar medium (ISM) is characterized by a turbulent multi-phase structure, showing a broad range of densities, temperatures and chemical compositions. In some regions, gravitational forces and pressure compress the gas to sufficiently high densities, so that the formation of molecules such as H_2 and CO is allowed.     

CO maps have revealed that Giant Molecular Clouds (GMCs) contain a very rich internal structure featuring filaments and clumps <cit.>. The typical sizes of the detected clumps range from 1 to 10 pc. Temperature and density of the gas can be estimated by combining line intensities with radiative transfer calculations. Such studies yield kinetic temperatures in the range T=15-200 K, associated with H_2 densities of n=10^3-4 cm^-3 <cit.>. A correlation between clump temperature and Hα flux suggests that denser clumps are warmer because of a larger UV radiation intensity, likely provided by external sources or internal star-formation activity.

Dense molecular clumps have also been detected within the Photo Dissociation Regions (PDRs) of OB stars, through observations in the infrared and millimiter bands <cit.>. Detections of fine-structure lines of [CI] and [CII], high-J CO rotational lines, and J=3-2 lines of HCN and HCO, show that PDRs are made of a low-density, more diffuse component (n≃ 10^2-4 cm^-3), and high-density structures (n≃ 10^6-7 cm^-3), such as in M17SW <cit.>, and in the Orion bar <cit.>. These clumps must have sizes as small as one tenth or a hundredth of pc, often showing elongated shapes. The presence of such clumps affects significantly the emission spectrum of stellar PDRs.

According to recent observations <cit.>, molecular clumps are also detected in outflowing gas around quasars. The radiation pressure drives a powerful wind (v∼ 0.1c with c speed of light) which collides with the ISM, so that a shock propagates forward into the ISM and a reverse shock propagates back into the wind <cit.>. The outflowing gas is heated by the shock to very high temperatures (T∼ 10^7 K), so that it is expected to be completely ionized. Nevertheless, detections of the CO, OH and H_2O lines <cit.> show that the outflow is in molecular form up to a radius of 1-10 kpc. To reach such distances, the molecular gas has to be structured in clumps, able to provide sufficient self-shielding against the strong quasar radiation field. 

The structure of a molecular clump is significantly determined by the presence of an ionizing/photo-dissociating radiation field, since incident photons with different wavelengths alter the chemical composition of the gas and its physical properties. Far ultraviolet (FUV) radiation (6 eV< h ν< 13.6 eV) is responsible for the dissociation of molecules, determining the formation of a PDR <cit.> at the surface of the clump itself. Furthermore, radiation above the Lyman limit (hν > 13.6  eV) ionizes neutral atoms, and it is completely absorbed within a shallow layer. 

The goal of this paper is to understand the evolution of radius and the density profile of molecular clumps exposed to a UV radiation field including hydrogen-ionizing photons produced by massive stars or quasars. The key point is that an ionized shell and an atomic shell form at the edge of the clump. The dynamics of this layered structure is determined by the fact that each layer is at a different temperature and pressure. For a clump with initial density n_0≃10^5 cm^-3, typical temperatures deep into the clump are T≃ 10-100 K, while an atomic (ionized) region can be heated up to around T≃ 10^3 K (T≃ 10^4 K). We denote this type of ionization/photodissociation regions as iPDR.

In particular, we apply our model to two scenarios. 

	
  * Stellar case: a molecular clump is in pressure equilibrium within a GMC in the proximity of an OB star, and it is suddenly engulfed by the expanding HII region.
	
  * Quasar case: a clump forms as a result of thermal instabilities within the outflow, finds itself embedded in the ionized wind and exposed to the quasar radiation.  

Previous works in the literature concentrated mostly on the effects of non-ionizing photons on photoevaporation of clouds <cit.>. The evolution in their case is simplified by the fact that the clump is exposed only to radiation below the Lyman limit. As a result, the clump develops a single shell structure.

The paper is organized as follows. In Sec. <ref> we describe the model adopted for the structure of molecular clumps, together with the physics of the dynamical and thermal processes involved. In Sec. <ref> we present the results for the evolution of the clump radius, and compute the evaporation time of a clump for different parameters of the system. The results are presented separately for clumps located near stellar or quasar sources. In Sec. <ref> we summarize our results. 




§ MODEL







 §.§ Gaseous environment





The gaseous environment where clumps are located plays a crucial role in determining properties such as temperature, density and confining pressure. We now describe the interclump medium (ICM) which surrounds clumps in the stellar and quasar case. The ICM properties are summarized in Tab. <ref>.

We assume that clumps in a stellar surrounding are in pressure equilibrium with atomic gas, whose temperature depends on the distance and luminosity from the stellar source. We take n_at≃ 10^3 as a typical density for this surrounding gas, and we compute the corresponding temperature according to the FUV flux, using the same method used to compute the temperature of the atomic phase in the clump (see Sec. <ref>). Such temperature ranges between 10^2 K and 10^3 K. The clump is exposed to the radiation of the massive star when it is engulfed in the growing HII region, for which the density is taken to be n_≃ 10.

In the case of quasars, molecular clumps likely form from thermal instabilities within the outflow. Clumps detach from the hot phase at the discontinuity between the fast wind and the ISM <cit.>, starting from the distance at which the outflow has become energy-driven. This critical radius has been estimated by <cit.>:

    R_c = 520  σ_200 M_8^1/2 v_0.1 pc  ,

where σ_200 is the velocity dispersion in the host galaxy in units of 200 km s^-1, M_8 is the mass of the SMBH in units of  10^8 and v_0.1 is the wind velocity in units of 0.1 times the speed of light. In the outflow, the gas is heated to T_≃ 2.2× 10^7 K and has a typical density n_≃ 60 cm^-3 <cit.>. The outflow fragments <cit.> because of thermal instabilties, so that one component cools to a low temperature. The existence of an equilibrium between a 10^4 K and a 10^7 K phase has been studied by <cit.>, while <cit.> show that an atomic clump requires a very short time to cool and turn to molecular form. Molecule formation can occur in the overdensities generated via thermal instabilities, since radiation is efficiently self-shielded and the gas deep into the clump is allowed to cool[This conclusion has been re-examined by <cit.>, who pointed out that molecule formation is problematic due to the efficient dust destruction by the outflow shock.]. When a clump starts to cool, we assume that it maintains pressure balance with the ICM until its temperature is T 10^4 K. Below such temperature the cooling time-scale is very short, and the evaporation process detailed in the next Sections happens before the clump can readjust to the external pressure. The final temperature of the molecular gas is about 100 K, in agreement with detections with CO and water vapour line emission <cit.>.




 §.§ Radiation field
 

Radiation affects the structure of a clump according to the shape of the emitted spectrum. In particular we are interested in ionizing (energy hν≥ 13.6 eV) and FUV photons (6 eV < h ν< 13.6 eV), whose flux G_0 is measured in units of the Habing flux[The Habing flux (1.6 × 10^-3  erg  s^-1 cm^-2) is the average interstellar radiation field of our Galaxy in the range [6 eV, 13.6 eV] <cit.>.].

For stellar sources, we use black body spectra with different effective temperatures T_eff. In terms of solar luminosity (), OB stars have typical luminosities ranging between 10^3 and 10^5. Then, the effective temperature is given by

    T_eff = ( L4π R_⋆^2 σ_sb)^1/4  ,

where σ_sb is the Stefan-Boltzmann constant, R_⋆ is the star radius and L is the bolometric luminosity. We compute R_⋆ through the mass-luminosity and radius-luminosity relations by <cit.>, which for an OB star give

    R_⋆ = 1.33  ( L1.02  )^0.142  ,

where  is the solar radius. Integrating the black body spectrum in the FUV band, typical values of the FUV flux are G_0=10^2 - 10^4 for gas at 0.3 pc from sources with luminosities in the range L=10^3 - 10^5. In the same way, we integrate the spectrum for hν≥ 13.6 eV to obtain the ionizing flux.
	
In the case of quasars, the fundamental difference is the wide extension of the spectrum to the X-rays, so that ionizing radiation is much more intense in this case. An analytical expression for the ionizing flux can be found with the same approach as in <cit.>, thus obtaining that the specific (ionizing) luminosity for ν > ν_L≃ 3.3× 10^15 Hz is:

    L_ν = 6.2× 10^-17(νν_L)^α-2 (Lerg  s^-1)	  erg  s^-1  Hz^-1  ,

where α= 0.5 for a radio-quiet quasar <cit.>. We assume eq. <ref> to be valid for energies below the cut-off value E_c=300 keV <cit.>. Furthermore, we can easily infer a relation between the bolometric and FUV luminosity, setting ν = ν_L in eq. <ref>. The spectrum is almost flat in the FUV band, with typical values G_0≃ 10^3 - 5 at 1 kpc, for L=10^45-47  erg  s^-1.

We investigate the evolution of clumps irradiated by stars or quasars, with L as the only free parameter determining the flux in the bands we are interested in.




 §.§ Clump structure
 




	




In this paper a clump is modeled as a Bonnor-Ebert (BE) sphere <cit.>, which is isotropically affected by an external impinging radiation field. Given the clump mass, the clump temperature and the confining pressure, the BE sphere model allows to compute the radial density profile inside the clump and its radius. In Fig. <ref> we show the density profile for clumps of different mass, with the same temperature T=10 K and a confining pressure P=10^-11 erg cm^-3 (corresponding to a confining medium with T_ = 100 K and n_ = 10^3). The clumps have different radii and same outer density, set by the pressure equilibrium between the clump and the surrounding gas. The density increases towards the centre, with a steeper profile for larger values of the mass. A clump undergoes a collapse if its mass is larger than the BE mass:

    M_be≃ 1.18 c_s^4√(P_0 G^3)  ,

where P_0 is the confining pressure, c_s is the isothermal sound speed and G is the gravitational constant. We are considering only the thermal contribution to pressure, not accounting for turbulent and magnetic pressure. We underline that, apart from the use of a BE density profile, gravity is not included in the hydrodynamical equations for the clump evolution presented in Sec. <ref>.

Given the clump density profile, we assume that the impinging radiation induces a shell-like structure, before any dynamical response of the gas to the photo-heating occurs (sudden heating approximation, see Fig. <ref>). The FUV radiation is responsible for the formation of an atomic layer (HI shell). The more energetic part of the spectrum partially or totally ionizes an outer shell (HII shell), depending on the intensity of the source. This sets up the initial condition for the subsequent hydrodynamical evolution of the clump.

The sudden heating approximation means that radiation dissociates and ionizes particles and heats the gas to its final temperature instantaneously, while the clump shape is unaltered. This situation is often referred to as a R-type ionization front <cit.>. To justify this assumption, we compare the sound-crossing time-scale t_cross with the ionization timescale t_i and the heating time-scale t_h. The Strömgren theory adapted for plane geometry (which can be assumed when the radius of the clump is much smaller than the distance from the source) allows to compute the HII shell thickness, 

    δ_ (t)=δ_ (1-e^-nα_bt)  ,

where n is the gas number density, and α_b is the case B recombination coefficient  <cit.>. Thus, we have t_i = 1/nα_b. On the other hand, the sound-crossing time-scale is t_cross = r_c/c_s where r_c is the clump radius and c_s∼√(k_bT/m_p) is the sound speed. Plugging in typical values, it is easily seen that the condition t_i≪ t_cross is always satisfied for physically reasonable values of r_c (0.01- 1 pc) and c_s (0.1 1 km s^-1 in the cold phase). Regarding the gas heating time-scale (t_h), a simple estimate gives us:

    t_h = k_B T_f Γ(T_f)  ,

where k_B is the Boltzmann constant, T_f is the final gas temperature and Γ(T_f) is the heating rate (in erg/s), mainly due to photoionization. To give some examples, Γ / n varies from ∼ 10^-23 erg cm^3 s^-1 at a distance of 1 pc from a star with L=10^3, to ∼ 10^-16 erg cm^3 s^-1 at 0.5 kpc from a quasar with L=10^47 erg s^-1 <cit.>. The result is that for any value of r_c and c_s of interest, the condition t_h≪ t_cross holds.
	
In what follows, we discuss how we compute the thickness and the temperature of each shell in the clump. The thickness of the HI shell (δ_) is defined as the depth at which hydrogen is found in molecular form. <cit.> find how such depth (expressed as hydrogen column density) scales with the gas density and the FUV flux:

    N_H∝ n^-4/3G_0^4/3  ,

and that the thickness for n=2.3 × 10^5 cm^-3 and G_0=10^5 is δ_ =1.5× 10^16 cm. Then from eq. <ref> we find for the shell thickness

    δ_  = 0.034  (n10^5  cm^-3)^-7/3 ( G_010^5)^4/3  pc  .

<cit.> outline that the PDR temperature is rather constant before it drops to the low values of the molecular core. They plot the temperature for different values of the density and the FUV flux, and we use a fit of their model to estimate the temperature of the HI shell. 

The outer shell presents a partial or total ionization, depending on its density and the intensity of the impinging radiation field. We compute the equilibrium temperature as a function of depth into the shell by balancing photoionization heating and  recombination cooling, line cooling and bremsstrahlung. The presence of a radiation field alters the heating and cooling rate: 1) the ionized fraction of each species is modified and thus the cooling rate by line emission is changed accordingly; 2) injection of photoionized electrons in the gas provides an extra heating term. Approximate heating and cooling functions, assuming collisional equilibrium but non-zero radiation field, are provided by <cit.>. We assume a fiducial value for metallicity, i.e. the mass fraction of elements heavier than helium, of Z=0.02 <cit.>, noticing that metals are important for the energetics of the gas, but their contribution to its dynamics (determined by gas pressure) is negligible. Moreover, we account for Compton heating, which we expect to be important for hard radiation fields:

    H_C = σ_T Fm_e c^2 (⟨ hν⟩ - 4k_bT)  ,

where σ_T is the Thomson cross section, m_e is the electron mass, F is the total flux and ⟨ hν⟩ is the average photon energy beyond the Lyman limit.

Once the temperature profile is computed, we obtain the ionization profile by balancing photoionization, collisional ionization and recombination:

    γ(T) n_e n_p + n_h∫_ν_L^∞F_νhνe^-τ_νa_ν(T) dν = α_b(T) n_e n_p   ,

where n_e, n_p and n_h are the electron, proton and neutral hydrogen density, respectively; F_ν is the specific flux from the source, τ_ν is the optical depth, a_ν and γ are the photoionization cross section and the collision ionization coefficient, respectively <cit.>. In Fig. <ref> we plot the ionization fraction x=n_e/n=n_p/n, as a function of the depth into the clump, for different source luminosities. The ionization profile varies smoothly throughout the HII shell, between the edge of the clump and the PDR, both in the stellar and quasar case. Nevertheless, the region where x is varying is of the order of 10^-5 pc, which is negligible with respect to typical clump radii (0.01 - 1 pc). Then we adopt a reference value for the HII shell thickness δ_, computed with the approximation of a sharp boundary between ionized and phase, in the same way as done for the Strömgren radius for a stellar HII region. Further assuming that the clump radius is much smaller than the distance from the source and that it is illuminated isotropically, the HII shell depth is

    δ_ = 1x^2_maxn^2 α_b(T)∫_ν_L^∞F_νh ν dν  ,

where x_max is the maximum ionization fraction, at the edge of the clump. The thickness δ_ is shown with dashed lines in Fig. <ref>.

We restrict our analysis to clumps where the shell thickness is much smaller than the molecular core radius. This allows us to determine the densities of the HI and the HII shells by using the outer density of the BE sphere. Furthermore, it simplifies our calculations, because we can consider separately the evolution of the shells and the core, since the dynamical time-scale of the former  is much shorter than the core one. In Fig. <ref> we plot the ratio of the HI shell thickness (left panels) and the HII shell thickness (right panels) to the total clump radius r_c as a function of the clump mass, for stars and quasars with different luminosities. The atomic shell is always thicker than the ionized shell, showing a large self-shielding effect of the ionized gas. The typical masses of clumps are different for stars and quasars, and we consider only masses smaller than the BE mass for collapse. The distances of clumps from the source are fixed in the two scenarios, and are reasonable for molecular gas engulfed by an expanding stellar HII regions (the distance scale is given by the Strömgren radius) and clumps forming in quasar outflows (critical radius given in eq. <ref>).

Clumps presenting a ratio δ_ /r_c=1 are completely dissociated on a time-scale t_i, and the analysis of this paper restricts to clumps where δ_≪ r_c. From the plot, we see that such condition is usually verified and breaks only in the stellar case for small clumps (M <
 0.01) and very intense sources (L≃ 10^5 ). 





 §.§ Shock dynamics inside the clump
 






Having set the initial conditions on a clump, i.e. a core-double shell structure, now we can study its dynamical evolution for t>0. The different layers in the clumps have different pressures, so that a shock or rarefaction waves originate, enforcing a continuous value of pressure and velocity across the contact discontinuity between two layers. 

The cooling time-scale of a gas at temperature T is

    t_cool = k_B TΛ(T)  ,

with Λ being the cooling function given by <cit.> for molecular gas  and by <cit.> for PDRs. For the range of temperatures and n, r_c and c_s values of interest here, t_cool≪ t_cross. Thus the fluid motion and the propagation of any disturbance in the gas (as shock and rarefaction waves) can be safely considered as isothermal processes. 

A qualitative diagram of the possible outcomes at an arbitrary discontinuity is shown in Fig. <ref>. In the situation considered in the upper inset, a gas has a pressure P_r to the right of an interface, and a pressure P_l to the left (with P_l>P_r). The velocity to the right is v_r=0, while we consider different values v_l,1, v_l,2 and v_l,3 for the velocity to the left. The solid lines connect to the initial state all the possible final states of the gas, when it is crossed by a rarefaction wave (RW) or a shock wave (SW). For example, the points on the blue line represent the possible final states of the gas to the right when it is crossed by a rarefaction wave.  The solution of the discontinuity problem is obtained when the lines departing from the two initial states of the gas to the left and to the right intersect, since the final values of P and v must be the same across the discontinuity. This also determines which kind of wave is required, i.e. a SW or a RW. 

The solution of the problem for given values of the initial pressure, density and velocity across the discontinuity is obtained numerically, imposing the final pressure and velocity to be continuous. The post-shock values are obtained solving the isothermal Rankine-Hugoniot conditions <cit.>


    ρ_0v_0   =ρ_1v_1
    ρ_0v_0^2+P_0   =ρ_1v_1^2+P_1
    
    T_0    = T_1


where the subscript 0 is used for pre-shock values and the subscript 1 for post-shock values, with v velocity in the shock front frame. Rearranging the equations <ref>, it is possible to write the following relations


    ρ_1   =ρ_0ℳ^2 
    
    P_1   =P_0ℳ^2 
    
    v_1   =v_0/ℳ^2


with ℳ = v_0/c_s being the shock Mach number. 

On the other hand, rarefaction waves are not discontinuities and values of flow variables across such waves are obtained following <cit.> and adapting the calculations to the isothermal case. Consider a wave originating at x=0 and propagating toward x>0, such that the final velocity after the wave has completely passed is v_f=-U. The profile between the wave head, moving at the initial sound speed c_s,0 in the gas, and the wave tail, moving at speed v_tail=c_s,0 - (γ+1)U/2, is 



    v(x)		   =-(c_s,0-x/t )
    ρ(x)	   =ρ_0exp(x/c_s,0t-1 )
    
    P(x)		   =P_0exp(x/c_s,0t-1 )


where ρ_0 and P_0 are the values of density and pressure before the rarefaction has passed.

Since shock waves are discontinuities, an interaction between two shocks can be treated as an arbitrary discontinuity between post-shock values of flow variables. To simplify our analysis, we also consider interactions involving rarefactions as discontinuities, by accounting only for the post-rarefaction values of flow variables. In App. <ref> we compare this approach with a numerical solution of the fluid dynamics equations, showing that the two results differ negligibly.

To compute the shock speed inside the clump, we have to account for the spherical geometry and for the density gradient given by the BE profile. Following <cit.>, the flow equations can be written as



    ∂_t ρ + ∂_r(ρ v) + ρ v ∂_r A(r)A(r)   = 0 
    ∂_t v + v ∂_r v + 1ρ∂_r P - 1ρ_0(r)∂_r P_0(r)    =0


where r is the radial coordinate, A(r)=4π r^2 in the spherical case, ρ_0(r) and P_0(r) are the initial density and pressure profiles for a BE sphere. Eq. <ref> and eq. <ref> can be combined to give the equivalent equation valid along the curves dr/dt = r+c_s in the (r,t) plane (called the C_+ characteristics):

    dP + ρ c_s d v + ρ c_s^2 vv+cA'(r)A(r) - ρ c_sv + c_s1ρ_0(r) P'_0(r)=0  ,

where the prime denotes the derivative with respect to r. 
According to <cit.> the shock trajectory in the (r,t) plane is approximately a C_+ characteristic, so that eq. <ref> can be applied along the shock. Then we can write eq. <ref> as a function of the Mach number ℳ, substituting the post-shock values P, ρ, v from eq. <ref>:

    dℳdr = - 12ℳ^2ℳ^2 - 1A'(r)/A(r) + 12ℳ^3ℳ+1P_0'(r)P_0(r)  ,

which is a differential equation for ℳ as a function of r. 

In Fig. <ref>, the solid lines show the the numerical solution of eq. <ref> for a molecular clump with mass M=0.1 and radius r_c=0.025  pc, for different values of the initial Mach numbers ℳ_0 of the shock at the edge of the clump, assuming an isothermal shock. For comparison, the dashed line is the classical analytical solution obtained by <cit.> for a homogeneous density distribution, and in the limit of a strong adiabatic shock:

    ℳ(r) = ℳ_0 (r_0r)^n(γ)  ,

where r_0 is the radius of the bubble, and n(γ) is an exponent depending on the adiabatic coefficient γ (e.g. n(5/3)≃ 0.543 for monoatomic gas and n(7/5)≃ 0.394 for diatomic gas). As opposed to Guderley solution, the isothermal shock speeds up considerably only at a smaller radius. After the shock wave has reached the centre, a reflected shock will travel outwards. The velocity as a function of radius has the same profile of the focusing shock.









§ RESULTS





In our model, a clump exposed to UV radiation develops an ionized PDR (iPDR) at its surface. First, we inspect the qualitative behaviour of the structure simply using the arbitrary discontinuity criterion. Fig. <ref> shows a diagram of shock and rarefaction waves propagating inside the clump because of the pressure difference between adjacent layers.

The HII shell pressure (P_) is higher both than the pressure of the HI shell (P_) and the pressure of the confining ICM (P_), hence two rarefaction waves cross the HII shell, originating from its edges. Since δ_≪ r_c, the evolution of the HII shell has a much shorter time-scale than the clump evaporation time. The rarefaction waves which propagate into it interact and reflect at its edges, determining a complex density profile. Nevertheless, the global effect is that the HII shell expands decreasing its density, eventually becoming completely transparent to the ionizing radiation (i.e. the mean free path of photons is much larger than the shell thickness) . 

On the other hand, P_ is lower than P_, but higher than P_h_2. Thus, a shock is driven from the HII shell into the HI shell, and a rarefaction wave propagates from the discontinuity with the core (see Fig. <ref>). Once the shock has crossed the HI shell, it reaches the core surface and speeds up its contraction. As a result, the inner boundary of the HI shell moves faster than the outer boundary, and the HI shell is also expanding and becoming transparent on a time-scale shorter than the core evolution time-scale.
 
The cold (T_2 ≃ 10- 100 K) molecular core is compressed because of the shock wave originating at the discontinuity with the atomic shell and propagating towards the centre. In addition, the shock wave originating at the HII/HI boundary reaches the core surface and catches up with the shock already propagating in the core, resulting in a single stronger converging shock wave. 

The shock wave is reflected at the centre of the clump, and eventually gets back to the core edge. The contraction is almost halted, and since the core has a much higher density than the surrounding medium, it starts to expand. The expansion velocity v_exp is computed considering the discontinuity between the core compressed by the reflected shock wave and the ICM at rest, using the arbitrary discontinuity algorithm.

We have explicitly verified that the core is so dense (n≃ 10^5-6 cm^-3) during the contraction phase that the FUV radiation penetrates to a negligible depth with respect to its radius. Thus we can ignore photoevaporation during the contraction phase. When the clump starts expanding, we have computed for each time t the thickness δ_ (t) of an HI shell (see eq. <ref>) for the corresponding core gas density. We get the core radius at t by subtracting δ_ (t) to the radius R(t) = R_0 + v_expt (R_0 is the core radius at the end of the shock-contraction phase). 



 §.§ Stellar case
 



As a first application of our analysis, we consider molecular clumps photoevaporating because of stellar radiation. We consider a cold clump (T_2 = 10 K) located in the surrounding of a star, embedded in an atomic region with density n_at = 10^3 cm^-3. Then we assume that the expanding HII region of the massive star engulfs the clump (the density of the HII region gas is n_ = 10, see Tab. <ref>), and we apply the machinery we developed in Sec. <ref>. We consider stars with bolometric luminosities L=1× 10^3, L=1× 10^4, L=5× 10^4 and L=1× 10^5. We assume the clump is located 0.3 pc from the source, since this distance is smaller than the Strömgren radius for every star in our set (for the fainter star the Strömgren radius is R_Str≃ 0.75 pc for a gas density n=10). The BE masses for the collapse of molecular clumps at 10 K are around few tenths of solar masses, and for each luminosity we consider only clumps with mass below that limit.

The time evolution of the molecular core radius is shown in Fig. <ref> (left panel), where a clump of initial mass M=0.1 is exposed to the stellar radiation field for the different luminosities considered. The radius has a similar evolution for the different luminosities, with a shorter time-scales for larger luminosities. Consider for example the (brown) curve for L=10^4. A clump with mass 0.1 at the distance of 0.3 pc has an initial radius of 0.023 pc when it is in pressure equilibrium with the ICM. In the shock contraction phase, the radius reduces to 6× 10^-4 pc in about 6000 yr because of the shock waves driven by the heated HI and HII shells. Then the expansion phase follows, and the core expands allowing the impinging radiation to penetrate and dissociate the molecules. This occurs significantly after the radius reaches its maximum value r≃ 0.025 pc. 

While the contraction phase has almost the same duration for the tracks of the three more intense sources, we see that it takes more time for the 10^3 star. In fact for this source the ionization fraction is low (see Fig. <ref>), since the temperature of the HII shell is only about 900 K. The shock driven by the HII shell is weak for this star, and needs more time to reach the centre of the core.

The lifetime of a clump (t_c) is defined as the time when the core radius goes to zero. In Fig. <ref> (right panel) we show t_c as a function of the clump mass and the source luminosity, at the same distance to the source (0.3 pc). 

<cit.> compute the lifetime of clumps located in a stellar PDR, in the absence of ionizing radiation. They account for photoevaporation by assuming that the clump continuously loses mass at a rate

    d M d t = - 4πρ_c r^2_c(t)c_  ,

where ρ_c is the mean mass density of the clump, and c_ the sound speed in the PDR of the clump. This implies that the clump loses mass also in the shock-compression phase. 02cit do not account for the shock reflection at the centre of the clump, so that the core does not expand after the compression phase. On the other hand, in our treatment, photoevaporation is negligible while the clump is being compressed, and the shock reflection allows for the following expansion of the core. As a result, radiation is allowed to penetrate and dissociate the molecules only when the gas is sufficiently expanded and diluted. Furthermore, 02cit find 
that under certain initial conditions[According to <cit.>, a clump undergoes a shock-compression only if its mean column density is 
    nr_c < 2.7×  10^21  cm^-2 (c_/c_2)^3  .
] there is no shock-compression, since the shock suddenly stalls just after its formation, and the clump directly expands and photoevaporates. In our treatment, we do not recover this scenario, since we always allow the shock to reach the centre of the core. Magnetic and turbulent contribution to pressure are included by 02cit, but not in this work. 

In Fig. <ref> we compare our predictions for the photoevaporation time (without ionizing radiation) with those from the 02cit model. A range of clump masses between 0.01 and the BE mass is considered, at a distance of 0.3 pc, for three different source luminosities, A modification in our code for iPDR is required, since the HII shell is not present when ionizing radiation is absent. In the 02cit model, the 10^4 and the 10^5 sources induce a shock-compression in the clump, while instead they predict an initial expansion for the 10^3 case. The evaporation time-scale for the low-L case differs by an order of magnitude with respect to the other two. In our model we do not find such dichotomy, and the evaporation time smoothly increases with L. However, the lifetimes are in agreement within a factor of 2 with those found in 02cit.

Finally, we make a comparison between the evaporation times obtained with our full iPDR model (see Fig. <ref>) and our model without ionizing radiation, i.e. with no HII shell (see Fig. <ref>). Clump lifetimes are always shorter when we consider the ionizing part of the spectrum, generally by a factor between 2 and 4 depending on clump mass and luminosity. This behaviour is expected, since the outer shell of the clump is heated to an high temperature and a stronger and faster shock propagates into the clump, decreasing its evolution time-scale. 





 §.§ Quasar case
 







We now describe the evolution of a clumps forming in the ionized outflows of quasars. We choose 1 kpc as a typical distance of a molecular clump from the source, that is of the order of the critical radius (eq. <ref>). The mass of clumps has been estimated by <cit.> to be around 8600, thus we consider masses up to 10^4. The ICM is the hot ionized medium of the wind, with temperature T_ ¸≃ 2.2× 10^7 K, and density n_≃ 60 cm^-3 at 1 kpc from the source. As explained in Sec. <ref>, we assume that the clump is in pressure equilibrium with the ICM until its temperature reaches 10^4 K (when the clump gas is still in atomic form). Afterwards, the gas turns into molecular form (T_2=100 K) in a very short time-scale, so that the density profile remains unchanged with respect to the 10^4 K BE sphere. We apply our model to quasars with bolometric luminosities L=1× 10^45  erg  s^-1, L=1× 10^46  erg  s^-1, L=5× 10^46  erg  s^-1 and L=1× 10^47  erg  s^-1, with a spectral energy distribution given by eq. <ref>.

The evolution of the molecular core radius of the clump is shown in the left panel of Fig. <ref>, while the right panel shows the lifetime of a clump at a distance of 1 kpc from a quasar, as a function of clump mass and source luminosity. Similarly to clumps around stars, the core radius presents a contraction phase followed by an expansion phase, where the core is dissociated and then ionized. Notice that the evaporation time is about ten times longer than in the stellar case, although the clump mass considered is about 10^4 times larger. This is because both the FUV and ionizing radiation field are much more intense since the quasar spectrum  extends to very high energies (see Table <ref> for reference values of the fluxes). This implies a higher temperature of the HII and HI shells, and a stronger compression ratio of the shock waves originating at the discontinuities.

Observations have detected molecular gas only up to a maximum distance of few kpc from quasars <cit.>. Photoevaporation has sometimes been invoked as an explanation for such limited extension. We have slightly modified our code to  account for the fact that radiation intensity decreases as the clump moves away from the source, being carried by the outflow.  have shown that molecular clumps forming at the base of the adiabatic outflow are ablated in a short time because of the friction by hot flowing gas. Therefore, we analyse the alternative scenario in which clumps form within the outflow, so that they are at rest with the outflow and they are not subject to a strong acceleration. In this way our model is able to predict the distance travelled by a clump during its lifetime, and we can compare this length with observations of molecular outflow extensions. 

We consider a subset of active galactic nuclei (AGNs) listed in 14cit (Table <ref>). In Fig. <ref> we plot the outflow extension obtained with our photoevaporation model as a function of source luminosity. We study clumps with 90% of their BE mass for collapse (i.e. the most massive clumps that do not collapse), forming at the contact discontinuity (CD) between the quasar wind and the surrounding ISM. According to <cit.> model, the initial position of the CD coincides with the critical radius R_c given in eq. <ref>, and it moves at a speed

    v_cd = 875  σ_200^2/3  km s^-1

in the energy-driven phase. In Fig. <ref> the dotted lines correspond to the position of the CD at different times, while the shaded regions represent the maximum distance that clumps with a range of masses (0.1-0.9 times their Bonnor-Ebert mass) can travel before photoevaporating. The observed extension of the outflows in the considered quasar sample exceed the maximum distance travelled by clumps before they are photoevaporated, if they form at R_c. This implies that there is no mechanism more efficient than photoevaporation in destroying molecular clumps. On the other hand, the existence of outflows with large extensions (up to 1 kpc) suggests that clumps continue to form within the outflow, when the CD moves outwards from its initial position. This can bee see from the other two cases shown in Fig. <ref> where we consider also 
clumps formed at a later times when the CD has moved to a radius R_c +v_cdΔ t, with Δ t=0.3 ,0.6 Myr. It appears that such delayed formation via thermal instabilities in the outflow can match the observed extensions. 






§ CONCLUSIONS





We have studied the evolution of molecular clumps exposed to radiation having both a far ultraviolet (FUV) and an ionizing component, determining the formation of an ionization/photodissociation region (iPDR) at the surface of clumps. The cases of a clump forming in the surroundings of an OB stars and a clump forming in the fast outflow of a quasar are studied separately. The clump is assumed to be an isothermal Bonnor-Ebert sphere with a mass lower than the critical mass for collapse. We assume a sudden heating scenario, inducing a shell structure in the clump, and then we analyse the evolution of its radius and density profile as a function of time, finally computing the clump lifetime (i.e. the time at which the molecular gas in the clump is completely dissociated). The clump evolution is solely determined by two parameters: its mass, M, and the bolometric luminosity L of the source. 

We show that the pressure difference between adjacent layers causes the propagation of shock and rarefaction waves into the clump. The core shrinks until the shock wave hits the centre and reflects back, while the external layers expand and become eventually transparent to radiation. The dense core is thus surrounded by a diluted medium and it starts an expansion phase. As a result, the core density decreases and the radiation propagates in the interior, progressively evaporating the whole core. In this analysis we have not included gravity effects which could limit the expansion following the shock-contraction phase of clumps. Gravity may also play a role for the clumps that become gravitationally unstable during the contraction phase, possibly triggering star formation <cit.>.

In the stellar case, we find that a higher luminosity speeds up considerably the shock-contraction phase: clumps of 0.1  at 0.3 pc from the source evaporate in 0.01 Myr for the brightest star considered (10^5), while it takes 0.06 Myr in the case of the 10^3 star. Indeed, the radiation from the fainter star is not able to completely ionize the surface layer of the clump, resulting in a lower pressure of the HII shell and a weaker shock-induce contraction phase. 

Our model agrees within a factor of 2 with the <cit.> model, in the case of clumps embedded in the PDR of a massive star and in absence of ionizing radiation. The main difference between the two models is the evaporation channel. In 02cit evaporation is due to a constant mass flow from the clump surface; in our model the clump evaporates as a consequence of the expansion and dilution driven by the reflected shock wave. We also notice that, in the absence of ionizing radiation,  evaporation times are always longer by a factor 2-4 with respect to the full iPDR model including both FUV and ionizing radiation. Therefore, considering ionizing radiation is important, since the evolution history of clumps is significantly modified.

In the context of high-redshift galaxies, this is significant for far infrared (FIR) emission, as [CII]. Indeed, most of the [CII] emission from high-redshift galaxy seems to be due to molecular clumps <cit.> and because of the high radiation field observed in such galaxies <cit.>, photoevaporation can play an important role. While <cit.> analyses the effect of clump photoevaporation using a time evolution based on 02cit, we argue that shorter photoevaporation time-scales obtained with our iPDR model could further affect the detectability of  high-redshift galaxies. However, we underline that other effects are also important: the contrast with cosmic microwave background (CMB) attenuates the observed FIR emission for redshift z≳ 5 <cit.>, which is relevant for low density gas <cit.>, while CO destruction by cosmic rays may enhance [CI] and [CII] emission <cit.>.  

The evolution in the quasar context is characterized by a similar behaviour. The duration of the contraction phase is roughly constant for different L, since all the quasars in the set are able to completely ionize and heat to about 10^4 K the outer shell of the clump. We obtain evaporation times of 0.21 Myr for the 10^47  erg  s^-1 quasar and 0.51 Myr for the 10^45  erg  s^-1 quasar. With comparison to the stellar case, the evaporation times are longer only by a factor ∼ 10, even though the clumps in the quasar case are ≃ 10^4 times more massive. This is consistent with the higher UV fluxes produced by quasars in spite of the larger spatial scales of the problem (see Tab. <ref>). 

Applying our algorithm to clumps embedded in quasar outflows, we have been able to predict the outflow extension. This is set by the maximum distance travelled by clumps before photoevaporating, assuming that they form at the contact discontinuity (CD) between the quasar wind and the ISM. We find that the observed molecular outflow extensions are always larger than the distance travelled by clumps forming at the initial position of CD, but they are compatible with clumps forming at the CD with a time delay Δ t≃ 0 - 0.6 Myr after the outflow has entered the energy-driven phase. Therefore, we argue that:


  * photoevaporation must be a crucial mechanism involved in the evolution of molecular gas structures in quasars, since none of the observed outflows has a smaller extension than what predicted with our photoevaporation model;

  * clumps need to form continuously within outflows, when the CD has moved farther from the quasar, in order to explain the most extended outflows.

A more comprehensive analysis of quasar outflows should consider a distribution of clump masses, the contribution of scattered light in a clumpy medium and the possible occurrence of star formation within the outflow <cit.>.




style/mnras






§ COLLISION OF RAREFACTION WAVES





The pressure and density profile of a fluid crossed by a rarefaction wave is given by eq. <ref> (assuming an isothermal  process). Since the profile is not a discontinuity, as for shock waves, collisions between rarefaction waves cannot be studied analytically. The MacCormack method <cit.> is a discretization algorithm for solving numerically hyperbolic differential. We apply this method to find a solution of fluid dynamics equation, where the initial condition is set by two approaching rarefaction wave, travelling in the opposite direction. In Fig. <ref> we have simulated the collision of the two rarefaction waves generating at the edges of the HII shell of a clump exposed to the UV radiation of a quasar (see Fig. <ref> and description in Sec. <ref>). The two rarefaction waves are set so that they travel 3× 10^-4 pc before colliding, a distance of the order of the HII shell thickness (see Fig. <ref>). The lines in figure represent the pressure and density profile at different intervals of time. The total length considered (3× 10^-2 pc) has been divided in a 10^5-cell grid, while the timesteps are chosen according to the Courant condition <cit.>. The result of the collision is two rarefaction waves travelling in the opposite direction, with an equilibrium pressure in the central region 

    P_mc≃ 1.0× 10^-7  erg  cm^-3  .
 
The same prediction for the outcome can be obtained with the arbitrary discontinuity algorithm discussed in Sec. <ref>. We approximate the rarefaction waves to jumps in the flow variables (P, ρ, v), i.e. we consider a contact discontinuity between the post-rarefaction states of the fluid. We obtain for the equilibrium pressure a value 

    P_ad≃ 9.8× 10^-8  erg  cm^-3  ,

in agreement with the value found previously within a few percent. Therefore, the discontinuity approximation introduces a small error, but it is more convenient from a computational point of view.



