© The Authors 2024

1 Introduction

Protoplanetary disks (PPDs) serve as cosmic nurseries where dust grains evolve into pebbles, planetesimals, comets, asteroids, and eventually planets. Given that comets and asteroids are less evolved than planets, analyzing their composition offers insights into the thermodynamic conditions present in the proto-solar nebula (PSN) and the temporal and radial evolution of its chemical and isotopic composition. Ices have been identified in molecular clouds and circumstellar disks (van Dishoeck 2004; Boogert et al. 2015) undergoing photochemistry-induced molecular diversity via various energetic processes (Danger et al. 2013; Öberg 2016). Complex organic molecules (COMs) are molecules composed of multiple atoms of carbon, hydrogen, and oxygen, such as CH₃CH₂OH, CH₃OCH₃, or C₂H₅OH (Belloche et al. 2009; Tenelanda-Osorio et al. 2022), and have been discovered in meteorites, comets, and star-forming regions (Briggs et al. 1992; Bisschop et al. 2007; Bergantini et al. 2017).

Laboratory experiments have demonstrated that COMs form through the alteration of astrophysical ices via UV irradiation, cosmic ray bombardment (Ehrenfreund & Charnley 2000; Shen et al. 2004), and electron (Maity et al. 2015) or ion (Moore et al. 1991; Hudson & Moore 1999) bombardment in PSN conditions. Among these ices, methanol is particularly abundant (Pontoppidan et al. 2003; Gibb et al. 2004; van Dishoeck 2014) and serves as an ideal molecule for generating significant molecular diversity (Öberg et al. 2009; Cuppen et al. 2011): Lyman-α irradiation of pure methanol ice at 20 K results in the detection of over 35 gas-phase molecules with 2–6 carbon atoms (Abou Mrad et al. 2014). The generation of this diverse array of COMs primarily depends on the dose received by the ice and the temperature at which the ice undergoes processing (Theulé et al. 2013; Tenelanda-Osorio et al. 2022).

While this study primarily investigates ice processes within the PSN, it is important to note that alternative experiments suggest COMs could have also formed in the interstellar medium (ISM), where the local temperature is ~10 K in dark molecular clouds (Combes et al. 1995; Qasim et al. 2020; Carraro 2021; Fedoseev et al. 2022). However, due to the complexity of PPDs involving various physical processes (turbulence, diffusion, advection, UV irradiation, etc.), the formation conditions of COMs have not been extensively explored in these environments. Irradiation levels notably fluctuate across the disk, ranging from higher rates at the surface to near absence in the midplane region (Ciesla & Sandford 2012; Mousis et al. 2018). Consequently, the trajectory of icy particles in PPDs can significantly influence their photochemical evolution.

In this study we computed the irradiation dose received by grains and icy pebbles during their radial and vertical transport using a time-dependent PPD model derived from Aguichine et al. (2020) and Schneeberger et al. (2023), following a Lagrangian scheme (Ciesla 2010, 2011; Ronnet et al. 2017). The results were then compared with experimental data derived by Tenelanda-Osorio et al. (2022), who irradiated methanol ices with a UV photon flux in a vacuum chamber to track COM formation and composition. This comparison allowed us to evaluate whether the formation conditions of COMs from methanol ice identified in the laboratory could occur in certain regions of the PSN.

Section 2 describes the laboratory experiments used to benchmark our transport and irradiation model. Section 3 details the different modules of our model, while Sect. 4 presents the comparisons between our numerical results and the experiments. Finally, Sect. 5 is dedicated to a discussion and conclusions.

2 Experimental frame

The experimental data utilized in our paper were derived from the research conducted by Tenelanda-Osorio et al. (2022). A pure methanol ice underwent irradiation at Lyman-α under varying doses with a lamp flux of 10¹³ photons cm⁻² s⁻¹ (Gerakines et al. 1996). Irradiation times of up to 24 hours were tested at both 20 and 80 K, resulting in an accumulation of up to 8.64 × 10¹⁷ photons/cm² on the icy sample. The icy sample consists of a thin layer of methanol, containing approximately 1.2 × 10¹⁷ molecules per cm², resulting in a photon-to-molecule ratio of 7.2. After irradiation, the sample was returned to room temperature, and the volatiles were released into the gas phase. Volatile organic compounds were then characterized using a gas chromatograph coupled with a mass spectrometer. Twenty-one different organic molecules containing up to five carbon atoms (including alcohols, aldehydes, ketones, esters, and ethers) were identified. The results demonstrate that molecular diversity increases with the irradiation dose from only MeOH. Interestingly, the abundances of C₁ and C₂ molecules tend to decrease or stabilize with increasing irradiation dose, while heavier molecules emerge. However, lower abundances of heavier molecules are present, consistent with radical reactivity. Among all the chemical families observed, ketones and esters contain the molecules with the largest number of carbon atoms.

Experiments conducted at 80 K require a higher irradiation dose to achieve the same molecular complexity observed at 20 K. In both cases, the data analyses indicate that 24 hours of irradiation suffice to enable COM formation at both temperatures, corresponding to a dose of 8.64 × 10¹⁷ photons cm⁻².

Our modeling aims to ascertain the conditions under which COMs form from irradiated methanol icy pebbles or grains during their migration in the PSN. To achieve this, CH₃OH-rich particles are released in PSN regions where the temperature is initially set to 20 or 80 K, the two temperatures considered in the experimental setup. Irradiation accumulates along the particles’ drift until it reaches the necessary doses for COM formation (8.64 × 10¹⁷ photons cm⁻²). Successful particles are those that reach the irradiation dose identified in experiments before crossing the methanol iceline, at approximately 124 K. Once the particles cross this iceline, they are assumed to sublimate.

3 Methodology

This section describes the modules employed for computing the evolution of the disk, as well as the transport and irradiation of the particles.

3.1 Structure of the disk

Our time-dependent accretion disk model is computed in the framework of the α-turbulent viscosity formalism, in which the viscosity is defined as $v=\alpha c_{\text{s}}^2/{\text{Ω}}_{\text{K}}$ (Shakura & Sunyaev 1973). Here, ${\text{Ω}}_{\text{K}}=\sqrt{G{M}_{\star }/r^3}$ is the Keplerian frequency, M_★ is the mass of the star, r is the distance to the star, and G is the gravitational constant. $c_{\text{s}}=\sqrt{R_{\text{g}}{T}_{\text{m}}/\mu }$ is the isothermal sound speed, R_g is the ideal gas constant, T_m is the midplane temperature of the disk, and µ = 2.31 g mol⁻¹ is the mean molar mass of the gas. α is the viscosity parameter, usually set in the 10⁻⁴–10⁻² range (Hartmann et al. 1998; Hueso & Guillot 2005; Desch et al. 2017), which we set to 10⁻³ in our model.

Our model computes the evolution of the gas surface density, Σ_g, at hydrostatic equilibrium via the following differential equation (Lynden-Bell & Pringle 1974), $$\frac{\partial {\text{Σ}}_{\text{g}}}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{\frac{1}{2}}\frac{\partial }{\partial r}\left(r^{\frac{1}{2}}{\text{Σ}}_{g}v\right)\right],$$(1)

which can be rewritten as $$\begin{array}{c}\frac{\partial {\text{Σ}}_{\text{g}}}{\partial t}=\frac{1}{2\pi r}\frac{\partial \dot{M}}{\partial r}\\ \dot{M}=3\pi {\text{Σ}}_{\text{g}}v\left(1+2\frac{\text{d}\ln \left({\text{Σ}}_{\text{g}}v\right)}{\text{d}\ln r}\right),\end{array}$$(2)

where Ṁ is the mass accretion rate.

Our code resolves this equation iteratively, using the following initial conditions (Lynden-Bell & Pringle 1974; Aguichine et al. 2020): $$\{\begin{array}{l}{\text{Σ}}_{\text{g},0}=\frac{{\dot{M}}_{\text{acc},0}}{3\pi \nu }{e}^{-{\left(\frac{r}{r_c}\right)}^{0.5}}\hfill \\ {\dot{M}}_0={\dot{M}}_{\text{acc},0}\left(1-{\left(\frac{r}{r_{\text{c}}}\right)}^{0.5}\right)e^{-{\left(\frac{r}{r_{\text{c}}}\right)}^{0.5}},\hfill \end{array}$$(3)

where r_c is the centrifugal radius (computed here as 1.2 AU at initialization), and Ṁ_acc,0, the accretion rate, is assumed to be 10^−7.6 M_⊙ yr⁻¹.

The temperature in the midplane is given by the sum of all the heating sources (Hueso & Guillot 2005): $$T_{\text{m}}^4=\frac{1}{2{\sigma }_{\text{sb}}}\left(\frac{3}{8}{\tau }_{\text{R}}+\frac{1}{2{\tau }_{\text{P}}}\right){\text{Σ}}_{\text{g}}v{\text{Ω}}_{\text{K}}^2+T_{\text{amb}}^4,$$(4)

where σ_sb is the Boltzmann constant. ${\tau }_{\text{R}}=\frac{{\kappa }_{\text{R}}{\text{Σ}}_{\text{g}}}{2}$ and τ_P = 2.4τ_R are the Rosseland and Planck optical depth, respectively. The mean Rosseland opacity, ${\kappa }_{\text{R}}={\kappa }_0{\rho }_{\text{m}}^a{T}^b$, is a function of the gas density and temperature in the midplane, with κ₀, a, and b determined by observations (see Bell & Lin 1994 and our Appendix A). Irradiation from the young star is neglected due to the assumption of shadowing in the outer part of the disk (Ohno & Ueda 2021). This shadowing is attributed to dust pileup at the water snowline and other icelines, which are inherent in our model (Aguichine et al. 2020; Schneeberger et al. 2023). Finally, T_amb = 10 K is the ISM temperature in dark molecular clouds.

Assuming the disk is vertically isothermal (a common assumption, as the adiabatic disk is isothermal except in the innerregions; Ciesla 2010; Ronnet et al. 2017), its vertical density, temperature, and pressure profiles are derived as follows: $${\rho }_{\text{g}}(r,z)={\rho }_{\text{m}}(r)e^{-\frac{z^2}{2H^2}},T(z,r)=T_{\text{m}}(r),P(z,r)=P_{\text{m}}(r)e^{-\frac{z^2}{2H^2}},$$(5)

with ρ_m, T_m, and P_m corresponding to the midplane density, temperature, and pressure respectively, and z is the vertical position. From those quantities, we define the gas disk scale height H = c_s/Ω_K. We refer the reader to Aguichine et al. (2020) for a full description of our disk model.

The temperature and pressure at which methanol condenses in the disk is computed using the Antoine’s law: $${\log }_{10}\left(P_{\text{vap\hspace{0.17em}}}\right)=A-\frac{B}{T+C}$$(6)

where A = 5.20409, B = 1581.341, C = −33.5, as provided by Ambrose & Sprake (1970) and obtained from the NIST Chemistry Webbook¹, and P_vap is the vapor pressure in bar. The methanol iceline is defined as the location in the disk where the partial pressure of methanol is equal to its equilibrium pressure. The methanol gas phase abundance assumed in the disk is 1.4018 × 10⁻⁶ with respect to H₂, based on an estimate of the CH₃OH/H₂O ratio in Comet 67P/Churyumov-Gerasimenko (Le Roy et al. 2015), and on a plausible water abundance with respect to H₂ in the PSN (Schneeberger et al. 2023).

3.2 Irradiation of particles

The irradiation of the disk was computed via the following relation for a flux F(r, z) of UV photons (Ciesla 2010; Ciesla & Sandford 2012): $$F(r,z)=F_0{e}^{-\tau (r,z)}$$(7)

with $$\tau (r,z)={{{\displaystyle \int }}^{\text{}}}_{\left|z\right|}^{\infty }{\rho }_{\text{g}}(r,z)\kappa \text{d}z,$$(8)

where τ(r, z) represents the optical depth of suspended matter above the particle, calculated at a given altitude and distance from the star, and κ denotes the opacity, which is assumed to be equal to the frequency averaged mean Rosseland opacity derived from experimental data (Bell & Lin 1994; Aguichine et al. 2020; Schneeberger et al. 2023). F₀ is the incident UV interstellar flux equal to 10⁸ photons cm⁻² s⁻¹ (= 1 G₀). 1 G₀ (in Habing units) is equal to 10¹² photons m⁻² s⁻¹, and corresponds to the quantity characterizing the average interstellar radiation field with an energy hv bracketed between 6 and 13.6 eV (between 91.2 and 206.7 nm). Here, v is the frequency and h is the Planck constant. G₀ can be expressed as follows (Yeghikyan 2009): $$G_0=\frac{4\pi {{{\displaystyle \int }}^{\text{}}}_{v(206.7\text{\hspace{0.17em}\hspace{0.17em}nm})}^{v(91.2\text{\hspace{0.17em}\hspace{0.17em}nm})}{J}_{\nu }\text{d}v}{1.6\cdot {10}^{-3}},$$(9)

with J_v the intensity of the radiation in erg cm⁻² s⁻¹ sr⁻¹ Hz⁻¹. We used G₀ = 9.986 · 10⁸ eV cm⁻² s⁻¹ for the UV interstellar (isotropic) photon flux (Habing 1968; Parravano et al. 2000; Yeghikyan 2009; Ciesla & Sandford 2012). However, depending on the source, the incident interstellar source can be a factor of 10 higher if it comes from the spiral arms of a galaxy (Yeghikyan 2009), and from 10 to 3× 10⁴ higher if it comes from clusters, where the median flux value experienced by a forming stellar system is 900 G₀ (Adams et al. 2006). We assumed in our calculations that the incident flux F₀ arrives perpendicularly to the disk, with the same intensity at each point of its surface, and that the irradiation from the central star is neglected (Ohno & Ueda 2021).

3.3 Vertical and radial transport of particles

For this study, we used particles ranging in size from 1 µm to 1 cm, with a density of 1 g/cm³. We assumed that the simulated particles do not collide with each other and that they remain constant in size. Given their size and density, the particles evolve according to the Epstein regime. In the PSN environment, particles are subject to specific drag regimes: the Stokes regime for larger particles and the Epstein regime for smaller particles. Particles are in the Epstein regime if they are smaller than the mean free path of the gas, $\lambda =\frac{m_{H_2}}{\pi {\rho }_g{r}_{H_2}^2}$, times 2.25 (Stepinski & Valageas 1996; Supulver & Lin 2000), where $m_{H_2}$ and $r_{H_2}$ are the mass and radius of H₂ molecules, and ρ_𝑔 is the density of the H₂ gas. In our model, particles up to 1 cm in size are in the Epstein regime when they are located beyond 0.5 AU from the star. Furthermore, considering their Stokes number, only the larger particles, approximately 1 cm in diameter, are subject to radial drift, while the smaller particles remain predominantly coupled to the gas. This is supported by the particles’ Stokes numbers, which range from 10⁻² to 10⁻⁷ (see Appendix B).

3.3.1 Vertical transport

The vertical transport of particles is depicted following a Lagrangian approach (Ciesla 2010). Their vertical evolution is depicted via the following advection diffusion equation (Dubrulle et al. 1995; Gail 2001; Ciesla 2010): $$\frac{\partial {\rho }_i}{\partial t}=\frac{\partial }{\partial z}\left({\rho }_{\text{g}}D\frac{\partial \frac{{\rho }_i}{{\rho }_{\text{g}}}}{\partial z}\right)-\frac{\partial }{\partial z}\left({\rho }_i{v}_z\right),$$(10)

where ρ_i is the material density, ρ_g the gas density, $D=\frac{v}{1+S{t}^2}$ the diffusivity, St the Stokes number, and υ_z the dust vertical velocity due to gravitational forces.

The vertical motion of an individual particle is computed using a Monte Carlo approach considering the random motion due to the viscous effects from eddies in the turbulent fluid (Ciesla 2010): $$z_i=z_{i-1}+v_{\text{adv},z}\delta t+R_1{\left[\frac{2}{{\sigma }^2}D(z)\delta t\right]}^{\frac{1}{2}},$$(11)

where z_i−1 is the initial vertical position of the particle. z is its position after a given time step $\delta t=\frac{1}{{\text{Ω}}_K}\cdot R_1\in [-1;1]$ is a random number, and σ² is its distribution variance (= 1/3 for uniform distribution) representing the random movement due to turbulence (Visser 1997). Hence, two particles beginning their trajectories at the same location will move along distinct paths.

The advection velocity term has three components coming from the advective term analogy of the advection diffusion equation (Ciesla 2010; Ronnet et al. 2017; Mousis et al. 2018): $${\upsilon }_{\text{adv},z}={\upsilon }_z+{\upsilon }_{gas,z}+\frac{\partial D}{\partial z}.$$(12)

The υ_z is then given by the terminal velocities of the particles, with the drag forces taken into account (Cuzzi & Weidenschilling 2006; Ciesla 2010): $${\upsilon }_z=-t_S{\text{Ω}}_K^{2}z.$$(13)

Here, t_s is the stopping time, which is the timescale required for the dust, assumed to consist of spherical grains, to transfer all its angular momentum to the gas. In the Epstein regime, t_s is given by (Perets & Murray-Clay 2011; Mousis et al. 2018) $$t_{\text{s}}=\frac{{\rho }_{\text{s}}}{{\rho }_{\text{g}}}\frac{R_{\text{s}}}{{\upsilon }_{\text{th}}},$$(14)

with ${\upsilon }_{th}=\sqrt{8/\pi }{c}_{\text{s}}$ the thermal sound velocity, R_s the solid dust grain radius and ρ_s its density (here 1 g cm⁻³).

The second term υ_gas,z of Eq. (12) is derived from the advection term of Eq. (10), and comes from the vertical variation of density of nebular gas, which is the medium through which the species is diffusing (Ciesla 2010; Ronnet et al. 2017; Mousis et al. 2018): $${\upsilon }_{\text{gas},z}=\frac{D}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial z}.$$(15)

Combining Eqs. (5) and (15), we derive ${\upsilon }_{gas,z}=-D\frac{z}{H^2}$.

The final term in Eq. (12) considers the vertical variation of the diffusion coefficient, with its value assumed to be zero under the assumption of a vertically isothermal disk.

3.3.2 Radial transport

A similar method, based on the one-dimensional radial advection–diffusion equation, was used to compute the radial transport of particles (Gail 2001; Ciesla 2011): $$\frac{\partial {\rho }_i}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{\rho }_{\text{g}}D\frac{\partial \frac{{\rho }_i}{{\rho }_{\text{g}}}}{\partial r}\right)-\frac{1}{r}\frac{\partial }{\partial r}\left(r{\rho }_i{\upsilon }_{\text{r}}\right).$$(16)

In Cartesian coordinates, the equation becomes (Ciesla 2011) $$\frac{\partial {\rho }_i}{\partial t}=\frac{\partial }{\partial x}\left({\rho }_{\text{g}}D\frac{\partial \frac{{\rho }_i}{{\rho }_{\text{g}}}}{\partial x}\right)-\frac{\partial }{\partial x}\left({\rho }_i{\upsilon }_x\right)+\frac{\partial }{\partial y}\left({\rho }_{\text{g}}D\frac{\partial \frac{{\rho }_i}{{\rho }_{\text{g}}}}{\partial y}\right)-\frac{\partial }{\partial y}\left({\rho }_i{\upsilon }_y\right),$$(17)

where υ_x and υ_y are the components of υ_r. This allows the description of the evolution of particles along the x and y axes in a way similar as the one provided for their motion along the z axis (Ciesla 2011): $$\begin{array}{c}{x}_i=x_{i-1}+{\upsilon }_{\text{adv},x}\delta t+R_1{\left[\frac{2}{{\sigma }^2}D(x)\delta t\right]}^{\frac{1}{2}}\\ y_i=y_{i-1}+{\upsilon }_{\text{adv},y}\delta t+R_1{\left[\frac{2}{{\sigma }^2}D(y)\delta t\right]}^{\frac{1}{2}},\end{array}$$(18)

where x_i−1 (resp. y_i−1) is the initial horizontal x (resp. y) position of the particle. x_i (resp. y_i) is its horizontal x (resp. y) position after a given time step $\delta t=\frac{1}{{\text{Ω}}_K}\cdot R_1\in [-1;1]$ is a random number, and σ² is its distribution variance (= 1/3 for uniform distribution) representing the random movement due to turbulence (Visser 1997).

The advection velocity along the x axis becomes the sum of the three terms: $${\upsilon }_{\text{adv},x}={\upsilon }_x+{\upsilon }_{\text{gas},x}+\frac{\partial D}{\partial r}\frac{x_{i-1}}{r_{i-1}}.$$(19)

The first term is the velocity of the material due to gas drag: $${\upsilon }_x\left(x_{i-1},y_{i-1}\right)={\upsilon }_r\frac{x_{i-1}}{r_{i-1}},$$(20)

where $${\upsilon }_r=-\frac{1}{1+\text{S}{\text{t}}^2}\frac{3}{2}\frac{v}{r}+\frac{2\text{St}}{1+\text{S}{\text{t}}^2}\frac{c_s^2}{r{\text{Ω}}_K}\frac{\text{d}\ln P}{\text{d}\ln r}.$$(21)

The first term of υ_r represents the radial velocity due to the gas velocity for a steady disk with a mass profile in r⁻¹ (Lynden-Bell & Pringle 1974; Clarke & Pringle 1988; Gail 2001; Cuzzi et al. 2003; Ciesla 2011; Guillot et al. 2014), while the second term accounts for the radial drift (Weidenschilling 1977; Nakagawa et al. 1986; Birnstiel et al. 2012; Aguichine et al. 2020).

The second term in Eq. (19) is the same as its equivalent along the z axis and accounts for the density gradient in the midplane: $${\upsilon }_{\text{gas},x}=\frac{D}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial r}\frac{x_{i-1}}{r_{i-1}}.$$(22)

The third term of Eq. (19), $\frac{\partial D}{\partial t}\frac{x_{i-1}}{r_{i-1}}$, is identical to the last term of Eq. (12), and corresponds to the contribution due to the diffusivity gradient along the x-axis. Equations along the y-axis can be developed in the same way as for the x-axis since axial symmetry is assumed.

4 Results

We conducted our computations with the goal of replicating the experimental conditions outlined in Sect. 2. This methodology enabled us to explore the conditions under which COMs form within the PSN. Each particle trajectory presented was averaged over a total of 500 independent particles released in the midplane under identical conditions, ensuring a statistically significant sample.

4.1 Map of disk irradiation

Figure 1 shows a two-dimensional map of F_acc, namely the UV irradiation accumulated by the disk over its time evolution. As the disk’s density shrinks and flattens over time, radiation progressively penetrates deeper into the midplane. After 10 kyr of disk evolution, the black zone in the figure, where the accumulated irradiation is lower than the most important dose obtained during the experiments of interest (8.64 × 10¹⁷ photons cm⁻²), extends radially up to ~ 16.5 AU, and up to a height of ~0.75 AU. After 100 kyr of irradiation, it reduces radially down to ~15.1 AU, and slightly vertically down to ~0.73 AU. After 1000 kyr of disk evolution, the black zone continues its shrinking radially down to ~12.7 AU, and vertically down to 0.31 AU. At any epoch of the disk evolution, its outer regions reach accumulated irradiation levels exceeding by far the most important dose obtained during the experiments.

4.2 Drift of irradiated particles

Figure 2 illustrates the positions of particles in the disk tracked over 300 kyr of PSN evolution, with each particle assumed to have a density of 1 g/cm⁻³. The particle trajectory shown in Fig. 2 represents the mean trajectory of 500 simulated particles. The insets on the right display the trajectory of an individual particle. Unlike the mean trajectories, which remain confined near the midplane, the individual particle trajectories exhibit greater vertical variation. However, the average vertical trajectory remains close to ɀ = 0 AU, due to the disk’s symmetry about the midplane. The top three panels depict the averaged trajectories of particles sized at 1 cm, 100 µm, and 1 µm, respectively, starting their migration from the midplane region of the disk, where the local temperature is 20 K (approximately 11 AU from the star). The bottom three panels exhibit the averaged trajectories of particles with the same size distribution, but originating from the midplane of the disk where the local temperature is 80 K (approximately 6.1 AU from the star).

On average, particles initiating their trajectories in the 20 K midplane region (panels a, b, c) require over 1 Myr of PSN evolution to cross the methanol iceline. The vertical trajectory of a particle is highly dependent on particle size. Specifically, 1 cm particles remain more confined to the midplane and never reach heights greater than 0.63 AU, whereas 1 µm particles are more likely to attain higher elevations, reaching up to 1.53 AU (as shown in the insets at the right of Fig. 2). On average, particles ranging from 1 cm to 1 µm cross the 10.2–10.8 AU region after 300 kyr of PSN evolution.

On average, particles that begin their trajectory at a local temperature of 80 K in the PSN (see panels d, e, and f in Fig. 2) require more than 300 kyr to cross the methanol iceline and sublimate: approximately 530 kyr for 1 cm particles, and over 1 Myr for 100 µm and 1 µm particles. After 300 kyr of PSN evolution, the mean radial trajectories of particles ranging from 100 µm to 1 µm intersect the ~5.7 AU region, while the mean radial trajectory of 1 cm particles reaches the ~4.4 AU region. The vertical trajectories also exhibit a distribution based on particle size (as shown in the insets at the right of Fig. 2): 1 cm particles are confined to the midplane and never reach altitudes above approximately 0.53 AU, whereas 1 µm particles can attain higher altitudes, up to about 1.23 AU.

Figure 3 illustrates the average temperature and pressure of particles sized at 1 cm, 100 µm, and 1 µm along their trajectories in the PSN, as well as their average accumulated irradiation, assuming the same initial conditions as those displayed in Fig. 2. The average accumulated irradiation is determined by averaging the irradiation accumulated by each individual particle. In the following, terms such as accumulated irradiation, temperature, or pressure denote the average accumulated irradiation, average temperature, and average pressure of the 500 simulated particles.

On average, smaller particles measuring 1 and 100 µm, originating at a temperature of 20 K, receive the maximum accumulated irradiation in the experiments (8.64 × 10¹⁷ photons cm⁻²) after approximately 15.2–21.4 kyr of PSN evolution. In contrast, larger particles with a size of 1 cm would, on average, require approximately 911 kyr to accumulate the equivalent laboratory dose.

After 300 kyr of PSN evolution, smaller particles accumulate approximately 2.1–3.2 ×10¹⁹ photons cm⁻², whereas larger pebbles accumulate about 6.2 ×10¹⁶ photons cm⁻², on average. Additionally, at the end of this period, particles attain temperatures and pressures in the range of approximately 21–25 K and 2.5–4 mPa, on average, respectively.

Particles that start at a temperature of 80 K reach their maximum accumulated irradiation in the experiments after an average of 141 kyr and 359 kyr of PSN evolution for the 1 µm and 100 µm sized particles, respectively. In contrast, 1 cm sized particles do not reach the irradiation dose even after 1 Myr. By the end of their respective time spans, the 1 µm and 100 µm particles attain temperatures and pressures in the ranges of approximately 77–78 K and 21–23 mPa on average, respectively. In comparison, 1 cm sized particles reach temperatures and pressures of 127 K and 405 mPa on average, respectively.

Fig. 1 Two-dimensional map of the accumulated irradiation after 10, 100, and 1000 kyr of the PSN’s evolution. The black area around the midplane corresponds to the region where the accumulated irradiation is lower than the most important dose obtained during the experiments of interest (8.64 × 1017 photons cm−2).

5 Discussion and conclusion

We investigated the influence of disk parameters on the mean trajectory and the average accumulated irradiation of the particles.

The opacity used in this model is the mean Rosseland opacity, derived from Bell & Lin (1994), which is averaged across all frequencies and may inadequately represent UV-range photon irradiation. To account for potential changes in opacity, we employed an opacity of 48 cm²/g following the model of Birnstiel et al. (2018) for particles up to 10 cm in size, averaged over the UV wavelength range of 122–160 nm used in the experiments by Tenelanda-Osorio et al. (2022). We conducted new simulations for 1 µm particles, which are more likely to receive higher irradiation doses.

After 300 kyr of PSN evolution, the 1 µm particles released at 11 AU (~20 K) receive an average accumulated irradiation dose of about ~3.2×10¹⁶ photons/cm⁻² with this new opacity, which is several orders of magnitude lower than the average accumulated irradiation dose of ~3.2×10¹⁹ photons/cm⁻² in our nominal case. Similarly, 1 µm particles released at 6.1 AU (~80 K) receive an average accumulated irradiation dose of about ~2.3×10¹⁶ photons/cm⁻² with the new opacity, compared to the ~2.6× 10¹⁸ photons/cm⁻² in our nominal case. This new opacity of 48 cm²/g results in a lower irradiation of the particles in our model, and on average, the particles do not receive the experimental irradiation dose. However, it should be noted that although this opacity is more suitable for UV radiation, it remains constant throughout the disk and does not vary with temperature or density, unlike our nominal opacity. Due to this significant limitation, we used the standard opacity data from Bell & Lin (1994).

An increase in the viscosity parameter, α, generally extends the particle trajectory in the PSN and tends to increase their irradiation exposure. For instance, with α = 1 × 10⁻² (10 times larger than the nominal value), the mean trajectory of 500 particles sized at 1 cm and released at 20 K begins at approximately 14 AU, which is about 3 AU farther than in our nominal computation. The particles drift inward in the PSN at a velocity approximately 3.5 times higher. They do not cross the methanoliceline even after 1 Myr of PSN evolution. Furthermore, the irradiation accumulated by the particles matches the laboratory dose after only approximately 16 kyr of PSN evolution, more than 50 times faster than in our nominal case. This acceleration happens because the outer region of the disk is optically thinner than the inner region, and a more viscous disk allows particles to reach greater altitudes. Particles measuring 1 cm and released at a local temperature of 80 K commence their journey at approximately 6.6 AU (slightly farther than in our nominal case) and drift approximately 4 times faster than in our nominal case. However, they do not cross the methanol iceline, even after 1 Myr of PSN evolution. Furthermore, they still do not receive the experimental irradiation dose after this duration. On the other hand, 1 µm sized particles released at a local temperature of 80 K receive the laboratory irradiation dose in approximately 36 kyr of PSN evolution (roughly 4 times faster than in our nominal case) and drift approximately 5 times faster than in our nominal case. However, they still do not cross the methanol iceline during the first million years of PSN evolution.

Conversely, assuming α = 1 × 10⁻⁴ (10 times less than the nominal value), 1 cm particles released at a local temperature of 20 K begin their migration at a distance of approximately 7.3 AU (about 3.7 AU closer to the star than in our nominal case) and drift inward in the PSN at a velocity approximately 3 times lower. They do not reach the irradiation threshold (accumulating approximately 10¹¹ times less irradiation than in our nominal case) and do not cross the methanol iceline, even after 1 Myr of PSN evolution. Similarly, particles measuring 1 cm released at a local temperature of 80 K start their drift at a radius of approximately 4.6 AU, very close to the methanol iceline, which is located at around 3.9 AU. They drift inward approximately 3 times more slowly than in our nominal case and do not cross the methanol iceline within 1 Myr of PSN evolution. Additionally, they do not receive the irradiation dose on this timescale. Particles measuring 1 µm released at a local temperature of 80 K also start at 4.6 AU and drift approximately 7 times more slowly than in our nominal case. They also fail to cross the iceline within 1 Myr of PSN evolution and experience irradiation at a rate approximately 7 times lower than in our nominal case. Therefore, the particles also fail to receive the irradiation dose during this simulation time frame.

Another key parameter is the initial disk’s mass accretion rate. A variation by a factor of 2 does not lead to significant differences. For instance, 1 cm sized particles starting from a local temperature of 20 K evolving in a disk with Ṁ₀ = 5 × 10⁻⁷ M_⊙/year (2 times higher than in our nominal case) receive approximately 1.6 times less irradiation and drift around 1.1 times faster. With Ṁ₀ = 1 × 10⁻⁷ M_⊙/year (2.5 times less than our nominal value), those particles receive roughly 1.5 times more irradiation and drift approximately 1.1 times more slowly. The 1 cm sized particles released at a local temperature of 80 K with Ṁ₀ = 5 × 10⁻⁷ M_⊙/year accumulate 10¹⁴ more irradiation than in our nominal case, and drift ~1.7 times more slowly. However, with Ṁ₀ = 1 × 10⁻⁷ M_⊙/year, those particles accumulate 10¹² times more irradiation than in our nominal case and drift 1.9 times more slowly. 1 µm sized particles released at a local temperature of 80 K with Ṁ₀ = 5 × 10⁻⁷ M_⊙/year drift inward approximately 1.1 times faster than in our nominal case and receive the irradiation dose about 120 kyr later, at 261 kyr of PSN evolution. With Ṁ₀ = 1 × 10⁻⁷ M_⊙/year, those particles drift approximately 1.2 times more slowly than in our nominal case and receive the irradiation dose after 144 kyr of PSN evolution, similar to our nominal case. This suggests that reducing Ṁ₀ cools the disk, leading to lower radial velocities for particles; however, even an initial mass accretion rate that is half of the nominal value does not significantly alter the results.

Another parameter influencing the trajectory of particles is their density. 1 cm sized particles at 20 K with twice the nominal density (2 instead of 1 g cm⁻³) receive approximately 20 times less irradiation (2.2 × 10¹⁴ instead of 4.4 × 10¹⁵ photons cm⁻² in 100 kyr of PSN evolution) because they remain more confined in the midplane, and they drift 1.6 times faster. Conversely, with a density half that of our nominal density (0.5 instead of 1 g cm⁻³), the particles are likely to reach higher altitudes, experience nearly 10 times more irradiation (4.3 × 10¹⁶ photons cm⁻² over 100 kyr of PSN evolution), and drift 1.5 times more slowly. In contrast, 1 µm sized particles released at 20 K with a density of 0.5 g cm⁻³ exhibit no significant differences from our nominal case in terms of trajectories or irradiation (the less dense particles accumulate 8.5 × 10¹⁸ instead of 8.44 × 10¹⁸ photons cm⁻² after 100 kyr of PSN evolution). Similar conclusions can be drawn with particles twice as dense, which are irradiated with 8.35 × 10¹⁸ photons.cm⁻² after 100 kyr of PSN evolution. Finally, 1 cm particles released at a local temperature of 80 K show similar results when we vary their density: Particles with a density twice as high as our nominal case (2 g cm⁻³ instead of 1 g cm⁻³) drift 1.7 times faster and receive less irradiation (less than 1 instead of 1.1 × 10⁶ photons cm⁻² over 100 kyr of PSN evolution). Conversely, particles with a density lower than in our nominal case (0.5 g cm⁻³ instead of 1 g cm⁻³) drift 1.5 times more slowly and receive more irradiation (3.3 × 10⁹ photons cm⁻² over 100 kyr of PSN evolution). The variation in particle density affects their radial evolution only slightly but significantly impacts their vertical trajectories and, therefore, their irradiation, particularly for larger particles measuring 1 cm.

Our findings indicate that 1 cm particles experience lower levels of irradiation due to gravitational forces floating them in the midplane. Conversely, smaller particles measuring 100 and 1 µm are prone to ascend to higher altitudes and thus experience increased irradiation. When released at a local temperature of 80 K (approximately 6.1 AU from the star), these particles require 359 and 141 kyr of PSN evolution, respectively, to accumulate sufficient irradiation for substantial molecular diversity to emerge (Tenelanda-Osorio et al. 2022). However, the irradiation timescale drastically reduces to less than 22 kyr for the formation of COMs when particles are released at a local temperature of 20 K (around 11 AU from the star). This notably brief timescale, relative to disk evolution, suggests that methanol-rich ices transported by small particles in the outer disk regions (particularly those beyond approximately 11 AU) efficiently convert into COMs before consolidating into larger solids. The growth timescale of particles is approximately 2.5 kyr at 11 AU after 10 kyr of PSN evolution, is shorter at greater distances, and is longer closer to the star (approximately 3.5 kyr at 6.1 AU). This growth timescale increases over time, reaching 80 kyr at 11 AU after 1 Myr of PSN evolution (Birnstiel et al. 2012; Aguichine et al. 2020).

Consequently, COMs could be prevalent in the outer regions of the PSN. Furthermore, particles smaller than 1 cm will not cross the iceline and therefore will not release molecules in the vapor phase during the first million years of PSN evolution if they are released at distances greater than ~6 AU. To cross the methanol iceline, particles must be either similar to or larger than the 1 cm particles simulated here or be released closer to the iceline. Additionally, Misener et al. (2019) has demonstrated that the growth and fragmentation of particles strongly influence their trajectories. Integrating the evolution of grain growth with their trajectories, and consequently their irradiation, presents an interesting avenue for future research.

Furthermore, it is crucial to consider that the experimental photon flux exceeds that in our simulations. In the experiments, an irradiation of 10¹³ photons s⁻¹ cm⁻² results in an experimental ratio of 7.2 photons per molecule, compared to a maximum irradiation in the disk of 10⁸ photons s⁻¹ cm⁻² in our simulations. The higher photon dose rate in laboratory experiments, relative to PSN conditions, likely enhances reactivity within the sample. This increased reactivity can promote the formation of more diverse and complex organic products, as recently formed radicals are exposed to additional photons. Conversely, the lower dose rate typical in PSN conditions allows for the recombination of radicals into their precursors. Hence, the same total irradiation dose would likely lead to the formation of fewer COMs in the PSN compared to laboratory experiments.

Rocky particles coated with an organic mantle (OMGs) exhibit greater adhesiveness compared to particles composed solely of organics or rocks (Kouchi et al. 2002; Homma et al. 2019; Bischoff et al. 2020). The enhanced fragmentation threshold velocity of OMGs within the 1.1–3.4 AU region (corresponding to a local temperature range of approximately 200–400 K) facilitates the formation of planetesimals (Kouchi et al. 2002; Homma et al. 2019). Nevertheless, our simulations suggest that particles released at local temperatures of 20 K or 80 K typically never reach this zone, even after 1 Myr of PSN evolution. Future studies should prioritize examining the evolution of newly formed COMs to explore their potential condensation as mantles on rocky particles within the region delineated by Kouchi et al. (2002). Another intriguing avenue of research would involve investigating the formation of COMs from molecules other than methanol, which may sublime at higher temperatures or require greater irradiation to develop such molecular diversity. To investigate this, integrating a chemical photoreaction model, such as those proposed by Takehara et al. (2022) or Ochiai et al. (2024), with the particle transport model could be highly beneficial. This approach could potentially facilitate the development of organic mantles around pebbles in regions characterized by the stickiest conditions in simulations.

Our primary focus has been on interstellar cosmic irradiation, which has led us to overlook other significant processes such as cosmic ray interactions with hydrogen gas. These interactions generate secondary UV radiation, albeit with low irradiation values of 1.4 × 10³ photons cm⁻² s⁻¹, in contrast to our nominal values (Prasad & Tarafdar 1983). Additionally, the effect of electron irradiation, though often disregarded, remains an important consideration. Furthermore, potential hydrogenation processes, known for their ability to destruct molecules (Linnartz et al. 2011), have not been factored into our analysis. These processes, despite appearing negligible, can contribute to irradiation in opaque regions of the disk, such as the midplane. Future studies should thus attempt to accurately quantify their impacts.

Planets, moons, and asteroids have grown from material originating from the PSN. Thus, our model offers vital insights into the drift ranges and thermodynamic conditions under which pebbles can foster molecular diversity from methanol ice. Should the distance range for COMs formation, as delineated by our model, coincide with the formation zones of planets and other celestial bodies, these bodies might have, at least partially, evolved from these COMs. Given the analogous properties between PPDs and circumplanetary disks, the particle transport module utilized in our study holds promise for application in circumplanetary disk models, thereby shedding light on the role of COMs in moon formation.

Fig. 2 Two-dimensional temperature map of the disk after 300 kyr of evolution, showing the methanol iceline at t = 0 (dashed blue line) and at t = 300 kyr (solid blue line). The white curve represents the average trajectory of 500 simulated particles. To provide a clearer illustration of individual particle behavior, the insets on the right display the trajectory of a single particle within the disk. The panels illustrate the trajectories of particles with sizes of 1 cm, 100 µm, and 1 µm, released at disk temperatures of 20 K (panels a, b, and c) and 80 K (panels d, e, and f). Initial and final positions are marked by a green triangle and a green cross, respectively.

Fig. 3 Averaged accumulated irradiation, temperature, and pressure for particles sized 1 cm, 100 µm, and 1 µm, released at disk temperatures of 20 K (panels a, b, and c) and 80 K (panels d, e, and f). Left: averaged accumulated irradiation of the particles over time (blue line) compared to the maximum irradiation accumulated by methanol ice during the experiments (orange line). Right: averaged temperature (red line) and averaged pressure (blue line) experienced by the particles over time.

Acknowledgements

O. Mousis acknowledges support from CNES and the Programme National de Planétologie (PNP) of CNRS/INSU. The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille Université-A*Midex, a French “Investissements d’Avenir program” AMX-21-IET-018. This research holds as part of the project FACOM (ANR-22-CE49-0005-01_ACT) and has benefited from a funding provided by l’Agence Nationale de la Recherche (ANR) under the Generic Call for Proposals 2022. This material is based upon work supported by NASA’S Interdisciplinary Consortia for Astrobiology Research (NNH19ZDA001N-ICAR) under award number 19-ICAR19_2-0041. AB acknowledges support from CNES. This work was supported by the Programme National de Planétologie (PNP) of CNRS-INSU cofunded by CNES. A.B. acknowledges support from the French government under the France 2030 investment plan, as part of the Initiative d’Excellence d’Aix-Marseille Université–A*MIDEX AMX-21-PEP-032. We would like to thank the anonymous referees for their helpful comments, which significantly improved the quality of this paper.

Appendix A Opacity values

The mean Rosseland opacity, denoted as κ_R, is a frequency-averaged quantity that depends on the gas density ρ_m and the temperature T_m in the midplane. It is expressed as ${\kappa }_{\text{R}}={\kappa }_0{\rho }_{\text{m}}^a{T}_m^b$, where κ₀, a, and b are parameters determined through observations by Bell & Lin (1994). Detailed values for these parameters are provided in Table A.1.

Appendix B Stokes number

In the PPD, particles can either be coupled to the gas or experience radial drift due to decoupling from the gas. The Stokes number governs the radial velocity component through Equation 21. A lower Stokes number indicates stronger coupling to the gas, while the influence of radial drift peaks when the Stokes number is close to 1 (Johansen et al. 2014; Drążkowska et al. 2023). The Stokes number is defined as St=t_sΩ_K (with $t_s=\frac{{\rho }_s{R}_s}{{\rho }_g\sqrt{8/\pi }{c}_s},c_s=\sqrt{R_g{T}_m/\mu }$ and ${\text{Ω}}_K=\sqrt{G{M}_{\star }/r^3}$). Thus, it depends on the disk density, temperature, and the particles’ position and size. Specifically, larger particles have a higher Stokes number (as shown in Fig. B.1), and denser disks result in a lower Stokes number (with the density profile displayed in Fig. B.2). Moreover, the stopping time varies depending on whether the particles are in the Stokes or Epstein regime (Johansen et al. 2014; Drążkowska et al. 2023). Particles with a size smaller than 9λ/4 are in the Epstein regime, as described in Sect. 3.3. In our nominal cases, we are always in the Epstein regime (cf. Fig. B.3, which shows the 9λ/4 map of the disk).

Figure B.4 shows the Stokes number for 1 cm sized particles in the PPD, along with the median trajectory of these particles. Figure B.1 presents the Stokes number for particles of various sizes, calculated using the mean density and mean temperature encountered by the particles and their mean trajectories, in each of our nominal cases. It is evident that particle size is the primary factor influencing the Stokes number. Specifically, 1 cm sized particles have a Stokes number ranging from 1 × 10⁻³ to 3 × 10⁻³, 100 µm sized particles range from 1 ×10⁻⁵ to 3 ×10⁻⁵, and 1 µm sized particles range from 1 ×10⁻⁷ to 4 ×10⁻⁷.

Fig. B.1 Stokes number computed for 500 simulated particles in each size category (1 cm, 100 µm, and 1 µm) and each based on the average temperature and density of the disk they encounter and their mean trajectories. Particles were released at temperatures of 20 K and 80 K.

Fig. B.2 Two-dimensional density profile of the disk after 100 and 300 kyr of evolution.

Fig. B.3 Mean free path of the gas, λ, multiplied by 2.25 after 100 and 300 kyr of PSN evolution. Particles smaller than 9λ/4 are in the Epstein regime, while larger particles fall into the Stokes regime. Blue and white lines represent the median trajectories of 500 simulated 1 cm particles released at 80 K and 20 K, respectively. The gray hatching highlights the region where 1 cm and smaller particles transition from the Epstein regime to the Stokes regime.

Fig. B.4 Stokes number in the disk after 100 and 300 kyr of PSN evolution, calculated for 1 cm sized particles. The Stokes number is lower for smaller particles. As in Fig. B.3, the median trajectories of 500 simulated 1 cm particles, released at 80 K and 20 K, are represented by blue and white lines, respectively.

References

All Tables

All Figures

	Fig. 1 Two-dimensional map of the accumulated irradiation after 10, 100, and 1000 kyr of the PSN’s evolution. The black area around the midplane corresponds to the region where the accumulated irradiation is lower than the most important dose obtained during the experiments of interest (8.64 × 1017 photons cm−2).
In the text

Fig. 2 Two-dimensional temperature map of the disk after 300 kyr of evolution, showing the methanol iceline at t = 0 (dashed blue line) and at t = 300 kyr (solid blue line). The white curve represents the average trajectory of 500 simulated particles. To provide a clearer illustration of individual particle behavior, the insets on the right display the trajectory of a single particle within the disk. The panels illustrate the trajectories of particles with sizes of 1 cm, 100 µm, and 1 µm, released at disk temperatures of 20 K (panels a, b, and c) and 80 K (panels d, e, and f). Initial and final positions are marked by a green triangle and a green cross, respectively.

In the text

Fig. 3 Averaged accumulated irradiation, temperature, and pressure for particles sized 1 cm, 100 µm, and 1 µm, released at disk temperatures of 20 K (panels a, b, and c) and 80 K (panels d, e, and f). Left: averaged accumulated irradiation of the particles over time (blue line) compared to the maximum irradiation accumulated by methanol ice during the experiments (orange line). Right: averaged temperature (red line) and averaged pressure (blue line) experienced by the particles over time.

In the text

	Fig. B.1 Stokes number computed for 500 simulated particles in each size category (1 cm, 100 µm, and 1 µm) and each based on the average temperature and density of the disk they encounter and their mean trajectories. Particles were released at temperatures of 20 K and 80 K.
In the text

	Fig. B.2 Two-dimensional density profile of the disk after 100 and 300 kyr of evolution.
In the text

Fig. B.3 Mean free path of the gas, λ, multiplied by 2.25 after 100 and 300 kyr of PSN evolution. Particles smaller than 9λ/4 are in the Epstein regime, while larger particles fall into the Stokes regime. Blue and white lines represent the median trajectories of 500 simulated 1 cm particles released at 80 K and 20 K, respectively. The gray hatching highlights the region where 1 cm and smaller particles transition from the Epstein regime to the Stokes regime.

In the text

	Fig. B.4 Stokes number in the disk after 100 and 300 kyr of PSN evolution, calculated for 1 cm sized particles. The Stokes number is lower for smaller particles. As in Fig. B.3, the median trajectories of 500 simulated 1 cm particles, released at 80 K and 20 K, are represented by blue and white lines, respectively.
In the text