The distance to CRL 618 through its radio expansion parallax

CRL 618 is a post-AGB star that has started to ionize its ejecta. Its central HII region has been observed over the last 40 years and has steadily increased in flux density at radio wavelengths. In this paper, we present data that we obtained with the Very Large Array in its highest frequency band (43 GHz) in 2011 and compare these with archival data in the same frequency band from 1998. By applying the so-called expansion-parallax method, we are able to estimate an expansion rate of 4.0$\pm$0.4 mas yr$^{-1}$ along the major axis of the nebula and derive a distance of 1.1$\pm$0.2 kpc. Within errors, this distance estimation is in good agreement with the value of ~900 pc derived from the expansion of the optical lobes.


INTRODUCTION
The fate of a star with a main sequence mass in the range from 0.8 to 8 M ⊙ is well established: it goes through the Asymptotic Giant Branch (AGB) phase, then into the Planetary Nebula (PN) phase and eventually it ends its evolution as a white dwarf (Kwok et al. 1978).However, the details of each evolutionary stage, such as the formation and the early evolution of PNe are not clear.For example, it is still uncertain which mechanism/s operate to produce the often complex morphology observed in a fraction of PNe and post-AGB stars, while the circumstellar envelopes (CSEs) around their progenitors (AGB stars) show regular and symmetric structures.
It is commonly accepted that the circumstellar shells of PNe are the product of the interaction between the residual, slowly expanding (∼20 km s −1 ) AGB envelope and a subsequent, rapid (100-1000 km s −1 ) shaping agent.High-speed collimated outflows, or jets, that operate during the early post-AGB evolutionary phase, or even towards the end of the AGB, have been proposed as the primary agent in shaping CSEs, the origin of such winds being linked to binarity (Soker 1998;De Marco & Izzard 2017).
To investigate this matter, many authors have tried to identify very young PNe or proto-Planetary Nebulae (PPNe), but this has turned out to be quite difficult, due to the short duration (≈1000 yr) of the transition from the AGB to the PN stage.One of the few objects that can be used for such studies is CRL 618, a PPN that started its post-AGB stage about 200 years ago and is rapidly evolving towards the PN stage (Westbrook et al. 1975).The rapid evolution of this object is explained by the nature of its central object, which has been found to ★ E-mail: luciano.cerrigone@alma.clbe an active symbiotic, whose companion star displays a WC8-type spectrum (Balick et al. 2014).This source provides us with a unique opportunity to study the physical processes taking place immediately before the birth of a planetary nebula.Multi-wavelength, multi-epoch observations (Sánchez Contreras et al. 2002;Sánchez Contreras & Sahai 2004;Sánchez Contreras et al. 2004;Pardo et al. 2004;Soria-Ruiz et al. 2013) have provided us with a complex picture of the source, consisting of: i) a large (∼20 ′′ ) molecular envelope made of the remnant AGB wind expanding at low velocity (∼17 km s −1 ); ii) multiple optical lobes, where shocked gas expands at high velocity (∼200 km/s); iii) a fast (up to ∼340 km s −1 ) bipolar molecular outflow along the polar axis; iv) a dense compact (∼1.5 ′′ ) molecular core surrounding the star and slowly (≤12 km s −1 ) expanding; v) a compact (about 0.2 × 0.4 ′′ ) HII region close to the central object, indicating the onset of ionisation in the envelope.
When the central star of a PPN becomes hot enough to ionize its CSE (T eff ≥20 000-30 000 K), radio continuum emission from the ionized gas can be detected in the centimetric range and this has been used to search for hot post-AGB stars where the ionisation has recently started (Umana et al. 2004;Cerrigone et al. 2011Cerrigone et al. , 2017)).CRL 618 was first detected at radio wavelengths by Wynn-Williams (1977) and was soon found to display increasing flux density over time, which was interpreted as due to the expansion of the ionized region (Kwok & Feldman 1981).Since then, several works have addressed the variability of the radio flux from this object, which displays a steady increase in the optically thick centimetric range and a more erratic behavior in the millimeter range, where the emission is optically thin (Martin-Pintado et al. 1988;Sánchez Contreras & Sahai 2004;Planck Collaboration et al. 2015;Sánchez Contreras et al. 2017).While the increase in flux density at frequencies where the emission is optically thick is interpreted as the expansion of the emitting surface, the erratic variability at frequencies of optically thin emission has been seen as an indication of on-going activity from a stellar post-AGB wind (Sánchez Contreras et al. 2017).
CRL 618 is also the only proto-PN studied in millimeter radio recombination lines (RRLs), which has allowed estimating a remarkable mass-loss rate of ∼8.4×10 −6 M ⊙ yr −1 for its post-AGB wind (Martin-Pintado et al. 1988;Sánchez Contreras et al. 2017).The existence of this still active ionized wind from the central star has also been invoked by Tafoya et al. (2013), who analyzed radio data of CRL 618 spanning years 1982-2007 and concluded that the ionisation of the circumstellar material started around 1971.Tafoya et al. (2013) also confirmed that the nebula is ionisation-bounded in the direction of its minor axis, indicating a much larger density in the equatorial direction than along the polar axis.
The distance to CRL 618 is still somewhat uncertain.Schmidt & Cohen (1981) estimated a distance of 1.8 kpc, assuming a main sequence luminosity of 2 × 10 4 L ⊙ for the B0 central star.Goodrich (1991) argued instead for a smaller luminosity of 10 4 L ⊙ , hence a distance of 0.9 kpc, although based on the radial velocity of the star and the Galactic rotation curve, they also mentioned a possible value of 3.1 kpc.Knapp et al. (1993) estimated a bolometric flux of 2 × 10 −7 erg cm −2 s −1 and quoted a distance of 1.1 kpc, for an assumed luminosity of 10 4 L ⊙ .Finally, Sánchez Contreras & Sahai (2004) measured proper motions in the high-velocity bipolar clumps that support the value of 0.9 kpc indicated by Goodrich (1991).
In the context of estimating the distance to a star, geometrical methods are the most reliable ones, when possible.One of these is the so-called expansion parallax, which allows deriving the distance by measuring the angular expansion of a structure such as a ring or shell, if its linear expansion velocity is known.Since the development of a PN implies the existence and expansion of an ionized shell, this method can be applied to this kind of sources even in their early phases, if their angular expansion and linear velocity can be measured.A version of this method was developed by Masson (1986) based on the comparison of radio interferometric visibility data sets and was later applied to several PNe (Guzmán-Ramírez et al. 2011, and references therein).

OBSERVATIONS AND DATA REDUCTION
We observed CRL 618 with the NSF Karl G. Jansky Very Large Array (VLA) telescope operated by the National Radio Astronomy Observatory (NRAO), while the array was in configuration A on 2011 June 12, at a frequency of ∼43 GHz (7 mm).The data were acquired within project 11A-171 (PI: G. Umana).The flux calibrator was 3C 138 and the complex-gain calibrator was JVAS J0414+3418.Two spectral windows were set up, spanning a total of 64 MHz in full polarization.The duration of the observations was 2 hours.
The data were reduced in CASA 5.6.1-8 (McMullin et al. 2007) with the VLA pipeline delivered with the software package, after correcting the intents 1 of the astrophysical sources in it: 3C 138 turned out not to have sufficient signal-to-noise ratio per channel to be used as bandpass calibrator, hence it was used only as flux calibrator, while 0414 + 3418 was used as both bandpass and complex-gain calibrator.After a first run of the pipeline, the data were inspected and some necessary manual flags were identified.The whole calibration 1 Intents are used by the pipeline to identify calibrators (in our case, pointing, bandpass, flux, and complex-gain calibrators) and science targets in an observation.was then repeated re-starting from the raw data and including the manual flags.Given the small band width of the data set, no correction was performed for the spectral slope introduced using the complexgain calibrator to calibrate the bandpass.
In this paper, we also make use of the data taken with the VLA on 1998-05-02 within project AW485 (PI: J. Wrobel).The data were acquired in continuum mode and span 100 MHz in frequency, centered around 43.3 GHz in full polarization.The array was in A configuration, but only a subarray of 12 antennas was used at 43 GHz, while a simultaneous subarray was used to observe in different frequency bands.The duration of the observation was about 8 hours.
This second data set was calibrated setting explicitly the flux scale of the phase calibrator JVAS 0443+3441 to the same values reported by Tafoya et al. (2013) 2 , since the flux calibrator in the original data is now known to be variable, hence not reliable.It must be noticed that the critical parameter for the analysis carried out in this work is not the absolute flux calibration in each data set, but the size of the emitting area, which is identified by a signal-to-noise ratio larger than 5, therefore it does not depend on the absolute calibration.Opacity and gain-curve correction were applied to the data, then a phase-only calibration table was generated, which was applied to generate the final amplitude and phase calibration.The reduction of the AW485 data was performed with the same CASA version used for the 11A-171 data set.
Both the 1998 and the 2011 data sets were self-calibrated only in phase, after the initial phase-reference calibration.This led to final flux densities of 0.98 ± 0.05 Jy beam −1 and 1.20 ± 0.06 Jy beam −1 in 1998 and 2011, respectively.The rms noise in each of the maps was instead of 0.4 mJy beam −1 in 11A-171 and 0.8 mJy beam −1 in AW485.
In Figure 1, we display the final  coverage for CRL 618 in the two data sets at 43 GHz.It is evident that the sampling is less dense in AW485 due to the smaller number of antennas, but at the same time, the longer integration time allows for a more uniform distribution over the plane.AW485 achieves about the same baseline lengths in both  and , which is expected to turn into a rounder synthesized beam.The 11A-171 data set achieves about the same baseline lengths in  but shorter in u.Although the  sampling is clearly different in the two data sets, these differences do not seem to be so substantial as to bias the final imaging products and their analysis, especially for values of || and || smaller than 3.5 M, which corresponds to a beam size of ∼0.06 ′′ .The differences can then be well mitigated through data weighting and beam convolution in the imaging step.In Figure 2, we display the dirty beams of the two data sets.
The 80th percentiles for the shortest and longest baselines in the two final data sets are respectively 3744 m and 15407 m in 11A-171 and 5377 m and 20408 m in AW485.The sampling in AW485 is therefore shifted towards longer baselines in comparison to 11A-171, due to the longer integration and different antenna distribution.To overcome the different angular resolution, we convolved the final products to a common beam with size matching both data sets.

NEBULAR EXPANSION
We performed our analysis of the expansion of the circumstellar ionized nebula in CRL 618 with the AW485 and 11A-171 data sets at 43 GHz following the method described by Masson (1986), to obtain a difference image from two radio interferometric data sets.
As explained above, the data sets were initially calibrated independently, then a round of self-calibration in amplitude and phase was executed on 11A-171, assuming as a model the image from AW485.As required by the radio parallax method (Masson 1986), this was done to align the two images on the same amplitude and phase reference, removing any offsets.
After this step, the AW485 data were subtracted in the visibility plane from 11A-171, using the CASA task uvsub.At this point, imaging of the difference data set was performed with tclean in CASA 5.6.1-8.In Figure 3, we display the images of the source at the two epochs and the difference obtained.The width of the convolution beam was set to 0.06 ′′ for both data sets.
The difference image (the bottom plane in Figure 3) displays a clear ellipse of emission.Along the direction of the major axis of the ellipse, two drops in emission at the eastern-most and westernmost locations can be seen (indicated by two arrows in the figure).This is compatible with the known difference in density between the material along the polar axis, which has been mostly swept away by the fast wind, and the material in the perpendicular direction, which is still an almost unaltered residual from the AGB phase (Sánchez Contreras et al. 2002).The inner region has negative emission (i.e., the older image has larger flux density in that area), as expected due to the expansion.

Expansion-parallax distance
After detecting the expansion in the difference image, we estimated the distance to our target as described by Guzmán-Ramírez et al. (2011).
We started from the image of CRL 618 obtained with the AW485 data set, then created a set of expanded images, each of them with an expansion factor of 1 + , where  ranged from 0.1 to 0.3 in steps of 0.003.From every expanded image we then subtracted the AW485 one, resulting into a grid of model difference images.Every model difference image thus obtained was compared to the real difference image pixel by pixel, calculating a  2 value where  ij is the value of the pixel with coordinates  and  in the model difference and  ij in the real difference, and  is the standard deviation of the  ij −  ij values over the whole map.
For a direct comparison, we display in Figure 4 the contours of the real difference image in black and those of the best-model difference image in red.
The values of  2 obtained are plotted in Figure 5 as a function of .The minimum was found by fitting a cubic curve in the range where 0.12 ≤  ≤ 0.2 (Figure 6) and then calculating its analytical minimum. 2 was then found to be minimum at  = 0.160 ± 0.003, where the error is the statistical error on the fit to the curve of  2 values.
The distance to the object can be calculated with the following   5,5,7,9,11,13,15) and rms of about 0.8 mJy beam −1 and 1 mJy beam −1 , respectively.Negative levels are displayed as dashed lines and positive ones as solid lines.formula (Guzmán-Ramírez et al. 2011): where  exp is the expansion velocity and  is the angular expansion rate in milliarcsec (mas) per year, which is equivalent to   mas 13.1205 yr (3) where  is the expansion factor that we have determined,  the angular radius at the time of the AW485 observation in mas, and 13.1205 yr is the time between the AW485 and 11A-171 observations.
To obtain an estimation of the size of the nebula in 1998, we first rotated the AW485 image by 2.5 degrees in the direction North to East, to align its main axis with the E-W direction, then slices of the nebula were taken in the N-S and E-W directions separately, measuring the length of every slice (amount of pixels above 5 multiplied by pixel length) and summing the flux density values along each of them.Finally, a flux-weighted mean size was calculated as         where  is the length of every slice and  the total flux from it.This returned a size of 0.34 ′′ × 0.65 ′′ , with a relative error of about 10%, considering the error of 5% on the flux density measurements.We can then use the length of the semi-axes to derive the angular expansion along the polar direction and orthogonally to it.This translates according to Eq.3 into  = 2.1 ± 0.2 mas yr −1 along the minor axis and 4.0 ± 0.4 mas yr −1 along the major axis, between years 1998 and 2011.These compare well with the values of 2.3 ± 0.6 mas yr −1 and 4.7 ± 1.1 mas yr −1 found by Tafoya et al. (2013) and averaged over time until 2007.
In the derivation of the distance from Eq. 2, a critical value is the velocity.It has been shown that expansion-parallax distances need to be corrected, if the velocity is obtained from optical spectral lines, due to the different velocities of the material traced by the radio continuum and the lines of ionized elements (Schönberner et al. 2018).To avoid this complication, we make use of the velocity derived from radio data by Martin-Pintado et al. (1988) and Sánchez Contreras et al. (2017), who estimate that the HII region expands at ∼20 km s −1 by modelling its radio recombination lines (namely, H30, H35, and H41), thus tracing the overall motion of the ionized gas, without abundance or line-excitation biases.As the expansion occurs mainly along the major axis, we associate this velocity (with a 10% error) to the expansion rate in this direction and thus obtain a distance of 1.1 ± 0.2 kpc.Though slightly larger, this value compares well with that of ∼900 pc derived from the expansion in the lobes, seen in the optical images of CRL 618 (Sánchez Contreras & Sahai 2004).

CONCLUSIONS
We have analyzed data of CRL 618 obtained in 1998 and 2011 at 43 GHz (7 mm) with the VLA.The 2011 data are presented here for the first time.A size increase of the nebula is immediately evident by eye in the images from the different epochs.The expansion was estimated by imaging the data after subtraction in the visibility plane and comparing this to a grid of model difference images obtained expanding the image from 1998.This assumes that the expansion is self-similar, which is not strictly true for CRL 618, because its major and minor axes have been found to expand at different rates.
Nevertheless, our analysis indicates that at first approximation this can be neglected and the method still returns reliable results.The expansion rates that we find are in fact compatible within errors with previous independent estimations: we find 2.1 ± 0.2 mas yr −1 along the minor axis and 4.0 ± 0.4 mas yr −1 along the major axis, leading to a distance to CRL 618 of 1.1±0.2kpc.While previous estimations of the distance to CRL 618 were derived from an assumed intrinsic luminosity of the source, this is its first direct measurement.The present value matches within errors with that of 0.9 kpc, which has been widely adopted in the last twenty years, and coincides with that given by Knapp et al. (1993) for a luminosity of 10 4 L ⊙ .

Figure 1 .
Figure 1. coverage of data sets 11A-171 and AW485 after calibration and flagging.

Figure 3 .
Figure 3. Images of CRL 618 in Q band from projects AW485 (top) and 11A-171 (middle); contours start at 5×rms and increase by steps of 5×rms (−5×rms is also plotted as a dashed line); the convolution beam of 0.06 ′′ is plotted in the bottom left; the rms noise is 0.8 and 0.4 mJy beam −1 for the AW485 and 11A-171 data, respectively.The bottom image displays the difference of the two data sets, linearly stretched between −5×rms and its maximum, with two arrows pointing at drops in flux density along the shell.

Figure 5 .
Figure 5. Distribution of  2 (as defined in Eq. 1) as a function of  (the expansion factor minus one).

Figure 6 .
Figure 6.Cubic curve fitted to the  2 values in the range 0.2 ≤  ≤ 0.225.