The Population of Galactic Planetary Nebulae: a Study of Distance Scales and Central Stars Based on the Second Gaia Release

We matched the astrometry of central stars (CSs) of spectroscopically confirmed Galactic planetary nebulae (PNe) with DR2 Gaia parallaxes (p), finding 430 targets in common with p > 0 and ∣ σ p / p ∣ < 1 . A catalog of PNe whose CSs have DR2 Gaia parallaxes is presented in Table 1. We compared DR2 parallaxes with those in the literature, finding good correlation between the two samples. We used PNe parallaxes to calibrate the Galactic PN distance scale. Restricting the sample to objects with 20% parallax accuracy, we derive the distance scale log ( R pc ) = − ( 0.226 ± 0.0155 ) × log ( S H β ) − ( 3.920 ± 0.215 ) , which represents a notable improvement with respect to previous ones. We found that the ionized mass versus optical thickness distance scale for Galactic PNe is not as constrained by the Gaia calibrators, but gives important insight into the nature of the PNe, and is essential to define the domain for our distance scale application. We placed the CSs whose distance has been determined directly by parallax on the HR diagram, and found that their location on the post-asymptotic-giant-branch (AGB) H-burning evolutionary tracks is typical for post-AGB stars.


Introduction
Planetary nebulae are essential probes of stellar evolution, being the gaseous and dusty remnant of asymptotic giant branch (AGB) evolution. When studied in populations, they are effective tracers of kinematic and chemical changes in the parent galaxy. Studies of PNe in the Milky Way have always been hampered by the difficulty of measuring their distances. Individually, only a handful of PNe have well-determined distances (see Section 2). The rest of the Galactic PN distances have always been estimated through distance scales (e.g., Daub 1982; Cahn et al. 1992, hereafter CKS;Stanghellini et al. 2008, hereafter SSV;Frew et al. 2016, hereafter F16), where a few PNe with known distances have been used to constrain a physical relation between pairs of PN physical parameters, one of which to be distance-dependent. We have studied distances to the CSs of Galactic PN using TGAS parallaxes after the first Gaia release (DR1), but the PN sample was very limited in DR1, and so was the distance scale study (Stanghellini et al. 2017, hereafter Paper I).
The release of the second Gaia catalog (DR2) has prompted us to study the parallaxes p of CSs of spectroscopically confirmed PNe, and to review the most widely used PN distance scales using Gaia parallaxes as calibrators. Planetary nebula catalogs and Gaia parallaxes have been matched by other authors since the DR2 release: the Gaia Collaboration et al. (2018), in the context of showing HR diagrams with Gaia targets, has examined a selection of a few CSs of PNe on the stellar evolutionary diagram, filling the evolutionary gap between AGB and WD stars. Kimeswenger & Barría (2018) have compared distances derived from DR2 to other distances and found reasonable agreement between the sets for short distances. González-Santamaría et al. (2019) have used Gaia parallaxes to estimate physical radii and expansion ages of PNe.
The main goal of the present study is to use DR2 central star parallaxes to calibrate the Galactic PN distance scale.
In Section 2 we summarize our search for Gaia matches to the CSs of Galactic PN. In this section we also compare Gaia DR2 parallaxes with other existing, reliable parallaxes and distances available to date, including some recent distances that had not been available at the time of previous distance-scale calibration, such as the expansion distances from HST images by Schönberner et al. (2018, hereafter SBJ).
In Section 3 we calibrate the commonly used PN distance scales with the Gaia parallax measurements. In Section 4 we present a limited study of the CSs and their location on the HR diagram, again using Gaia parallaxes. In Section 4 we discuss the results of our study and future directions.

The Galactic PN Sample with Reliable Parallaxes from Gaia
We initially matched the J2000.0 coordinates of spectroscopically confirmed Galactic PNe (Acker et al. 1992;Kerber et al. 2003) with Gaia DR2 catalog and found 655 unique Gaia targets that match the CSs with very high confidence. In this sample there are 497 stars with s < p 1 p | | , while the sample with s < p 1 p | | and p>0 consists of 430 CSs. Table 1 (published electronically in its entirety) includes all Galactic PNe whose CSs have a corresponding DR2 parallax. The columns containing the following information: column(1): PN name as in Simbad (the PN G number); column(2): Gaia DR2 ID of the coordinate match; column(3): DR2 parallax and its uncertainty, in milliarcseconds; column(4): nebular morphology; column(5): angular radius in arcseconds; columns(6) and (7): logarithmic Hβ flux and extinction constant and their uncertainties in cgs units; column (8): linear radius of the nebula in parsecs, calculated from the angular radius and parallax; column(9): references for the ancillary data, as explained in the table note.
Obviously, the sample of Galactic PNe with Gaia parallaxes represents an incomplete sample, which markedly declines beyond ∼2.5 kpc. This result should be taken into account when using this sample to derive population-broad results.

Parallaxes of Galactic PNe from Gaia Compared to Other Independent Parallaxes and Distances
There are 25 PNe in DR2 whose parallaxes or distances had been previously measured with other reliable methods. In Figure 1 (upper panel) we compare directly DR2 parallaxes with other trigonometric parallaxes (Harris et al. 2007, hereafter H07), and also with parallaxes derived by inverting other independent distances. Parallaxes used in Figure 1 are given in Table 2. Distances by SBJ have been measured via multi-epoch expansion parallax and Hubble Space Telescope (HST) images, but suffer from indeterminate expansion velocities of the PN shells. Distances via spectroscopic parallaxes, such as those of Ciardullo et al. (1999, hereafter C99), suffer from model dependency.
The correlation between the DR2 and independent parallaxes is tight, with linear correlation index R xy =0.98. In the lower panel of Figure 2 we show the residuals, whose average is 0.22 for the whole sample. This comparison represents a positive assessment of Gaia parallaxes for CS PNe, making them ideal probes to calibrate the distance scale.

Galactic PN Distance Scale Calibrated with Gaia Parallaxes
We aim to use the DR2 parallaxes to study and calibrate the statistical distance scales for Galactic PNe. Statistical distance scales for Galactic PNe are based on the physical relation between distance-dependent and distance-independent parameters that are commonly measured or derived in PN studies. Such relations, once calibrated on PNe with known distances, are used to infer the statistical distances to all other PNe, whose distance is not known. As seen in the previous section, before the release of Gaia DR2 there were only 25 PNe whose distances were reasonably well known, half of whose became available only recently, thus the calibrator sample was really limited. The DR1 Gaia release included a handful of calibrators, which we examined in Paper I, with the conclusion that the physical radius versus surface brightness statistical distance scale was promising. The DR2 Gaia release gave us high motivation to re-explore PN distance scales. The next section includes a study on the recalibration of the Galactic PN distance scales based on DR2 parallaxes. As in Paper I, we  Notes. a Morphological class codes are: 1: Round, 2: Elliptical; 3: Bipolar Core; 4: Bipolar; 5: Pointsymmetric. b The three number sequence corresponds to the references for the angular radius θ, the flux at b H , and the extinction constant, respectively. An asterisk ( * ) marks the missing reference. The numbers correspond to the following references: (1) this study, (2) Stanghellini et al.  studied the distance scales relating (i) the physical radius to surface brightness and (ii) the ionized mass to optical thickness.

The Physical Radius versus Surface Brightness Distance Scale
This scale is realized by the assumed linear relation between log(R pc ) and log(Sb Hβ ), representing the distance-dependent and distance-independent parameters, respectively. The physical radius in parsecs, R pc , is computed as θ/(206265×p), where θ is the angular radius and p the DR2 parallax, both in arcseconds; the Hβ surface brightness is computed as The sample of Figure 2, consisting of 112 targets, shows the distribution in this plane of PNe whose parallax is available from DR2 with accuracy better than 20%, and whose physical parameters used for the distance scale are also available (see Table 1). The error bars in Figure 2 reflect the propagation of uncertainties of all parameters. With this pruning, we retain a quite sizable distribution of calibrators whose parallax uncertainties compare favorably with the ones expected from the cosmic scatter of the statistical scale (see Buckley & Schneider 1995;Schneider & Buckley 1996), making them physically more revealing.
By limiting s p p | | we could, in principle, introduce a Lutz-Kelker-type bias that would principally affect the upper left end where larger distances play a role. However, in the present analysis we did not take it into consideration, given that our calibration is also affected by the intrinsic scatter of the ( )-( ) plane. All these effects could be better investigated using future Gaia data releases.
A simple weighted least-squares fit in presence of asymmetric error bars would not suffice for this distance scale. The issue here is that we measure the parallax (of the CSs), the PN radius, and the PN line flux and extinction for the calibrators, and uncertainties, but we do not directly measure distances, physical radii, and surface brightness. Thus we approached the problem by utilizing only the measured quantities and their original uncertainties, to calibrate the distance scale. An optimal way to perform this calibration is a maximum likelihood approach, which is illustrated in the Appendix, and used to estimate the parameters of the PN statistical scale (see Kelly 2007).
In Figure 2 we show the maximum likelihood result as a thick solid line, corresponding to the following fit: The distribution of PNe around this linear scale presents considerable scatter. This parameter distribution-not observed in the TGAS sample (Paper I) due to the paucity of objectsbesides depending to some extent on the scatter in the measured PN parameters, is empirical evidence that not all evolving PNe follow a unique physical correlation on this plane.
To gain insight into the actual line (or strip) that PNe follow in their evolution, we needed to analyze parameter dependency. In Figure 3 we show the same data distribution of Figure 2, but we examine specific groups of PNe. We differentiated PNe with small  (1999). G86: extinction distances, Gathier et al. (1986). , lower right, see the next section for a definition). These are the critical PN groups for which the errors may be higher than is reported in the literature. We found that by eliminating these extreme groups of PNe, the correlation of the distance scale is tighter. The parameter limits selected for Figure 3 are quite arbitrary, and were used to recognize that the PN distances derived from the physical radius versus surface brightness relation (Equation (1)) may have larger uncertainty for PNe observed very early in their evolution (smaller radii), or large extinction, or low Hβ fluxes. The distribution indicates that the relation is poorly followed by PNe with very low ionized masses, like the ones that are still in an ionization-bound state.
It is worth noting that the angular radii available in the literature have been measured in different ways. The best way to measure angular radii is by photometry, defining the photometric radii as the one encompassing 85% of the total nebular flux. This method has been used almost exclusively when HST images of the PNe are available. To date, the only homogeneous set of photometric radii of Galactic PNe published so far being the one by Stanghellini et al. (2016), which is limited to small PNe, unresolved from the ground. Given the importance of a larger data set with photometric measured data, we measured the PN sizes of an additional sample of PNe whose HST images are available on the data archive. We have performed aperture photometry on these images at various distances from the center, which corresponds to the location of the central star or to the geometrical center of the PN if a central star is not seen in the image. We then measure the total flux of the PN with broad aperture photometry, and find a radius that encircles 85% of total flux. The technique used here is identical to that of Stanghellini et al. (2016). The resulting radii are given in Table 1 with reference code 1.
In Figure 4 we plot the physical radius versus (Hβ) surface brightness locus for the s < p 0.2 p | | sample, i.e., the same data set of Figure 2, where we show the sample of extended PNe with photometric radii measured from HST images, which exhibit a tight correlation on this plane (R xy =−0.92). The PN sample with photometrically measured radii is too small to derive accurate coefficients through maximum likelihood analysis, and thus we will not give its fit, but we can see from Figure 4 that these eight targets define a much tighter correlation than the general sample, showing the importance of measuring photometric radii. Note that these PNe have different morphologies. With the final Gaia release, and with future HST observations, there may be more of these calibrators available for the distance scale.
We explored the distance scale with calibrator PNe of different morphological types. From the data set of Figure 2 we have selected PNe that have been morphologically classified, and plotted them in Figure 5 on the same plane as in previous figures, omitting error bars for clarity. In this paper we used the classification scheme by Stanghellini et al. (1993) with Round, Elliptical, Bipolar Core, Bipolar, and Pointsymmetric PNe. Morphological classifications are given in Table 1, based on images published by Manchado et al. (1996), Schwarz et al. (1992), and Balick (2007). The sample of Figure 5 includes 16 Round PNe, with linear correlation coefficient R xy =−0.97 in the scale plane. This was expected from the result of Figure 4, since the photometric and geometric radii correspond in Round PNe.
With the maximum likelihood method we determine the fit to Round PNe with s Our study assumes that σ θ /θ=0.2. Errors on radii measurements are generally not published, and we used this conservative estimate in our likelihood analysis of the scale parameters. Biases in these measurements can arise if the PN is not spherically symmetric. In fact, our Figure 5 shows that Round PNe are much less dispersed about the scale relation than the rest of them. It would be difficult, with the sample at hand, to quantify the influence of such possible biases; the differences between the fits of Equations (1) and (2) could be attributed partly to such an effect, though they are still compatible with the stated statistical uncertainties on theta.
We evaluate the accuracy of the scales represented by Equations (1) and (2) with the method presented by Smith (2015, hereafter S15). We calculate the distance ratios K 3 by multiplying distances from the distance scales-excluding in all cases PNe with m <log 2 ( ) , which do not follow our physical scale-to DR2 parallaxes. We used Equation (1(a)) in S15 to estimate the variance of the distance ratios, s K 2 , and Equation (5) for the typical relative errors, α s . We used the scale distances d s and DR2 parallaxes p±σ p directly in Equation 1(a), while we estimated s s 2 for the scale distances by propagating the formal errors of Equations (1) and (2) and also accounting for correlations in the observed parameters.
In Table 3 we give the results of our accuracy evaluation for Equations (1) and (2), both compared to DR2 parallaxes, for m >log 2 ( ) . We also calculate these statistics for Round PNe exclusively, for both scales, to determine the best one to use in each case. In the first three columns of Table 3 we give the where we indicate with filled symbols the extended PNe with photometric radii from this study (see Table 1), and with crossed symbols the PNe with m <log 2 ( ) . The straight line represents Equation (1).
Equation ( scale used for distances, the sample studied, and the sample population. Column(4) gives the typical error of the scale, α s ; columns (5) and (6) contain the average distance ratio, á ñ = á´ñ K d p s , and the average variance of the distance ratio, s á ñ K . In this analysis, the best distance scale is with K as close as possible to unity, and with the smallest σ K , since σ K represents the fractional uncertainty in the distances from the studied scale.
If we compare the scales with parallaxes without considering the nebular morphology, we find that Equation (1) gives a slightly short scale (−2%) with respect to Gaia parallaxes, while Equation (2) gives a slightly long scale (+7%). It makes more sense to assess the scales with the best parallaxes, i.e., s < p 0.2; p | | in these cases, s á ñ~0.25 K . By limiting the comparison to Round PNe both distance scales are very accurate.
The best results in Table 3 are those obtained from comparing Equations (1) and (2)  . Also, it is worth noting that Equation (1) works as well for the sample used to define it as for Round PNe with secure parallaxes. Finally, in principle, the mean of Equation (1) and Equation (2) applied to Round PNe with secure parallaxes will be almost entirely free of bias.
In terms of average distance ratios, i.e., by comparing the distances from a statistical scale with DR2 parallaxes, both Equations (1) and (2)  We also performed a direct comparison between our scales and other widely used distance scales by using the = S S k d d i i y,i x,i estimator by Phillips (2002); see also Equation (2) in S15; note that k is different from the above K ), where d y,i is the distances of the ith object in the old distance scale, and d x,i is the distance of the same object in our scale. By comparing F16 distances with our scales we found k=1.243 and 1.165, respectively, for Equations (1) and (2). Both comparisons are based on the 103 PNe in common. The same comparisons with the SSV distances yielded k=1.136 and k=1.068. The scale of Equation (1) should be used when DR2 Gaia parallaxes are not known.

The Ionized Mass versus Optical Thickness Distance Scale
In Figure 6 we show the relation between the ionized mass and optical thickness with PNe with low parallax uncertainties (i.e., s < p 0.2 p | | ). A similar plot has provided the basis for the distance scales of Galactic PNe by CKS, and by SSV, where the former has been calibrated with Galactic PNe, and the latter with Magellanic Cloud PNe. Note that the domain of calibrators in those published scales was limited to the parameter space defined by the small rectangle in the figure.
The concept of this scale is that PNe evolve from optically thick to optically thin, and once they reach optically thin (or density bounded) status, their ionized mass is constant. This scale includes the strong assumption that all PNe have identical ionized mass, and that they evolve the same way from optically thick to optically thin.
The abscissa, t q =´F log 4 2 ( ) , is a measure of the inverse optical thickness of the nebula, where θ is the angular radius and F is the 5 GHz flux of the PN. The ordinate, m q =´-p F 2.266 10 21 5 3 0.5 ( ), the ionized mass, is the distance-dependent term that has been calculated with p (the Gaia parallax in arcsec) and the 5 GHz fluxes for the data points of Figure 6 (5 GHz fluxes from CKS; Acker et al. 1992).
The optically thick sequence in the figure is represented by the sloping segment of the solid and broken lines, while the optically thin sequences are the horizontal lines. Neither sequence is well defined by the Gaia calibrators (i.e., when using p in the calculation of the ionized mass).
The optically thin sequence is reasonably reproduced by Round PNe, identified in the figure as filled symbols. The average dispersions of log(μ) are −0.931±0.72 and −0.81±0.39, respectively, for the whole sample and the Round PN sample. This finding agrees with the spherical assumption of the ionized mass term in CKS's calibration and shows that this scale may not be accurate for other PN morphologies. The data set is too small to make a meaningful calibration of this distance scale based on photometric radii.
We note that PNe with m <log 2 ( ) seem to lie in a parallel sequence where the ionized mass is underestimated, or they simply evolve differently from the others. If we plot the PNe with m <log 2 ( ) on the R log pc ( ) versus log(Sb Hβ ) scale we find that they are all located far from Equation (1) (see Figure 3, lower right panel).
PNe with such low ionized mass are a minority, representing ∼20% of the whole population. We explored the parameter space for these PNe. We showed in Figure 3 that several of them have a combination of small apparent radii, low fluxes, and high extinction, although not all high-extinction PNe have low ionized mass. Filling factors for these PNe cannot be calculated, due to the lack of electron density and temperature information in the literature. The distance scale of Equation (1) Figure 6. Ionized mass vs. optical thickness relation. The solid line is the CKS distance scale, while the broken line is the SSV scale, calibrated on Magellanic Cloud PNe. The box represents the extent to which the relation was populated when the CKS and SSV scales were derived. The solid circles are the PNe with Round morphology. calculated exclusively for PNe whose m >log 2 ( ) is within the 1σ of Equation (1) (i.e., within the shaded area of Figure 2).
Since optical thickness varies with metallicity (see SSV) we also inspect their metallicity relative to solar through their oxygen abundances (from Stanghellini & Haywood 2018), and found that most of these targets are oxygen-rich (supersolar). Given that the low-ionized mass PNe are not exclusively supersolar, since several targets with m >log 2 ( ) are also supersolar, we could not conclude this to be the only factor causing the targets to stray from the distance-scale correlations. Additional analysis of these targets is in order for more insight.
Finally, we noted that low-extinction PNe exclusively populate the optically thin sequence. Note that a similar result was found by Kimeswenger & Barría (2018); in fact, by comparing distances from Gaia DR2 parallaxes and those by SSV they found that reddened objects are comparatively worse than blue objects for targets with distances <4 kpc.

Populating the HR Diagram with CSs of Galactic PNe
The Gaia DR2 allowed, for the first time, determination of the CS distances to a considerable sample of Galactic PNe. In order to study these CSs on an HR diagram, we selected published Johnson V magnitudes and effective temperatures of the CSs from the literature, matching the search with PNe with known p. For this initial study, we limited the CSs to those with V magnitudes taken from Tylenda et al. (1991, hereafter T91). We also selected only He II Zanstra temperatures from a group of references that calculated them with the same method (Kaler 1983;hereafter K83;Shaw & Kaler 1989;hereafter SK89;and Stanghellini et al. 2002;hereafter S02). If the magnitude was not available in T91 we used the ones given by K83 or SK89.
In Figure 7 we show the CSs of Galactic PNe on the T L L log log eff ( )-( )  plane. In order to gauge their evolutionary stage and their mass, we superimpose the evolutionary tracks are from Vassiliadis & Wood's (1994) hydrogen-burning post-AGB tracks for solar metallicity (Z=0.016) and CS masses of 0.57, 0.6, 0.63, 0.68, 0.75, and 0.9 M e , corresponding to turnoff masses of 1.0, 1.5, 2.0, 2.5, 3.5, 5.0 M e .
In Figure 7 we also indicate the locus of the HR track by Miller Bertolami (2016), with an initial mass of 1.0 M e , the same as the lowest initial mass in the Vassiliadis & Wood (1994) sample. The slight difference in the 1.0 M e tracks is due to different CS (final) mass, which is 0.53 M e in the Miller Bertolami case, as opposed to 0.57 M e in Vassiliadis and Wood's track. This difference is due to differences in the massloss treatment. Since identical post-AGB mass tracks do not exist, a perfect comparison of different track sets is not possible, but we do not expect that this would change any of our conclusions for this section.
We only plot those stars whose parallax and effective temperature uncertainties are smaller than 20%. The input effective temperatures (column 3), magnitudes and their quality, if available, as in T91 (column 4), luminosities (column 5), and derived CS mass (column 7) of the plotted stars are listed in Table 4. In column (8) of the same table we give the reference code, with the first code referring to the effective temperature and the second code referring to the magnitude.
The sample of 39 CSs represented in Figure 7 is not a complete sample, and it is skewed toward nearby PNe, both for the choice of using Gaia parallaxes and for the selection of available magnitudes. This study is an assessment of the realm of CSs and their evolution, to confirm that Gaia distances produce a reasonable distribution of stars in the evolutionary stage that we think PN nuclei should be at, compared to the post-AGB evolutionary tracks. Without any further selection or other assumptions, we found that most of these CSs are correctly encompassed by the classical evolutionary tracks for H-burning post-AGB stars.
It is possible that after the final Gaia release and a careful analysis of their temperatures and luminosities, there will be enough parallaxes to perform a similar study based on a CS sample that is complete within a representative distance (e.g., all PNe within 5 kpc), or for a luminosity-limited sample. At this time, this is just an assessment of the quality of the distances, and a confirmation of the actual post-AGB nature of the CSs.

Discussion and Outlook
We presented a parallax catalog of spectroscopically confirmed, Galactic PNe (Table 1), where we matched their CSs with the DR2 Gaia catalog. This catalog also includes published and newly measured parameters for the PNe with Gaia parallax counterparts, including morphology, angular diameters, H β fluxes, extinction constants, and linear radii.
We compared the parallaxes to all reliable individual PN parallax and distances available in the literature, to determine the accuracy of other methods. We found good agreement between the two sets of parallaxes, including those derived by inverting other accurate distances.
We then used the Gaia parallaxes to calibrate the most commonly used distance scales and examine under which conditions these calibrations would hold. The statistical method on which our calibration is based uses the Gaia parallaxes and their uncertainties directly, with their assumed Gaussian distribution, naturally incorporating negative parallaxes.
We determined a distance scale (Equation (1)), based on the correlation between the nebular physical radius and its surface brightness, that yields to very accurate distances for PNe whose ionized mass is not extremely low. The accuracy in terms of average distance ratio K, which is the average of distances from the scale multiplied by Gaia DR2 parallaxes, with K=1 for scales that reproduce the DR2 parallaxes exactly, and σ K representing the fractional uncertainty in the distances from the studied scale, is s á ñ  á ñ =  K 0.9478 0.2483 K for the s < p 0.2 p | | sample. The scatter of Round PNe on the log(R pc ) versus log(Sb Hβ ) plane is smaller than that for the general sample. This is important for assessing the distance scale, since it proves that the basic relation is sound from a physical point of view for PNe that are not compact, faint, or have high extinction and low ionized mass. The scale calibrated on Round PNe exclusively is excellent, but does not improve on the Equation (1) scale.
We also found that the ionized mass versus optical thickness distance scale is not well constrained by Gaia parallaxes; we determined that Round PNe in the optically thin sequence of this scale define the sequence fairly well. The scale has provided guidance for parameter limits of the surface brightness-physical radius scale, but ultimately does not work as well as a distance scale for Galactic PNe.
As a working conclusion, we found that PN distances whose parallaxes are not available from Gaia or from other parallax measurements, should preferably be derived with the physical radius versus (Hβ) surface brightness statistical distance scale of Equation (1). If the PN is not spherically symmetric, an effort should be made to derive the PN radius photometrically on the emission-line image. The scale of Equation (1) represents a considerable improvement over previously used statistical distance scales for Galactic PNe.
We located the CSs of the PNe on the HR diagram using Gaia parallaxes, and effective temperatures and magnitudes from the literature. By comparing their location to post-AGB evolutionary tracks, we determined their masses, as shown in Table 4. The average mass recovered from this incomplete and inhomogeneous sample is ∼0.62 M e . PNe and their CSs with  Note. Reference codes are as follows: K83: Kaler (1983); T91: Tylenda et al. (1991); S02: Stanghellini et al. (2002); SK89: Shaw & Kaler (1989). We performed a numerical solution to this problem, for a grid of a and b values to build the bivariate posterior probability function, p(a, b), from which we determined the maximum likelihood estimates of a and b and their confidence intervals. As an example, Figure 8 shows the computed likelihood plot for the case of s < p 0.2 p | | .