© The Authors 2023

1 Introduction

The signatures of polycyclic aromatic hydrocarbons (PAHs) have been identified in space through the detection of prominent mid-infrared (MIR) emission features located at 3.3, 6.2, 7.7, 8.6, and 11.2 µm (Allamandola et al. 1985, 1989). In protoplanetary discs, given their small size, PAHs are well coupled to the gas and are therefore present in the disc’s photosphere. There, PAHs efficiently absorb UV radiation and convert the energy to MIR radiation, which makes them an excellent tracer of the photosphere of protoplanetary discs. A famous example is the protoplanetary disc HD 97048, whose flaring has been shown by VLT Imager and Spectrometer for mid-InfraRed (VISIR) imaging of the PAH emission at 8.6 µm (Lagage et al. 2006; Doucet et al. 2007).

Although many protoplanetary discs show PAH features, there are several discs in which the IR PAH signatures are absent (Meeus et al. 2001; van Boekel et al. 2004; Acke & van den Ancker 2004). A study by Acke et al. (2010) shows that PAHs are detected in about 60–70% of Herbig Ae/Be stars. Valegård et al. (2021) show that 27% of their sample discs around intermediate mass TTauri stars (IMMT, defined as 1.5 M_⊙ ≤ M_* ≤ 5 M_⊙ and with a spectral type between F and K3) contain PAHs, where they adopted the criterion that at least two distinct PAH features mustbe present foracleardetection. If only one feature is considered, the fraction of IMMT discs with PAHs increases to 44%. For classical TTauri stars, Geers et al. (2006) find a lower detection limit of 8%, excluding 14 tentative detections of PAHs in their total sample of 38 T Tauri discs.

Despite a very general correlation between PAH detection and stellar temperature due to the emitted UV field, it is not clear why the PAH emission strength also varies strongly within a stellar type. To explain this phenomenon, it is necessary to look at protoplanetary discs individually and relate their morphology, environment, and properties to the PAH abundance and emission strength. In the past several authors proposed mechanisms for the variation of PAH strength in discs. For the lack of PAH emission in T Tauri stars Siebenmorgen & Krügel (2010) and Siebenmorgen & Heymann (2012) modelled the destruction of PAHs through absorption of extreme ultraviolet (EUV) and X-ray photons. The authors show that for typical hard photon luminosities, PAHs can be sufficiently destroyed, and hence can explain the lack of PAH emission in these disc. Other authors such as Geers et al. (2009) provide different possible explanations. Gas-phase PAHs can interact with icy dust grains and freeze out on the grain surface. When these dust grains are then exposed to UV radiation, the PAHs can be ionised, triggering ion-mediated reactions on the grain surface that deplete the PAHs (e.g. Ehrenfreund et al. 2006; Bouwman et al. 2010). Alternatively, in dense regions of the discs PAH molecules can aggregate to PAH clusters and form stack-like structures, as described by Rapacioli et al. (2005, 2006). Due to their larger heat capacity, PAH clusters have a lower excitation temperature after absorption of a UV photon and consequently the strength of short-wavelength emission bands is reduced (Bakes et al. 2001). Furthermore, the emission profile of PAH clusters is altered compared to the emission of individual PAH molecules since additional intermolecular modes occur in the MIR (Rapacioli & Spiegelman 2009) that may be responsible for the broad underlying emission bands observed in the interstellar medium (ISM; Allamandola et al. 1989).

In our previous study (Lange et al. 2021) we investigated the stability of PAH clusters that are fully exposed to the stellar UV radiation. We found that in the inner regions of the protoplanetary disc, PAH clusters dissociate within a very short time regardless of their size. In the outer regions of the disc, however, clusters can no longer be photodissociated once they exceed a critical size as their peak temperature is too low to drive for dissociation. As a result, we expect that PAH clustering in protoplanetary discs is a realistic scenario to surpress the emission of the short-wavelength PAH features, provided that the critical cluster size can be reached so that clusters can be stable in the disc photosphere.

In this study we want to investigate the freeze-out of PAHs on dust grains and PAH clustering simultaneously as they are competing processes that depend on the local gas-phase abundance of PAHs and grain density. We consider both processes as key processes to explain the absence of detectable PAH emission in some protoplanetary discs and focus our model on a Herbig Ae/Be star disc. We want to quantify how the presence of dust grains in the protoplanetary disc affects the evolution and growth of PAH clusters and explore the possibility of using PAHs as an independent tool to obtain information on the dust grain population supplementing methods based on submillimetre emission (e.g. with ALMA; Birnstiel et al. 2018) and scattered light observations (e.g. Tazaki et al. 2019). In our models the gas-phase abundance of PAHs after its evolution in the disc is proportional to the total dust grain surface area. Therefore, if the PAH abundance of a protoplanetary disc is inferred from near-infrared (NIR) and MIR observations of the PAH features, our model is able to check whether the dust population is consistent with ALMA/SPHERE characterisations of the dust grains. Additionally, we model the desorption of PAH monomers and PAH clusters from the surface of dust grains to investigate their ability to enter the gas-phase again. As the James Webb Space telescope (JWST) was launched in December 2021, the NIR and MIR wavelength range will be more accessible than ever through the Mid-Infrared Instrument (MIRI) and the Near Infrared Camera (NIRcam) with high sensitivity and significant spatial resolution. Hence, modelling work is crucial in order to understand the relevant processes and to analyse and interpret the upcoming observations.

This work is organised in the following way. In Sect. 2 we describe the coagulation and desorption models of the PAH clusters that were used. In Sect. 3, we present the results of our coagulation, desorption, and turbulent mixing model and investigate their implications on the PAH abundance in certain regions of the protoplanetary disc. Then in Sect. 4 we discuss the consistency of our results to previous studies, discuss our results in the context of sputtering, and heterogeneous PAH clusters.

2 Methods

First, we describe the relevant processes that influence the PAH abundance in our Herbig Ae/Be protoplanetary disc and describe where they take place in our model. Here we neglect possible effects by hard photons (EUV, X-ray) as they are not dominant in Herbig star discs. For this purpose, we divide the protoplanetary disc into layers depending on which process is important (Fig. 1). The upper layer of the disc is the photosphere where PAHs and PAH clusters are exposed to the full stellar radiation field and are able to desorb from dust grains. There, clusters are rapidly dissociated into their constituent monomers. Below the photosphere, the UV field is partially shielded so that most clusters cannot be dissociated. However, it is strong enough so that dimers can be dissociated more quickly than they can grow from monomers inhibiting the growth of new clusters. Hence, we call it the quiet layer as PAHs keep their current state. Near the midplane, the coagulation layer is the region where PAHs rapidly coagulate into clusters and adsorb on dust grains. We note that the coagulation layer includes not only the classical cold midplane, but also the warmer layers above that are shielded sufficiently from UV radiation.

The coagulation layer is relatively dense and PAH clustering is very efficient. Protected from stellar UV radiation, photodissociation of PAH clusters can be neglected there, and hence the coagulation layer provides an ideal place for the growth of PAH clusters. Because dust particles are also present, PAH clusters will collide with them and freeze out to the dust particles. As the adsorbed PAH is then coupled to the thermal reservoir of the dust particle, the clusters quickly adapt to the temperature of the grains and thermal evaporation of the clusters becomes impossible in most regions of the disc. Therefore, we expect that in the coagulation layer PAHs are found exclusively as clusters on dust grains and that no PAHs are present in the gas phase.

In contrast, the photosphere is dominated by UV radiation from the central star. Based on our previous study (Lange et al. 2021), the growth of PAHs to clusters is negligible as photodissociation is much faster than the first step of clustering, the formation of dimers through collisions of monomers. Already existing clusters, however, can be dissociated and individual PAHs adsorbed on dust grains can be released from the grain surface when a UV photon is absorbed. Further, adsorbed PAHs receive sufficient UV radiation such that desorption from the grain is possible. As the photodissociation and photodesorption processes in the photosphere are strong, we do not expect cluster formation and freeze-out to be relevant there.

For each disc region we developed a specific model taking into account the processes relevant there. To understand the clustering and adsorption of PAHs on dust grains in the coagulation layer, we developed a coagulation code, which we describe in Sect. 2.1. We then present our model for the desorption of clusters and monomers in the photosphere (Sect. 2.2), which we later apply to the results of the coagulation code to determine what fraction of the original PAHs can be recovered from the dust grains when they are mixed from the coagulation layer into the photosphere.

Fig. 1 Sketch of PAH processing cycle in the protoplanetary disc driven by vertical mixing. PAH monomers will eventually mix into the UV depleted coagulation layer where PAH collisions lead to sticking. The cluster growth continues until the PAH clusters freeze out on dust grains. If the cluster carrying grains are mixed into the UV-rich photosphere, PAH monomers and small clusters can desorb while large clusters remain on the grain. Every cycle of coagulation, adsorption, desorption, and dissociation leads to a depletion of gas-phase PAHs as with every cycle large clusters form and remain on the grains.

2.1 PAH coagulation model

The growth of PAHs into clusters is a mechanism that reduces the feature-to-continuum ratio of PAH bands (Geers et al. 2009) because temperature fluctuations are reduced compared to individual monomers. Furthermore, the additional vibrational modes of PAH clusters (Rapacioli et al. 2007) are responsible for underlying broad bands below the PAH features (Allamandola et al. 1989) which further reduce the observable feature-to-continuum ratio. So far, we have only investigated photodissociation of clusters (Lange et al. 2021) and estimated the cluster growth timescale based on isolated gas-phase PAH collision rates. In this work we want to determine how fast PAH clusters can grow in the coagulation layer given that dust grains are also able to interact with PAHs through collisions. On the one hand, the deposition of PAHs onto dust grains is another mechanism for reducing the strength of the MIR PAH bands, as the adsorbed PAHs are thermally coupled to the dust grains, and therefore have lower peak temperatures. On the other hand, the presence of dust grains allows PAHs to stick to the surface of the grains so that these PAHs are no longer available for PAH cluster formation in the gas-phase. Since the coagulation layer is shielded from UV radiation, we neglect photodissociation of forming PAH clusters in this model. Furthermore, the expected temperature in the coagulation layer is lower than the required 400 K (for the smallest simulated PAH coronene) above which thermal desorption and thermal dissociation become relevant, and therefore we also neglect these processes. In our coagulation model we consider a local box filled with PAHs, dust, and gas without transport of PAHs and dust in the radial or vertical direction because the coagulation is much faster than the transport processes, as we show below.

We describe the size of PAH clusters by the number of carbon atoms N_C (not monomer molecules) it contains. Then the cluster growth can be described by the classical Smoluchowski equation (Smoluchowski 1916) $$\begin{array}{l}\frac{df\left(N_{\text{C}}\right)}{\text{d}t}={\displaystyle {\int }_0^{N_{\text{C}}/2}f\left({N^{\prime }}_{\text{C}}\right)f\left(N_{\text{C}}-{N^{\prime }}_{\text{C}}\right)k\left({N^{\prime }}_{\text{C}},N_{\text{C}}-{N^{\prime }}_{\text{C}}\right)\text{d}{N^{\prime }}_{\text{C}}}\hfill \\ -{\displaystyle {\int }_0^{\infty }f\left(N_{\text{C}}\right)f\left({N^{\prime }}_{\text{C}}\right)k\left(N_{\text{C}},{N^{\prime }}_{\text{C}}\right)\text{d}}{N^{\prime }}_{\text{C}}\hfill \\ -Q\_\left(N_{\text{C}}\right),\hfill \end{array}$$(1)

where f (N_C, r, θ) is the number density of a given PAH cluster size at a given radial position r and height θ; k describes the coagulation kernel; and Q_ is a sink term to account for the adsorption of PAHs from the gas phase to grains. We do not consider the case that PAHs can grow into clusters on grains here. The first term in Eq. (1) describes the formation of a cluster with N_C carbon atoms from the collision of two smaller clusters; the second term describes the loss of clusters with size NC through collisions with other clusters. At van der Waals binding energies of 1 eV between PAH monomers and the cluster, very high temperatures of more than 13 500 K are needed for destructive monomer-cluster collision outcomes (Rapacioli et al. 2006). With the relevant thermal energies in the protoplanetary disc, we can safely assume perfect sticking for PAHs monomer-monomer, monomer-cluster, cluster-cluster, and all PAH-grain collisions.

For PAHs in the gas phase, we only consider Brownian motion and neglect turbulent processes for the coagulation model as sources for relative velocities as these do not contribute, due to the strong coupling of PAHs to the gas. The coagulation kernel can then be determined by the relative velocity ∆υ_i,j and the collision cross-section σ_i,j, which both depend on the size of the collision partners i and j involved: $$k={\sigma }_{i,j}\text{Δ}{\upsilon }_{i,j}.$$(2)

The relative velocities in Brownian motion are given by $$\text{Δ}{\upsilon }_{i,j}=\sqrt{\frac{8k_{\text{B}}T}{\pi {\mu }_{i,j}}},$$(3)

with µ_i,j being the reduced mass of the collision partners, k_B the Boltzmann constant, and T the kinetic temperature of the gas that we set to T = 100 K. The collision cross-section σ_i,j between PAH clusters is assumed to be the maximum geometrical cross-section with radius r $${\sigma }_{i,j}=\pi {\left(r_i+r_j\right)}^2.$$(4)

PAH cluster sizes and adsorption

Based on the studies of Rapacioli et al. (2005), small PAH clusters tend to prefer stack-like structures, while larger clusters favour extended 3D structures made from smaller stacks. However, since the temperature of a cluster is never zero, other configurations are energetically possible, and even small clusters can form modified stack structures or 3D structures (Rapacioli et al. 2007). Consideration of all the different allowed cluster structures would unnecessarily complicates the model. Hence, we decided to simplify the model and divide the effective radius into three size ranges: monomers, which are a planar molecules; small clusters, which grow linearly to account for the formation of stacks until the diameter of a monomer is equal to the length of a stack; and large clusters, which grow spherically as compact dust grains. The detailed formulas for the cluster radius r_i over the entire size range can be found in Appendix A.

As we did for the growth of PAH clusters, we model adsorption on dust grains with Eqs. (2) and (3), where the collision partners here are a PAH cluster i and a dust grain l with grain size a_l. We take this into account in a semi-analytical approach by summing over all grain sizes l, $$Q\_\left(N_{\text{C}}\right)=f\left(N_{\text{C}}\right){\displaystyle {\int }_{a_{\min }}^{a_{\max }}kn\left(a\right)\text{\hspace{0.17em}d}a,}$$(5)

where n_l is the number density of the dust grains. Consequently, the adsorption of PAHs depends on both the local dust grain abundance and the grain size distribution. Therefore, to characterise the dust population with a single parameter, we use the total collisional cross-section per unit volume of the dust as the relevant quantity defined by $${\sigma }_{\text{dust}}{\displaystyle {\int }_{a_{\min }}^{a_{\max }}\pi n\left(a\right)a^2\text{d}a.}$$(6)

We note that this approach is limited to the case where the PAH clusters are much smaller than the dust grains (r_i· ≪ a), which is unrestrictedly valid in our set-up.

To solve the coagulation Eq. (1), we follow the principles described in the appendix of Dullemond & Dominik (2005). To ensure mass conservation, we use their conservation method, but apply it only to the two integral terms where numerical mass loss may inadvertently occur. We add for the time step determination a safety factor of 0.3 and stop the simulation as soon as less than 10⁻⁴% of the PAH mass is left in the gas phase. We start our simulation with PAH monomers in the gas phase only, where N_C,0 describes the number of carbon atoms in a monomer. To ensure numerical stability, we use a Gaussian of width σ = 5C atoms centred around N_C,0 as our initial distribution.

2.2 Photodesorption of PAH clusters from dust grains

After presenting the coagulation model, we look at PAHs frozen onto dust grains and how they can be brought back into the gas phase. Given that cluster formation and freeze-out require a strong UV shielding in the coagulation layer to be faster than their counter processes, desorption must occur higher up in the disc. Additionally, PAH charging through photo-ionisation and pick-up of free-electrons (i.e. from C->C+) lead to a difficult PAH charge distribution (Thi et al. 2019) that can prevent PAH sticking higher up in the disc. Assuming a turbulent disc, the PAH-bearing dust grains need to be transported vertically through the disc by turbulence until they eventually reach the UV-rich photosphere to desorb from the grains. When the equilibrium temperature of the grains is above the evaporation temperature of the PAHs, PAHs can thermally evaporate. In addition, we consider that adsorbed PAH clusters absorb UV photons that increase the internal excitation of the PAH vibrational modes, including those that bind the PAHs to the grain. This internal excitation will be lost through coupling to the grain phonon modes. Therefore, the competition between sublimation and energy transfer to the grain must be studied to determine the sublimation rate. Due to the high heat capacity of the dust grains, the temperature fluctuations after absorption of a photon by the dust grains are too small to sublimate PAHs. Only direct photon absorptions by the PAHs can lead to sublimation.

Like the dissociation of PAH clusters by UV radiation, the sublimation of PAHs can be simulated by a Monte Carlo model where we consider photon absorption, cooling, and desorption of the entire cluster as possible events.

2.2.1 Photon absorption

We model the absorbed energy flux of a PAH Φ_a,E at a given photon energy E by $${\text{Φ}}_{\text{a,E}}=F_{\text{E}}{\sigma }_{\text{E}},$$(7)

where σ_E is the photon absorption cross-section of a PAH at photon energy E, and F_E is the stellar flux density. We estimate the stellar flux in the photosphere with a black-body spectrum with stellar temperature and luminosity, as given in Table 1. To determine the flux deeper in the disc, we moderate the flux with an optical depth through the radiative transfer code RADMC-3D¹ and the opacity calculator OpTool (Dominik et al. 2021) using the standard DIANA opacities (Woitke et al. 2016). For σ_E we use photon absorption cross-sections from time dependent density functional theory (DFT) calculated by Malloci et al. (2007) for individual PAH species ranging from coronene (C₂₄H₁₂) to circumovalene (C₆₆H₂₀). We divide the energy space between 0.1 and 13.6 eV into 100 bins with average energy E_m so that we can derive an absorption rate of $$r_{\text{a,m}}=\frac{1}{E_m}{\displaystyle {\int }_{E_{\text{i,min}}}^{E_{\text{i,max}}}{\text{Φ}}_{\text{a,E}}\text{d}E.}$$(8)

We do not consider photon energies above 13.6 eV. Our assumed black-body spectrum for a typical Herbig star does not include a significant amount of EUV or X-ray photons with energies of a few 100 eV or more. Therefore, the ionisation of hydrogen increases the optical depth around 13.6 eV significantly so that the relevant photons with slightly higher energies cannot penetrate deeper into the disc (Ryter 1996), and therefore do not contribute to the evolution of PAHs.

2.2.2 Cooling through thermal coupling to grains

Cooling of PAH clusters can happen by IR photon emission and by energy transfer to the dust grain. However, since cooling by IR photons happens on a timescale of seconds (Bakes et al. 2001) and is thus much slower than energy transfer to the dust grain, we can neglect IR cooling. We further assume that the temperature of the dust grains is constant at the radiative equilibrium value because the grains are large enough to neglect temperature fluctuations (Tielens 2005). Since we focus on the optically thin photosphere, we use Eq. (23) to estimate the equilibrium temperature of the grains.

For the van der Waals bond between the grain and the PAHs (resp. PAH clusters), we assume that there is a flat contact surface between the grain-facing PAH and the dust grain so that the potential energy is minimised. Then we can assume that the energy transfer can be modelled as heat transfer between two graphene sheets through the c-plane. We expect the PAH and grain temperature to be in the range of 100–1000 K, and hence we assume a thermal conductivity of κ = 0.04 Wm⁻¹ K⁻¹ = 4000 erg s⁻¹ cm⁻¹ K⁻¹ based on the theoretical calculations for graphene by Alofi & Srivastava (2014). The effective contact surface is determined by the size of the PAH monomer species through $$A=\pi {\left(\sqrt{N_{\text{C,0}}}\cdot 0.9\text{Å}\right)}^2=N_{\text{C,0}}\cdot 2.5\times {10}^{-16}{\text{cm}}^2,$$(9)

where N_C,0 is the number of C atoms of a PAH monomer (Tielens 2008). The energy transfer rate from the PAH cluster to the grain is then determined as $$\frac{\text{d}Q}{\text{d}t}=\frac{\kappa A\text{Δ}T}{\text{Δ}d},$$(10)

where ∆d is the bond length between the surface molecule and the grain surface. We assume the typical alternating graphite interlayer distance of ∆d = 3.34 Å (de Andres et al. 2008) for the bond length, which is comparable to the van der Waals bond distance of two plane-parallel PAH molecules of ≈3.5 Å (Rapacioli et al. 2005). The energy transport from the cluster to the dust grain is not a classical macroscopic heat transport.

Instead, energy is transferred quantum-wise with the natural frequency of the van der Waals mode of the bond that can be estimated by $${\upsilon }_z=\sqrt{\frac{2N_{\text{s}}{E}_{\text{A}}}{{\pi }^{2}m}},$$(11)

where N_s ≈ 2 × 10¹⁵ sites cm⁻² is the number of binding sites of the adsorbent, E_A is the binding energy, and m the mass of the adsorbate (Tielens 2005). For a coronene monomer with a binding energy of E_A = 1.4 eV, the corresponding frequency is υ_z = 1.4 × 10¹² s⁻¹. To match the vibrational frequency, we choose the energy of a heat quantum to be 1% of the remaining excitation energy of the cluster (∆E = 0.01(E − E_grain)) so that the cooling rate $$r_{\text{c}}=\frac{\text{d}Q}{\text{d}t}=\text{Δ}{E}^{-1}\approx {10}^{12}{\text{s}}^{-1}$$(12)

agrees with the estimated van der Waals frequency. We stop the cooling process as soon as the temperature of the cluster deviates less than 0.1 K from that of the grain and use the equilibrium temperature of the dust grain until the next absorption of a photon to improve the efficiency of the code.

2.2.3 Evaporation rate

The evaporation of adsorbates can be generally described by the Arrhenius equation, which can be derived from the Rice-Ramsperger-Kassel-Marcus (RRKM) theory (Tielens 2005). The classical expression is $$r_{\text{k}}=k_0\cdot \text{exp}\left(\frac{-E_{\text{A}}}{k_{\text{B}}{T}_{\text{e}}}\right),$$(13)

where k₀ is the evaporation rate constant, E_A the binding energy of the PAH cluster to the grain, and T_e is the effective temperature of the PAH cluster determined by $$T_{\text{e}}=T_{\text{m}}\left(1-0.2\frac{E_{\text{A}}}{E}\right).$$(14)

We assume k₀ = 2.5 × 10¹⁷ s⁻¹ for all PAH clusters similar to the dissociation of PAH clusters derived in Lange et al. (2021). We use the heat capacity of Bakes et al. (2001), and thereby describe the temperature-energy relation as (Tielens 2021) $$T_{\text{m}}=\{\begin{array}{ll}3750{\left(\frac{E\left(\text{eV}\right)}{3N-6}\right)}^{0.45}\text{\hspace{0.17em}K}\hfill & \text{if\hspace{0.17em}}{T}_{\text{m}}<1000\text{\hspace{0.17em}K}\hfill \\ 11000{\left(\frac{E\left(\text{eV}\right)}{3N-6}\right)}^{0.8}\text{\hspace{0.17em}K}\hfill & \text{if\hspace{0.17em}}1000\text{\hspace{0.17em}K}<T_{\text{m}}\hfill \end{array},$$(15)

where N is the total number of atoms in the cluster. We estimate the activation energy E_A for the evaporation of PAH clusters from dust grains by using an approximation from DFT calculations for the adsorption of PAHs on a graphene surface by Li et al. (2018) $$\frac{E_{\text{A}}}{N_{\text{C}}}=0.046+0.021\frac{N_{\text{H,0}}}{N_{\text{C,0}}},$$(16)

with N_H,0 and N_C,0 being respectively the number of hydrogen and carbon atoms in a PAH monomer. This approximation matches the findings of the experimental study by Zacharia et al. (2004). The calculated binding energies can be found in Table 2. Our derived evaporation rates agree with Montillaud & Joblin (2014), who interpolated evaporation rates for clusters calibrated through molecular dynamics studies.

2.2.4 Monte Carlo scheme

Finally, the desorption rate is determined by evaluating the energy fluctuations of the PAH clusters after absorption of UV photons with the desorption rates. The energy probability distribution G(E) of the adsorbed PAH cluster is determined by a Monte Carlo simulation in which photon absorption, heat transport, and desorption events are considered. We do this using a Monte Carlo scheme that follows the fundamental principles explained in Zsom & Dullemond (2008). The total rate of events for a simulated PAH molecule is thus $$r={\displaystyle \sum _{\text{m}}{r}_{\text{a,m}}+r_{\text{c}}+r_{\text{k}}.}$$(17)

We note that although multi-photon absorption events are included in this equation, they are extremely rare because the energy transfer from the PAHs to the grain (τ ≈ 10⁻¹⁰ s) is very fast compared to the absorption rate (τ ≈ 1 s) of the UV photons. The energy probability distribution is therefore entirely dominated by single photon absorptions. One Monte Carlo time step δt is determined by a random draw from an exponential distribution with mean 1/r $$f\left(\delta t\right)=r\exp \left(-r\delta t\right).$$(18)

We determine which event has occurred by the set of all possible events with rates R = {r_a,m, r_c, r_k} with the probability of an event j with rate r_j ∊ R given by $$P\left(j\right)=\frac{r_j}{r}.$$(19)

The final effective desorption rate can be determined by integrating all possible energy states of G(E) with the individual desorption rates r_k: $$k={\displaystyle {\int }_0^{\infty }G\left(E\right)r_{\text{k}}\left(E\right)\text{d}E.}$$(20)

The disc model

The values of dust density and dust size distribution in the protoplanetary disc are crucial in our coagulation model, while temperature and radiation intensity are essential for the evaporation of PAHs from dust grain surfaces. Therefore, we need to set up a disc model to define all relevant quantities. Our disc has a gas surface density profile of $$\sum \left(r\right)=730\frac{\text{g}}{{\text{cm}}^2}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-1.5}\left(\frac{M_{*}}{M_{\odot }}\right)$$(21)

scaled to a minimum mass solar nebula (MMSN) power-law profile, as first described by Weidenschilling (1977). Our disc is in vertical hydrostatic equilibrium, $$\rho \left(z,r\right)={\rho }_{\text{mid}}\left(r\right)\exp \left(\frac{-z^2}{2H^2}\right),$$(22)

where the gas pressure scale height is determined by H = c_s/Ω_K. The sound speed c_s is calculated by $c_{\text{s}}=\sqrt{k_{\text{B}}T/\text{μ}{m}_{\text{p}}}$, where we assumed µ = 1.37 and an optically thin surface layer so that the temperature profile can be described through the approach of Hayashi (1981): $$T\left(r\right)=280\text{\hspace{0.17em}K}{\left(\frac{r}{\text{au}}\right)}^{-0.5}{\left(\frac{L_{*}}{L_{\odot }}\right)}^{0.25}.$$(23)

The midplane density is determined by ${\rho }_{\text{mid}}\left(r\right)=\sum \left(r\right)/\sqrt{2\pi }H\left(r\right)$.

For the dust grains, we assume a Mathis–Rumpl–Nordsieck (MRN) size distribution (Mathis et al. 1977). Then the number density of the dust grains of size a is given by $$n\left(a\right)\propto a^{-\gamma },$$(24)

with γ = 3.5. Consequently, the smallest dust grains provide the largest share of the available surface area, and are therefore the dominant Brownian-motion collision partners for PAH molecules. Hence, PAHs are mainly adsorbed on the smallest available dust grains. In our standard model these grains are also closely coupled to the gas. However, the large dust grains also contribute to the total dust surface area. We want to consider an evolved dust disc that is in equilibrium between settling and mixing. To do this we use the results of Fromang et al. (2007) who made magnetohydrodynamic (MHD) simulations and apply their fitted turbulent toy model as an estimate for our disc. In their model the dust grains do not exert any feedback on the gas, which is why we use the equations for the hydrostatic equilibrium for the gas density. The vertical density distribution of the dust grains can be described by $$\frac{\partial {\rho }_{\text{d}}}{\partial t}-\frac{\partial }{\partial z}\left(z{\text{Ω}}^2{\tau }_{\text{s}}{\rho }_{\text{d}}\right)=\frac{\partial }{\partial z}\left[D_{\rho }\frac{\partial }{\partial z}\left(\frac{{\rho }_{\text{d}}}{\rho }\right)\right],$$(25)

where the steady-state solution is given through $$\frac{\partial }{\partial z}\left(\ln \frac{{\rho }_{\text{d}}}{\rho }\right)=-\frac{{\text{Ω}}^2{\tau }_{\text{s}}}{D}z.$$(26)

The diffusion coefficient D of the dust grains is defined by the local velocity fluctuations δυ_z $$D=\delta {\upsilon }_{\text{z}}^2{\tau }_{\text{corr}},$$(27)

where the correlation time of the velocity fluctuations τ_corr is determined by $${\tau }_{\text{corr}}=0.15\frac{2\pi }{{\text{Ω}}_{\text{k}}}.$$(28)

From their hydrodynamical studies, Fromang et al. (2007) found that δυ_z can be best described by $$\delta {\upsilon }_{\text{z}}=\{\begin{array}{ll}\delta {\upsilon }_{\text{z,\hspace{0.17em}mid}}+\left[\delta {\upsilon }_{\text{z,up}}-\delta {\upsilon }_{\text{z,mid}}\right]{\left(\frac{\left|z\right|}{2H}\right)}^2\hfill & \text{if\hspace{0.17em}}\left|z\right|<2H\hfill \\ \delta {\upsilon }_{\text{z,\hspace{0.05em}up}}\hfill & \text{else}\hfill \end{array},$$(29)

with δυ_z,up = 0.15c_s and δυ_z,mid = 0.05c_s. Finally, the stopping time of the particles τ_s is defined as $${\tau }_{\text{s}}=\frac{{\rho }_{s}a}{\rho c_{\text{s}}}.$$(30)

We numerically solve Eq. (26) to obtain the local dust grain density ρ_d(a, r, θ). To ensure numerical stability we solve for ln(ρ_d/ρ) rather than ρ_d. The parameters used for the disc set-up are provided in Table 1. Additionally, a test case set-up to compare to Fromang et al. (2007) is provided in Appendix B.

Fig. 2 Ratio of total PAH collisional cross-section σPAH to the total dust collisional cross-section σdust for the standard Herbig disc model STD with ISM PAH abundance. As the surface area of the dust is dominated by the smallest grains, which are closely coupled to the gas, the ratio is almost constant σPAH/σdust ≈ 310 resp. log(σPAH/σdust) ≈ 2.5 at relevant coagulation altitudes (z < 3H). Only in the outer disc does this ratio change; however, the relevant coagulation height decreases as well (see Fig. 8). The gas pressure scale heights of z =1 H, 3 H, 5 H are indicated by dashed lines.

2.4 Vertical mixing

In order to estimate how likely the evaporation of PAH clusters is in a turbulent disc, we want to compare the evaporation rate with the mean-residence time t_mrt of an adsorbed PAH cluster in the photosphere given a typical α-viscosity disc (Shakura & Sunyaev 1973). For this purpose we use the standard turbulence model typically used in the models of Dubrulle et al. (1995), Dullemond & Dominik (2004), Schräpler & Henning (2004), and Fromang et al. (2007), among others. In this model, the turbulent velocities of the largest vortices are proportional to the speed of sound c_s. Then the turbulent turnover time can be estimated by $$t_{\text{ed}}=\frac{{\alpha }^{1-2q}}{{\text{Ω}}_{\text{k}}},$$(31)

with the preferred scaling of q = 0.5 so that the length scale of an eddy is given by $$l_{\text{ed}}={\alpha }^{0.5}H.$$(32)

We estimate that the probability of a grain coupling to a neighbouring eddy is given by the mass fluxes through the foot j_z and top $j_{z+l_{\text{ed}}}$ of the eddy via $$j\left(z\right)=\upsilon \rho \left(z\right)={\alpha }^{0.5}{c}_{\text{s}}\rho \left(z\right)$$(33)

evaluated at the foot and top of an eddy. The probabilities can then be evaluated through mass conservation respectively by $$p_{\text{up}}\left(z\right)=\frac{j\left(z+l_{\text{ed}}\right)}{j\left(z\right)+j\left(z+l_{\text{ed}}\right)}$$(34)

and $$p_{\text{down}}\left(z\right)=\frac{j\left(z\right)}{j\left(z\right)+j\left(z+l_{\text{ed}}\right)}.$$(35)

Using the α description model for the mass flux, the mixing probabilities become vertically dependent on the densities at the given height so that Eqs. (34) and (35) respectively become $$p_{\text{up}}\left(z\right)=\frac{\rho \left(z+l_{\text{ed}}\right)}{\rho \left(z\right)+\rho \left(z+l_{\text{ed}}\right)}$$(36)

and $$p_{\text{down}}\left(z\right)=\frac{\rho \left(z\right)}{\rho \left(z\right)+\rho \left(z+l_{\text{ed}}\right)},$$(37)

which under hydrostatic equilibrium reduce to a simple exponential form. By tracing the particle movement through the eddies, we can determine how long carried particles are fully exposed to the incident UV radiation in the photosphere and how often they mix between the photosphere and the disc midplane. We follow a particle through 10⁸ randomly sampled steps in the disc with an upper boundary of 4H_p where particles are forced to mix down.

3 Results

3.1 Coagulation and freeze-out of PAH clusters

In the following section, we present the results of the PAH coagulation model with and without the presence of dust grains to show the effect on the PAH population. We apply our Monte Carlo method to the standard disc model STD from Table 1 and simulate the coagulation of coronene (C₂₄H₁₂) with an initial ISM abundance (ratio of hydrogen atoms to C atoms locked up in PAHs; C_PAH:H = 1.5 × 10⁻⁵; Tielens 2008). We use r = 10 au and θ = 0 as representative of the entire coagulation layer since the ratio of PAH to dust collision cross-section is initially nearly constant (Fig. 2). Therefore, our coagulation simulations are universally applicable to the coagulation layer if the location is optically thick to the incident UV radiation to prevent dissociation of the clusters (τ_UV ≫ 1). Additionally, the PAHs need to be uncharged as otherwise the electric repulsion force inhibits the sticking of the PAHs. We describe our PAH population by the number of cluster members (N monomers = N_c/N_C,0) and give the corresponding PAH size distribution as $N_{\text{C}}^{2}f$, so that the surface area under the curve is equivalent to the total number of C atoms in this size range (resp. proportional to the PAH mass).

Figure 3 shows the coagulation of PAH clusters without the presence of any dust grains. Within the first few days after the start of the simulation, dimers, trimers, and quadrumers form whose respective peak width is given by the propagation of the initial conditions. The further growth of the PAH clusters is rapid. Within a year, clusters can form whose size corresponds to the size of very small carbonaceous grains (VSGs; a ≈ 10 Å resp. N_C ≈ 10⁴ C atoms), as defined in the models by Desert et al. (1990) and Draine & Li (2001), among others. This growth continues until all PAHs are eventually locked up in very large PAH clusters (with N ≥ 10⁵ monomers). Depending on the underlying PAH species, these very fast-forming clusters cannot be photodissociated at larger distances from Herbig stars (e.g. r ≥ 40 au for circumcoronene) when they are transported into the photosphere as they exceed the critical cluster size (Lange et al. 2021). Consequently, if the conditions are met for this scenario (lack of UV radiation and dust grains), we expect the formation of very large clusters that can hardly be broken up if transported to UV-rich environments, even for Herbig star discs.

However, a different picture emerges when adsorption on dust grains is considered with the same initial conditions as before (Fig. 4). The first phases of cluster growth are not disturbed by the dust grains as cluster growth is initially many times faster than adsorption. However, this changes when the PAHs have grown to larger cluster sizes as the absolute number of clusters has decreased so much that the probability of PAH-PAH collisions becomes similar to PAH-dust particle collisions. Consequently, clusters themselves start to stick on dust grains, which further reduces the number of possible PAH-PAH collision partners. From this point on, further cluster growth is mostly inhibited and thus the maximum cluster size is limited by the total available dust collision cross-section. Based on the adsorbed distribution, most of the mass of PAHs is adsorbed as large clusters. However, the most frequent adsorbed PAH size is still the monomer ƒ (Fig. 5), and the occurrence decreases with size. The adsorbed size distribution is identical for all dust particle sizes a since the size of the PAHs is negligible compared to the dust particles. Moreover, the number of PAHs on a specific grain size is proportional to the total surface area σ_dust(a) of the dust particle size range (see Eq. (5)). Generally, PAHs are more likely to freeze out on the smallest dust grains that are locally available.

Furthermore, we want to investigate how clustering is affected when the initial PAH abundance differs from the ISM abundance as we do not know the initial PAH content in the molecular cloud. For comparison with the other model, we keep the dust grain population the same. To do this, we vary the PAH abundance in the range 0.01–10 times the ISM abundance (C_PAH:H = 1.5 × 10⁻⁵; Tielens 2008) and compare the adsorbed PAH cluster distribution (Fig. 6). If the PAH abundance is very low (≤0.01 ISM abundance), the amount of adsorbed PAHs decreases with decreasing PAH abundance. Due to the low number of PAH collision partners, the PAH clusters cannot grow very quickly and mainly small PAH clusters are formed that adsorb on the dust grains. If the PAH abundance is sufficiently high (≥0.1 ISM abundance), the absolute number of adsorbed PAHs is limited for the smaller cluster sizes (N ≤ 10) and does not increase further, even though the initial PAH abundance has been increased. The largest cluster size at which this maximum adsorbate abundance is reached depends on the initial PAH abundance and shifts with it to larger cluster sizes.

In general, the existence of this maximum adsorbate abundance can be explained by comparison of the adsorption rate and the clustering rate. A necessity for this is that PAH-PAH collisions are much more likely than PAH-dust particle collisions so that the PAH cluster growth is initially unaffected by the presence of the grains. Then the rate at which PAH clusters collide with another PAH cluster r_PAH-PAH is proportional to f (N_C)² (Eq. (1)), while the rate of collisions with dust particles r_PAH-dust is proportional to f (N_C), (Eq. (5)). Consequently, the timescale τ on which all PAHs have grown is proportional to f (N_C)⁻¹ (τ = f(N_C)/r_PAH-PAH² ∝ f(N_C)/f(N_C)² =f(N_C)⁻¹. Accordingly, the absolute number of PAHs and clusters that can adsorb during this time is independent of the initial amount of PAH monomers f (N_C), (f_adsorbed (N_C) = τ · r_PAH-dust ∝ f (N_C)⁻¹ · f (N_C) = const). Most importantly, the maximum adsorbed amount for small clusters (N ≤ 10) is already reached when the PAH abundance is equal to the ISM abundance which becomes relevant when we analyse the PAH desorption in the next section.

In summary, the presence of PAHs and dust grains in a strongly UV-shielded region of the protoplanetary disc inevitably leads to clustering of PAHs with subsequent adsorption of PAH monomers and PAH clusters onto dust grains. Dust grains inhibit the further growth of PAH clusters, which otherwise would grow to carbonaceous nanograins. The adsorbed PAH clusters follow a size distribution whose number distribution is dominated by monomers and small clusters, but mass-wise is dominated by large clusters. Based on this, we estimate in the next section how difficult it is to recover these clusters and monomers from the dust grains depending on their size.

Fig. 3 PAH coagulation model without the presence of dust grains at r = 10 au and θ = 0 in the coagulation layer for an initial ISM PAH abundance of coronene monomers. The distribution is shown for selected time cuts after the start of the simulation. The growth of PAH monomers to very small carbonaceous grains is very fast, and within years only the largest clusters remain.

Fig. 4 Similar to Fig. 3 with the same time cuts, but in the presence of dust grains. The dashed lines show the size distribution of adsorbec PAH clusters. Initially, the density of PAHs is high enough that clustering is faster than adsorption. Over time, PAHs get adsorbed on dus grains until adsorption and clustering act on similar timescales. Ther growth is halted, and clusters will be stuck on grains. Most PAHs an adsorbed as large clusters.

Fig. 5 Similar to Fig. 4, but the particle size distribution is shown. While most of the mass is adsorbed at the end of the coagulation, the most frequent cluster size is the monomer, which can only adsorb during the early stages of coagulation.

Fig. 6 Size distribution of adsorbed PAH clusters as a function of cluster size for initial PAH abundances in the range 0.01–10 times the ISM abundance. For small clusters (N ≤ 10) a maximum abundance is reached with 0.1 × ISM abundance that cannot be increased by further increasing the initial PAH abundance.

3.2 Desorption of PAH clusters from grains

In this section, we investigate how fast the adsorbed clusters that formed in the coagulation layer can desorb again from the dust grains when they are brought into a UV-rich environment such as the photosphere. For this purpose we use the Monte Carlo model described in Sect. 2.2 to determine the energy probability distribution G(E) of adsorbed PAH clusters to obtain the thermal and UV-assisted desorption rate. We apply our model to clusters with a size of 1–16 monomers, and determine the desorption of the whole clusters from the grains as a function of the distance r from the central star. These calculations are performed for various PAH species whose monomers have between 24 (coronene) and 66 (circumovalene) C atoms and for which the absorption spectrum is available in the Cagliari theoretical spectral database of polycyclic aromatic hydrocarbons (Malloci et al. 2007).

Figure 7 shows the desorption rate of different cluster sizes for selected PAH species in the standard Herbig disc model STD from Table 1 assuming that all PAHs are fully exposed to the stellar UV radiation. The calculations for additional PAH species can be found in Fig. C.1. The total desorption rate is shown as solid lines and the thermal evaporation rate of clusters at the equilibrium temperature of the grains is shown as a dash-dotted line (using Eqs. (13), (14), and (23)). For all PAH species the desorption rate in the inner disc regions is completely dominated by the thermal evaporation of the clusters from the grain surface. Recalling the rapid dissociation of PAH clusters in the photosphere (Lange et al. 2021), this high evaporation rate implies that PAHs are exclusively present as gas-phase monomers in the inner photosphere. We note that the slightly lower thermal desorption rate for monomers and small clusters reflects the heat bath correction in Eq. (14).

Increasing the distance to the central star, the thermal desorption rate decreases due to the lower equilibrium temperature of the dust grains and desorption through single UV photon absorption by the PAH cluster becomes the main cause of cluster evaporation. The evaporation rate decreases approximately with the distance square r⁻² due to the reduction of the photon density, but only approximately since the dust temperature continues to decrease and the energy-temperature relation is not linear. Overall, the larger a PAH cluster is, the less likely it is to evaporate from the grain surface as the peak excitation temperature decreases due to the larger heat capacity of the cluster. The same effect occurs when a larger PAH species is considered, but the evaporation rate is further reduced by the fact that the grain-contact PAH monomer of the cluster is larger, resulting in a stronger van der Waals bond to the dust grain. Hence, a coronene dimer (C₂₄H₁₂)₂ detaches more quickly than a dicoronylene monomer (C₄₈H₂₀).

We compare the desorption rates of the clusters to the timescale these clusters would be exposed to the full stellar UV radiation in the protoplanetary disc turbulence model (see Sects. 2.4 and 3.3). The inverse of the mean-residence times in the photosphere are plotted as dotted lines in Fig. 7. Provided the cluster is large enough that the desorption rate is lower than the inverse mean residence time in the photosphere, the cluster can remain adsorbed on the grains. Therefore, we define the size at which this condition is met as the critical cluster size. We track the critical cluster size and use it in Sect. 3.5 to determine how many PAHs can be recovered from the grains. The derived critical cluster sizes at 10 au can be found in Table 2.

In particular, we find that the evaporation rate of circumcoronene monomers is much lower than the exposure rate except for the inner disc regions. Therefore, PAHs with more than 54 carbon atoms will stay frozen out on dust grains in most parts of the disc. In view of these results, we can conclude that the coagulation and freeze-out of PAHs has two effects if these processed PAHs become exposed to stellar UV radiation. The first effect is a permanent freeze-out in most parts of the disc as a large fraction of the PAH clusters will grow to sizes that cannot be recovered. Even though the effectiveness of this process depends on the PAH species, we expect that the gas-phase PAH abundance is smaller than the abundance in the ISM. The second is a size selection effect. Under the same conditions, smaller PAH species are more likely to be recovered than larger PAH species. Thus, we expect only small to medium gas-phase PAHs (N_C,0 ≤ 54) in a vertically mixing disc, and the photosphere is expected to be dominated by the small PAH molecules. We note that we only consider a single PAH species in this model. Clusters of mixed PAH species are more difficult to treat because any PAH species could be the bonding molecule to the grain, which will affect the desorption rates. Next, we want to quantify the strength of this depletion mechanism depending on the PAH abundance and dust grain population after one cycle of coagulation, freeze-out, and desorption.

Fig. 7 Desorption rates in the Herbig disc STD photosphere calculated with the Monte Carlo model for increasing cluster sizes, radii, and PAH species. Shown are thermal evaporation rates calculated from Eqs. (13) and (23) (dash-dotted lines) and turbulent mean-residence time in the photosphere calculated in Sect. 3.3 (dotted lines). For increasing PAH size and cluster size, desorption becomes slower until clusters are unlikely to desorb at all. PAH coagulation and adsorption on dust grains therefore lead to loss of gas-phase PAHs to the grains. More PAH species are available in Fig. C.1.

3.3 The structure of the PAH disc

We use the radiative transfer code RADMC-3D³ and supplementary tool radmc3dPy⁴ with our disc model STD to estimate the expected UV field between the midplane and the photosphere as this region is attenuated. For the opacities of the dust grains, we use the tool OpTool (Dominik et al. 2021) with the standard DIANA opacities (Woitke et al. 2016). We measure the UV field in the range 6–13.5 eV in units of the mean interstellar UV field G₀ = 5.33 × 10⁻¹⁴ erg cm⁻³ (Habing 1968). As desorption occurs only through single photon events, we can use a power law description similar to that in Lange et al. (2021) and scale the desorption rate by G₀ assuming that the attenuation from the dust grains does not change significantly between 6 and 13.5 eV. Figure 8 shows the strength of the UV field G₀ in the given disc model. The τ_UV = 1 line (black) where G₀ has dropped by 1/e and the coagulation front (red) where the adsorption and desorption rate of PAH monomers are equal. We find that the disk is vertically structured in three physically layers: the photosphere, the quiet layer, and the coagulation layer.

In the photosphere, the PAH molecules are fully exposed to stellar radiation. Small clusters can evaporate from grains on the timescales given in Sect. 2.1, and all clusters that desorb can also be dissociated on a faster timescale Lange et al. (2021). Therefore in this layer, PAHs exist mainly as gas-phase monomers and large adsorbed clusters on grains. As gas-phase PAHs move into the quiet layer, the UV field weakens and the density increases. However, as shown by Thi et al. (2019), PAHs are the major carrier of negative charges at these heights. Therefore, we do not expect cluster formation through collisions to be effective because of the electric repulsion force between two molecules. At the same time, small adsorbed clusters coming from the midplane stay adsorbed on the grains as the UV field is significantly weaker than in the photosphere. Hence, gas-phase monomers and small adsorbed cluster rarely desorb, cluster, or adsorb. For simplicity, we neglect all these processes here as a full model tracking the state of all PAHs is required. The upper boundary of the coagulation layer, the coagulation front, is located where the UV field equals the interstellar field (G₀ ≈ 1). At this location, the conditions match as PAH charging is not significant any longer, due to the lack of UV photons, and the densities are high enough that coagulation and adsorption are faster than their UV-driven counter processes.

3.4 Recoverable PAH abundance

From the results of the previous sections, we derive the number of PAHs that can re-enter the gas phase (hereafter recoverable PAHs) as these are the major emitters of the infrared bands through their strong temperature fluctuations. In order to be recoverable, a PAH cluster must be able to desorb from the grain and dissociate within the mean residence time in the photosphere. For simplicity, we approximate the total time to become a monomer from a cluster frozen on a grain by the desorption timescale as this is the longest timescale, otherwise a full model tracing the PAHs through the photosphere and quiet layer would be necessary to estimate the time for dissociation and desorption. To obtain a generally applicable result we parametrise the PAH and dust population by the total available collision cross-section per unit area in the PAHs σ_PAH compared to the total collision cross-section per unit area of the dust σ_dust. This is possible as the competition between clustering and freeze-out determines the fraction of PAHs below the critical cluster size that can desorb again.

To obtain the recoverable fraction of adsorbed PAH clusters, we run our coagulation model with varying initial PAH abundances and determine the adsorbed cluster size distribution ƒ(N_C)·According to the calculated critical cluster sizes, we then calculate what fraction of PAHs is recoverable from the dust grains by integration of all clusters that are smaller than the critical cluster size $$\xi =\frac{f_{\text{tot,recov}}}{f_{\text{tot}}}=\frac{{\displaystyle {\int }_0^{N_{\text{C,crit}}}f\text{d}{N}_{\text{C}}}}{{\displaystyle {\int }_0^{\infty }f\text{d}{N}_{\text{C}}}}.$$(38)

In our case, this ratio can be used generally throughout the disc as σ_PAH/σ_dust is constant (see Fig. B.1). However, when a different dust population is considered and even allowed to evolve, then σ_dust might not be suitable any more as the surface area can vary with height and in time.

Figure 9 shows the recoverable ratio σ_PAH/σ_dust and the lost fraction ξ at 10 au after one cycle of coagulation and freeze-out as a function of the initial surface ratio of PAHs to dust grains for different PAH species. The typical ISM PAH abundance is shown as a dotted line. If σ_PAH/σ_dust ≤ 10² (with the STD model this equals 0.1 times ISM abundance of PAHs), the recoverable σ_PAH/σ_dust ratio reaches a plateau, which is the maximum recoverable fraction from the dust grains. As explained in Sect. 3.1, this is due to the much faster growth of clusters compared to adsorption on dust particles, so that the absolute amount of PAHs that can be adsorbed is constant and independent of the PAH abundance. As a consequence after one cycle, any information on the initial σ_PAH/σ_dust is lost, as long as σ_PAH/σ_dust is large enough. Hence, no higher ratio than the maximum recoverable σ_PAH/σ_dust can be expected when PAHs have processed through one cycle.

In contrast, if the PAH abundance is very low at the beginning, almost all PAHs can be recovered due to slow cluster growth and only a small fraction of PAHs is lost in a coagulation-adsorption cycle. The lower σ_PAH/σ_dust becomes, the smaller the fraction that is lost per cycle as only small clusters form during the coagulation. Thus, when PAHs are processed in many cycles, the loss of PAHs over time slows down until the lost fraction is negligible and the PAH abundance is effectively constant. For this reason, we want to estimate how many PAHs are lost through adsorption within the typical protoplanetary disc life time.

Fig. 8 Intensity of the mean UV radiation field G0. The black line indicates the τUV = 1 line, and the red line shows where the monomer desorption rate is equal to the monomer adsorption. In the photosphere the small PAH clusters quickly desorb and dissociate into monomers. In the coagulation all gas-phase PAHs coagulate and adsorb on dust grains, while in the quiet layer monomers, clusters, and PAH-carrying grains can co-exist.

3.5 A simple turbulent PAH model

To determine the time evolution of the PAH abundance in a protoplanetary disc through persistent processing by clustering and adsorption, we consider a simple ID vertical model. We look at a vertical column and assume that gas-phase PAHs and the smallest dust grains can be transported up and down the disc by vertical eddies following the equations presented in Sect. 2.4. We then follow a coupled particle with a Monte Carlo approach through the disc to determine the mean residence time t_mrt in the disc photosphere and the mean cycle time i_mcyc. Starting with the entry into the coagulation layer, the mean cycle time is the time needed to re-enter the coagulation layer provided that the photosphere has been reached before. For the radial distance dependent boundaries of the photosphere and coagulation layer we use the results from the radiative transfer model (Fig. 8).

Figure 10 shows the distribution of t_cyc (upper panel) and t_rt(lower panel) for the α = 10⁻² case at 40 au as an example.

The mean cycle time can be estimated through an inverse-Gauss function, as the travel between photosphere and coagulation layer can be described as the first passage problem solved by Smoluchowsky (1915). However, a slight deviation is expected as the classical random walk uses a 50% probability for each step, while in the turbulence model the probabilities are slightly dependent on height. We find that both distributions are heavily skewed to the left and dominated by smaller values. Given that the number of cycles that can be completed in a typical protoplanetary disc life time of 2.5 Myr (Mamajek 2009) is low, the averages of t_cyc and t_mrt should be used with care; a detailed model tracking PAHs and its processing at each height would be more appropriate. However, for the sake of simplicity, we use the mean values in this model.

We calculate these timescales at all radial positions between r_min and r_max. Then we can estimate what the largest recoverable cluster size is (see Sect. 2.2) and determine the loss of unrecoverable PAHs deposited as clusters on grains. We assume that after one mean cycle time t_cyc all PAHs have undergone one complete cycle of mixing from the coagulation layer to the photosphere and back. Then we can estimate the average PAH loss in one cycle with Eq. (38) by obtaining the recoverable fraction of PAHs through the relative collision cross-section ratio σ_PAH/σ_dust. Figure 11 shows the depletion of gas phase PAHs through turbulent processes after 1 Myr and 5 Myr as a function of distance from the central star, turbulent strength, and PAH molecule. For a strong mixing disc (α ≈ 10⁻²), we expect a depletion of at least a factor 100 in the inner disc regions independent of the sampled PAH species as many cycles can be completed. With a lower α fewer cycles can be achieved in the same time span. For a very weakly mixing disc with α=10⁻⁴, no cycle can be achieved in 1 Myr. Beyond 50 au our disc model is not optically thick enough (we require G₀ ≈ 1) to allow for the formation of the coagulation layer. Hence, we do not expect turbulent processing in the outer parts of the disc as no new clusters can form. After 5 Myr, a depletion also occurs in the very inner regions (≈ 5 au), while in the α=10⁻² and α=10⁻³ cases the depletion is stronger and more extended.

Despite the clear α turbulence dependence, we note that the used dust model from Fromang et al. (2007) was constructed to fit a α=10⁻³ disc below four gas pressure scale heights, which we have used for all three α cases. Therefore, with a different α the vertical distribution of the grains will change (e.g. Dullemond & Dominik 2004) as well as the height of the coagulation front and photosphere, which are tightly coupled to the optical thickness caused by the dust grains. These effects were not considered in this study.

Fig. 9 Recoverable σPAH/σđust ratio and lost fraction of PAHs after one cycle of clustering for different PAH species and different initial σPAH using the standard dust model. The dotted black line assumes a standard ISM abundance of PAHs. For PAH abundances higher than σPAH/σdust ≥ 30 (≈ 0.1 × ISM) the maximum abundance is reached (Sect. 3.1) and information on the initial PAH abundance is lost. For σPAH/σdust < 1, almost all PAHs can be recovered because adsorption is much faster than coagulation and processing has no effect on the PAH abundance any longer as all PAHs are preserved.

Fig. 10 Distribution of measured cycle times tcyc and residence time in the photosphere trt with a Monte Carlo framework at 40 au with α = 10−2. The cycle time can be approximated with an inverse-Gaussian distribution. The mean-residence time is dominated by short timescales corresponding to an immediate mixing down into the quiet layer after reaching the photosphere.

4 Discussion

4.1 Initial conditions and disc infall

In our study we assume that the PAH abundance at the beginning of the disc stage is comparable to the value found in the ISM. However, we neglect the evolution of PAHs during the infall stage of the disc starting from the molecular cloud phase. As the local density increases during the infall, coagulation of PAHs and adsorption on pristine dust grains might occur. In order for the PAHs to not be dissociated or desorbed immediately, the local UV field must match the G₀ ≈ 1 condition. Whether this is possible is questionable as the local UV field is driven by the external UV field and the evolving and embedded protostar. In the case that these conditions can be achieved in the infall stage, then PAHs will already enter the disc stage as adsorbed clusters. As the first coagulation cycle already reduces the recoverable fraction to 1/30 of the initial value (e.g. for coronene), the disc would lose a substantial fraction of gas-phase PAHs. Consequently, if all PAHs were coagulated and adsorbed once already before the disc stage, then the outer regions of the disc (≥50 au) would show a depletion of PAHs and the additional depletion through turbulent processing in the inner disc would have a much smaller effect. Hence, the depletion of PAHs in the inner disc by continuous vertical mixing would be much less prominent or not even existent.

However, in our Herbig disc model we only take into account photons with less than 13.6 eV and do not consider the hard photons expected for T Tauri star discs originating from stellar accretion. The interaction of individual PAH molecules in a turbulent discs with such hard photons has been investigated by Siebenmorgen & Krügel (2010) and Siebenmorgen & Heymann (2012). The authors conclude that typical X-ray luminosities can be an efficient mechanism to destroy individual PAH molecules in the disc. However, as PAH clusters have a higher heat capacity, and frozen-out PAHs can transfer heat much faster than the IR-emission timescale considered in these works, the interaction of hard photons with frozen-out PAHs and PAH clusters should be considered in a future study.

Fig. 11 Turbulent depletion of PAHs through processing by clustering, adsorption, desorption, and photodissociation in the disc at 1 Myr and 5 Myr after formation. The calculations were performed for coronene C24H12, bisanthene C28H14, and ovalene C32H14. Strong (α = 10−2) and intermediate turbulence (α = 10−3) cause a depletion of recoverable PAHs by a factor of 50–500 compared to the initial ISM abundance in the model STD. At distances larger than 50 au we do not find a layer where adsorption and therefore a full cycle can happen. For larger species than ovalene, the depletion increases when the PAH species is larger.

4.2 Typical PAH size and abundance in protoplanetary discs

Our results from the desorption model strongly constrain the expected PAH size in protoplanetary discs. Even a single cycle of clustering and adsorption leads to a size selection effect, where the expected abundance of smaller PAHs (coronene) is a factor of 3 higher than the abundance of medium-sized PAHs (ovalene) if initially they had the same abundance. Additionally, PAHs that are larger than circumcoronene (C₅₄H₁₈) cannot be desorbed, even as a monomer. Therefore, according to our model, we expect a distribution of small PAHs with less than 66 C atoms to be the dominant gas-phase PAH species and, due to their low heat capacity, also the dominant emitters of the IR features as the temperature of the adsorbed PAHs is similar to the grain temperature. Our result agrees with the study of Seok & Li (2017), where the authors fit emission models to protoplanetary disc spectra. Most of their best fit models predict PAH sizes between 3.5 and 6 µm, which translate to PAH monomer sizes between 15 and 44 C atoms, where the upper size limit is a bit smaller than our largest desorbable PAH circumcoronene (C₅₄H₁₈). Further, Seok & Li (2017) report a small negative correlation between PAH size and stellar age, which the authors expect to be outgassing of comets or grinding of planetesimals. However, as the depletion per cycle for larger PAH species in our model is greater than for smaller PAH species and locally all PAHs experience the same mixing, we expect a correlation between PAH size and disc age by turbulent processing as well.

Typically, the common assumption is that the monomer size of a protoplanetary disk PAH is N_C,0 = 100 C atoms (Siebenmorgen & Krügel 2010; Maaskant et al. 2014). Similarly to the observationally constrained size by Seok & Li (2017), we favour PAH species with less than 66 C atoms as the typical disc PAH, as larger species can only desorb thermally in the disc. Therefore, we propose a significant size difference for PAHs between the interstellar medium and protoplanetary discs. The ISM estimates for astronomical PAHs are usually based on observations of the ISM (Allamandola et al. 1989) and nebulae like NGC7023 (Croiset et al. 2016), which emphasises the difference between the interstellar PAH population and disc PAHs. Unfortunately, the ISM estimates are derived from temperature fluctuations of the PAHs using the direct link between heat capacity and size. This size estimate is hardly possible in protoplanetary discs without spatially resolved images at multiple PAH emission wavelengths and a disc model to estimate emission heights and optical depth. The James Webb Space Telescope (JWST) with its instruments NIRcam and MIRI will be one step to understand the difference between ISM PAHs and disc PAHs. With its filters at 3.3 µm and 11.3 µm, we can start to analyse possible PAH structures in protoplanetary discs.

We also expect a depletion of PAHs in protoplanetary discs because of the processing of PAHs through coagulation, adsorption, desorption, and dissociation during their life times. Even one cycle can lower the observable PAH abundance by a factor of 30 in our standard model. Geers et al. (2009) find after an analysis of protoplanetary disc spectra a depletion factor of PAHs of 10–20 compared to the ISM. Maaskant et al. (2014) also reduce their PAH mass to obtain a PAH-to-dust mass ratio of 5 × 10⁻⁴ (corresponding to 1/5 ISM PAH abundance) in order to match their radiative transfer model to the observed disc spectrum of HD 97048. These results might indicate that those discs experienced no cycles or only a few cycles of coagulation so that most PAH were preserved in the gas-phase. In contrast, discs without detectable PAH emission might have experienced more cycles of PAH processing so that most PAHs in these discs are frozen out on grains and are unable to enter the gas-phase again.

4.3 Dust population and vertical mixing

Our evolution model is parametrised by σ_PAH/σ_dust making it possible to adopt it to other dust populations than our standard model (STD). Therefore, the dust population is a free parameter that allows for a variation in the depletion of the gas-phase PAHs. The first PAH processing cycle is especially important for the evolution of the PAHs as the largest fraction is lost if σ_PAH/σ_dust is limited by the maximum recoverable fraction. As this value is achieved with the ISM abundance of PAHs, the PAH depletion is strongly dependent on the collisional cross-section σ_dust of the dust population. More precisely, in the case of rapid coagulation (σ_PAH/σ_dust > 30), the PAH loss during the first cycle is proportional to σ_dust as the abundance of PAHs after one cycle is independent of the initial σ_PAH/σ_dust ratio.

As a consequence, two similar discs with similar host stars but different dust populations will have a different depletion of PAHs. The larger the dust grains and smaller the total grain population cross-section, the more PAHs will stay frozen out on the dust and the larger the depletion. For this purpose, we introduce three other dust populations: a model with only small grains (SGs), a model with dominantly large grains (LGs), and a model with very large grains (VLGs). For each of them we can calculate the initial σ_PAH/σ_dust ratio and then the PAH depletion by one cycle of coagulation, adsorption, and desorption. Using the initial σ_PAH /σ_dust ratios in Table 1 we find a PAH depletion factor of ≈300 for the standard model (STD), a depletion factor of ≈2 in the small grain model (SG), a depletion factor of ≈20 000 for the large grain model (LG), and an extreme depletion of ≈ 100 000 for the very large grain population (VLG).

A Herbig star disc with a sufficient UV field but no visible PAH signatures must therefore have depleted its PAHs in most parts of the disc. For the inner part (<50 au), the disc must have large grains present where PAHs coagulate. Furthermore, if the disc is weakly turbulent (α ≤ 10⁻⁴), the disc must be sufficiently old to have enough time to allow for at least a few PAH processing cycles. If instead the disc is intermediate or strongly turbulent (α > 10⁻⁴), the disc can be young as a PAH processing cycle only takes 1 Myr. Finally it must destroy its PAHs or not have any PAHs present at all. However for the depletion of PAHs in the outer disc regions, we require from our model that the processing of PAHs occurs before the disc stage, such as the infall; that there is a dense disc structure that induces G₀ ≈ 1 regions to also allow for adsorption; and that the PAHs are destroyed or that there are no PAHs at all. Furthermore, other effects that are beyond the scope of this work, such as disc structures, partial shadowing of the disc, inclination of the disc, or dust evolution will likely affect the observability of PAH features as well.

4.4 PAH desorption through high-velocity grain collisions

In this section we discuss the importance of dust grain collisions for the PAH evolution. So far, we have only considered desorption by photon-induced processes and neglected other desorption mechanisms. Therefore, we want to consider the relevance of mechanical processes such as the abrasion of PAHs by dust grain collisions. For this to be relevant, there must be a sufficient number of high-velocity impacts that mechanically erode the surface of the dust grains. We estimate the minimum velocity for these events by considering the necessary collision velocity required for the abrasion of graphitic surfaces. Jones et al. (1996) reports for graphitic surfaces that the minimum critical impact velocity for the occurrence of splinters is ≈1 km/s. These required velocities are supported by the laboratory studies of Zamith et al. (2020), who report a required collisional energy of 0.7–1 eV (0.8–1 km s⁻¹) for the collisional dissociation of pyrene clusters. Such high velocities cannot be achieved by Brownian motion of dust grains because extreme temperatures are required. However, turbulent collisions between dust grains can reach much higher collision velocities than pure thermal motion.

Cuzzi & Hogan (2003) and Ormel & Cuzzi (2007) give close-form expressions to calculate the relative particle velocities between two particles with different Stokes numbers in a turbulent disc. In their model the highest relative velocities are reached when a particle with a Stokes number St = t_stop/t_edd = 1 (the ratio of the particle stopping time to the eddy turnover time) is involved in a collision. Then the relative velocity is comparable to the local turbulent velocity of the gas: $$\text{Δ}\upsilon \approx {\upsilon }_{\text{g}}=\sqrt{\alpha }{c}_{\text{s}}.$$(39)

This means that even for an 1000 K gas an unrealistically high turbulence parameter of α ≈ 1 is required to achieve collision speeds of 1 km/s. These conditions are typically not met in protoplanetary discs despite the lack of St = 1 grains close to the photosphere. Therefore, we consider the mechanical removal negligible for the general evolution of adsorbed PAHs on dust grains.

4.5 Heterogeneous PAH clusters

One of our simplifying assumptions is that initially all PAHs are present as one species, and therefore we do not consider heterogenous clusters made from different PAH species. In reality, PAHs will be present as a mixture of different species, possibly a selection of the grandPAHs, a few stable and compact PAH molecules that seem to dominate the astronomical PAHs (Tielens 2013; Andrews et al. 2015), as extracted from photodissociation region (PDR) spectra in the ISM. The heat capacity of a PAH cluster only depends on the number of carbon atoms; therefore, the temperature fluctuations of a mixed cluster will be very similar to a homogenous cluster with the same number of carbon atoms. However, the size of the grain-touching PAH that binds to the dust grains plays an important role as the binding energy of the PAH species determines the detachability of the overall cluster. A cluster with a small PAH binding to the grain can be detachable, while the same cluster with a large PAH connecting to the grain might not be detachable. Thus, the number of carbon atoms in a heterogenous cluster is not sufficient to characterise the cluster and determine its chance of evaporation. In the case that a small PAH is binding to the grain, more large PAH species can be recovered than estimated through the homogenous model. The whole cluster will desorb and then the large PAHs will dissociate from the cluster as dissociation can also be achieved through multi-photon events. In the contrasting case that a large PAH is binding to the grain, then the cluster will be very unlikely to desorb. However, as smaller PAHs in the cluster have a lower binding energy to the cluster than the surface PAH to the grain, evaporation of small PAHs can be possible leaving only the medium to large PAH species in the adsorbed cluster behind. We therefore estimate that our observed size selection effect for homogenous clusters favouring small PAHs for protoplanetary disc conditions is also present if heterogenous clusters are considered, but it is not as effective.

In this regard we note that heterogenous clusters are able to rearrange their molecules. Bowal et al. (2019) use molecular dynamics (MD) to study the morphology of heterogeneous clusters at high temperatures (800 K and 1600 K), where the individual monomers can rearrange within the cluster. Regardless of the initial arrangement of the heterogeneous clusters, the authors find that, on average, the larger PAHs are found closer to the cluster core than the smaller PAHs, which preferably accumulate on the surface of the cluster. However, the clusters in these studies are significantly larger than the clusters we can detach in our models (16 monomers). In addition, the clusters we simulated have a much lower equilibrium temperature, so that rearrangement within the cluster can only take place partially under adsorption of UV photons. Separate MD studies of gas-phase clusters and adsorbed clusters are needed to clarify whether systematic rearrangements of small PAHs towards the surface are possible and relevant.

5 Conclusions

We have modelled the coagulation of PAHs and their freeze-out on dust grains in a protoplanetary disc environment around a typical Herbig Ae/Be star. We find that under disc conditions coagulation is very fast. Considering only PAHs, clusters can grow to sizes with more than 10⁵ cluster members within a year. However, in the presence of dust grains, cluster growth is suppressed as clusters start to freeze out on grains. In this set-up, all PAHs adsorb on a timescale of a year, but the largest cluster sizes that can grow are a factor of 10–100 smaller than in the case without dust grains. Modelling the UV-driven desorption of these clusters, we find that a critical cluster size exists where PAH clusters cannot desorb any longer. The larger the cluster building PAH species, the smaller the critical cluster size. This number is on the order of 10–20 cluster members for small PAHs, such as coronene (C₂₄H₁₂), and decreases until monomers as large as circumovalene (C₆₆H₂₀) cannot be desorbed at all through UV photons. Based on these results, we expect a PAH depletion of at least a factor of 10 in disc environments compared to the interstellar medium if only one cycle of coagulation, adsorption, desorption, and dissociation can happen either during the disc stage or the infall stage. Given that larger clusters of smaller PAH species can desorb more easily compared to clusters made from large PAH species, we expect a size selection effect with which an initial PAH species distribution is expected to shift to smaller PAHs with every cycle.

Furthermore, we have used these results in a vertical mixing model where the turbulent turnover time of the largest eddies is the dynamical time. With a radiative transfer model, we determined that the inner disc (<50 au) provides enough UV shielding (G₀) so that PAH clusters can freeze out on dust grains. However, the outer disc does not provide enough shielding to meet these conditions so that a continuous processing of PAHs is unlikely as PAHs simply cannot adsorb on dust grains. We find that through turbulent processing, the inner disc can have multiple coagulation-adsorption cycles over the lifetime of the disc. There the depletion of PAHs depends on the turbulent parameter α, but also on the PAH species and age. However, we find that a weak mixing disc α < 10⁻³ is barely able to process PAHs in less than 5 Myr. Therefore, we expect a similar PAH abundance in the inner and outer disc given that the disc is sufficiently young (1 Myr< t < 5 Myr).

Nevertheless, the abundance difference depends on the initial PAH content of the disc stage and if the PAHs have been processed in the infall stage to adsorbed clusters already. The James Webb Space Telescope with its instruments NIRcam and MIRI will provide crucial information to further investigate the life of PAHs in protoplanetary discs.

Acknowledgements

The authors thank the referee Ralf Siebenmorgen for his comments to improve the quality and comprehension of the manuscript. K.L. acknowledges funding from the Nederlandse Onderzoekschool Voor Astronomie (NOVA) project number R.2320.0130. C.D. acknowledges funding from the Netherlands Organisation for Scientific Research (NWO) TOP-1 grant as part of the research program “Herbig Ae/Be stars, Rosetta stones for understanding the formation of planetary systems”, project number 614.001.552. Studies of interstellar PAHs at Leiden Observatory are supported through a Spinoza award from the Dutch research council, NWO.

Appendix A Calculation of the cluster size

Fig. A.1 PAH radii given the model in Appendix A. For large clusters, the clusters grow in three dimensions, like a dust grain. Assuming PAHs grow initially in stacks, we linearly interpolate between the monomer and the cluster size where a stack has similar width and length.

Given the growth cases in section 2.1, we need to calculate the collisional cross-section of all cluster sizes. Therefore, we define the effective cluster radius as $$r=\{\begin{array}{ll}{r}_0=0.9\text{Å}\sqrt{N_{\text{C,0}}}\hfill & \text{if\hspace{0.17em}}N=1\hfill \\ r_1=0.9\text{Å}(\left(N_{\text{C}}-N_{\text{C,0}}\right)\frac{r_2\left(m{N}_{\text{C,0}}\right)-r_0}{m{N}_{\text{C}}-N_{\text{C,0}}}+r_0)\hfill & \text{if\hspace{0.17em}}1\text{\hspace{0.17em}<\hspace{0.17em}N\hspace{0.17em}<\hspace{0.17em}m}\hfill \\ r_2=09\text{Å}{\left(N_{\text{C}}\right)}^{1/3}\hfill & \text{if\hspace{0.17em}}N\ge m\hfill \end{array},$$(A.1)

where N is the number of molecules in the cluster. The parameter r₀ defines the size of a monomer, which is proportional to the carbon atoms in the PAH $\sqrt{N_{\text{C,0}}}$ since a PAH monomer grows in a plane. We expect the largest cluster to grow as dust grains, in three dimensions without structure, therefore $r\propto N_{\text{C}}^{1/3}$. In our model this occurs once the length of a stack (Rapacioli et al. 2005) of PAHs exceeds the diameter of a monomer. This occurs at a cluster with m monomers: $$m=\frac{2\cdot 0.9\text{Å}\sqrt{N_{\text{C,0}}}}{1.42\text{Å}}+1.$$(A.2)

Here we assume that the distance between two stack-bound monomers is similar to the interlayer-distance of graphene (1.42Å). Between the monomer and grain-like case, PAHs will likely grow as single stacks, even though multi-stack structures are also allowed depending on the temperature of the cluster. As this introduces an orientation dependent collisional cross-section, we choose to simplify this problem by linear interpolation between r₀ and r₂.

Appendix B Test case for dust density distribution

In order to test the vertical dust density distribution we want to compare our implementation with figure 13 in Fromang et al. (2007). In their toy model, they use the following coding parameters: gravitational parameter GM = 1, sound speed at R₀ is c₀ = 1, midplane gas density ρ₀ = 1, reference radius R₀ = 1, aspect ratio H/R = 0.1, and Stokes number in midplane at reference radius (Ωτ_s)₀ = 0.001. Without specifying a disc model, this can be achieved by setting the grain size to a = 0.0001 and the bulk grain density ρ_s = 1. Because Fromang et al. (2007) evaluate their disc model not at R₀ but average between 3 ≤ R ≤ 5, we evaluate the vertical dust grain distribution at R = 4 rather than averaging. As a consequence, the sound speed c_s, the Kepler frequency Ω, the scale height H, and the midplane gas density ρ need to be calculated at R = 4 according to their disc model. This can be done using $c_{\text{s}}=c_0\sqrt{R_0/R}$, $\text{Ω=\hspace{0.17em}}\sqrt{R^{-3}}$, H = 0.1R, and ρ = ρ₀ (R₀/R)^1.5. These values are used with the given expressions in equations (26) – (30) to solve the differential equation (26) with a fourth-order Runge-Kutta solver.

Fig. B.1 Diagnostic plot to test the vertical grain distribution (same as figure 13 in Fromang et al. 2007). Left: Variation of the midplane velocity fluctuations. Dotted: δυz,mid/cs = 0.025, dashed: δυz,mid/cs = 0.05, dash-dotted: δυz,mid/cs = 0.075. Right: Variation of the upper layer velocity fluctuations. Dotted: δυz,up/cs = 0.075, dashed: δυz,up/cs = 0.15, dash-dotted: δυz,up/cs = 0.3.

Figure B.1 shows the resulting density distribution with varied velocity fluctuations as in figure 13 in Fromang et al. (2007). This figure closely agrees with figure 13.

Appendix C Desorption rates for additional PAHs

Figure C.1 shows additional PAH species for which we have performed desorption calculations. Bisanthene, tetraben-zocoronene, and dicoronylene are species that are between the smallest and largest shown species by size in the results section. Circumovalene is the smallest PAH species that we cannot desorb through UV photons from the grains any more.

Fig. C.1 Desorption rates in the Herbig disc STD photosphere calculated with the Monte Carlo model for increasing cluster sizes, radii, and PAH species. Shown are thermal evaporation rates calculated from eq. (13) and eq. (23), (dash-dotted lines) and turbulent turn-over rates calculated from eq. (31), (dotted lines). For increasing PAH size and cluster size, desorption becomes slower until clusters are unlikely to desorb. PAH coagulation and adsorption on dust grains therefore lead to loss of gas-phase PAHs to the grains.

References

All Tables

All Figures

Fig. 1 Sketch of PAH processing cycle in the protoplanetary disc driven by vertical mixing. PAH monomers will eventually mix into the UV depleted coagulation layer where PAH collisions lead to sticking. The cluster growth continues until the PAH clusters freeze out on dust grains. If the cluster carrying grains are mixed into the UV-rich photosphere, PAH monomers and small clusters can desorb while large clusters remain on the grain. Every cycle of coagulation, adsorption, desorption, and dissociation leads to a depletion of gas-phase PAHs as with every cycle large clusters form and remain on the grains.

In the text

Fig. 2 Ratio of total PAH collisional cross-section σPAH to the total dust collisional cross-section σdust for the standard Herbig disc model STD with ISM PAH abundance. As the surface area of the dust is dominated by the smallest grains, which are closely coupled to the gas, the ratio is almost constant σPAH/σdust ≈ 310 resp. log(σPAH/σdust) ≈ 2.5 at relevant coagulation altitudes (z < 3H). Only in the outer disc does this ratio change; however, the relevant coagulation height decreases as well (see Fig. 8). The gas pressure scale heights of z =1 H, 3 H, 5 H are indicated by dashed lines.

In the text

	Fig. 3 PAH coagulation model without the presence of dust grains at r = 10 au and θ = 0 in the coagulation layer for an initial ISM PAH abundance of coronene monomers. The distribution is shown for selected time cuts after the start of the simulation. The growth of PAH monomers to very small carbonaceous grains is very fast, and within years only the largest clusters remain.
In the text

Fig. 4 Similar to Fig. 3 with the same time cuts, but in the presence of dust grains. The dashed lines show the size distribution of adsorbec PAH clusters. Initially, the density of PAHs is high enough that clustering is faster than adsorption. Over time, PAHs get adsorbed on dus grains until adsorption and clustering act on similar timescales. Ther growth is halted, and clusters will be stuck on grains. Most PAHs an adsorbed as large clusters.

In the text

	Fig. 5 Similar to Fig. 4, but the particle size distribution is shown. While most of the mass is adsorbed at the end of the coagulation, the most frequent cluster size is the monomer, which can only adsorb during the early stages of coagulation.
In the text

	Fig. 6 Size distribution of adsorbed PAH clusters as a function of cluster size for initial PAH abundances in the range 0.01–10 times the ISM abundance. For small clusters (N ≤ 10) a maximum abundance is reached with 0.1 × ISM abundance that cannot be increased by further increasing the initial PAH abundance.
In the text

Fig. 7 Desorption rates in the Herbig disc STD photosphere calculated with the Monte Carlo model for increasing cluster sizes, radii, and PAH species. Shown are thermal evaporation rates calculated from Eqs. (13) and (23) (dash-dotted lines) and turbulent mean-residence time in the photosphere calculated in Sect. 3.3 (dotted lines). For increasing PAH size and cluster size, desorption becomes slower until clusters are unlikely to desorb at all. PAH coagulation and adsorption on dust grains therefore lead to loss of gas-phase PAHs to the grains. More PAH species are available in Fig. C.1.

In the text

Fig. 8 Intensity of the mean UV radiation field G0. The black line indicates the τUV = 1 line, and the red line shows where the monomer desorption rate is equal to the monomer adsorption. In the photosphere the small PAH clusters quickly desorb and dissociate into monomers. In the coagulation all gas-phase PAHs coagulate and adsorb on dust grains, while in the quiet layer monomers, clusters, and PAH-carrying grains can co-exist.

In the text

Fig. 9 Recoverable σPAH/σđust ratio and lost fraction of PAHs after one cycle of clustering for different PAH species and different initial σPAH using the standard dust model. The dotted black line assumes a standard ISM abundance of PAHs. For PAH abundances higher than σPAH/σdust ≥ 30 (≈ 0.1 × ISM) the maximum abundance is reached (Sect. 3.1) and information on the initial PAH abundance is lost. For σPAH/σdust < 1, almost all PAHs can be recovered because adsorption is much faster than coagulation and processing has no effect on the PAH abundance any longer as all PAHs are preserved.

In the text

	Fig. 10 Distribution of measured cycle times tcyc and residence time in the photosphere trt with a Monte Carlo framework at 40 au with α = 10−2. The cycle time can be approximated with an inverse-Gaussian distribution. The mean-residence time is dominated by short timescales corresponding to an immediate mixing down into the quiet layer after reaching the photosphere.
In the text

Fig. 11 Turbulent depletion of PAHs through processing by clustering, adsorption, desorption, and photodissociation in the disc at 1 Myr and 5 Myr after formation. The calculations were performed for coronene C24H12, bisanthene C28H14, and ovalene C32H14. Strong (α = 10−2) and intermediate turbulence (α = 10−3) cause a depletion of recoverable PAHs by a factor of 50–500 compared to the initial ISM abundance in the model STD. At distances larger than 50 au we do not find a layer where adsorption and therefore a full cycle can happen. For larger species than ovalene, the depletion increases when the PAH species is larger.

In the text

	Fig. A.1 PAH radii given the model in Appendix A. For large clusters, the clusters grow in three dimensions, like a dust grain. Assuming PAHs grow initially in stacks, we linearly interpolate between the monomer and the cluster size where a stack has similar width and length.
In the text

	Fig. B.1 Diagnostic plot to test the vertical grain distribution (same as figure 13 in Fromang et al. 2007). Left: Variation of the midplane velocity fluctuations. Dotted: δυz,mid/cs = 0.025, dashed: δυz,mid/cs = 0.05, dash-dotted: δυz,mid/cs = 0.075. Right: Variation of the upper layer velocity fluctuations. Dotted: δυz,up/cs = 0.075, dashed: δυz,up/cs = 0.15, dash-dotted: δυz,up/cs = 0.3.
In the text

Fig. C.1 Desorption rates in the Herbig disc STD photosphere calculated with the Monte Carlo model for increasing cluster sizes, radii, and PAH species. Shown are thermal evaporation rates calculated from eq. (13) and eq. (23), (dash-dotted lines) and turbulent turn-over rates calculated from eq. (31), (dotted lines). For increasing PAH size and cluster size, desorption becomes slower until clusters are unlikely to desorb. PAH coagulation and adsorption on dust grains therefore lead to loss of gas-phase PAHs to the grains.

In the text