© The Authors 2025

1. Introduction

The chemistry of the molecular clouds plays a vital role in our understanding of the star formation process. Observationally, the dynamics of molecular clouds are constrained through the Doppler shift of molecular line emission. Chemistry also plays an important role in the thermodynamics of molecular clouds through atomic and molecular line cooling (Draine 2011; Tielens 2013). Therefore, astrochemical models of molecular clouds are central in both theoretical and observational studies of star formation and the physics of molecular clouds.

In astrochemical modeling, the effective attenuating column density for the external radiation field is a primary input (Wakelam et al. 2013; Bovino & Grassi 2024). This effective column density, N_eff (the column density of hydrogen nuclei), or the effective extinction, A_V, eff (magnitude), gives the local far-ultraviolet (FUV) flux, χ/χ₀ ∝ e^{−A_V, eff}, where χ₀ is the FUV flux at the cloud boundary. In the case of isotropic radiation, the effective column density at a point, x, is the weighted average (see below) of all directions of the integrated column density from point x to the cloud boundary. It is thus different from the observational column density, which is integrated through the entire cloud along an observed line of sight. Given the importance of this parameter to astrochemical models, estimating it is crucial for accurate chemical models of molecular clouds.

Astrochemical models have to make a choice for N_eff, such as assuming it is equal to the observed column density or leaving it as a free parameter in model fitting. In molecular cloud simulations, N_eff is often calculated through ray-tracing (e.g., Nelson & Langer 1997; Glover & Mac Low 2007; Glover et al. 2010; Van Loo et al. 2013; Safranek-Shrader et al. 2017; Seifried et al. 2017; Hu et al. 2021; Wünsch 2024), where the hydrogen nucleus column density is summed along each ray out of the domain from each cell using, for example, a HEALPIX-based method for the ray directions, or the simpler six-ray approach along each cardinal direction. The effective extinction is the average,

$$\begin{array}{c}\hfill A_{\mathrm{V},\mathrm{eff}}=-\frac{1}{\gamma }log\left(\frac{1}{N_{\mathrm{rays}}}{\displaystyle \sum _i^{N_{\mathrm{rays}}}}{e}^{-\gamma A_{V,i}}\right),\end{array}$$(1)

where the coefficient, γ, is a normalization constant chosen to be the photo-destruction rate coefficient for particular molecules (e.g., Glover et al. 2010; Hu et al. 2021) or is based on dust absorption (e.g., Nelson & Langer 1997). When no ray-tracing is used, N_eff is computed using an assumed physical length scale and the cell’s density (see the discussion in Safranek-Shrader et al. 2017). Recently, Bisbas et al. (2023, hereafter B23) provided a fit to several hydrodynamic simulations for A_V, eff − n_H,

$$\begin{array}{c}\hfill A_{\mathrm{V},\mathrm{eff}}=0.05exp\left[1.6{\left(\frac{n_{\mathrm{H}}}{{\mathrm{cm}}^{-3}}\right)}^{0.12}\right].\end{array}$$(2)

However, this fit is biased by the assumptions and initial conditions in the underlying simulations: the simulations were predominantly for Milky Way-like clouds, with both a resolution limit and a numerical column density floor. Safranek-Shrader et al. (2017) and Hu et al. (2021) also fitted a polynomial from their simulations and for n_H > 10 cm⁻³ found $A_V\propto n_{\mathrm{H}}^{0.40}$.

In this Letter we present a simple analytic model connecting the local gas density to the effective column density. The model is derived such that, with few assumptions about the cloud-scale turbulence and dynamics, an effective column density can be estimated that is in general agreement with the results from three-dimensional simulations. In Sect. 2 we present the underlying model components. In Sect. 3 we show the results of the model in comparison to the previous three-dimensional calculations. Finally, in Sect. 4 we discuss the potential applications for this model and its limitations, and provide some concluding remarks.

2. Methods

We built an analytic model to reproduce a characteristic N_eff(H)−n_H relationship from a simple, physically intuitive construction. We built up the model by considering three different characteristic length scales: (i) turbulent, (ii) gravitational, and (iii) numerical resolution. Figure 1 shows a chart of the scales of interest and the various parameters that enter into the model.

Fig. 1. Chart showing the different components and where parameters enter into the model.

2.1. Turbulence length scale

Within the molecular cloud, any particular parcel of gas is surrounded by an envelope of turbulent gas that attenuates the external FUV radiation. We assumed this column density is, on average, related to an effective length scale, ℓ_turb, outwards to the cloud boundary through the turbulent envelope multiplied by the average number density, n₀, in the turbulent envelope,

$$\begin{array}{c}\hfill N_{\mathrm{turb}}\approx {\ell }_{\mathrm{turb}}{n}_0.\end{array}$$(3)

If the cloud depth is known, for example from three-dimensional cloud maps that use observations from Gaia (Zucker et al. 2021; Rezaei Kh. et al. 2024; Edenhofer et al. 2024), and there is a measured column density distribution function, this can be estimated by dividing the mean of the log-normal component by the cloud depth.

We assumed that at a given number density, n_H, we can determine a characteristic length scale ℓ, given assumptions on the macroscopic cloud properties. We utilized the virial parameter, α_V, to calculate this effective size as a function of density, defined by

$$\begin{array}{c}\hfill {\alpha }_V=\frac{5\sigma {(\ell )}^2}{G{\mu }_H{m}_H{n}_{\mathrm{H}}{\ell }^2},\end{array}$$(4)

where σ(ℓ) is the turbulence velocity dispersion as a function of the length scale, G is the gravitational constant, μ_H is the mean molecular weight per hydrogen nucleus, and m_H is the mass of hydrogen. To get the above expression, we have assumed that the mass and density are connected via M ≈ ρℓ³. Relating a length scale to the density and virial parameter assumes a hierarchical density structure. The linewidth-size relationship is

$$\begin{array}{c}\hfill \sigma (\ell )={\sigma }_0{\left(\frac{\ell }{L_0}\right)}^{\beta },\end{array}$$(5)

where σ₀ and L₀ are the normalization velocity and scale (or driving scale), respectively, and β is related to the physics behind the turbulence cascade (Larson 1981; Heyer & Dame 2015). Recent surveys have shown that β typically ranges between 0.3 and 0.67 (Rice et al. 2016; Kauffmann et al. 2017; Spilker et al. 2022). Combining Eqs. (4) and (5) gives us the length scale,

$$\begin{array}{c}\hfill \ell ={\left[\frac{G{L}_0^{2\beta }{\mu }_H{m}_H{n}_{\mathrm{H}}}{5{\sigma }_0^2}{\alpha }_v\right]}^{1/[2(\beta -1)]}.\end{array}$$(6)

At low densities or high virial parameters, this length scale, ℓ, can exceed the cloud size, L_cl. Therefore, we computed a modified length to limit ℓ to be less than L_cl,

$$\begin{array}{c}\hfill \overline{\ell }=\ell \times S_{\ell }(z_{\ell }),\end{array}$$(7)

where z_ℓ = L_cl/ℓ. The dampening function requires that for ℓ < < L_cl, $\overline{\ell }\to \ell$, and for ℓ ≳ L_cl, $\overline{\ell }\to L_{\mathrm{cl}}$. We adopted $S_{\ell }(z_{\ell })=\mathrm{erf}(\frac{\sqrt{\pi }}{2}{z}_{\ell })$ since it closely approximates a step-function filter while remaining smooth at all scales. This choice primarily impacts the lowest densities. The length scale out of the cloud is then

$$\begin{array}{c}\hfill {\ell }_{\mathrm{turb}}=R_{\mathrm{cl}}-\frac{1}{2}\overline{\ell },\end{array}$$(8)

where $R_{\mathrm{cl}}=\frac{1}{2}{L}_{\mathrm{cl}}$ is the cloud radius. In low-density gas, more likely to be in the outer regions, ℓ_turb → 0, and so N_turb → 0. Conversely, in the dense gas, ℓ_turb → R_cl, and so N_turb → n₀R_cl. Since the effective column density for an isotropic external radiation field is most sensitive to the lowest column density to the cloud boundary, it can be readily seen for any approximately spherically symmetric cloud that the maximum ℓ_turb is R_cl.

2.2. Gravitational length scale

In the high-density limit, the effective column density is dominated by gravitationally bound gas in the immediate surroundings. The length scale of relevance here is the Jeans length, λ_J. We defined this component as

$$\begin{array}{c}\hfill N_J(n_{\mathrm{H}})\approx \frac{1}{2}{\lambda }_J{n}_{\mathrm{H}}=\frac{1}{2}{\left(\frac{\pi }{G{\mu }_H{m}_H}\right)}^{1/2}{c}_s{n}_{\mathrm{H}}^{1/2},\end{array}$$(9)

where $c_s=\sqrt{k_{B}T/{\mu }_p{m}_H}$ is the sound speed, T is the gas temperature, k_B is the Boltzmann constant, and μ_p is the mean molecular mass per free particle. Under isothermal conditions, $N_J\propto n_{\mathrm{H}}^{1/2}$.

We accounted for the fact that gravity is not dominant at all densities by introducing a threshold that determines when gravity becomes dynamically important. We denote the threshold function as S_n(n_H/n_T), where n_T is the transition density. A good parameterization of n_T equates the Jeans length and turbulence sonic scale, and following Burkhart & Mocz (2019) for isothermal turbulence,

$$\begin{array}{c}\hfill \frac{n_T}{n_0}=\frac{{\pi }^2}{15}{\alpha }_V{M}_s^2,\end{array}$$(10)

where M_s and α_V are calculated at L₀ using σ₀, which is similar to the expression in Padoan & Nordlund (2011).

We adopted a softened threshold, S_n(z_n) = z_n/(1 + z_n), where z_n = n_H/n_T. The choice of threshold impacts how rapidly or slowly the effective column density transitions between the turbulence- and gravity-dominated regimes. The total gravitational column density is

$$\begin{array}{c}\hfill N_{J,T}(n_{\mathrm{H}})=S_n(n_{\mathrm{H}}/n_T)\times N_J.\end{array}$$(11)

2.3. Resolution limits

In hydrodynamic simulations, the total attenuating column is constrained by the limited resolution (Glover et al. 2010; Van Loo et al. 2013). The resolution limit imposes a floor for the attenuation column density,

$$\begin{array}{c}\hfill N_{\mathrm{\Delta }x}\approx C_{\mathrm{\Delta }x}{n}_{\mathrm{H}}\mathrm{\Delta }x,\end{array}$$(12)

where the prefactor C_Δx depends on how the self-contribution is accounted for, for example a whole cell size or half a cell size.

2.4. Combined model

The final N_eff(H)−n_H relation without numerical resolution effects is

$$\begin{array}{c}\hfill N_{\mathrm{eff}}(n_{\mathrm{H}})=N_{\mathrm{turb}}(n_{\mathrm{H}})+N_{J,T}(n_{\mathrm{H}}).\end{array}$$(13)

If a restriction based on resolution is needed, then the expression is modified to

$$\begin{array}{c}\hfill N_{\mathrm{eff}}(n_{\mathrm{H}})=max(N_{\mathrm{turb}}(n_{\mathrm{H}})+N_{J,T}(n_{\mathrm{H}}),N_{\mathrm{\Delta }x}(n_{\mathrm{H}})).\end{array}$$(14)

To enable comparisons with B23, we also used a floor of 0.05 mag.

3. Results

For our fiducial model, we aim to model Milky Way-like molecular clouds. We adopted the fiducial parameters in Table 1, mean molecular weights, μ_H = 1.4 and μ_p = 2.33, a conversion factor between the hydrogen nucleus column density and visual extinction of A_V = 6.29 × 10⁻²²N(H) (Röllig et al. 2007), and a normalization scale for the linewidth-size relationship of L₀ = 1 pc. In Appendix A we investigate the impact of assuming a fixed c_s and μ_p and show that assumptions on μ_p have a minimal impact, while a variable c_s can increase N_eff at intermediate and low densities. However, these have less of an impact than varying key physical parameters.

Figure 2 shows the fiducial model in comparison to the B23 fit. At moderate densities, the model matches B23 within a factor of two, although toward low densities, a floor is needed to remain consistent with B23. The slope flattens around n_H ≈ 500 cm⁻³ due to the transition between the turbulent to gravitational components. At high densities, N_eff is dominated by the gravitational component with a slope flatter than 0.5 due to the contribution of the turbulent column, even with an isothermal assumption. The fiducial model slope between 10³ ≤ n_H ≤ 10⁵ is consistent with the $\propto n_{\mathrm{H}}^{0.42}$ and $\propto n_{\mathrm{H}}^{0.4}$ found in the simulations of Safranek-Shrader et al. (2017) and Hu et al. (2021), respectively. The inclusion of the relevant resolution limit and floor brings the model into excellent agreement with that of B23.

Fig. 2. Fiducial model results for the Neff(H)−nH relationship. The different components and totals are annotated. The red line shows the Bisbas et al. (2023) fit, and the band shows a factor of 2 deviation from it. The secondary y-axis shows the extinction using a scalar conversion factor as described in the text. “Total w/ biases” includes the resolution limit and floor.

Fig. 3. Parameter space results for the Neff(H)−nH relationship. In each panel the annotated parameter is varied but the others are fixed to the fiducial values given in Table 1.

We also investigated the impact of varying the major underlying parameters given in Table 1. Figure 3 shows the impact of varying six different parameters on the N_eff(H)−n_H relationship. For the parameter study, we only varied one parameter at a time, keeping the rest fixed at the fiducial value. At high densities, the only parameters that have a noticeable impact are the gas temperature, which increases the Jeans radius, and the average ambient turbulent density, which globally shifts the attenuating column.

As expected, increasing the cloud size or the average turbulent density increases the turbulent column and causes the slope to flatten around n(H)≈10³ cm⁻³. Further, decreasing the strength of the turbulence or the linewidth-size slope increases the turbulent column density. For isotropic hydrodynamic turbulence, the N_eff(H)−n_H relationship maintains a relatively constant slope at all densities, consistent with the N_eff ∝ n^0.42 found in the hydrodynamic simulations of Safranek-Shrader et al. (2017). The virial parameter does not have a significant impact on the results, since it impacts ℓ and n_T in opposite directions. In Appendix B we detail the calculation of the slopes of the N_eff(H)−n_H function, compare them to the fitted power-law slopes, and show that toward high densities the model is in good agreement with simulations.

4. Discussion and conclusions

An analytic expression for the N_eff(H)−n_H relationship is advantageous for both simulations and chemical models of molecular clouds. Astrochemical codes, in particular zero-dimensional models, rely on the input of a meaningful attenuating column density along with the density, temperature, and surface FUV radiation field. Further, the N_eff(H)−n_H relationship can be used to generate a one-dimensional density distribution for rapid one-dimensional photodissociation region models (see, e.g., Bisbas et al. 2021, 2023). For observations of unresolved molecular clouds, the use of N_eff(H)−n_H one-dimensional photodissociation region models can be used with the input physical parameters to gain insights into the physics of the cloud.

In hydrodynamic simulations that model molecular clouds, performing radiative transfer can be cost-prohibitive. In such instances, this analytic model can be utilized as a sub-grid model to rapidly estimate the attenuating column density. While this will not replicate the accuracy of a full radiative transfer simulation, it is worth highlighting that even in such simulations, the spread in A_v as a function of density decreases with increasing density (see e.g. Hu et al. 2021); therefore, toward higher densities, the impact of using this model versus full ray-tracing is reduced. Finally, attenuating column densities are not just used for FUV radiation, but also for cosmic-ray energy losses into clouds (see, e.g., Padovani et al. 2009; Gaches et al. 2022).

In this Letter we have presented an analytic model for estimating the effective column density, N_eff, which acts as the attenuating column density of an isotropic external radiation field, as a function of the local hydrogen-nucleus number density, n_H. We have demonstrated that our fiducial model, which aims to represent clouds similar to those of the solar neighborhood, qualitatively reproduces the function presented by B23, which is a fit to several three-dimensional simulations with radiative transfer, and the power-law behavior toward high densities found in other works, such as Glover et al. (2010), Van Loo et al. (2013), Safranek-Shrader et al. (2017), and Hu et al. (2021). The analytic model can enable the rapid estimation of effective column densities for astrochemical models and act as a sub-grid model in simulations of molecular clouds that do not include radiative transfer. Future work will combine this model into grids of photodissociation region codes to estimate macroscopic cloud parameters from molecular line observations.

Acknowledgments

We thank the anonymous referee for their useful suggestions, which improved this Letter, and Thomas Bisbas and Daniel Seifried for their thoughtful discussions. BALG is supported by the German Research Foundation (DFG) in the form of an Emmy Noether Research Group - DFG project #542802847 (GA 3170/3-1).

References

Appendix A: Impact of density-dependent temperature and μ_p

Our model presented in the main text only considered constant temperatures and mean molecular mass per free particle, resulting in constant sound speeds. This is clearly an oversimplification, so we investigated the impact of a density-dependent sound speed, c_s(n_H), and mean molecular mass per free particle, μ_p(n_H), on the calculated N_eff(H)−n_H  relationship. We used the 1D A_V − n_H density distribution from B23 (Eq. 2) with the public photo-dissociation code 3D-PDR (Bisbas et al. 2012), including a total H₂ cosmic-ray ionization, ζ₂ = 10⁻¹⁶ s⁻¹ and an external FUV field of χ = 1 in the unit of the Draine (1978) field. We used a reduced chemical network of 33 species and 331 reactions, with initial abundances to reproduce Milky Way-like molecular clouds (see Bisbas et al. 2012, for details).

This is not a self-consistent treatment, but since the N_eff(H)−n_H  function with constant c_s and μ_p is similar to the B23 fit, the use of the photodissociation region model provides a qualitative understanding of the impact of variable c_s and μ_p. We compute c_s(n_H) and μ_p(n_H) from the gas temperature and chemical abundances directly. Figure A.1 shows the computed μ_p(n_H) and c_s(n_H), where we have shown c_s calculated with a constant and a variable μ_p.

Fig. A.1. Variation in μp (left) and cs (right) with number density, nH, using one-dimensional models from 3D-PDR. The figure shows cs assuming both a constant μp and a variable μp.

Figure A.2 shows the impact of a density-dependent c_s and μ_p. Varying μ_p has a marginal impact: it enters weakly into the Jeans length and only deviates from a constant 2.33 in regimes where the gravitational component is anyway negligible. The impact of a variable sound speed is most noticeable at lower densities, increasing the gravitational component due to two reasons. First, as shown in Eq.10, increasing the sound speed will decrease the Mach number for a given σ₀ and hence decrease n_T. Second, the increased sound speed at low density increases the Jeans length (see Eq. 9). The result of both is that the N_eff(H)−n_H  is brought slightly closer to the B23 fit. The most self-consistent method to include variable c_s and μ_p is an iteration process between a photodissociation region model and the analytic relationship.

Fig. A.2. Neff(H)−nH  relationship modified by using a variable – cs(nH) (left), μp(nH) (middle), or both cs(nH) and μp(nH) (right).

Appendix B: Power-law index of the N_eff(H)−n_H  relationship

Previous works, such as Van Loo et al. (2013) and Safranek-Shrader et al. (2017), fitted the results from their three-dimensional simulations for a polynomial expression for N_eff(H)−n_H  in dense gas. We calculated the slope of the analytic expressions from the derivative d log N_eff/d log n_H to examine the behavior in the asymptotic limits and transition regimes using the NUMPY gradient function.

Figure B.1 shows the computed slopes across the parameter space in Table 1. All models, with the exception of high values of n₀, have their slopes asymptote to the isothermal Jeans-only limit of 0.5 at high density. The power-law indices at intermediate densities are highly dependent on the assumed parameters due to changes in where the transition between turbulence and gravitational regimes occurs. In models with lower σ₀, the slope reaches a minimum between 10³ ≤ n_H ≤ 10⁴, followed by a slight steepening of the slope toward 0.5. At low density, the model produces a range of super-linear power-law slopes from 1.5 to 2, with most of the solutions asymptoting around N_eff ∝ n^1.5, with n₀ and α_V being the most impactful. The average density, n₀, is also the most impactful at the high density asymptotic limit.

The high-density power-law index asymptotes to $N_{\mathrm{eff}}\propto n_{\mathrm{H}}^{1/2}$. This is close to the slopes computed from the simulation results discussed above, which found $N_{\mathrm{eff}}\propto n_{\mathrm{H}}^{0.4}$. This is explained in Fig. B.2, which compares the slope of the fiducial N_eff(H)−n_H  to a linear fit in log-log space restricted to high density. While the constant fitted slope is consistent with previous works, the slope as a function of density varies between 0,3 and 0.6. The slightly flatter slope of the power-law fit is due to the inclusion of the turbulent column density component biasing the power-law slope fit to a flatter slope than a pure isothermal Jeans component.

Fig. B.1. Same as Fig. 3 but showing the derivative, d log Neff/d log n. The horizontal lines show a slope of 0.4 and 0.5, which correspond to the results from hydrodynamic simulation fits (see text) and the isothermal Jeans-only limit.

Fig. B.2. Comparison of the gradient calculation and a power-law fit in log-log space. Left: Neff(H)−nH  relationship for the fiducial model (black), as shown in Fig. 2, and the polynomial fit for the dense gas (red). Right: Derivative, d log Neff/d log n, as a function of density (black). The power-law fit is shown with the horizontal line. The dense gas regime with nH > 103 is highlighted with a bold line. The power-law relationship is annotated at the top.

All Tables

All Figures

	Fig. 1. Chart showing the different components and where parameters enter into the model.
In the text

	Fig. 2. Fiducial model results for the Neff(H)−nH relationship. The different components and totals are annotated. The red line shows the Bisbas et al. (2023) fit, and the band shows a factor of 2 deviation from it. The secondary y-axis shows the extinction using a scalar conversion factor as described in the text. “Total w/ biases” includes the resolution limit and floor.
In the text

	Fig. 3. Parameter space results for the Neff(H)−nH relationship. In each panel the annotated parameter is varied but the others are fixed to the fiducial values given in Table 1.
In the text

	Fig. A.1. Variation in μp (left) and cs (right) with number density, nH, using one-dimensional models from 3D-PDR. The figure shows cs assuming both a constant μp and a variable μp.
In the text

	Fig. A.2. Neff(H)−nH  relationship modified by using a variable – cs(nH) (left), μp(nH) (middle), or both cs(nH) and μp(nH) (right).
In the text

	Fig. B.1. Same as Fig. 3 but showing the derivative, d log Neff/d log n. The horizontal lines show a slope of 0.4 and 0.5, which correspond to the results from hydrodynamic simulation fits (see text) and the isothermal Jeans-only limit.
In the text

Fig. B.2. Comparison of the gradient calculation and a power-law fit in log-log space. Left: Neff(H)−nH  relationship for the fiducial model (black), as shown in Fig. 2, and the polynomial fit for the dense gas (red). Right: Derivative, d log Neff/d log n, as a function of density (black). The power-law fit is shown with the horizontal line. The dense gas regime with nH > 103 is highlighted with a bold line. The power-law relationship is annotated at the top.

In the text