Molecules in the carbon-rich protoplanetary nebula CRL 2688

We present observations of the carbon-rich protoplanetary nebula (PPN) CRL 2688 made with the Institut de Radioastronomie Millimetrique (IRAM) 30 m telescope in the 3mm and 2mm bands. In total, 196 transition lines belonging to 38 molecular species and isotopologues are detected, among which, to our best knowledge, 153 transition lines and 13 species are the first report for this object. Additionally, in order to contribute to future research, we have collected observational data on the molecular lines of CRL 2688 from the literature and compiled them into a single unified catalog. We find that the molecular abundance of CRL 2688 cannot be explained by the standard model of a circumstellar envelope. The implications of metal-bearing molecules on circumstellar chemistry are discussed.


INTRODUCTION
Mass loss of asymptotic giant branch (AGB) stars fuels the recycling of materials in the galaxy (see, e.g., Herwig 2005;Ziurys 2006;Höfner & Olofsson 2018, for reviews). Suffering from complex dynamic and chemical processes, the resultant circumstellar envelopes (CSEs) are active sites for the formation of substantial amount of dust and molecular gas, making them bright infrared (IR) and radio sources. Along the evolution beyond AGB, the CSEs gradually expand, resulting in decreasing density, and the central stars would become visible at a certain moment (Bloecker 1995). With further evolution, the gas in CSEs are partly or fully dissociated and ionized by intense ultraviolet (UV) radiation from the central stars, and emits strong ionic emission lines at optical wavelengths, marking the beginning of a planetary nebula (PN) stage. Passing through the PN evolution, the circumstellar matters finally disperse into and contaminate the interstellar medium (ISM). Recently, there are increasing observational evidences showing that some molecules are able to survive in the PN evolution and then seed the chemistry in diffuse and molecular clouds (Tenenbaum et al. 2009;Schmidt & Ziurys 2016, 2018Zhang 2017).
The variations of molecular composition during the CSE evolution is a hot topic in astrochemistry (Cernicharo et al. 2011). So far, more than 70 kinds of molecular species have been detected in circumstellar environments (McGuire 2018(McGuire , 2021, and references therein). Many of these detections greatly benefited from molecular line surveys performed with advanced singledish telescopes. However, spectral line surveys have been overwhelmingly concentrated on the brightest nearby carbon star IRC+10216 (Cernicharo et al. 2000(Cernicharo et al. , 2010He et al. 2008;Patel et al. 2011;Gong et al. 2015;Zhang et al. 2017, etc.). In comparison, line surveys of other CSEs, although existing, are far sparser, which would have the risk of leading to biased views of circumstellar chemistry. It is evidently necessary to extend such observations to CSEs at different evolutionary stages.
Protoplanetary nebulae (PPNs) represent a transient phase before the PN stage when the AGB winds cease (Kwok 1993), during which the physical conditions of the CSEs change drastically on spatial and temporal scales. Given their relative short evolutionary timescale (∼10 3 yr), the PPNs that are suitable for line-survey observations are rare, and thus the chemical processes in this transition phase are far from fully understood. CRL 2688 and CRL 618 are two ideal PPNs for line-survey observations (Pardo et al. 2007;Park et al. 2008;Zhang et al. 2013). Both are carbon rich, but CRL 618 appears to be closer to the PN stage (Bujarrabal et al. 1988). The spectra of the two targets exhibit dramatically different molecular species, indicating that quick The antenna temperature (T A ) is converted to the main beam temperature (T mb ) with the expression of T mb = T A × F eff B eff , where the forward efficiency (F eff ) and main-beam efficiency (B eff ) of the telescope at each band are listed in Table 2. The data are reduced using the CLASS software of the GILDAS package 1 . The spectra at the same frequency coverage were co-added to improve the signal-to-noise ratio according to the weights proportional to 1/σ 2 , where σ is the root-mean-square uncertainty in the continuum. The identification of each transition line was made by using the information from the NIST database recommended rest frequencies for observed interstellar molecular microwave transitions 2 (Lovas 2004), the Cologne Database for Molecular Spectroscopy catalog 3 (Müller et al. 2001(Müller et al. , 2005, and splatalogue database for astronomical spectroscopy 4 . To validate the detection and identification of a given line, we checked the other transition lines from the same species that are available in the current observations and those in literature. In order to get a more complete view of the chemistry in CRL 2688, this work was complemented with observational data reported in the literature. We have collected as many results as possible from previous molecular radio observations of CRL 2688 using single-dish telescopes, except for a few old observations where the line intensity was not given correctly. For the past reports that did not explicitly give the half-power bandwidths (HPBW) of used telescopes, we take the values derived by the expression of HPBW = 75531. ′′ 2 ν s −1 and from −68 to −65 km s −1 , respectively. The NMA observations show that they are associated with the regions located at about 2 ′′ from the center to the northwest.

Silicon monosulfide -SiS
SiS in CRL 2688 was firstly discovered with the Nobeyama Radio Observatory (NRO) 45 m telescope by Fukasaku et al. (1994). So far, the SiS J = 5 → 4, 6 → 5, 13 → 12, and 14 → 13 lines and the Si 33 S J = 8 → 7 line have been detected in this object (Fukasaku et al. 1994;Zhang et al. 2013). The four SiS and isotopic lines, SiS J = 9 → 8 and 10 → 9, 29 SiS J = 10 → 9, and Si 34 S J = 10 → 9, lie within our frequency coverage. All are clearly detected, as shown in Fig. 3. The double-peaked profiles suggest extended (>14 ′′ ) and optically thin natures. There are intriguing similarities between the profiles of SiS and CS lines; both exhibits a broad component distributed from v LSR = −70 to 0 km s −1 and the strengthened red peak, probably suggesting that sulfur-bearing molecules are formed through a common chemical pathway.

Nitrogen sulfide -NS
We detected two NS lines with >3σ. Figure 3 shows the fitting result to the NS 2 Π 1/2 J = 7/2 → 5/2 l = e line. The other tentatively detected lines of NS is blended with other lines and cannot be fitted well. The NS line is too faint to accurately retrieve its profile, but it seems to be asymmetric probably due to the blend of hyperfine components. The NS 2 Π 3/2 J = 7/2 → 5/2 line at 162.67 GHz is also within our frequency coverage, but is not clearly detected due presumably to its intrinsic weakness.

Phosphorous mononitride -PN
The PN J = 2 → 1 line is detected at a 4.5σ level, as shown in Fig. 3. This is the only PN line lying within our frequency coverage. This molecule has been previously discovered by Milam et al. (2008) in this object, IRC+10216, and an oxygen-rich supergiant star VY Canis Majoris, using the Arizona Radio Observatory (ARO) 12 m telescope.

Aluminum fluoride -AlF
AlF was firstly discovered in CRL 2688 with the IRMA 30 m telescope by Highberger et al. (2001) through the J = 4 → 3, 5 → 4, and 8 → 7 lines. Two AlF lines lie within our frequency coverage, one of which (J = 5 → 4) is detected, as shown in Fig. 3. This line exhibits an approximately flat-topped profile, probably suggesting that AlF is located within a central compact region. The J = 3 → 2 line at 98.9 GHz is unfortunately located at the edge of the spectrometer where the signal to noise is relatively low.
3.1.6. Sodium cyanide -NaCN NaCN in CRL 2688 was firstly detected with the ARO 12 m telescope by Highberger et al. (2001) at the 2 mm band. Nine transition lines of NaCN are clearly detected in our observations, as shown in Fig. 3. Eight of them are newly detected. The NaCN J K a ,K c = 10 1,11 → 10 1,10 line has been detected by , which shows a consistent intensity and profile with ours. The NaCN J K a ,K c = 5 0,5 → 4 0,4 line is partly blended with the HC 7 N J = 69 → 68 line (see Fig. 3), and thus two stellar-shell components are used to decompose its profile. Most of the transition lines show a double-peaked profile, indicating that NaCN emission is optically thin and the emission region is resolved by the telescope beam.
3.1.8.  The detection of SiC 2 in CRL 2688 was firstly reported with the IRMA 30 m telescope by Bachiller et al. (1997) for the J K a ,K C = 6 0,6 → 5 0,5 line. We detected nine transition lines from SiC 2 and one isotopic transition from Si 13 CC, as shown in Fig. 3. Except for two transition lines at 165.5 GHz that are heavily blended with the CH 3 CN J = 9 → 8 line (see Fig. 3), all have a relatively high signal to noise ratio and exhibit a symmetric double-peaked profile. Moreover, four 29 SiC 2 isotopic transition lines are marginally detected, as shown in Fig. 2, but cannot be well fitted using the stellar-shell model assuming spherical symmetry.
3.1.9. Dicarbon sulfide -C 2 S Two weak C 2 S transition lines are detected, as shown in Fig. 3. The (N, J) = (7, 8) → (6, 7) line is partly blended with one vibrationally excited transition line of C 4 H. Park et al. (2008) reported a detection of the C 2 S (N, J) = (8, 8) → (7, 7) line with a flux density of 3.5 Jy. However, this line was not revealed in the more sensitive spectrum of Zhang et al. (2013). Based on our measurements and the intrinsic strength ratio between the (N, J) = (7, 8) → (6, 7) and (8, 8) → (7, 7) lines (1 : 1.16), we estimated the flux density of the latter line to be less than 0.012 Jy, confirming the speculation of Zhang et al. (2013) that this line was unlikely detected by Park et al. (2008).
3.1.10. Ethynyl radical -13 CCH Huggins et al. (1984) first reported the detection of C 2 H through the (N, J) = (1, 3/2) → (0, 1/2) line in CRL 2688 using the NRAO 11 m telescope. Although our spectral range does not cover any C 2 H lines, we detected two 13 CCH isotopic transition lines with hyperfine components, as shown in Fig. 3. They are obviously optically thin.
3.1.11. Hydrogen cyanide -HCN, hydrogen isocyanide -HNC, and formyllum cation -HCO + The J = 2 → 1 rotational transition lines of HCN, H 13 C 15 N, HNC, and HCO + and two vibrationally excited transition lines of HCN at the ν 2 = 1 state are detected, as shown in Fig. 3. The HCN J = 2 → 1 line is blended with the much weaker HCN ν 2 = 1 J = 2 → 1 l = 1e line. Zhang et al. (2013) reported a tentative detection of the HCO + J = 1 → 0 line. The existence of a small amount of HCO + in CRL 2688 is confirmed by our detection of the weak HCO + J = 2 → 1 line, although it is partly blended with the HC 5 N J = 67 → 66 line.
Similar to those of sulfur-bearing molecules, the HCN and HNC lines exhibit a broad and asymmetric component with a stronger blue wing. A prominent absorption feature, ranging in velocity from −50 to −60 km s −1 , is revealed in the broad wing of the HCN J = 2 → 1 line. The absorption ranging from −50 to −55 km s −1 has been detected in the CO J = 1 → 0, 3 → 2, and 4 → 3 , HCN J = 4 → 3, H 13 CN J = 4 → 3, and CS J = 7 → 6 lines (Kawabe et al. 1987;Jaminet et al. 1992;Young et al. 1992;Cox et al. 2000). Our observations reveal that this feature is composed of two components, ranging from −50 to −55 km s −1 and −55 to −60 km s −1 , respectively. Based on a Voigt fitting (Ginsburg & Mirocha 2011), we obtained the column densities of the two HCN absorption components to be 5.39 × 10 13 and 1.51 × 10 13 cm −2 , respectively.
3.1.12. Cyanoethynyl radical -C 3 N The first detection of C 3 N in CRL 2688 was presented by Lucas et al. (1986). We detected the N = 8 → 7, 10 → 9, 17 → 16, and 18 → 17 lines, as shown in Fig. 3. The C 3 N N = 8 → 7 and 10 → 9 lines show an optical thin double-peak profile, while the a parabolic profile is revealed for the N = 17 → 16 lines, suggesting that the higher-N lines might arise from the inner compact region that cannot be resolved by the telescope beam. The two N = 18 → 17 lines are heavily blended with each other, and cannot be well fitted using the stellar-shell model.
3.1.13. Propynylidyne radical -l-C 3 H and Cyclopropenylidene radical -c-C 3 H l-C 3 H has been discovered in CRL 2688 (Zhang et al. 2013) . In this work, we detected eight vibrationally excited transition lines of l-C 3 H, as shown in Fig. 3. Each line splits into two hyperfine structure. All show an optical thin double-peak profile.
The c-C 3 H N K −1 ,K +1 = 2 1,2 → 1 1,1 line, including two hyperfine structure components at 91.5 GHz is detected with a low signal to noise ratio (see Fig. 2). c-C 3 H has no other lines within our frequency coverage, except for the line at 91.7 GHz. However, the spectrum around 91.7 GHz is too noisy to allow us to confirm the secure detection. We thus infer that the abundance of c-C 3 H is much lower than its linear counterpart.
3.1.14. Cyanoacetylene -HC 3 N HC 3 N is detected through the J = 9 → 8, 18 → 17, and 20 → 19 lines, as shown in Fig. 3. These lines consistently exhibit three peaks and a broad wing. We also detect six ν 7 = 1 vibrationally excited lines of HC 3 N and several isotopic lines from H 13 CCCN, HC 13 CCN, and HCC 13 CN (see Fig. 3). The transition lines of HC 13 CCN and HCC 13 CN from the same energy levels are partly blended with each other. We are not able to adequately fit the profile of the J = 18 → 17 and 20 → 19 lines of HC 13 CCN and HCC 13 CN using the stellar-shell model.
The first detection of c-C 3 H 2 in a PPN (CRL 618) was made by Bujarrabal et al. (1988), who also observed CRL 2688. But they did not report any c-C 3 H 2 line in CRL 2688 although some lines are obviously lying in their spectral coverage. We detected three c-C 3 H 2 lines, as shown in Fig. 3. The double-peaked profile suggests that they are likely to be optically thin.
3.1.16. Butadiynyl radical -C 4 H C 4 H in CRL 2688 was detected for the first time by Lucas et al. (1986) through two N = 10 → 9 transition lines. But it was only tentatively seen in CRL 618 by Bujarrabal et al. (1988). All the eight rotational lines lying within our spectral coverage are detected with a prominent double-peaked profile, as shown in Fig. 3. In addition, we detected six and two vibrationally excited transition lines within the 2 Π 3/2 and 2 Π 1/2 states, respectively, as shown in Fig. 3. The C 4 H ν 7 = 1 2 Π 3/2 J = 33/2 → 31/2 line is significantly blended with the C 13 CCCH N = 17 → 16 J = 35/2 → 33/2 line (see Fig. 3), and cannot be well fitted with the stellar-shell model. Due to the low signal to noise ratio, it is hard to tell whether the profile of these vibrationally excited lines is double-peaked or flat-topped.
3.1.17. Pentynylidyne Radical -l-C 5 H The first detection of l-C 5 H in CRL 2688 was made by Highberger et al. (2001) through the 2 Π 1/2 J = 65/2 → 63/2 l = e & f line. Several vibrationally excited transition lines of l-C 5 H at 2 Π 3/2 and 2 Π 1/2 states were subsequently detected Zhang et al. 2013). We detected six vibrationally excited lines of l-C 5 H, as shown in Fig. 3. Each line is split into l = e and f hyperfine structure components. Other six vibrationally excited lines of l-C 5 H are lying within our spectral coverage but not being detected due presumably to their intrinsic weakness; their intensity upper limits are estimated to be 10 mK.
3.1.20. Hexatriynyl radical -C 6 H C 6 H was marginally detected in CRL 2688 with the ARO 12 m telescope by Highberger et al. (2001), and then the detection of the line at 102 GHz was confirmed  using the IRMA 30 m telescope. Six and five vibrationally excited transition lines of C 6 H, respectively in the 2 Π 3/2 and 2 Π 1/2 states, are detected in our observations as shown in Fig. 3. The profiles of the eight strongest C 6 H lines are well fitted by using the stellar-shell model, suggesting an optically thin nature.

Unidentified lines
Four U-lines that cannot be assigned to known transition lines of any molecular species are newly detected, increasing the number of U-lines in CRL 2688 to 13. It is clear that the frequency spans of these U-lines cannot be accounted for by a single rotation constant, suggesting that they are unlikely to arise from the same species. Probably they are associated with unknown circumstellar molecules or internal rotations of polymer molecules, which call for further experimental investigation. For comparison, we examine the U-lines in the spectrum of IRC+10216 ranging from 129.0 to 172.5 GHz presented by Cernicharo et al. (2000), but do not find any match. We cannot rule out the possibility that some of them are artificial features.

A comparison with previous line survey observations of CRL 2688
Our spectral range has overlaps with those of Park et al. (2008) and Zhang et al. (2013) in frequency coverages, allowing us to make a comparison between these measurements. Within the overlapped frequency range from 91.38-99.16 GHz, the CS J = 2 → 1 line is the only line discovered by Park et al. (2008), while we unambiguously detect 35 lines. All the 23 lines observed by Zhang et al. (2013) in the overlapped frequency range from 75. 70-83.48 GHz are visible in our spectra, in which we totally detect 102 lines. The HPBWs of the IRAM 30 m and ARO 12 m telescopes are different (69. ′′ 4-64. ′′ 0 and 26. ′′ 9-24. ′′ 8, respectively, for the overlapped frequency range), resulting in different main beam temperatures for the the same line between ours and that measured by Zhang et al. (2013). Through a comparison between these measurements, we can roughly estimate the size of the region where the molecular lines originate. The same technique has been used by Takano et al. (2019) and Qiu et al. (2020) to investigate the molecular emission regions in extragalactic objects.
With the assumption of a Gaussian intensity distribution, the velocity-integrated intensities measured by different telescopes should satisfy the relationship I mb (θ represent the intensities and beam sizes corresponding to the IRAM 30 m and ARO 12 m telescopes, respectively, and θ s represents the Gaussian full width at halfmaximum (FWHM) of the source size. To estimate θ s , we define a parameter for the ith line where σ i is the uncertainty of I mb,i /I ′ mb,i . Then an F factor is introduced by F = i f i . Here, f i is the residual of the ratio between the velocity-integrated intensities measured by IRAM 30 m and ARO 12 m telescopes considering the beam dilution weighted by 1/σ i . The source size θ s can be obtained by minimizing the F factor. We find that the F factor reaches the minimum value if the source size is 15. ′′ 5, as shown in the upper panel of Fig. 4. The CO mapping of CRL 2688 in the J = 1 → 0 line shows that the emission is dominantly from the central core with a radius of 10 ′′ -13 ′′ (Kawabe et al. 1987), while the 13 CO image of the J = 1 → 0 line reveals a high-brightness central core, a high-velocity component near the center, and a less bright extended envelope with a diameter of about 15 ′′ (Yamamura et al. 1995). These observations are in reasonable agreement with our size estimation based on single-dish data. The F factor cannot reach zero (see Fig. 4), probably suggesting that the molecular gas has a stratified structure and/or the brightness distributions deviate from the Gaussian assumption. To further examine this, we plot in the lower panel of Fig. 4 ] ratios for all the lines lying within the overlapped frequency range, where θ s = 15. ′′ 5. The deviations from one in this plot is an indication that different molecules may have different spatial distributions, and thus the obtained source size only represent the average value.
Alternatively, we could assume that the source has a disk geometry, where the brightness distribution is constant within a diameter of θ d . Then we have . Following a similar way as described above, but defining f i by we find θ d = 24. ′′ 7 when F approaches the minimum (see the upper panel of Fig. 4). This means that we would obtain the same beam-filling factors if assuming θ d = 24. ′′ 7 for a disk distribution or assuming θ s = 15. ′′ 5 for a Gaussian distribution. Therefore, b + θ 2 s )] ratios for θ s = 15. ′′ 5, as shown in the lower panel of Fig. 4. It is natural that θ s is smaller than θ d since the former is the FWHM of the surface brightness profile while the latter represents the extent of the source.

Rotation diagram, column Densities, and fractional Abundances
The column densities (N) and excitation temperatures (T ex ) of 22 molecules and isotopologues, including CS, 13 CS, C 33 S, C 34 S, HCN, HNC, HCO + , SiS, C 2 S, SiC 2 , l−C 3 H, C 3 N, HC 3 N, H 13 CCCN, HC 13 CCN, HCC 13 CN, C 4 H, c-C 3 H 2 , C 6 H, CH 3 CN, HC 5 N, and HC 7 N, are estimated by means of the rotation diagram. Under the assumption of local thermal equilibrium (LTE) and negligible optical depth, the level populations can be expressed with where k and h are the Boltzmann constant and Planck constant, ν is the rest frequency in Hz, Q(T ex ), A ul , g u are the partition function, spontaneous emission coefficient, and upper state degeneracy, respectively. Under the assumption of a Gaussian distribution of the surface brightness, the source brightness temperature T S is derived after correcting for beam dilution, T S = T mb (θ 2 b + θ 2 s )/θ 2 s , where θ b is the antenna HPBW and θ s is the source size. T S dν is the velocity-integrated intensity. Q(T ex ), A ul , g u , and E u /k are taken from the CDMS catalog. Following the previous works (Park et al. 2008;Zhang et al. 2013), we simply assume a source size of 20 ′′ for the calculations of all the molecules. The resultant rotation diagrams are presented in Fig. 5, and the results are listed in Table 4.
If the optical depth (τ) poses an effect, the dots in the rotation diagram might deviate from a linear distribution. If this is the case, an optical depth correction factor of C τ = τ/(1 − e −τ ) needs to be introduced to construct the rotation diagrams (Goldsmith & Langer 1999). However, such a deviation cannot be clearly seen in Fig. 5. Using the population diagram method developed by Goldsmith & Langer (1999), we estimate ln C τ < 0.03, suggesting that the effect of optical depth is negligible.
We use the expression given by Olofsson (1997) to derive the molecular abundances with respect to molecular hydrogen ( f X ), where v e is the expansion velocity in km s −1 , D is the distance in pc (assumed to be 1 kpc hereafter),Ṁ H 2 is the mass-loss rate in M ⊙ yr −1 , ν is the rest frequency in GHz, E l is the energy of the lower level, T mb dv is the velocity-integrated intensity in K km s −1 , x e,i = R e,i /(θ b D) with R e and R i for the outer and inner radii of the shell, and θ b , Q(T ex ), g u , A ul , and k are the same as in Equation 1. To simplify the calculations, we assume the values of x e and x i to be 1 and 0 for all the species, respectively. For the molecules that have more than one transition line detected, we adopt the average abundances weighted by their velocity-integrated intensities. As mentioned in Zhang et al. (2009), the abundances is essentially insensitive to T ex in the case of E l ≪ kT ex . The results are listed in Table 4. It should be noted that if the transition lines are optically thick, the derived f X values should be considered as lower limits.

The fractional abundances derived from radiative transfer method
In order to determine the fractional abundances ( f X in Table 4), we have assumed that all the lines originating from a given molecule are characterized by a single excitation temperature, which is generally valid if the molecule is located within a narrow layer. However, CRL 2688 exhibits a moderate radial temperature gradient (Truong-Bach et al. 1990). For the abundance computation of molecules lying within a region with wide temperature variation, a more sophisticated approach is to use a radiative transfer model. As such, we redetermine the fractional abundances using the three-dimensional non-local-thermodynamical (NLTE) radiative transfer code LIME (Brinch & Hogerheijde 2010). We limited the calculations to a few strong lines, aiming to examine the systematic uncertainties in the determination of f X .
For the modelling, we assume a spherically symmetric structure with a velocity law of v(r) = v exp (1 − r 0 /r) 0.5 , where r 0 is the inner radius of the envelope and is taken to be 10 15 cm (Truong-Bach et al. 1990). Such a velocity law is the solution of the equations of the stellar wind driven by a force decreasing with r −2 . The terminal velocity v exp is taken from Table 3. In most parts of the envelope (r > 10 16 cm), v(r) is very close to v exp . We thus simply assume a H 2 number density law of n(r) =Ṁ/(4πm H 2 r 2 v exp ), where m H 2 is the mass of a H 2 molecule, and the mass-loss rateṀ is taken to be 1.7 × 10 −4 M ⊙ yr −1 (Truong-Bach et al. 1990) . Even though the assumed mass-loss rate may be slightly larger than the current value of CRL 2688, we consider that it would be acceptable as a first approximation. Moreover, we apply a radial distribution of the kinetic temperature obtained by Truong-Bach et al. (1990) through high-resolution mapping of the CO J = 1 → 0 and 2 → 1 lines. The molecular data, including the energy states, Einstein coefficients, and collisional rates, are taken from the LAMDA database 5 (Schöier et al. 2005). The output spectral resolution is set to be 1 km s −1 .
Fixing the distance between CRL 2688 and earth at 1 kpc and allowing the inner radius (R i ) and outer radius (R e ) of individual molecules to vary, we determine the fractional abundance by fitting the LIME model to the observations. The fittings for the CS J = 2 → 1, SiS J = 10 → 9, PN J = 2 → 1, SiC 2 J K a ,K c = 4 0,4 → 3 0,3 , HCN J = 2 → 1, HNC J = 2 → 1, HCO + J = 2 → 1, HC 3 N J = 9 → 8, c-C 3 H 2 J K a ,K c = 2 0,2 → 1 1,1 , and CH 3 CN J K = 5 3 → 4 3 lines are shown in Fig 3. The LIME models appear to match the main components of the observed line profiles reasonably well, but fail to reproduce the broad wings of CS, SiS, HCN, and HNC lines. An interferometric observation has revealed that HCN and HC 3 N emission primarily traces dense gas in polar and equatorial bipolar outflows (Dinh-V-Trung et al. 2020). Therefore, the mismatches between the models and the observations may suggest a deviation from a spherical geometry. In order to reproduce the present observations more closely, it would be necessary to assume a more realistic geometry, which is possibly determined by high-angular resolution interferometry. However, such the observation is beyond the scope of the present research. Table 5 summarizes our fitting results, showing that the agreement between the molecular abundances derived from the LIME models and f X is better than one order of magnitude.
3.5. Isotopic abundance ratios 5 http://www.strw.leidenuniv.nl/ moldata Isotopic abundance ratios reflect the nucleosynthesis in the stellar interior and mixing processes (Herwig 2005). We determine the isotopic abundance ratios of C, S, Si, and Mg elements in CRL 2688 based on the intensities of the corresponding molecular transition lines. The results, along with those of the Sun and IRC+10216 (Lodders 2003;Kahane et al. 2000), are given in Table 6. The 12 C/ 13 C ratio is believed to be one of the most useful tracers of the relative degree of primary to secondary processing in stars, which drops during the first dredge-up and is enhanced by the third dredge-up. Our measurements suggest a 12 C/ 13 C ratio of 15-40 in CRL 2688, slightly lower than that in IRC+10216 and significantly lower than the solar value. This can be attributed to the more evolved nature of CRL 2688. Strikingly, the 12 C/ 13 C ratios derived from HC 3 N lines tend to increase with increasing J (see Table 6). Chemical fractionation may alter the isotopic ratio. The reaction, 13 C + + CN → 13 CN + 12 C + , preferentially takes place under low temperature (Watson et al. 1976;Langer 1992;Milam et al. 2005). As CN is the key parent molecule of cyanopolyynes, the enhanced 13 CN would leads to a higher abundance of the 13 C isotopologues of HC 3 N in the cool region. If the gas temperature of the CSE has a negative gradient form the inner side to the outer side, the low-J lines preferentially arise from the outer low-temperature region, respective to the high-J lines. This offers a seemingly plausible explanation for the observations that the HC 3 N low-J lines result in a lower 12 C/ 13 C ratio. However, chemical fractionation is efficient only at an extremely low kinetic temperature. Through theoretical calculations, Milam et al. (2005) found that for a molecular cloud with a kinetic temperature ranging from 10-100 K, the effect of fractionation causes negligible alteration of the 12 C/ 13 C ratio during the typical life timescale (10 6 yr). Although PPNs are denser than molecular clouds, tending to decrease the reaction timescale, their lifetime is much shorter (∼10 3 yr). Moreover, PPNs are energetic sources with a kinetic temperature of much higher than 10 K. Therefore, chemical fractionation is unlikely to affect the 12 C/ 13 C ratio in CRL 2688. Isotope-selected photodissociation could also affect the 12 C/ 13 C isotopologue ratio. For instance, 12 CO is more abundant than 13 CO, and thus is more self-shielding from the interstellar UV field, resulting in an increasing 12 C/ 13 C ratio. HC 3 N is the product of photochemistry of HCN and CN, which thus would result in a radial 12 C/ 13 C gradient for HC 3 N if their isotopologues have different self-shielding factors. However, unlike CO, HCN and CN are dissociated mainly in the continuum rather than bands, suggesting that their isotopologues are equally affected by UV radiation (Saberi et al. 2017). As a result, we can rule out the isotope-selected photodissociation as the cause of the 12 C/ 13 C gradient observed for HC 3 N. Giesen et al. (2020) found that the C 13 CC/ 13 CCC isotopic ratio in the massive star-forming region Sgr B2(M) is higher than the statistically expected value, which is probably due to the lower zero energy of C 13 CC relative to 13 CCC making the positionexchange reaction of converting 13 CCC to C 13 CC favorable. We find that the column densities of HC 13 CCN, H 13 CCCN, and HCC 13 CN are quite consistent with each other, suggesting that such an effect is not significant for HC 3 N isotopolgues.
Presumably, the isotopic abundance ratios of very heavy elements are not altered during the AGB evolution. We do not find a significant difference of the S, Si, and Mg isotopic abundance ratios between CRL 2688, IRC+10216, and the Sun, being consistent with the prediction of stellar theory. The only exception is the CS J = 2 → 1 line which gives a relatively low 32 S/ 34 S ratio. But considering that this line is likely to be optically thick, the resultant S isotopic ratio only represents a lower limit.

Temporal variations
Evidences of temporal variations of molecular lines in IRC+10216 were presented by Cernicharo et al. (2014), who found that most of the abundant molecules, except SiC 2 , exhibit strong variations in line intensities and the variations are more pronounced for high-J lines. This can be attributed to the variations in the IR pumping rate. In this scenario, optically thick lines, such as low-J lines of SiC 2 , are insensitive to the IR pumping. Such studies are relevant to the reliability of the intensities of molecular lines in AGB envelopes as the flux calibrator and the tracer of mass loss rate.
CRL 2688 is a commonly used standard source for flux calibration. The brightness of its lobes has been found to increase by −0.22 V magnitude during 1994-2008 (Hrivnak et al. 2010). Klochkova et al. (2000) found that the radial velocity of the envelope varies with time. Presumably, the rapid variable stellar brightness and circumstellar dynamics may lead to temporal variations of the intensities of some molecular lines in the CSE. The C 4 H J = 21 → 19, C 4 H J = 19 → 17, SiC 2 J K a ,K c = 4 2,2 → 3 2,1 , and H 13 CCCN J = 11 → 10 lines in CRL 2688 were also observed in May 1985 by Lucas et al. (1986) using the same observational configurations as ours. We investigated the intensity variations using the data sets spanning over 32 years. The velocity-integrated intensity ratios of the four lines between the measurements of Lucas et al. (1986) and ours are 1.6 ± 0.1, 1.7 ± 0.1, 0.6 ± 0.2, and 1.2 ± 0.1, respectively. It is notable that the two C 4 H lines show a consistently increasing velocity-integrated intensity ratio, whereas no significant variation is found for the SiC 2 and H 13 CCCN lines. Although we cannot conclusively attribute this to the luminosity variations, it would be interesting to study the different variational behaviors of these molecular lines in the future.

A comparison with previous line survey observations of IRC+10216
A λ = 2 mm line survey of IRC+10216, the prototype of carbon-rich CSEs, has been presented by Cernicharo et al. (2000) with the IRAM 30 m telescope. An overlap of their and our spectra occurs in the spectral range from 160.8 to 168.6 GHz. To investigate the alteration of physical and chemical environments during the AGB-PPN evolution, we make a comparison of the velocity-integrated intensities of the lines in the overlapping range between IRC+10216 and CRL 2688. The line intensities, normalized such that I(HC 3 N J=18→17) = 1, are shown in Fig. 6, which reveals systematic differences. The 13 C isotopologues of HC 3 N shows slightly strengthened intensities in CRL 2688 with an unknown reason (see Section 3.5). CRL 2688 also exhibits enriched cyanides, such as HC 5 N and CH 3 CN, which can be ascribed to the enhancement of CN fragments resulted from the photolysis of HCN during the AGB-PPN transition. The vibrational excitation line of HC 3 N is largely strengthened in CRL 2688, suggesting a higher temperature in the inner regions of the PPN where this line arises. In CRL 2688, the lines of C-rich molecules, such as C 3 N, l-C 3 H, and C 4 H, and of Si-bearing molecules, such as SiS and SiC 2 , appear to be consistently weaker than those in IRC+10216, This can be interpreted if IRC+10216 is intrinsically richer in carbon abundance and the gas-phase refractoryelement molecules in CRL 2688 have been heavily depleted onto grains during the PPN expansion. The depletion of SiC 2 of CRL 2688 relative to IRC+10216 is consistent with the 3 mm observations of Lucas et al. (1986). However, Lucas et al. (1986) suggested an enhanced C 4 H in CRL 2688, contrary to our result and that of Zhang et al. (2013). This might be due to the effect of excitation conditions and spatial distribution of this molecules, as pointed out by Zhang et al. (2013). We find that the column density of l-C 3 H is higher than that of c-C 3 H by a factor of 3 in CRL 2688, where the l-C 3 H/c-C 3 H ratio is close to that of IRC+10216 (∼2; Cernicharo et al. 2011). The higher l-C 3 H/c-C 3 H ratio in IRC+10216 relative to that in the dark cloud TMC-1 has been regarded as an evidence of atom-neutral reactions (Kaiser et al. 1996).

Modelling and circumstellar chemistry
Chemical modelling is a powerful tool for deepening our understanding of circumstellar chemistry (e.g. Millar & Herbst 1994; Agúndez et al. 2020). The abundance patterns of many AGB molecules can be satisfactorily reproduced by state-of-the-art astrochemistry models, which, however, seemingly perform less well at explaining the formation of complex molecules in PPNs. New observations continue to show unexpected results, prompting the development of models. To investigate the shortcomings of the models in studying PPN chemistry, we compare the observations with model predictions. For the modelling, we use the package presented by McElroy et al. (2013), which is based on an updated UMIST Database for Astrochemistry (RATE12). The reaction network consists of 6173 gas-phase reactions and 195 reactions on grain surfaces involving 467 atomic and molecular species. We assume an uniform loss-rate of 3.0 × 10 −5 M ⊙ yr −1 (Lo & Bechis 1976) and an spherically symmetrical expansion velocity of 16.5 km s −1 (Crampton et al. 1975). The gas density radially decreases as 1/r 2 with an inner radius of 1.5 × 10 16 cm (∼1 ′′ at D = 1 kpc) where the parent species are injected into the CSE and an outer radius of 1.5 × 10 18 cm (∼100 ′′ at D = 1 kpc) where the majority of molecules are dissociated by a standard interstellar ultraviolet radiation field. The radiation field of the central star is ignored.
The modelling results are listed in Table 4 and are compared with the observations in Fig. 7. We find that for 19 out of the 26 detected molecular species, the predicted and observed column densities show agreement within one order of magnitude. An excellent agreement is shown for the abundance patterns of cyanopolyynes (HC 2n−1 N: n = 1-4); the only exception is HCN, which exhibits a much lower observed abundance and could be attributed to the effect of optical depth. For other C-rich species, C n H (n = 2-6) and C n S (n = 1-3), the models can reproduce the observations to a reasonable accuracy; there are only slight disagreements for C 5 H and and C 6 H. Based on the 20 GHz band data obtained with the 100 m Green Bank Telescope (GBT), Gupta et al. (2009) derived the C 6 H column density of 4.5 × 10 12 cm −2 , which is about six times lower than our result, and even worse when comparing with the modelling result. Thus the disagreements are likely to result from the uncertainties of the model. The other molecules that cannot be modelled well are PN, H 2 CS, H 2 CO, and SiO, which show the observed-to-calculated abundance ratios of 826, 42, 151, and 0.001, respectively. Figure 8 shows the modelled abundances of the molecules in CRL 2688 as a function of radius as well as the observed values. It is obvious that the observed abundances of PN, H 2 CS, and H 2 CO are higher than the modelled results at any radius.
Apart from the uncertainties introduced by observations and data analysis, the possible causes of the disagreements include poor knowledge of the physical conditions, missing routes in the modelling, and inaccurate reaction rates. To decouple the affects of uncertain physical conditions, we compare the observed-to-calculated abundance ratios of the molecules in IRC+10216, the progenitor object of CRL 2688. The results are shown in Fig. 7. The modelled and observed data of IRC+10216 are taken from McElroy et al. (2013) and the references therein. Generally, the agreement of the molecules in IRC+10216 is slightly better than those in CRL 2688. However, it is interesting to note that the models overestimate the abundances PN, H 2 CS, and H 2 CO in CRL 2688 and IRC+10216 by about the same degree, suggesting that significant reaction routes for forming these molecules might been neglected or undervalued in the modelling.
The enhancement of SiO has been commonly observed in C-rich AGB stars, and can be attributed to shock-induced chemistry (Schöier et al. 2006) that has not been taken into account in this model. CRL 2688 shows rich shock features caused by fast outflows. However, the largely overestimated SiO abundances by the model are against with the predictions of shock chemistry. We infer that the model might substantially underestimate the freeze-out of SiO in the PPN stage. Alternatively, because the fast outflows in CRL 2688 are highly collimated, the shock chemistry processes only in rather small regions and does not play an important role on the enhancement of gas-phase SiO. Severe depletion of SiO and SiS has been also found in another PPN IRAS 22272 + 5435 (Zhang 2020), which exhibits strong 21 µm feature. The CO observations suggest that the material ejection of IRAS 22272 + 5435 is dominated by the spherical wind, but a jet is developing and interacting with ambient materials (Nakashima et al. 2012). CRL 2688 has well developed jets and is thought to just metamorphose from a '21 µm phase' (Geballe et al. 1992). Therefore, we reasonably conjecture that the depletion of Si-bearing molecules and the development of jets may be associated to the phenomenon of the 21 µm feature in PPNs. Sensitive observations of a large PPN sample are required to draw firm conclusion.

Metal-bearing molecules
CSEs surrounding evolved stars are characterized by unambiguous identification of metal-bearing (Na-, Mg-, Al-, K-, Ca-, and Fe-bearing) molecules. These molecules have never been discovered in other types of sources, even in the rich molecular cloud TMC-1, which might be attributed to severe depletion onto dust grains (Turner et al. 2005 (Cernicharo et al. 2019;Pardo et al. 2021, and references therein), indicating a rich gas-phase metal chemistry in C-rich circumstellar environments. The first identification of metal halides was made by Cernicharo & Guélin (1987), who also reported a tentative detection of AlF, later confirmed by Ziurys et al. (1994). According to the chemical equilibrium calculations (Tsuji 1973), these stable closed-shell diatomic molecules are mostly formed in the warn inner shell, where the LTE condition holds. Kawaguchi et al. (1993) identified MgNC for the first time in the spectrum of IRC+10216 presented by Guélin et al. (1986). The metal cyanide has an open-shell structure and should be formed in the outer envelope along with other radicals, facilitating the successful detection of NaCN in IRC+10216 (Guélin et al. 1993;Turner et al. 1994). The different spatial distributions of metal halides and metal cyanides are revealed by line profiles and high spatial resolution maps of IRC+10216 (Ziurys 2006, and references therein). When the CSE expands, materials move from the inner to outer regions. As a result, one could expect that the AlF/MgNC abundance ratio decreases through the transition from the AGB to the PPN phase. This can be tested by measuring the abundance of metal-bearing molecules in CRL 2688.
For comparison, Figure 9 and Table 7 show the results of Highberger et al. (2001) who used the same transitions for AlF, the N = 11 → 10, 12 → 11, and 13 → 12 transitions for MgNC, and the J K a ,K c = 9 0,9 → 8 0,8 , 10 0,10 → 9 0,9 , 10 3,8 → 9 3,7 , 10 3,7 → 9 3,6 , and 10 2,8 → 9 2,7 transitions for NaCN. We find that the results of AlF are consistent with those of Highberger et al. (2001), whereas there a slight and significant disagreement for MgNC and NaCN, respectively. Under the assumption of a source size of 5 ′′ , we obtain the column density of NaCN an order of magnitude lower than that by Highberger et al. (2001). However, the difference can be alleviated if NaCN arises from a more extended shell. Based on observed line profiles,  suggested that NaCN in CRL 2688 should be sited much farther away from the central star than that in IRC+10216. In Fig. 10, we show the column densities calculated by using different source sizes for the beam-dilution corrections. We can see that if the source sizes of AlF and MgNC are 5 ′′ and 20 ′′ , respectively, our results are perfectly consistent with those of Highberger et al. (2001). However, even for source size of 30 ′′ , the column density of NaCN derived by Highberger et al. (2001) is still 6 times larger than our value, suggesting that the disagreement cannot be completely ascribed to the uncertainty of the source size.
The fractional abundances are determined following the approach described in Highberger et al. (2001) and Highberger & Ziurys (2003). The input parameters and the results are listed in Table 8. The assumed geometries are the same as those used by Highberger et al. (2001) and Highberger & Ziurys (2003); AlF is sited in a spherical envelope while NaCN and MgNC are confined on spherical shells. We derive a AlF/MgNC abundance ratio of ∼0.7, significantly lower than that in IRC+10216 (∼3), supporting the theoretical expectation of the CSE expansion. Further support for the chemical evolution comes from the observations of CRL 618, the presumable descendant of CRL 2688, in which MgNC is the only discovered metal-bearing molecule (Highberger & Ziurys 2003).
The fractional abundances of these metal-bearing molecules can provide a lower limit of the gas-phase elemental abundance. Our calculations show that magnesium locked in MgNC has a Mg/H ratio of 2 × 10 −8 (Table 8). According to the element depletion parameters given by Jenkins (2009), one can estimate that the gas-phase Mg/H ratio in the ISM is larger than 2 × 10 −6 . If assuming that solar abundances can well reflect the total elemental abundance and the depletion factor in the PPN is similar to that of the ISM, we could conclude that less than 1% gas-phase magnesium is in the form of MgNC, and thus MgNC may not be the dominant Mg-bearing molecule in CRL 2688.
Fluorine has a cosmic abundance much lower than aluminium, and thus is the key element of AlF. The astrophysical origin of fluorine remains a widely debated problem (see, e.g., Indelicato et al. 2017;Spitoni et al. 2018, and references therein). Because the only stable fluorine isotope, 19 F, is fragile in hot stellar interiors, the most possible explanation for cosmic fluorine abundance is that most of the F escapes from the harsh stellar interiors immediately after its production. The AGB star is an ideal source of F, where 19 F is synthesized and brought to stellar surface during the He-burning thermal pulses and the third dredge-up (Forestini et al. 1992;Jorissen et al. 1992). The F enhancement in AGB stars has been confirmed by the observations of F ions in PNs (e.g. Zhang & Liu 2005). The most abundant interstellar F-bearing molecule is HF, which, however, is inaccessible from ground telescopes. CF + has been found to be the second abundant one in photodissociation regions (Neufeld et al. 2006). Our research group has systematically searched for CF + in CSEs, but failed to see any positive detection (Zhang et al. 2014). Therefore, AlF provides an unique opportunity to constrain the F abundance in CSEs. A discussion on this aspect has been made by Highberger et al. (2001). We reexamine this problem using our data of CRL 2688. The derived AlF/H 2 abundance ratio is 2.7 × 10 −8 (Table 8), suggesting a F/H abundance of >1.4 × 10 −8 . The obtained lower limit of F/H is only 2.4 times lower than the solar value of Lodders (2003), and is presumably far below the fluorine abundance in CRL 2688 since AlF is not the dominant F-bearing molecule in circumstellar environments. Although the α element Al in the PPN can be enhanced through the third drudge-up, given the fact that available hydrogen is much richer than aluminium in circumstellar environments (solar Al/H =3 × 10 −6 ; Lodders 2003), the abundance of HF should be larger than AlF by orders of magnitude. As a result, we infer that CRL 2688 has a significant F enhancement. Combined with the C-rich nature of CRL 2688, our results further confirm the third drudge-up in AGB stars as the source of F production.

SUMMARY
CRL 2688 is a prototypical source suitable for studying the chemical alteration of CSEs during the transition from AGB to PPN phase. It might hold key clues to investigate the puzzle of the 21 µm feature. This paper reports a new observation towards this PPN at the 3 mm and 2 mm bands, which results in clear detections of 196 transition lines arising from 38 molecules. Although no new circumstellar molecule is discovered, 153 transition lines and 13 molecules are new detections in this object. We compile a comprehensive list of the gas-phase molecules discovered in CRL 2688 so far, providing a valuable reference for the community of circumstellar chemistry. By comparing the line intensities observed by different telescopes, we estimate a Gaussian brightness FWHM of 15. ′′ 5 or a source size of 27. ′′ 7 for a disk geometry, although the exact sizes are different from species to species. The standard chemical modelling is unable to reproduce the observed abundances of PN, H 2 CS, and H 2 CO, suggesting that some key reaction routes may have been missed in the chemical network. Through a comparison between the molecular abundances in IRC+10216, CRL 2688, and CRL 618, we conclude that the differences can be partially attributed to the effect of CSE evolution. The three CSEs are the only sources where metal-bearing molecules are detected. The results of the analysis support the theoretical expectations that AlF is confined in the inner shell while MgNC and NaCN are located on extended regions. An enhancement of F in the PPN is revealed by the fractional abundance of AlF relative to H 2 .
8 Derived from the archive date of Gaia TGAS, which gives a parallax of 2.78 ± 0.26 mas.
9 Derived from the archive date of Gaia DR2, which gives a parallax of 0.31 ± 0.41 mas.
10 Obtained by assuming a distance of 1 kpc (Crampton et al. 1975) and a isotropic IR                                Notes.
The temperature scale used in this table is given in term of main beam temperature. Columns 1, 2, and 3 give the rest frequency, the species, and the transition, respectively. Columns 4 and 5 give the expansion velocity and Horn-to-Center ratio of the the line profile, respectively, obtained from the stellar-shell model fitting. Columns 6 and 7 give the central temperature of the line profile and velocity-integrated intensity, respectively. Columns 8, 9, and 10 give the rms noise of the base line, the corresponding velocity resolution, and the FWHM of each line, respectively. The telescopes used to perform the observations and the corresponding references are listed in Columns 11 and 12, respectively. Column 13 provides comments for the details of lines. Boldface denotes newly detected lines in this work. The rows are ordered by increasing carbon-atom number. Symbol: ? means the values that can not be derived from our detection or are not listed in the archive data.
× means the values that can not be obtained because of unsolved hyperfine structures or two blended spin-rotation components. : means uncertain measurements. a Blended fine-structure group (J = 3/2 → 1/2).
b Unsolved hyperfine structure lines.
f The 13 CS J = 2 → 1 and C 3 S J = 16 → 15 lines are blended with each other. The former is stronger.
g Blended line.
i The AlF J = 4 → 3 and SiC 4 J = 43 → 42 lines are blended with each other. The former is stronger.
j The HC 3 N J = 29 → 28 and AlF J = 8 → 7 lines are blended with each other. The former is stronger.
m Two blended spin-rotation components.
o The line intensity was obtained by integrating the main beam temperature over the velocity coverage from -80 to -20 km s −1 .
q The line intensity was obtained by integrating the main beam temperature over the velocity coverage from -110 to -20 km s −1 .
t The line intensity was obtained by integrating the main beam temperature over the velocity coverage from -60 to -10 km s −1 .
v The HCN J = 2 → 1 and HCN ν 2 = 1 J = 2 → 1 l = e lines are blended with each other. The former is stronger.
w The HCO + J = 2 → 1 and HC 5 N J = 67 − 66 lines are blended with each other. The former is stronger.
x The line intensity was obtained by integrating the main beam temperature over the velocity coverage from -70 to 0 km s −1 .
z The line intensity was obtained by integrating the main beam temperature over the velocity coverage from -90 to 10 km s −1 .
α The HC 3 H J = 12 → 11 and HC 5 H J = 41 → 40 lines are blended with each other. The former is stronger.
β Blended with the corresponding transition lines from HCC 13 CN.
γ The HC 13 CCN and HCC 13 CN lines arising from the same J transition lines are blended.
c Taken from .
d This work.