© ESO 2018

1 Introduction

Protoplanetary disks have complex structures due to the many physical and chemical processes that shape their environment. This includes but is not limited to radiative heating from the central pre-main sequence (PMS) star, viscous heating, molecular line cooling, photodissociation and ionization, dust grain growth and radial drift, and the freeze-out of molecular species in cold disk regions (for a review of disk structure and evolution, see Williams & Cieza 2011). These processes culminate in the formation of terrestrial and giant planets, often before the gas disk is dispersed. The study of protoplanetary disk structure while the disks still retain their large gas reservoirs is important to understand the environments in which planets will form.

Several gas-rich protoplanetary disks have been imaged at high spatial resolution with the Atacama Large Millimeter/submillimeter Array (ALMA), revealing physical structures that may indicate the presence of low-mass companions or planets in the disk, for example, IRS 48, HD 142527, HD 100546, HL Tau, TW Hya, HD 97048, HD 163296, HD 169142, and AS 209 (van der Marel et al. 2013; Casassus et al. 2013; Walsh et al. 2014a, 2016; ALMA Partnership et al. 2015; Andrews et al. 2016; Isella et al. 2016; Fedele et al. 2017, 2018). Disks have a stratified vertical structure with an atomic upper layer, a chemically active warm intermediate layer, and a dense, cold midplane. The environment most conducive to planet formation is at the disk midplane, where various molecular species such as H₂O, CO₂ , CO, and N₂ freeze out onto dust grains, creating an icy mantle that enhances grain sticking efficiency (Bergin & Tafalla 2007; Blum & Wurm 2008). Molecular ices can then be more easily incorporated into the bulk of planetary bodies than their gas-phase counterparts. The location at which ~50% of a given molecule has condensed into ice is called a snow line or ice line.

Probing the conditions of the midplane of the disk is difficult. Dust and molecular line opacities can obscure lower layers of the disk, particularly at radii close to the central star. Molecular snow lines can reside too near to the star to be directly observed, as is the case for H₂ O (Zhang et al. 2013; Piso et al. 2015; Banzatti et al. 2015). They can also be obscured by opacity effects in the outer disk, as is the case for the farther out CO snow line, where ¹²CO, ¹³ CO, and even C¹⁸O can remain optically thick at large radii (Qi et al. 2015; Fedele et al. 2017). Direct determination of the CO snow line can be done using emission from the rarest CO isotopolgues (e.g., Yu et al. 2016; Zhang et al. 2017), but only for the closest objects. To characterize the disk midplane environment, less abundant species must be observed which trace chemical processes occurring deep in the disk, such as molecular freeze-out. The DCO⁺ molecule has been suggested as an optically thin molecular tracer of the midplane regions around the CO snow line and as a simultaneous tracer of ionization occurring in the intermediate layers of the disk due to its formation via cold ( ≲30 K) and warm ( ≲100 K) deuterium fractionation pathways available in protoplanetary disks (Mathews et al. 2013; Favre et al. 2015; Huang et al. 2017).

The disk around HD 169142 makes an excellent test bed in which to explore the chemistry of DCO⁺ in protoplanetary disks and its usefulness as a tracer of disk midplane conditions. HD 169142 is one of a handful of disks found to have millimeter dust rings, and which also exhibits CO emission that extends beyond the edge of the millimeter grains (ALMA Partnership et al. 2015; Andrews et al. 2016; Isella et al. 2016; Walsh et al. 2016; Fedele et al. 2017). HD 169142 is an isolated system with a Herbig Ae spectral type A8 Ve star and stellar mass M_* = 1.65 M_⊙ (Grady et al. 2007; Blondel & Djie 2006). Recent distance measurements by Gaia put the system at a distance of d = 117 ± 4 pc (Gaia Collaboration 2016). The new distance results in a revised luminosity that is lower by a factor of ~0.65, which places the age of the system closer to ~10 Myr (Pohl et al. 2017). The new age estimate is older than previous estimates of 6${}_{-3}^{+6}$ Myr (Grady et al. 2007), but within the errors. With disk inclination i = 13° and position angle PA = 5° (Raman et al. 2006; Panić et al. 2008), the system is viewed close to face-on, allowing for accurate characterization of the radial distribution of the continuum and molecular line emission. With an estimated total gas mass of 1.9 ×10⁻² M_⊙ and ¹² CO extending out to ~200 AU (Fedele et al. 2017), the HD 169142 disk has a high concentration of gas. There is already known substructure in the dust around HD 169142. A hot inner ring of dust at ~0.2 AU was detected (Wagner et al. 2015) within a central dust cavity, and two dust rings at ~25 AU and ~60 AU are clearly visible in the 1.3 mm continuum with ALMA (Fedele et al. 2017) and in scattered light with the Gemini Planet Imager and the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument on the Very Large Telescope (VLT/SPHERE; Monnier et al. 2017; Pohl et al. 2017). The gap carved out between the rings may be indicative of ongoing planet formation. An outer gap at ~85 AU just beyond the edge of the 1.3 mm continuum emission was also detected in scattered light (Pohl et al. 2017) and in 7 and 9 mm emission with the Very Large Array (Macías et al. 2017). While the millimeter grains terminate at ~85 AU, the micron-sized grains are present throughout radial extent of the gaseous disk. Determining the midplane conditions of this disk would provide insight into the cold disk environment during the planet-building epoch.

This paper presents ALMA observations of the J = 3 − 2 transition of DCO⁺ toward HD 169142 and characterizes its distribution throughout the disk. Section 2 describes the observations and data reduction. The detection and distribution of DCO⁺ throughout the disk is detailed in Sect. 3. Modeling of the disk structure and DCO⁺ emission is explained in Sect. 4. Section 5 discusses the relationship between DCO⁺ and the disk environment, followed by the conclusions in Sect. 6.

2 Observations and reduction

The star HD 169142 (J2000: RA = 18^h 24^m29.776^s, Dec =–29°46′50.000′′) was observed with ALMA in band 6 (211–275 GHz) with 35 antennas on 2015 August 30 at a spatial resolution of ~0.3′′. The project code is ADS/JAO.ALMA#2013.1.00592.S. The data used in this work were reduced in the same manner as Fedele et al. (2017). Their paper provides further details on calibration, self-calibration, and continuum subtraction. Data reduction was performed with version 4.3.1 of the Common Astronomy Software Applications (CASA; McMullin et al. 2007). Images were created using the CASA task CLEAN, with natural weighting for the lines to enhance sensitivity.

The full data set contained observations of the 1.3 mm continuum and the molecular lines ¹² CO J = 2 − 1, ¹³ CO J = 2 − 1, C¹⁸ O J = 2 − 1, and DCO⁺ J = 3 − 2. The DCO⁺ J = 3 − 2 line at 216.1128 GHz was observed in the lower sideband and had a frequency (velocity) resolution of 61.0 kHz (0.084 km s⁻¹ ). This work focuses on the analysis of the DCO⁺ J = 3 − 2 data and makes use of the 1.3 mm continuum and C¹⁸O J = 2 − 1 images. A previous analysis of the continuum and the three CO isotopologue lines was reported in Fedele et al. (2017). Macías et al. (2017) presented a brief analysis of the C¹⁸O J = 2 − 1 and DCO⁺ J = 3 − 2 data. In their comparison of the molecular radial intensity profiles, the authors use a uv taper to increase the signal-to-noise ratio (S/N), resulting in a lower resolution DCO⁺ image. We instead use a Keplerian mask to improve the signal-to-noise ratio of the DCO⁺ integrated intensity image, thus retaining the high spatial resolution. We also present extensive modeling of the DCO⁺ emission to explore the sensitivity of the emission to the disk physical conditions. Table 1 summarizes the observational parameters for the DCO⁺ J = 3 − 2 emission in this work.

Fig. 1 Top-left: integrated intensity map of the DCO+ J = 3 − 2 line from 5.4–8.8 km s−1 after applying a Keplerian mask to the image cube. The synthesized beam and AU scale are shown in the lower corners. Top-right: DCO+ J = 3 − 2 integrated intensity map overlaid with white marking the model inner dust cavity, the inner dust ring (R1 , hatched), the dust gap, and the outer dust ring (R2, hatched). The synthesized beam and AU scale are shown in the lower corners. Bottom-left: velocity-weighted coordinate map of the DCO+ J = 3 − 2 line, clipped at 3.5σ. Solid black contours show the 233 GHz/1.3 mm emission at 7.0 × 10−5 Jy beam−1 (1σ) × [5, 50, 200]. The synthesized beam and AU scale are shown in the lower corners. Bottom-right: disk-integrated spectrum of the DCO+ J = 3 − 2 line before Keplerian masking, Hanning smoothed to 0.17 km s−1 channels. The horizontal dashed black line indicates the continuum-subtracted spectral baseline. The vertical red line shows the systemic velocity at 6.9 km s−1 .

3 Results

The DCO⁺ J = 3 − 2 line in the disk around HD 169142 was readily detected and imaged at 0.37′′ × 0.23′′ [43× 27 AU at 117 pc] spatial resolution, with beam PA = –74.8°. The systemic velocity is 6.9 km s⁻¹ (Fedele et al. 2017). The spectrum shown in Fig. 1 was extracted from the original self-calibrated, continuum-subtracted CLEAN image. The right ascension and declination axes of the image cube are collapsed over a circular region with radius 1.75′′ centered on the source position.

To enhance the S/N of the DCO⁺ emission maps and radial profile, a mask in right ascension, declination, and velocity was applied to the original image cube data, following Carney et al. (2017) and Salinas et al. (2017). The mask is based on the velocity profile of a rotating disk, which is assumed to be Keplerian around a central stellar mass of M = 1.65 M_⊙ (Blondel & Djie 2006). A subset of pixels are identified in each velocity channel where the pixel Keplerian velocity matches the Doppler-shifted line velocity. Pixels with velocities that do not match the Keplerian rotational profile criteria are masked. Appendix B shows the DCO⁺ J = 3 − 2 channel maps with the Keplerian mask outline visible as the blue contours. To obtain the integrated line flux for DCO⁺ J = 3 − 2 reported in Table 1, the spectrum was extracted from the image cube after applying the Keplerian mask and integrated over the velocity range 5.4–8.8 km s⁻¹.

Fig. 2 Azimuthally averaged radial intensity profiles of DCO+ (top), C18O (top-middle), the 1.3 mm continuum (bottom-middle), and the J band ( ~1.2 μm) polarized light (bottom). Shaded regions represent 1σ errors on the intensity.

3.1 Radial distribution of DCO⁺

The DCO⁺ emission has a ring-like morphology in this disk, with the majority of emission originating in a region between 0.4′′ –1.4′′ [47–164 AU at 117 pc], based on the velocity-weighted coordinate (first-order moment) map in Fig. 1 obtained from applying a 3.5σ clip to the emission in the Keplerian-masked DCO⁺ image cube. The ring extends significantly beyond the outer edge of the 1.3 mm continuum, similar to the ¹² CO, ¹³ CO, and C¹⁸ O molecular lines (Fedele et al. 2017). Figure 1 also shows the integrated intensity (zero-order moment) map from velocity channels 5.4–8.8 km s⁻¹ . Applying the Keplerian mask improved the S/N of the integrated intensity image by a factor of three, from 5 to 15. The radial profile for DCO⁺ in Fig. 2 is obtained by taking the mean intensity in azimuthally averaged elliptical annuli projected to an inclination i = 13° and position angle PA = 5°. The radial bin size was set to 0.1′′ [11.7 AU at 117 pc]. Errors are calculated as the standard deviation of the pixel intensity contained within each annulus divided by the square root of the number of beams.

With the increased S/N of the averaged annuli, it is clear from the radial intensity profile that DCO⁺ extends out to ~230 AU and peaks at a radius of ~100 AU. Within 100 AU, there is a gap between ~30–60 AU where the intensitydrops, with some DCO⁺ emission returning at radii ≲30 AU. Given the errors on the curve, the actual drop in emission in the 30–60 AU region may be small (see Fig. 2). At ~150 AU there is a knee in the radial profile, and at ~200 AU there is a distinct bump.

The intensity at r ≲ 30 AU already suggests that there is a warm component to the DCO⁺ emission, as temperatures in this region of the disk are too high to allow the cold deuterium fractionation pathway to be active. As seen in Fig. 2, the dip in DCO⁺ intensity from r = 30–60 AU corresponds well to the gap between the two dust rings, indicated by the filled regions. With less dust and gas in the gap, the overall surface density profile falls dramatically, causing a corresponding dip in the DCO⁺ radial profile, more prominently than seen in C¹⁸O.

The DCO⁺ intensity increases significantly within the outer dust ring, which is to be expected if DCO⁺ is forming near the midplane where the dust temperature is sufficiently low for some degree of CO freeze-out. Interestingly, the peak in the radial profile at ~100 AU is beyond the outer edge of the second dust ring, and emission is present throughout the outer disk. This suggests that beyond the 1.3 mm continuum the disk remains cold due to the presence of micron-sized dust grains, as observed by Quanz et al. (2013), Monnier et al. (2017), and Pohl et al. (2017).

In addition, Fig. 2 compares the radial profiles of DCO⁺, C¹⁸ O, the 1.3 mm continuum, and the J band (1.2 μm) polarized intensity. The polarized intensity data is from VLT/SPHERE, with the radial profile obtained after azimuthally averaging the deprojection of the r²-scaled J band Q_ϕ image and normalizing to the maximum brightness of the inner ring (Pohl et al. 2017). The dust rings R1 and R2 are clearly visible in both the millimeter emission from ALMA and the micron emission from VLT/SPHERE. The J band profile shows small grains existing throughout the extent of the gaseous disk out to ~200 AU. The radial profiles show that the DCO⁺ emission is highly sensitive to changes in the disk structure, whereas C¹⁸ O is less affected. The feature at ~200 AU reveals that there is some mechanism in the disk causing more DCO⁺ emission than would be expected for a smoothly decreasing abundance. There may be an accompanying slope change of the C¹⁸ O at ~190 AU, but it is difficult to discern in Fig. 2. A feature at large radii in the C¹⁸ O is more apparent in the radial slice along the major axis of the disk shown in Fig. 3 of Fedele et al. (2017), where a small bump can be seen at ~1.7′′ [200 AU at117 pc], hinting at outer disk structure in C¹⁸ O. At these radii the DCO⁺ is tracing midplane substructure in the disk that is not as apparent in the more abundant, optically thick CO isotopologues.

3.2 Column density and disk-averaged abundance in LTE

We estimated the disk-averaged abundance of the observed DCO⁺ based on the total integrated line flux, an assumed excitation temperature, and the total disk mass. Following the formula used by Remijan et al. (2003) and Miao et al. (1995) for optically thin emission in local thermodynamic equilibrium (LTE), we can estimate thecolumn density as $$N=2.04\frac{{\displaystyle \int I_{\nu }}\text{d}v}{{\theta }_{\text{a}}{\theta }_{\text{b}}}\frac{Q_{\text{rot}}\text{exp}(E_{\text{u}}/T_{\text{ex}})}{{\nu }^3\langle S_{\text{ij}}{\mu }^2\rangle }\times {10}^{20}{\text{cm}}^{-2},$$(1)

where ∫ I_νdv is the integrated line flux in Jy beam⁻¹ km s⁻¹, θ_a and θ_b correspond to the semi-major and semi-minor axes of the synthesized beam in arcseconds, T_ex is the excitation temperaturein K, and ν is the rest frequency of the transition in GHz. The partition function (Q_rot), upper energy level (E_u, in K), and the temperature-independent transition strength and dipole moment (S_ij μ², in debye²) for the DCO⁺ molecule are taken from the Cologne Database for Molecular Spectroscopy (Müller et al. 2005).

DCO⁺ is expected to form primarily in the midplane close to the CO freeze-out temperature, where gas densities are typically higher ( ~10⁹ cm⁻³; Walsh et al. 2014b) than the critical density of the J = 3 − 2 transition at 20–30 K ( ~2 × 10⁶ cm⁻³; Flower 1999). Under these conditions, LTE is a reasonable assumption. Furthermore, the density of the H₂ gas taken from the Fedele et al. (2017) model (see Fig. 4) is greater than the critical density of DCO⁺ for z∕r < 0.3. Therefore, only if DCO⁺ is present solely in the diffuse upper disk layers where the gas and dust temperature have decoupled would LTE be an unreasonable assumption. Currently, formation routes for DCO⁺ place the molecule in significant abundance only in the intermediate disk layers ( ≲100 K) and near the midplane, further justifying the use of LTE.

We explore excitation temperatures of 25, 50, and 75 K, which cover the range of expected DCO⁺ emitting regions (Mathews et al. 2013; Favre et al. 2015). The total integrated line flux and excitation temperature are used to calculate a disk-averaged DCO⁺ column density. Assuming optically thin emission, the disk-averaged column density was then used to estimate the total number of DCO⁺ molecules in the disk, N(DCO⁺) = N_avg × (a × b), where (a × b) is the total emitting area of DCO⁺. Assuming the total disk mass is primarily molecular hydrogen, we can estimate the total number of H₂ molecules, $N({\text{H}}_2)=M_{\text{disk}}/m_{{\text{H}}_2}$, where $m_{{\text{H}}_2}$ is the molecular hydrogen mass. The emitting area is set to a = b = 3′′ based on the diameter of emission in the integrated intensity map, and the total disk mass is 1.9 × 10⁻² M_⊙. Table 2 shows the disk-averaged column density and abundance for T_ex = 25, 50, and 75 K, which are consistent to within a factor of two over the temperature range. DCO⁺ column densities of order 10¹¹−10¹² are similar to the values reported for HD 163296 (Mathews et al. 2013; Salinas et al. 2017), TW Hya (Qi et al. 2008), and DM Tau (Teague et al. 2015).

4 Modeling DCO⁺ emission

The aim of modeling the DCO⁺ emission in this disk is to determine the midplane conditions which create sufficient production of DCO⁺ in the outerdisk, and to estimate the contribution of cold and warm formation routes to the overall DCO⁺ abundance. The initial physical structure of the HD 169142 disk is adopted from Fedele et al. (2017), who constrained the density and temperaturestructure by simultaneously fitting the radial distribution of the 1.3 mm continuum and three CO isotopologues: ¹² CO J = 2 − 1, ¹³ CO J = 2 − 1, and C¹⁸ O J = 2 − 1. The disk structure is then optimized to include a small grain dust population throughout the disk that was absent in the original model. With the optimized disk structure, we then reproduce the DCO⁺ radial intensity profile in a parameterized way with a simple deuterium chemical network.

4.1 Fiducial physical structure

We use the thermo-chemical code Dust And LInes (DALI; Bruderer et al. 2012; Bruderer 2013) to obtain the physical disk structure. Input for DALI consists of a blackbody radiation field with T_eff = 8400 K to estimate the stellar photosphere and a power-law gas surface density with an exponential drop-off $${\Sigma }_{\text{gas}}={\Sigma }_{\text{c}}{\left(\frac{R}{R_{\text{c}}}\right)}^{-\gamma }\left[\text{exp}-{\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }\right],$$(2)

where R_c (100 AU) is the critical radius, Σ_c (6.5 g cm⁻²) is the value of the gas surface density at the critical radius and γ (1.0) is the power-law exponent. The initial dust surface density is extrapolated from the gas surface density by assuming a gas-to-dust ratio (Δ_gd = 80) such that Σ_dust = Σ_gas∕Δ_gd. The vertical gas density is described by a Gaussian distribution with a scale height $h=h_{\text{c}}{(R/R_{\text{c}})}^{\psi }$ that depends on the disk radius and the flaring exponent ψ (0.0) with a critical scale height, h_c (0.07), defined at the critical radius.

Dust settling is approximated in DALI by considering two different populations of dust grains following the power-law description from D’Alessio et al. (2006), with a power-law exponent p = 3.5. The small grains (0.005–1 μm) have a scale height h while large grains (0.005–1000 μm) have a scale height hχ, where the settling parameter χ (0.2) is in the range 0–1. The fractional distribution between the two populations of dust grains is set by the parameter f_large (0.85), which results in dust surface densities of Σ_dustf_large for large grains and Σ_dust(1 − f_large) for the small grains.

The DALI code solves for dust temperatures and radiation field strength in each grid cell using two-dimensional radiative transfer, then determines the heating-cooling balance of the gas, molecular excitation, and chemical abundances based on an input chemical network. The DALI model described in this work uses the ISO chemical network, which includes CO freeze-out and CO isotope-selective photodissociation (Miotello et al. 2014, 2016).

Fedele et al. (2017) modified the surface density profile from our Eq. (2) to include gas depletion in the inner disk and two millimeter dust rings: R₁ from 20–35 AU and R₂ from 56–83 AU. These modifications result in a radially variable gas-to-dust ratio throughout the disk, with Δ_gd = 80 valid only for the R₂ outer dust ring from 56–83 AU. A full description and table of values including the parameter ranges and best-fit parameters of their fiducial DALI model can be found in Sect. 5.2, Fig. 5, and Table 2 of their paper. Their fiducial model only included a small grain dust population within the millimeter rings since they fit only the millimeter emission. We expand on their fiducial model by including small grains in other regions of the disk and optimize the model parameters to keep the fit to the spectral energy distribution, 1.3 mm emission, and CO isotopologues. See Appendix A for details.

Fig. 3 One-dimensional, radial structure in the HD 169142 disk. Vertical brown shaded regions represent the dust rings. Left: DCO+ radial column density (dotted blue) calculated from Eq. (1) using the DCO+ the radial intensity profile (see Fig. 2) and the midplane gas temperature of the optimized DALI model (dot-dash red) as the Tex profile. The gas surface density profile (dashed green) was used to derive the DCO+ radial abundance. Right: radial abundance structure of DCO+ with respect to H2. The blue shaded region represents 1σ errors on the abundance.

4.2 Vertically averaged radial abundance profile in LTE

A radial abundance profile for the observed DCO⁺ emission can be calculated using the method outlined in Sect. 3.2. Rather than estimate a global excitation temperature for DCO⁺, with the model physical structure outlined in the previous section we can now obtain a radial description of the excitation temperature. To estimate T_ex , which will be equivalent to the kinetic gas temperature assuming LTE, the midplane gas temperature in each radial bin was taken from the optimized DALI model. The line intensity was extracted from the integrated intensity map in radial bins as in Sect. 3.1, with a bin size of 0.1′′ [11.7 AU]. Equation 1 was then used to determine the radial column density of DCO⁺. Assuming the gas is composed primarily of H₂, the H₂ mass can beused to convert the gas surface density profile of the DALI model into a gas column density profile. Dividing the DCO⁺ radial column density by the gas radial column density gives a radial abundance profile for DCO⁺.

Figure 3 shows the DCO⁺ radial abundance. The profile has an inner radius of 13 AU and an outer radius at 180 AU, corresponding to the inner and outer radii of the model gas surface density. Beyond r = 56 AU, the abundance increases with radius by a factor of about 5 with values ranging from 1− 5× 10⁻¹², which are comparable values to the DCO⁺ radial abundance estimates for HD 163296 (Salinas et al. 2017). Similar trends of DCO⁺ abundance increasing with radius have been observed in TW Hya and DM Tau (Qi et al. 2008; Teague et al. 2015). The sharp increase in abundance at r = 56 AU is due to the δ_{gas, gap} = 0.025 depletion factor in the surface density of the gas for r < R_{gap out} (56 AU). Within the errors, the abundance profile remains relatively flat at radii less than R_{gap out}. Due to the proximity to the central star, DCO⁺ in this region is likely formed via the warm deuterium fractionation pathway.

4.3 Parameterized models

We move from the one-dimensional derivation of the DCO⁺ radial profiles to a two-dimensional DCO⁺ structure to explore variations in abundance for different radii and heights in the disk. For this we parameterize the HD 169142 physical disk structure obtained from the optimized DALI model, shown in Fig. 4. The DALI code calculates the dust temperature, the local radiation field, heating and cooling rates, molecular abundances, and gas temperature self-consistently, making it difficult to isolate and explore individual parameters that may affect DCO⁺ emission. We examine the effect of alterations to the disk gas temperature and CO abundance on DCO⁺ production by employing a simple, parameterized modeling technique using the steady-state, analytic chemical code (hereafter DCO+ CHEMNET) from Murillo et al. (2015). This time-independent chemical model is preferred because the chemical timescales for gas-phase reactions are sufficiently fast that a steady state is achieved at times much shorter than the estimated lifetime of the HD 169142 disk ( ~10 Myr). Murillo et al. (2015) have already shown that the DCO+ CHEMNET code reproduces the trends of full, time-dependent chemical models for protostellar envelopes. We extend this treatment to the protoplanetary disk environment. Section 4.1, Table 1, and subsect. 4.2.3 in their paper give a full description of the chemical network, the simplified set of DCO⁺ formation and destruction reactions via the cold deuterium fractionation pathway, and a comparison to a full chemical network.

The DCO+ CHEMNET code takes as input the gas density, gas temperature, and CO abundance structure. The input profiles are taken from the output of the optimized DALI model described in Appendix A. Because DCO+ CHEMNET is a single-point code, each grid cell in (r, z) of the DALI disk grid was run separately to capture the full two-dimensional model structure. The gas density of the parameterized models is assumed to be equal to the molecular hydrogen density, $n_{{\text{H}}_2}$, and the HD abundance is constant at X_HD = n_HD∕n_H = 10⁻⁵ throughout the disk. The cosmic ray ionization rate of H₂ is set to ζ_cr = 1.26 × 10⁻¹⁷ s⁻¹ .

The parameterized models are used to probe the conditions near the midplane of the outer disk, a region to which previously observed CO isotopologue tracers are not sensitive. The CO observations do not probe all the way down to the midplane in the outer disk and only a small amount of CO is absent from the total column density, while DCO⁺ production is highly sensitive to the midplane conditions (Mathews et al. 2013; Qi et al. 2015; Favre et al. 2015; Huang et al. 2017). The CO lines become optically thick beyond R_{dust out} (83 AU) near the midplane of the optimized DALI model (z∕r < 0.2 and < 0.1 for ¹²CO and C¹⁸O, respectively) so the physical structure is not constrained in this region. Therefore, in all parameterized models, only the region r > 83 AU and z∕r < 0.1 was altered since this is the τ_mm = 1 surface of C¹⁸O J = 2 − 1 from the DALI model. The region z∕r < 0.1 quickly becomes highly optically thick for C¹⁸O J = 2 − 1 with optical depths of τ_mm = 10 near the midplane. For z∕r > 0.1, the C¹⁸O abundance of the model would be sensitive to changes in the gas temperature. Thus, where z∕r < 0.1, we can alter the disk structure without affecting the intensity profiles of the CO isotopologues.

Five parameterized model scenarios are tested: a disk with high midplane CO abundance, a disk with low midplane CO abundance, a cold disk, a shadowed cold disk, and a shadowed cold disk including CO depletion. The aim is to initially test the CO abundance and the gas temperature separately to investigate which parameter has a stronger influence on the production of DCO⁺. These two parameters are intrinsically interlinked (i.e., more CO freeze-out will occur in lower temperature environments), and the parameterized models allow us to explore these effects in isolation. In the parameterized models we also include an additional DCO⁺ constant abundance region to act as a proxy for warm DCO⁺ formation.

Fits to the radial intensity profile of the data are used to evaluate the model parameters. To obtain model intensity profiles, synthetic DCO⁺ image cubes are created using the two-dimensional gas density, gas temperature, and DCO⁺ abundance structure from DCO+ CHEMNET as input to the LIne Modeling Engine (LIME; Brinch & Hogerheijde 2010) radiative transfer code. The LIME code was run in LTE with 30 000 grid points to create synthetic images of the DCO⁺ J = 3 − 2 transition. The images are continuum-subtracted and sampled in the uv plane using the PYTHON vis_sample¹ routine, which reads the uv coordinates directly from our observed ALMA measurement set and creates synthetic visibilities of the model. The model visibilities are imaged in CASA using CLEAN with natural weighting, and an integrated intensity map is created over the same velocity range as the data (5.4–8.8 km s⁻¹). Azimuthally averaged elliptical annuli projected to the disk inclination and position angle are used to extract the integrated intensity of the model with the same radial bins as the data.

Fig. 4 Two-dimensional physical structure of the HD 169142 disk from the optimized DALI model ( Appendix A). Top: gas density structure. The gas density contours (dashed black) are shown as log(ngas). Middle: gas temperature structure below 100 K. Temperature contours are shown in dashed blue. Bottom: CO abundance structure with respect to H2.

4.3.1 DCO⁺ from deuterium fractionation

Formation of DCO⁺ occurs as a result of deuterium fractionation, which is an enhancement in the D/H ratio observed in certain deuterium-bearing molecules. Typically, deuterium fractionation occurs in colder environments such as pre-stellar cores and below the surface layers of protoplanetary disks because of the lower zero-point energies of the deuterated molecular ions (Brown & Rice 1986; Millar et al. 1989). Deuterium fractionation in the low-temperature regime occurs via the reaction (Wootten 1987) $$\text{HD}\text{+}{\text{H}}_{\text{3}}^{\text{+}}\leftrightarrow {\text{H}}_{\text{2}}\text{+}{\text{H}}_{\text{2}}{\text{D}}^{\text{+}}\text{+}\Delta E,$$(3)

where ΔE = 220 K (Roberts & Millar 2000; Gerlich et al. 2002; Albertsson et al. 2013). This deuterium fractionation pathway is typically efficient at temperatures below ~30 K due to the energy barrier ΔE for the back-reaction. The DCO⁺ molecule is then formed by the following reaction $${\text{H}}_{\text{2}}{\text{D}}^{\text{+}}\text{+}\text{CO}\leftrightarrow {\text{H}}_{\text{2}}\text{+}{\text{DCO}}^{\text{+}}.$$(4)

Gas-phase CO is needed to produce DCO⁺, but CO will also rapidly combine with H${}_3^{+}$, quenching the production of H₂ D⁺. Thus there is a balance where CO must be sufficiently depleted for H₂ D⁺ to remain abundant yet enough gas-phase CO must be present so that DCO⁺ may form. The simple chemical network from Murillo et al. (2015) is used in this work to model DCO⁺ produced via the cold deuterium fractionation pathway. The model shown in Fig. 4 suggests that DCO⁺ is expected throughout much of the disk midplane, as temperatures are cold enough for CO freeze-out and the production of H₂ D⁺.

Deuterium fractionation in disks can also occur when HD and CH${}_3^{+}$ combine to create H₂ and CH₂D⁺ (Millar et al. 1989). The molecule DCO⁺ is then formed via CH₂D⁺ reactions directly or via one of its products, CH₄D⁺. The energy barrier for the back-reaction of this deuterium fractionation pathway was recently revised (Δ E = 654 K; Roueff et al. 2013) and suggests that CH₂D⁺, and therefore DCO⁺, could be formed efficiently at higher temperatures. In recent models of the TW Hya disk with a full deuterium chemical network, Favre et al. (2015) find that T ≥ 71 K is sufficient to switch off the production of DCO⁺ formed via CH₂D⁺. The warmer CH₂D⁺ fractionation route dominates over the cold H₂ D⁺ fractionation route for temperatures greater than ~30 K because of the higher energy barrier and the fact that H₂ will readily destroy H₂ D⁺ above ~30 K.

4.3.2 CO abundance versus gas temperature

Table 3 shows the parameters used for each model scenario. The high CO and low CO models respectively increase and decrease the CO abundance in the region of interest (r > 83 AU and z∕r < 0.1) in order to determine how changes in the availability of gas-phase CO influence the DCO⁺ emission. The cold disk model tests the influence of the gas temperature on the DCO⁺ emission by moderately decreasing the temperature profile in the region of interest. The shadowed cold disk model expands on the cold disk model by creating a secondary colder region just outside of the edge of the millimeter grains.

In this work we do not include a chemical network describing the formation of DCO⁺ via the warm deuterium fractionation pathway. The complexity of hydrocarbon cation chemistry on which the warm deuterium fractionation pathway depends introduces large uncertainties in the results of even basic chemical networks. As a first-order approximation, we instead adopt a region of constant DCO⁺ abundance for 30 K  ≤ T ≤ 70 K, with the lower temperature limit based on the energy barrier of the back reaction for Eq. (3) and the upper temperature limit based on the results of Favre et al. (2015), where they find that formation of DCO⁺ via the warm deuterium fractionation pathway is effectively switched off at 71 K. In this way we introduce a proxy in the model that is representative of DCO⁺ production in the high-temperature regime, with the reasonable expectation that the warmer pathway will contribute little to the emission of the outer disk (Öberg et al. 2015), which is the focus of this work. Further investigation on the detailed contribution of the warm deuterium fractionation pathway to the overall DCO⁺ production, with particular attention to the inner regions of the disk at radii ≲50 AU, will be the focus of future work.

The constant abundance component from 30–70 K is tuned so that the model radial profile matches the intensity of the observed DCO⁺ emission from 40–70 AU, where warm DCO⁺ is the primary contributor before peaking and turning over at 70–80 AU. This gives $X({\text{DCO}}_{\text{warm}}^{+})$ = 2.0 × 10⁻¹², which is consistent with the DCO⁺ abundance found by Willacy & Woods (2009) in models of a protoplanetary inner disk at radii less than 30 AU, and consistent with the abundance between 30–70 K in more recent work by Öberg et al. (2015) in their model of IM Lup. It is roughly one to two orders of magnitude lower than the abundance found by Favre et al. (2015), but their warm DCO⁺ was confined to a thin layer spanning only 1–2 AU at radii less than 60 AU. In HD 169142 the 30–70 K layer spans roughly 5–15 AU, depending on the radius, and results in a DCO⁺ column density on the order of 10¹² cm⁻², which is consistent with Favre et al. (2015). With this treatment, the warm deuterium fractionation pathway contributes <20% to the DCO⁺ radial intensity profile for r > 83 AU in all models considered. For r < 83 AU, the warm component is the primary contributor to the overall DCO⁺, producing >80–95% of the emission, depending on the model. For r ≲ 30 AU, no warm component of DCO⁺ exists in the model because the midplane gas temperature reaches 70 K at ~30 AU. Because the warm component is set to the same abundance value for all parameterized models, we consider the outer disk at r > 83 AU for comparison between the data and the model.

Figure 5 shows the gas temperature map and ¹² CO abundance map used as input, the DCO⁺ abundance map calculated by DCO+ CHEMNET and including the constant abundance warm component, and the DCO⁺ radial intensity profiles derived from the LIME synthetic images with only the cold deuterium fractionation pathway active (C only) and with the constant abundance warm component included (C + W). The no modifications model in the first row of Fig. 5 uses the input from our optimized DALI model with no alterations to the midplane disk structure. It is already clear that the DCO⁺ emission is under-produced by the DCO+ CHEMNET code for the disk structure outlined in Appendix A.

In the low CO abundance case, an order of magnitude decrease in midplane gas-phase CO abundance results in a factor of two to three decrease in DCO⁺ emission because there is no longer enough gas-phase CO available to efficiently form DCO⁺. The high CO abundance model with an one order of magnitude increase in gas-phase CO abundance results in only a factor of two increase in DCO⁺ emission, highlighting the non-linear nature of the chemistry. CO abundances in this model become much higher than the canonical value of ~10⁻⁴ with respectto H₂ in the midplane of disk where CO depletion is expected; hence, this model is physically unrealistic. Order of magnitude variations to the CO abundance do not have a significant influence on the formation of DCO⁺ near the outer disk midplane for the physical structure given by the optimized DALI model.

The second and third rows of Fig. 5 illustrate increasingly colder midplane scenarios. In the cold disk model case, a decrease of 2 K in the gas temperature provides an improved fit to the DCO⁺ radial profile but still fails to fully recover much of the DCO⁺ between 80–140 AU. It matches the observed intensity profile well beyond 150 AU. The shadowed cold disk model expands on the cold disk model to invoke a secondary cold region from r = 83−120 AU with T = T_gas − 8 K. For the secondary cold region, the outer boundary and temperature drop parameters are explored from 100–140 AU in steps of 10 AU and from T_gas − 3 K to T_gas − 10 K in steps of 1 K, respectively. The outer boundary at 120 AU and the temperature drop of T_gas − 8 K was found to provide the best improvement on the fit to the DCO⁺ profile. More DCO⁺ is produced in the shadowed cold disk model between 80–150 AU than in the cold disk model case, but it is still not enough to capture the 100 AU peak. The final model invokes CO depletion in the secondary cold region of the shadowed cold disk model, which would be expected for significantly colder disk regions. The CO abundance is reduced by a factor of five between r = 83−120 AU. CO depletion factors of two, five, and ten are tested, with a factor of five resulting in the right amount of DCO⁺ emission. In this scenario, the DCO⁺ radial intensity profile, including the 100 AU peak, is reproduced well.

The efficiency of the cold deuterium fractionation pathway, and thus the production of DCO⁺, is also affected by cosmic ray ionization, ζ_cr, and the ortho-to-para (o/p) ratio of H₂. Ionization of H₂ in the disk by cosmic rays will affect the number of H${}_3^{+}$ ions and the number of free electrons. The DCO+ CHEMNET code was rerun for the model with no modifications to test changes to the cosmic ray ionization rate, initially set to ζ_cr = 1.26 × 10⁻¹⁷ s⁻¹ . There is evidence that the local ISM in the HD 169142 region may reach values of ζ_cr = 1−5 × 10⁻¹⁶ s⁻¹ (Indriolo et al. 2007; Neufeld & Wolfire 2017). An increase in ζ_cr of one order of magnitude results in DCO⁺ emission comparable to the cold disk model, meaning that the amount of DCO⁺ observed at r> 120 AU may be a consequence of a higher local cosmic ray ionization rate or may be due to a moderate decrease in the gas temperature of the outer disk. Ionization of the outer disk by UV radiation may influence the cold deuterium fractionation pathway and free electron population differently, but the effect is not modeled here. Such ionization would influence primarily the disk upper layers, as UV photons will not reach the disk midplane with a sufficient flux to outpace ionization due to cosmic rays.

The ortho-to-para (o/p) ratio of H₂ influences the survival of H₂ D⁺ in cold regions of the disk because the back-reaction in Eq. (3) has a lower energy barrier for o-H₂ (Δ E = 61 K) than for p-H₂ (Δ E = 232 K; Walmsley et al. 2004). We apply a thermal treatment of the o/p ratio of H₂ in LTE following Murillo et al. (2015) with o/p = 9 ×exp(−170∕T), where T is the gas temperature, and again rerun DCO+ CHEMNET code for the no modification model including the H₂ o/p ratio, which was not considered in the original models. In this case, too much DCO⁺ emission is produced in the outer disk for r > 150 AU, but DCO⁺ is under-produced for r < 150 AU. However, this requires an efficient equilibration of the H₂ spin temperature with the gas temperature, rather than the 3:1 obtained from grain surface formation of H₂ . The degree of cosmic ray ionization and the precise distribution of the o/p ratio of H₂ will affect the exact values for the temperature drop and CO depletion necessary to obtain a fit to the DCO⁺ radial intensity profile, but our overall conclusions remain unchanged.

Based on the best-fit shadowed and depleted cold disk model, a significantly colder outerdisk midplane with increased CO freeze-out is the most likely scenario for DCO⁺ production beyond the millimeter grains. The DCO⁺ emission at r > 120 AU could be reproduced with a small drop in the gas temperature or an increase in the cosmic ray ionization rate. The modeling efforts presented here show that DCO⁺ emission can reveal structure in low-temperature regions in the midplane of disks that are not apparent in the CO isotopologues.

5 Discussion

The interpretation of DCO⁺ emission regarding its chemical origins and location within the disk is complicated by the multiple deuterium fractionation pathways available. The two main pathways (cold, via H₂ D⁺ fractionation and warm, via CH₂D⁺ fractionation) are efficient over different temperature ranges, therefore it is useful to consider DCO⁺ emission from distinct regions of the disk where the conditions are expected to be more favorable for one fractionation pathway over the other. The inner disk provides warmer temperatures that can switch off the H₂ D⁺ pathway, while the outer disk hosts a cold midplane that allows the H₂ D⁺ pathway to operate efficiently.

5.1 Inner disk DCO⁺

In models of TW Hya including a full deuterium chemical network, Favre et al. (2015) found that DCO⁺ observed in the inner tens of AU is not primarily formed by the cold deuterium fractionation pathway via H₂ D⁺ because of the warm temperatures of the inner disk. The physical structure in Fig. 4 shows that the disk around HD 169142 is far too warm in the inner 50 AU for the H₂ D⁺ fractionation pathway to be the main contributor. Instead, DCO⁺ in this region is likely formed by the warmer CH₂D⁺ fractionation pathway. The disk temperature at the midplane is greater than 70 K at r < 30 AU, therefore even the warm component is switched off. In this disk, DCO⁺ formed via the CH₂D⁺ fractionation pathway may continue to be active at temperatures greater than 70 K. Alternatively, we may be missing a cool inner component, such as inner disk dust rings that could keep temperatures low enough for the warm deuteration fractionation pathway to remain active. Andrews et al. (2016) observed optically thick millimeter dust rings on the order of a few AU in the disk around TW Hya. Recent work by Ligi et al. (2018) presented a tentative detection of another dust ring in the HD 169142 disk located at ~0.1′′ [12 AU at 117 pc] using VLT/SPHERE radial differential imaging with the Infrared Dual-band Imager and Spectrograph and the Integral Field Spectrograph instruments.

The molecule DCN provides another avenue to probe the warm component of DCO⁺ emission as it is also formed via warm deuterium fractionation in disks (Millar et al. 1989). Co-spatial peaks in DCO⁺ and DCN would indicate that the warm deuterium fractionation pathway is a strong contributor to the production of DCO⁺. Recent observations of DCN and DCO⁺ in several T Tauri and Herbig Ae/Be sources show them peaking in different regions of the disk, but with some DCO⁺ present where the DCN peaks, indicating that the warm deuterium fractionation pathway contributed partially to the DCO⁺ emission(Qi et al. 2008; Öberg et al. 2012; Huang et al. 2017; Salinas et al. 2017).

Fig. 5 Parameterized modifications to the optimized DALI model described in Appendix A. Right and right-middleplots: gas temperature and CO abundance, respectively. The black dashed box shows the modified region. Left-middle plot: DCO+ abundance map calculated with DCO+ CHEMNET with the 30–70 K constant abundance warm component included. Left plot: DCO+ radial intensity profiles of the data (solid black with shaded gray 1σ errors) and of the model with the cold component only (C only; dashed blue) and with warm component included (C + W; dashed red). Top: optimized DALI model with no modifications made to the gas temperature or CO profile. Middle: cold disk model. The gas temperature is decreased by 2 K in the black dashed region. Bottom: shadowed and depleted cold disk model. This is an extension of the cold disk model to simulate a secondary colder region beyond the millimeter dust edge. The gas temperature is decreased by 8 K with CO depleted by a factor of five from r = 83−120 AU. The gas temperature is decreased by 2 K for r > 120 AU.

5.2 Outer disk DCO⁺

In order to recover the observed DCO⁺ radial intensity profile, it was necessary to modify the structure of the optimized DALI model from Fig. 4 to include a much colder, CO depleted region beyond the edge of the millimeter grains. The cause of the decrease in temperature could be due to the number and location of micron-sized grains in the outer disk, which have been observed out to at least 200 AU in scattered light imaging (Quanz et al. 2013; Momose et al. 2015; Monnier et al. 2017; Pohl et al. 2017). The exact, local distribution of the micron-sized grains in the HD 169142 disk might be different than the distribution approximated in the optimized DALI model, leading to lower gas temperatures in the outer disk.

The fact that the peak in the DCO⁺ radial intensity profile occurs at 100 AU is already evidence that the colder H₂ D⁺ deuterium fractionation pathway is responsible for the majority of DCO⁺ production in this disk. As mentioned in Sect. 5.1, the warm component is not significant beyond the edge of the millimeter grains (r = 83 AU), with DCO⁺ formed via the cold network in DCO+ CHEMNET contributing >80% to the radial intensity profile. In this context, it is probable that the 100 AU peak observed in the DCO⁺ profile is tracing the CO snow line as CO begins to freeze out onto dust grains, but is still sufficiently abundant in the gas-phase to allow for the formation of DCO⁺. Macías et al. (2017) previously compared the C¹⁸O and DCO⁺ radial intensity profiles and similarly concluded that a co-incident slope change in the C¹⁸ O and peak in the DCO⁺ is evidence that DCO⁺ is tracing the CO snow line, which they locate at ~100 AU. The modeling done in this work shows that invoking a cold, CO depleted region in the outer disk midplane from ~80–120 AU is necessary to recover the DCO⁺ radial intensity profile in the parameterized models, indicating that the bulk of the DCO⁺ does indeed trace regions of CO freeze-out in the HD 169142 disk.

There is tentative evidence for an outer DCO⁺ ring, seen as a bump in the DCO⁺ radial profile at the edge of the gaseous disk around 200 AU. Similar DCO⁺ structure has been seen in other disks, e.g., IM Lup, LkCa 15, and HD 163296 (Öberg et al. 2015; Huang et al. 2017; Salinas et al. 2017; Flaherty et al. 2017). An outer ring could indicate that some CO ice is returned to the gas phase in the outer disk via thermal desorption caused by a temperature inversion or via photodesorption by UV radiation (Huang et al. 2016; Cleeves 2016; Cleeves et al. 2016; Facchini et al. 2017). The observed bump is weaker than the outer DCO⁺ ring in other disks, suggesting that the HD 169142 disk remains cold out to large radii or that the degree of CO freeze-out in this disk may not be as extreme as in other sources.

Follow-up observations of N₂ H⁺ in this diskcould provide further insight into the relationship between DCO⁺ and CO freeze-out. Formation of N₂ H⁺ requires low temperatures and significant CO depletion, and therefore is expected to be abundant only in the cold midplane where there is a high degree of CO frozen out onto the surface of dust grains (Walsh et al. 2012; Qi et al. 2015; van’t Hoff et al. 2017). Co-spatial emission of N₂ H⁺ and DCO⁺ in the outer disk would place strong constraints on the extent to which DCO⁺ is directly tracing CO freeze-out.

6 Conclusions

In this paper we present ~0.3′′ resolution ALMA observations of DCO⁺ J = 3 − 2 in the protoplanetary disk around HD 169142. We update the fiducial DALI model from Fedele et al. (2017) to include small dust grains throughout the disk and employ a simple deuterium chemical network to investigate the production of DCO⁺ formed by the cold deuterium fractionation pathway via H₂ D⁺. The CO abundance and gas temperature structure of the DALI disk model is adapted in a parameterized way using the DCO⁺ CHEMNET code to recover the observed DCO⁺ radial intensity profile. The warm deuterium fractionation pathway via CH₂D⁺ is approximated with a constant abundance between 30–70 K in the parameterized models. The conclusions of this work are the following:

• DCO⁺ has a broad, ring-like morphology over radii ~50–230 AU in the disk around HD 169142, with most of the emission located outside of the millimeter continuum edge. There is an inner component to the DCO⁺ radial profile with emission returning in the central ~30 AU. The DCO⁺ radial intensity profile peaks at ~100 AU, witha tentative secondary bump at ~200 AU.

• Parameterized modeling of the HD 169142 disk shows that lowering the gas temperature in optically thick regions near the midplane of the outer disk by several Kelvin has a significant effect on the DCO⁺ profile. Order of magnitude changes to the CO abundance in the same region cause the DCO⁺ radial profile to increase or decrease only by a factor of 2–3, still under-producing the amount of observed DCO⁺. It was necessary to invoke both effects to successfully reproduce the full radial intensity profile.

• The best-fit parameterized model has a shadowed, cold region with CO depletion near the disk midplane just beyond the edge of the millimeter grains. For z∕r < 0.1 and r > 83 AU, the model recovers the radial intensity profile of DCO⁺ with a T = T_gas − 8 K region with a factor of five CO depletion from r = 83−120 AU, and a T = T_gas − 2 K region for r > 120 AU. The exact values for the drop in T_gas and CO abundance will also depend on the ionization degree and the ortho-to-para ratio of H₂. The fact that the added shadowed region is needed to recover the 100 AU radial intensity peak highlights the sensitivity of DCO⁺ to small changes in the gas temperature and CO abundance structure.

• The best-fit model suggests that the contribution to the overall DCO⁺ emission from the cold deuterium fractionation pathway via H₂D⁺ is >85% at r > 83 AU, while the contribution from the warm deuterium fractionation pathway via CH₂D⁺ is >80% at radii less than r < 83 AU using a constant abundance of $X({\text{DCO}}_{\text{warm}}^{+})=2.0\times {10}^{-12}$ from 30–70 K. The warm component does not recover the return of DCO⁺ emission within ~30 AU for the current disk structure.

• DCO⁺ is an optically thin molecular tracer that acts as a filter to detect disk substructure that is not observable in more abundant, and therefore more easily detectable, molecular tracers such as ¹²CO, ¹³CO, and C¹⁸O. In this work, DCO⁺ observations reveal the low-temperature ( ≲25 K) midplane structure in the disk around HD 169142.

This work shows that DCO⁺ can be used as a valuable tracer of protoplanetary disk midplane conditions with simple models and chemical networks. To fully characterize the complex chemical origins of DCO⁺ in the disk around HD 169142, full dust evolution models (e.g., Facchini et al. 2017) and an expanded deuterium chemical network such as the one used by Favre et al. (2015) would be required. Further constraints could be placed on the origins of DCO⁺ with future ALMA observations of additional chemical tracers, such as DCN and N₂ D⁺, that would correlate with DCO⁺ formed via the warm and cold deuterium fractionation pathways, respectively.

Adjustments to the dust distribution will affect the SED, and the accompanying opacity and temperature variations will affect the molecular abundances and thus radial profiles of ¹² CO, ¹³ CO, and C¹⁸ O. The new parameters for the small grain (0.005–1 μm) dust population were chosen such that the fit to the SED was maintained and the fits to the CO isotopologues were consistent to within a factor of about 10%.

The PAH structure was also modified. Originally, the fiducial model dust temperature was set by PAH thermal emission in regions with no dust grains, with a global PAH abundance set equal to the interstellar medium (ISM) abundance. PAH abundances have been observed to be low in disks with respect to the ISM (Li & Lunine 2003; Geers et al. 2006; Thi et al. 2014). However, to achieve the required opacity for a reasonble fit to the ¹² CO radial profile for r < 25 AU, it was necessary to keep a high PAH abundance in the inner regions. The PAH abundance in our optimized model was therefore set equal to the ISM abundance for r ≤ 83 AU and set to 1% with respect to the ISM abundance for r > 83 AU, which is comparable to current estimates for the HD 169142 disk (Seok & Li 2016).

Figure A.2 shows the fit to the SED and radial profile of the CO isotopologues from the optimized DALI model structure. Each point in the SED is well-fit by the DALI model with the small grain population modifications. The PAH feature is slightly overproduced because of the high abundance in the model at r ≤ 83 AU. Fits to the ¹²CO, ¹³ CO, and C¹⁸O radial intensity profiles remain consistent with the previous fits from Fedele et al. (2017) to within 10%.

Acknowledgements

The authors thank the anonymous referee for useful comments that helped to improve the paper. M.T.C. thanks A. Pohl for providing the VLT/SPHERE J band polarized intensity data. M.T.C. and M.R.H. acknowledge support from the Netherlands Organisation for Scientific Research (NWO) grant 614.001.352. C.W. acknowledges the NWO (grant 639.041.335) and the University of Leeds for financial support. D.F. and C.F. acknowledge support from the Italian Ministry of Education, Universities and Research, project SIR (RBSI14ZRH). This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.00592.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. A.M. acknowledges an ESO Fellowship.

Appendix A Model details

Modifications to the fiducial DALI model from Fedele et al. (2017) were necessary to more accurately incorporate the (sub)micron-sized grain population into the HD 169142 disk model and to obtain a better description of the dust and gas temperatureof the outer disk beyond the millimeter continuum edge. Here the modifications are discussed in more detail.

The Fedele et al. (2017) model fit the 1.3 mm dust emission, and dust in their models exists only within dust rings from r = 20–35 AU and from r = 56–83 AU. There is no dust of any kind in their models outside of these radii. Observations by Quanz et al. (2013), Monnier et al. (2017), and Pohl et al. (2017) have revealed the presence of a micron-sized dust population in the millimeter ring gap and in the outer disk beyond the millimeter-sized dust population. Monnier et al. (2017) discussed an outer disk in the context of micron-sized grains observed in the J band with the Gemini Planet Imager. The authors detected a dip in polarized light at ~55 AU which they model as a 40–70 AU gap, and they observed a flared outer ring that peaks at ~75 AU. The authors adopted a distance of 145 pc for HD 169142. Their radii scaled to the new Gaia distance of 117 pc gives a 32–56 AU gap and a 60 AU outer ring peak in the J band, consistent with the values for the millimeter gap and outer ring observed by Fedele et al. (2017).

Fig. A.1 Top: dust structure of the fiducial DALI model from Fedele et al. (2017) with small (0.005–1 μm) and large (0.005–1000 μm) grains present only in the millimeter rings. Bottom: dust structure of the optimized DALI model with small grains in the millimeter gap and in the outer disk. The hot inner dust ring is not visible due to the extremely small scale height at r = 0.2 AU.

The original dust structure of the Fedele et al. (2017) model was $$n_{\text{dust}}=\{\begin{array}{ll}0\hfill & \text{for}r<R_{\text{dust}\text{in}}\hfill \\ {\delta }_{\text{dust}}\times n_{\text{dust}}\hfill & \text{for}{R}_{\text{dust}\text{in}}<r<R_{\text{gap}\text{in}}\hfill \\ 0\hfill & \text{for}{R}_{\text{gap}\text{in}}<r<R_{\text{gap}\text{out}}\hfill \\ n_{\text{dust}}\hfill & \text{for}{R}_{\text{gap}\text{out}}<r<R_{\text{dust}\text{out}}\hfill \\ 0\hfill & \text{for}r>R_{\text{dust}\text{out}}.\hfill \end{array}$$(A.1)

However, the gap is not empty and the micron-sized grains are detected throughout the outer disk, thus their fiducial DALI model was neglecting the presence of micron-sized grains that would insulate the outer disk midplane and thus lower dust and gas temperatures.

We modified the dust structure of the Fedele et al. (2017) fiducial model to expand the small grain (n_dust,small; 0.005–1 μm) dust population to be present within the gap between the dust rings (r = 35–56 AU) and at radii beyond the dust rings (r > 83 AU). In order to fit the spectral energy distribution (SED), we included a thin inner region of hot dust at r = 0.2 AU, as observed by Wagner et al. (2015). The optimized dust structure of the model is $$n_{\text{dust}}=\{\begin{array}{ll}{10}^{-5}\times n_{\text{dust},\text{small}}\hfill & \text{for}0.2<r<0.4\text{AU}\hfill \\ 0\text{otherwise}\hfill & \text{for}r<R_{\text{dust}\text{in}}\hfill \\ {\delta }_{\text{dust}}\times n_{\text{dust}}\hfill & \text{for}{R}_{\text{dust}\text{in}}<r<R_{\text{gap}\text{in}}\hfill \\ {10}^{-2}\times n_{\text{dust},\text{small}}\hfill & \text{for}{R}_{\text{gap}\text{in}}<r<R_{\text{gap}\text{out}}\hfill \\ n_{\text{dust}}\hfill & \text{for}{R}_{\text{gap}\text{out}}<r<R_{\text{dust}\text{out}}\hfill \\ n_{\text{dust},\text{small}}\hfill & \text{for}r>R_{\text{dust}\text{out}}.\hfill \end{array}$$(A.2)

Fig. A.2 Fits to observables for the optimized DALI model. Unlike the azimuthally averaged radial profiles in Fig. 2, the radial profiles shown here are cuts along the disk semi-major axis (PA = 5°). Top left: near-infrared to submillimeter spectral energy distribution: the model is shown by the black solid curve and the data (red dots) are from Malfait et al. 1998, IRAS (Helou & Walker 1988), Herschel (Pascual et al. 2016), SCUBA (Sandell et al. 2011), and ALMA (Fedele et al. 2017). Top middle: radial intensity profile of the 1.3 mm continuum. Top right: radial intensity profile of 12 CO J = 2 − 1. Bottom middle: radial intensity profile of 13CO J = 2 − 1. Bottom right: radial intensity profile of C18O J = 2 − 1.

Figure A.1 shows the change in the dust structure between the fiducial model and our optimized model. The gas density structure of the optimized model remains unchanged $$n_{\text{gas}}=\{\begin{array}{ll}0\hfill & \text{for}r<R_{\text{gas}\text{in}}\hfill \\ {\delta }_{\text{gas}}\times n_{\text{gas}\text{cavity}}\hfill & \text{for}{R}_{\text{gas}\text{in}}<r<R_{\text{dust}\text{in}}\hfill \\ {\delta }_{\text{gas}}\times n_{\text{gas}\text{gap}}\hfill & \text{for}{R}_{\text{dust}\text{in}}<r<R_{\text{gap}\text{out}}\hfill \\ n_{\text{gas}}\hfill & \text{for}{R}_{\text{gap}\text{out}}<r<R_{\text{gas}\text{out}}\hfill \\ 0\hfill & \text{for}r>R_{\text{gas}\text{out}}.\hfill \end{array}$$(A.3)

Appendix B Channel maps

Figures B.1–B.3 show channel maps of the DCO⁺ J = 3 − 2 line observed with ALMA and from the modeling presented in this paper. Channel maps in Fig. B.1 are from the ALMA observations imaged with CLEAN in CASA. Channel maps in Fig. B.2 are from the best-fit model synthetic image cube created with LIME, sampled in the uv plane with the PYTHON vis_sample routine, and imaged with CLEAN in CASA. Channel maps in Fig. B.3 show the residuals of the model subtracted from the data in each velocity channel.

Fig. B.1 DCO+ J = 3 − 2 channel maps from 5.74–7.95 km s−1, Hanning smoothed to 0.17 km s−1 channels. Black contours show 5.5 mJy beam−1 (1σ) × [3, 6, 9]. The blue contours show the outline of the Keplerian mask. Channel velocity is shown in the upper-right corner. The synthesized beam and AU scale are shown in the lower-left panel.

Fig. B.2 DCO+ J = 3 − 2 best-fit model channel maps. Black contours show 5.5 mJy beam−1 (1σ) × [3]. Channelvelocity is shown in the upper-right corner. The synthesized beam and AU scale are shown in the lower-left panel.

Fig. B.3 DCO+ J = 3 − 2 residual (data-model) channel maps. Black contours show 5.5 mJy beam−1 (1σ) × [3]. Dashed contours are negative at the same intervals. Channel velocity is shown in the upper-right corner. The synthesized beam and AU scale are shown in the lower-left panel.

References

All Tables

All Figures

Fig. 1 Top-left: integrated intensity map of the DCO+ J = 3 − 2 line from 5.4–8.8 km s−1 after applying a Keplerian mask to the image cube. The synthesized beam and AU scale are shown in the lower corners. Top-right: DCO+ J = 3 − 2 integrated intensity map overlaid with white marking the model inner dust cavity, the inner dust ring (R1 , hatched), the dust gap, and the outer dust ring (R2, hatched). The synthesized beam and AU scale are shown in the lower corners. Bottom-left: velocity-weighted coordinate map of the DCO+ J = 3 − 2 line, clipped at 3.5σ. Solid black contours show the 233 GHz/1.3 mm emission at 7.0 × 10−5 Jy beam−1 (1σ) × [5, 50, 200]. The synthesized beam and AU scale are shown in the lower corners. Bottom-right: disk-integrated spectrum of the DCO+ J = 3 − 2 line before Keplerian masking, Hanning smoothed to 0.17 km s−1 channels. The horizontal dashed black line indicates the continuum-subtracted spectral baseline. The vertical red line shows the systemic velocity at 6.9 km s−1 .

In the text

	Fig. 2 Azimuthally averaged radial intensity profiles of DCO+ (top), C18O (top-middle), the 1.3 mm continuum (bottom-middle), and the J band ( ~1.2 μm) polarized light (bottom). Shaded regions represent 1σ errors on the intensity.
In the text

Fig. 3 One-dimensional, radial structure in the HD 169142 disk. Vertical brown shaded regions represent the dust rings. Left: DCO+ radial column density (dotted blue) calculated from Eq. (1) using the DCO+ the radial intensity profile (see Fig. 2) and the midplane gas temperature of the optimized DALI model (dot-dash red) as the Tex profile. The gas surface density profile (dashed green) was used to derive the DCO+ radial abundance. Right: radial abundance structure of DCO+ with respect to H2. The blue shaded region represents 1σ errors on the abundance.

In the text

	Fig. 4 Two-dimensional physical structure of the HD 169142 disk from the optimized DALI model ( Appendix A). Top: gas density structure. The gas density contours (dashed black) are shown as log(ngas). Middle: gas temperature structure below 100 K. Temperature contours are shown in dashed blue. Bottom: CO abundance structure with respect to H2.
In the text

Fig. 5 Parameterized modifications to the optimized DALI model described in Appendix A. Right and right-middleplots: gas temperature and CO abundance, respectively. The black dashed box shows the modified region. Left-middle plot: DCO+ abundance map calculated with DCO+ CHEMNET with the 30–70 K constant abundance warm component included. Left plot: DCO+ radial intensity profiles of the data (solid black with shaded gray 1σ errors) and of the model with the cold component only (C only; dashed blue) and with warm component included (C + W; dashed red). Top: optimized DALI model with no modifications made to the gas temperature or CO profile. Middle: cold disk model. The gas temperature is decreased by 2 K in the black dashed region. Bottom: shadowed and depleted cold disk model. This is an extension of the cold disk model to simulate a secondary colder region beyond the millimeter dust edge. The gas temperature is decreased by 8 K with CO depleted by a factor of five from r = 83−120 AU. The gas temperature is decreased by 2 K for r > 120 AU.

In the text

	Fig. A.1 Top: dust structure of the fiducial DALI model from Fedele et al. (2017) with small (0.005–1 μm) and large (0.005–1000 μm) grains present only in the millimeter rings. Bottom: dust structure of the optimized DALI model with small grains in the millimeter gap and in the outer disk. The hot inner dust ring is not visible due to the extremely small scale height at r = 0.2 AU.
In the text

Fig. A.2 Fits to observables for the optimized DALI model. Unlike the azimuthally averaged radial profiles in Fig. 2, the radial profiles shown here are cuts along the disk semi-major axis (PA = 5°). Top left: near-infrared to submillimeter spectral energy distribution: the model is shown by the black solid curve and the data (red dots) are from Malfait et al. 1998, IRAS (Helou & Walker 1988), Herschel (Pascual et al. 2016), SCUBA (Sandell et al. 2011), and ALMA (Fedele et al. 2017). Top middle: radial intensity profile of the 1.3 mm continuum. Top right: radial intensity profile of 12 CO J = 2 − 1. Bottom middle: radial intensity profile of 13CO J = 2 − 1. Bottom right: radial intensity profile of C18O J = 2 − 1.

In the text

	Fig. B.1 DCO+ J = 3 − 2 channel maps from 5.74–7.95 km s−1, Hanning smoothed to 0.17 km s−1 channels. Black contours show 5.5 mJy beam−1 (1σ) × [3, 6, 9]. The blue contours show the outline of the Keplerian mask. Channel velocity is shown in the upper-right corner. The synthesized beam and AU scale are shown in the lower-left panel.
In the text

	Fig. B.2 DCO+ J = 3 − 2 best-fit model channel maps. Black contours show 5.5 mJy beam−1 (1σ) × [3]. Channelvelocity is shown in the upper-right corner. The synthesized beam and AU scale are shown in the lower-left panel.
In the text

	Fig. B.3 DCO+ J = 3 − 2 residual (data-model) channel maps. Black contours show 5.5 mJy beam−1 (1σ) × [3]. Dashed contours are negative at the same intervals. Channel velocity is shown in the upper-right corner. The synthesized beam and AU scale are shown in the lower-left panel.
In the text