Astronomy & Astrophysics manuscript no. AApaper6 January 21, 2009 (DOI: will be inserted by hand later) Σ − D The relation for planetary nebulae D. Uroˇsevi´c1,5, B. Vukoti´c2, B. Arbutina1, D. Ili´c1,5, M. Filipovi´c3, I. Bojiˇci´c4, S. Sˇegan1, and S. Vidojevi´c1 9 1 DepartmentofAstronomy,FacultyofMathematics, Universityof Belgrade, Studentskitrg16, P.O.550, 11000 0 Belgrade, Serbia 0 2 Astronomical Observatory,Volgina 7, 11160 Belgrade 74, Serbia 2 3 University of Western Sydney,Locked Bag 1797, Penrith South DC, NSW 1797, Australia n 4 Department of Physics, Macquarie University,Sydney,NSW 2109, Australia a 5 Isaac Newton Instituteof Chile, Yugoslavia Branch J 5 1 Abstract. We present an extended analysis of the relation between radio surface brightness and diameter – the so-called Σ−D relation for planetary nebulae(PNe).Wereviseourpreviousderivation of thetheoretical Σ−D ] R relation for the evolution of bremsstrahlung surface brightness in order to include the influence of the fast wind fromthecentralstar.Differenttheoreticalformsarederived:Σ∝D−1 forthefirstandsecondphasesofevolution S andΣ∝D−3 forthefinalstageofevolution.Also,weanalyzedseveraldifferentGalacticPNsamples.Allsamples . h areinfluencedbysevereselectioneffects,butMalmquistbiasseemstobelessinfluentialherethaninthesupernova p remnant (SNR) samples. We derived empirical Σ−D relations for 27 sample sets using 6 updated PN papers - from whichanadditional21newsetswereextracted.Twentyfourofthesehaveatrivialformofβ≈2.However, o we obtain one empirical Σ−D relation that may be useful for determining distances to PNe. This relation is r t obtained by extracting a recent nearby (< 1 kpc) Galactic PN sample. s a [ Key words. planetary nebulae: general – radio continuum: ISM – methods: analytical – methods: statistical – ISM: supernovaremnants 1 v 4 3 1. Introduction Observational improvements, including the increased 3 use of radio interferometers, account for the high number 2 Planetarynebulae(PNe)areusuallyidentifiedbytheirop- ofresolvedradioPNeoverthepasttwodecades.TheΣ−D . tical spectra, consisting mainly of recombination lines of 1 relationfor PNe has only beensporadicallydiscussed(for 0 hydrogen, helium and collisionally excited light elements. example, Amnuel et al. 1984) until Uroˇsevi´cet al. (2007) 9 Theyarecomposedofgaseousshellsionizedbyahotcen- (hereafter Paper I). 0 tralstar,whichallowforthenecessaryconditionstocreate : The statistical analysis of the Galactic PN v theirspectra.PNewerefirstdiscoveredover200yearsago. distances was presented in Daub (1982) and in i Today,therearemorethan 2500confirmedGalacticPNe X Cahn, Kaler & Stanghellini (1992). Using radio data, withtotalnumberestimatesatleasttentimeshigher.PNe r the primary statistical method for determining distance emission from radio to X-ray have been detected overthe a was based on the correlation between PN radius and past 30 years. Several hundreds of these PNe have been brightness temperature, known as the R−T relation1. observedinradio-continuumalone.Theirnumberfarout- b Van de Steene & Zijlstra (1995) derived an empirical rank the known Galactic supernova remnant (SNR) pop- R−T relation for 131 Galactic PNe, felt to be located ulation (approximately 250). b near the Galactic center, at approximately the same The relation between radio surface brightness and di- distance. Additionally, a smaller sample of 23 (non- ameterforSNRs,theso-calledΣ−Drelation,hasbeenthe Galactic bulge) PNe were studied havingwell-determined subject of extensive discussions over the past forty years. distances.Thesamecalibrationrelationswereobtainedin This Σ−D relation has a power law form: both cases. Zhang (1995) studied the statistical distance Σ∝D−β, (1) scaleforPNe basedonanothersample of132PNehaving well-determined individual distances. He found that whereΣistheradiosurfacebrightness,D isdiameterand the derived distance scale was in good agreement with β the slope of the log-log plot. those derived by Van de Steene & Zijlstra (1995). More Send offprint requests to: Dejan Uroˇsevi´c, e-mail: [email protected] 1 ThisrelationisessentiallyequivalenttotheD−Σrelation. 2 D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae recently, Phillips (2002) (hereafter Ph02) derived R−T radiationmechanismusedforderivationofthetheoretical b relations for 44 nearby PNe whose distances are less than Σ −D relation is thermal bremsstrahlung, which is the 0.7 kpc. Thus, different samples of Galactic PNe with basic process of radio-continuum emission in Hii regions. knowndistancesweredefinedinthethreepreviouslycited We begin our derivation using the following notation: papers.All of these derivedR−T relations were applied b to PNe with unknown distances. We also note that R – shell radius, s a recent paper by Stanghellini, Shaw & Villaver (2008) (hereafter SSV) gives an updated reliable sample of R – outer radius of the shell, o individual Galactic PN distances. The statistical analysis of Galactic radio PNe is less R – inner radius of the shell, i developedthanthoseofGalacticSNRs2.Therefore,selec- tion effects are felt to greatly influence Galactic radio PN V – velocity of the shell, s samples causing the statistical determination of distance to be highly uncertain. V – AGB wind velocity, In the present paper, we use the methodology for the v – fast wind velocity, statistical analysis of SNRs to improve that for Galactic radio PNe. Primary objectives include the following: M˙ – mass loss in the AGB phase, (i) Development of an improvedform of the theoretical Σ−DrelationforPNe(comparedtotheonederived m˙ – fast wind mass loss, inPaperI)byincorporatingtheinteractionbetween asymptotic-giant-branch (AGB) star wind and the Ms – mass in the shell, fast wind from the central star. (ii) Derivation of improved empirical Σ − D relations ρs – density of the shell, using current and previous samples to determine if these are all under the influence of selection effects. t=0 – the moment when the AGB wind stops, (iii) Validation of Σ−D relations for the determination of radio PN distance. t=τ – the moment when the fast wind starts. 2. Dynamics of PNe 2.1. The momentum conserving phase In considering the dynamics of PNe, we adopt the in- If we assume momentum conservationduring the interac- teracting stellar winds model (ISW; Zhang & Kwok 1993, tion, then: Kwok 1994, and references therein). This model implies that a nebular shell results from the sweeping up of cir- P =M V =P +P , (2) cumstellar material moving within the AGB star wind by s s s AGB FW afastwindemanatingfromthePNcentralstar.Thespeed ofthefastwindis100timesthatoftheAGBwindandacts asa”snowplow”,pilingupmatterintoahighdensityshell. P = P +P s AGB FW This interaction produces two shocks. The inner shock is Ro v(t−τ) formed at the interaction between central star wind and = 4πr2ρ Vdr+ 4πr2ρ vdr Z AGB Z FW shell.TheouteroneissituatedintheAGBwindjustout- Vt Ri side the shell. Most of the volume interior of the shell is Ro M˙ v(t−τ) m˙ = 4πr2 Vdr+ 4πr2 vdr made up of shocked central-star wind at a high tempera- Z 4πr2V Z 4πr2v ture (millions of degrees). The outer shock is likely to be Vt Ri = M˙ R −Vt +m˙ v(t−τ)−R , (3) isothermal – similar to the shocks in the radiative phase o i (cid:0) (cid:1) (cid:0) (cid:1) of SNR evolution, where a dense shell is formed. while the mass in the shell is: The dependence between gas density and radius of a PN is necessary for deriving the theoretical Σ−D rela- tionforPNe.This dependencewasshowninPaperIfora M = M +M model consisting of one steady wind emanating from the s AGB FW Ro v(t−τ) central star. A trend of decreasing radio surface bright- = 4πr2ρ dr+ 4πr2ρ dr Z AGB Z FW ness with increasing diameter was shown in Paper I. The Vt Ri Ro M˙ v(t−τ) m˙ 2 The radio samples of Galactic PNe are richer in terms of = 4πr2 dr+ 4πr2 dr Z 4πr2V Z 4πr2v thenumberofobjects,butincompletenessisdrasticallyhigher Vt Ri than in SNR samples – less than 5% of total number of PNe M˙ m˙ can bedetected byradio surveys(Kwok 1994). = Ro−Vt + v(t−τ)−Ri . (4) V v (cid:0) (cid:1) (cid:0) (cid:1) D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae 3 If, additionally, the shell is thin: R → R = R then we 2.2. The energy conserving phase i o s may write: Themomentum-conservingphaseisperhaps,areasonable P = M˙ −m˙ R −M˙ Vt+m˙v(t−τ), (5) descriptionofearlyPNevolution.However,ifapartofthe s s energyofthefastwindistransformedintothermalenergy, (cid:0) (cid:1) M˙ m˙ then the pressure of the hot bubble (P) will provide ad- Ms = − Rs−M˙ t+m˙(t−τ) (6) ditional accelerationto the nebular shell. In this “energy- V v (cid:0) (cid:1) conserving”phase,the equationofmotionandthe energy (see Equation 6 in Zhang & Kwok 1993). equation is expressed as: If V = const., the shell is expanding following the law: s d2R M˙ M s =4πR2P − (V −V)2, (15) vVτ vτ s dt2 s V s R = +V t− , (7) s v−V s v−V (cid:0) (cid:1) d 1 4π 1 P · R3 = m˙v2−4πR2PV (16) which may be reduced to Equation (5) in dt(cid:16)γ−1 3 s(cid:17) 2 s s Zhang & Kwok (1993) if v ≫V: (Kwok 2000). R ≈Vτ +V t−τ . (8) Thesimilaritysolutionmethodcanbeusedfordefining s s (cid:0) (cid:1) the change of PN radius with time. The power-law form If t≫τ then, R ∝ tα, is expected for a PN in the energy conserving s phase of evolution. We can define dimensionless variable Rs ≈Vst, (9) in form: and the shell velocity can be expressed as: ξ ≡R tlρm E˙n, (17) s AGB V = Ps = M˙ −m˙ Rs+(m˙v−M˙ V)RVss. (10) where E˙ is the power of fast wind which injects energy s Ms (cid:0)(cid:0)MV˙ − m(cid:1)v˙ (cid:1)Rs+(m˙ −M˙ )RVss (iltahrriotyugahnatlhyesiisn,ntehresfhoollcokw)iningtoextphreehssoitonbuisbboblet.aAinfetder: sim- Solving this equation we obtain: E˙ 1 3 Rs =ξ 5t5. (18) M˙ −m˙ +(v−V) M˙m˙ (cid:0)ρAGB(cid:1) V = q vV , (11) s M˙ − m˙ We obtain a Rs ∝ t (α = 1) dependence assuming V v ρ ∝ R−2 and using Equation (18). Finally, a station- AGB s which is Equation (4) in Zhang & Kwok (1993). ary solution V = R˙ = const. can be used for further s s The average density inside the shell can be expressed derivations. as: The mass of the compressed AGB wind in the shell, using V =const., is: M s s ρ = , (12) s 43πfRs3 Ms =M˙ 1 − 1 Rs. (19) (cid:16)V Vs(cid:17) where f is the volume filling factor. From Equations (6) Ifforf,weassume∆=const.,thenEquation(12)implies and (9) we see that if t≫τ: ρ ∝R−1. s s M˙ m˙ m˙ −M˙ M = − + R . (13) s s h(cid:0)V v (cid:1) Vs i 2.3. The final phase If M˙ ≫m˙, and v ≫V, then M ≈(M˙ + m˙−M˙ )R and We only have AGB wind momentum and mass conserva- s V Vs s tion in the final phase of evolution when the pressure of M˙ the hot bubble, P →0: ρ ∝ . (14) s 4πRs2V P = Rs4πr2ρ Vdr+P′ s Z AGB Essentially, we have perturbed the AGB wind R′ (ρAGB = 4πRM2A˙GBV). From Equation (11) we es- = M˙ (cid:0)Rs−R′(cid:1)+P′, (20) timate V ≈ 15 km s−1 using typical values of s M˙ = 10−5 M⊙yr−1, m˙ = 10−8 M⊙yr−1, V = 10 km s−1and v = 2000 km s−1. We obtain the depen- M = Rs4πr2ρ dr+M′ dance ρs ∝Rs−1 using Equations (12) and (13) assuming s ZR′ AGB that the shell thickness ∆ = Ro − Ri = const. and M˙ f =1−(RRoi)3 =1−(1− R∆s)3 ≈3∆/Rs. = V (cid:0)Rs−R′)+M′, (21) 4 D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae where P′ and M′ are the momentum and mass acquired relationcouldbesteeperduetothetheseeffects,ascanbe by the time the shell radius reaches some value R′. For seenfromEquation(26).Therelationgivenhereisderived the shell velocity we now have: under the assumption of spherical symmetric expansion. However, we know that the large number of PNe show a V = Ps = M˙ Rs−R′ +P′ V → V, (22) high asymmetry and can have a very complex morphol- s Ms M˙ R(cid:0)s−R′ (cid:1)+M′V ogy. This may explain why the theoreticalΣ−D relation (cid:0) (cid:1) is not necessary connected to the empirical relations that as Rs →∞. we discuss inSection 4.Evenif wehad asphericallysym- Duringtheearlyphases,thePNshelltendstothecon- metricPNewithΣ∝D−1,theconstantofproportionality tact discontinuity, while in the last phase, it will be just wouldbe differentfordifferentphasesofevolution,ascan behind the isothermal forwardshock.If we apply isother- be inferred from Section 2. Moreover,parameters govern- mal shock jump conditions, we can estimate: ing the evolution (V, v, M˙ , m˙) are not the same for all PNe, so no unique empirical relation for PNe is feasible. ρ =M2ρ , (23) s AGB In the final stage of evolution, after the fast wind has ceased, the isothermal shock wave continuously perturbs whereM=V /c is the Machnumberandc ,the isother- s s s AGBwind.Nonetheless,thesurfacebrightnesssteeplyde- mal sound speed. The density jump is not limited to a maximum value of 4, as in an adiabatic shock. If V ≈ clines,asdemonstratedinPaperI,becausethemechanical s const., then ρ∝R−2 behind the shock. energyinput(acceleration)fromthecentralstarnolonger s exists. Using ρ ∝R−2, the Σ−D relation in the final AGB s stage of PN evolution can be expressed as: 3. Theoretical Σ−D relation for PNe Σ∝D−3. (27) Assumingathermalbremsstrahlungmechanismisrespon- sible for radiationofHii regionsatradiowavelengths,we 4. The empirical Σ−D relation for PNe may write the volume emissivity ε of a PN as: ν Themostimportantprerequisitefordefiningaproperem- ε [ergs s−1 cm−3 Hz−1]∝n2T−1/2, (24) ν pirical Σ−D relation is the extraction of a valid sample (Rohlfs & Wilson 1996)wherenisthevolumedensityand of PNe. This sample should consist of calibrators having T the electron temperature. well determined distances. The surface brightness can be expressed as: Thedistancestothecalibratorshavetobedetermined by reliable methods. These may include trigonometric or Σ ∝ Lν ∝ε D, (25) spectroscopic parallaxes of PN central stars and meth- ν D2 ν odsbasedonnebulaeexpansion.Althoughlimitedbycur- rent technology, trigonometric parallax is the most direct where L is luminosity and D =2R , the diameter of the ν s method of measuring distances. PN.WecombineEquations(24)and(25),usethe derived dependence from Section 2 such that n ∝ ρ ∝ D−1, and All extracted PN samples suffer from severe selection s assume that temperature is constant.3 This radio surface effects causedby limitations insurvey sensitivity andres- olution. For Galactic samples, Malmquist bias4 may be brightnesstodiameterrelationhasadifferentformforthe the most severe effect; this is similar to the situation for first and second phase of evolution (in comparison to the Galactic SNR samples (Uroˇsevi´c et al. 2005). earlier result given in Paper I): For the literature extractedPN samples in this paper, Σ∝n2T−1/2D ∝D−1. (26) weusedfluxdensitiesat5GHz.ThisisbecausemostPNe are expected to be optically thin at this frequency (see Generally, we may write n ∝ D−x, where x >∼ 1 (if ∆ = Section 5). Radio diameters would be preferable,however const.). for most PNe used here, only optical diameters are avail- We do not expect the temperature to be strictly con- able. stant throughout the nebula. There are temperature gra- We extracted 27 PN sample sets as shown in Tables 1 dients in PNe arising from radiation hardening. More en- and 2. While both detail fit parameters for Σ− D and ergetic photons will travel further and when they are ab- L−D,Table1doessobyauthor(18PNsets)andTable2 sorbed by the PN, they will impart a greater kinetic en- bythespecificdistancedeterminationmethodused(9PN ergy to the ions thereby producing a higher temperature. sets). This would only slightly increase the value for β (+0.1 Data for sample Nos. 1 – 4 in Table 1 are taken from approximately), although the Te −D dependence is not SSV and Cahn et al. (1992). We used Equation (7) from quite of the power-law form (see numerical model results Van de Steene & Zijlstra (1995)to convertthe brightness from Evans & Dopita 1985). temperatures listed in Ph02 to their 5 GHz flux. PN dis- Thevalueβ =1undertheaboveconditions(including tance ranges (kpc) are listed in Table 1 (e.g. “1.0< d≤ ∆ = const.) is only a theoretical lower limit; the Σ−D 4 IntrinsicallybrightPNearefavoredinanyfluxlimitedsur- 3 Hiiregions are approximately isothermal at T ∼104 K. vey because they are sampled from a larger spatial volume. D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae 5 2.0” indicates that all PNe in that set have distances be- correlation. The abundance of low rL−D values as shown tween 1 and 2 kpc). Three PNe (He 1-5, HDW 6 and in Tables 1 and 2, indicate that the dependence between HaWe 5) indicated by “†” were removed from some sets Σ and D is uncertain. Therefore, the validity of Σ−D because of inconsistencies with other data, as we will de- relations is diminished even if free from Malmquist bias. tail in Section 5. Aside from a few potential sample sets, including No. 1 Sampleset15istakenfromZhang (1995)andconsists listed in Table 1, most should not be used when calibrat- of 132 PNe with well-determined individual distances. ing the PN Σ−D relation. Set 16 contains 23 PNe with well-determined distances, Sample No. 1 in Table 1 is the only set with a re- while, set 17 consists of 132 Galactic bulge PNe. Both spectable correlation coefficient, rL−D = −0.67, yet hav- samplesweretakenfromVan de Steene & Zijlstra (1995). ingalowprobabilityofMalmquistbias6.ThisSSVsample Set 14 consists of 109 Galactic bulge PNe taken from consistsof13relativelynearbyGalacticPNe,withreliable Bensby & Lundstr¨om (2001). Set 18 (88 PNe) is actually individually determined distances less than 1 kpc. They the Van de Steene & Zijlstra (1995) Galactic bulge sam- are different in morphological form (e.g. round, elliptical, plewithstrongercriteriaappliedtoextractGalacticbulge bipolarcore,bipolar,etc.).Forthis sample,weobtainthe PNe as defined in Bensby & Lundstr¨om (2001). Σ−D relation: The first eight sample sets listed in Table 2 were ob- Σ=7.83+4.1710−22D−2.61±0.21, (30) tained from SSV and Ph02 data. Set 9, USNO-PN, con- −2.73 sists of PNe with well-measuredtrigonometricparallaxes, and a corresponding Σ −D and L −D diagrams are ν ν determined by United States Naval Observatory (USNO) shown in Fig. 2. ground based observations (Harris et al. 2007). Details of We define the fractional error as: this USNO-PN sample data are presented in Table 3. d −d SSV Σ The trigonometric parallaxes of PNe determined by fe= , (31) (cid:12) d (cid:12) the Hipparcos mission suffer from large measurement er- (cid:12) SSV (cid:12) (cid:12) (cid:12) rors(Harris et al. 2007)andthereforeareoflimitedprac- in or(cid:12)der to obta(cid:12)in an additional estimate of accuracy. tical use for our study. The averagemeasurementerrorin IndependentlyderiveddistancesdSSV aretakenfromSSV Hipparcos trigonometric parallaxes for some of our PNe and dΣ is the distance derived (for each individual PN) is 57 %. Thus, the USNO-PN sample set is more reliable fromthe Σ−D relation.The maximumandaveragefrac- andhavemoreaccuratedistancemeasurents5.TheΣ−D tional errorsare femax =1.08and f¯e=0.35,respectively. relationfor the USNO-PNsample has the following form: Accordingtoourcriteriaandtherelativelysmallaver- age fractional error,Equation (30) may be used as a cali- Σ=7.39+−42..570910−22D−2.38±0.56. (28) brationrelationfordistancedetermination.SincethePNe in the corresponding sample are nearby, the Malmquist Corresponding Σ −D and L −D diagrams are shown ν ν bias should not be significant. The L−D relation exists in Fig. 1. and the slope β is reasonable. However, this slope is still MalmquistbiastendstoincreasetheslopeoftheΣ−D most likely a ”mixture” of the evolutionary and selection relation for SNRs (Uroˇsevi´c et al. 2005). Similar implica- effects. Additionally, the sample number is small (only 13 tions shouldbe validfor PNe aswell.Slopes fromTable 1 PNe) and they do not belong to the same morphological suggest that some of our sample sets may suffer from a type. Therefore, all results obtained using Equation (30) Malmquist bias. Specifically, Σ−D slopes derived from should be taken with caution as we discuss below. the SSV and SSV+Ph02 are steeper for higher sampling distances (i.e. >2 kpc). 5. Discussion In spite of the possibility of Malmquist bias, the ma- jority of the Σ−D slopes listed Tables 1 and 2, are very In Section 2, we derived the theoretical dependence be- close to the so-called trivial Σ−D relation with β = 2 tween the gas density in the PN shell and its radius with (Arbutina et al. 2004). This can be understood mathe- respect to the interacting stellar winds model. Next, we matically from: showed that the slope of the theoretical Σ−D relation is S L(D)/d2 changed significantly by including these equations. Σ∝ ∝ ∝D−2L(D), (29) The theoretically derived slope (β ≈ 1) is shallower Ω D2/d2 than slopes from the empirical relations given in Tables whereS isthefluxdensity,Ωisthesolidangleoccupiedby 1 and 2. This discrepancy may possibly be explained by the source and d is the source distance. The usefulness of poor quality of Galactic PN samples or by incorrect as- theΣ−D relationcanbetestedwiththecorrelationcoef- sumptions used for the derivation of the theoretical rela- ficientrL−D ofthe L−D relation,whichshouldapproach tion. Due to variations in the power law density distri- –1.0 if the data can be explained by a linear relation. bution derived here and the approximately constant tem- It is clear from Equation (29) that the Σ− D rela- perature of the expanding PN envelope, the theoretical tion has a trivial Σ ∝ D−2 form in the absence of L−D 6 The correlation coefficients for sets 3 and 9 from Table 1 5 An average error of measured trigonometric parallaxes in are the steepest, but Malmquist bias is also likely. All other this sample is 17.2 %. Σ−D slopes are approximately trivial. 6 D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae Table 1. The results of the Σ−D and L−D fits, where the parameters of the fit (β and α, respectively) and the correlationcoefficientraregivenforeachsample.ThenumberofPNeinthesampleisgiveninthelastcolumn(N).SSV stands for Stanghellini, Shaw & Villaver (2008), B&L for Bensby & Lundstr¨om (2001) and Ph02 for Phillips (2002). No. Samplea β r α r N Σ−D Σ−D L−D L−D 01 SSVd<1.0 2.61±0.21 −0.97 0.61±0.21 −0.67 13 02 SSV1.0≤d<2.0 1.89±0.43 −0.78 –0.11±0.43 0.08 15 03 SSVd≥2.0 3.96±0.62 −0.91 1.97±0.61 −0.73 11 04 SSV 2.56±0.22 −0.89 0.56±0.22 −0.40 39 05 SSV+ Ph02 d≤0.4 1.47±0.40 −0.74 –0.53±0.40 0.38 13 06 SSV+ Ph02 0.4<d≤0.6 2.18±0.22 −0.91 0.18±0.22 −0.17 23 07 SSV+ Ph02 0.6<d≤1.0 2.37±0.20 −0.95 0.37±0.20 −0.43 18 08 SSV+ Ph02 1.0<d≤2.0 2.64±0.40 −0.88 0.64±0.40 −0.41 15 09 SSV+ Ph02 d>2.0 4.44±0.47 −0.96 2.44±0.47 −0.90 9 10 SSV+ Ph02 2.40±0.17 −0.86 0.40±0.17 −0.27 78 11 SSV+ Ph02 d≤0.4 † 2.31±0.38 −0.90 0.31±0.38 −0.26 11 12 SSV+ Ph02 0.4<d≤0.6 † 2.33±0.15 −0.96 0.33±0.15 −0.44 22 13 SSV+ Ph02 † 2.57±0.14 −0.92 0.57±0.14 −0.45 75 14 B&L 2.31±0.11 −0.90 0.31±0.11 −0.26 109 15 Zhang(1995) 2.17±0.14 −0.80 0.17±0.14 −0.10 132 16 Vande Steene& Zijlstra (1995) 2.41±0.34 −0.84 0.41±0.34 −0.26 23 17 Vande Steene& Zijlstra (1995) 2.24±0.09 −0.92 0.24±0.09 −0.24 132 18 Vande Steene& Zijlstra (1995) (stronger criteria byB&L) 2.21±0.09 −0.94 0.21±0.09 −0.25 88 a PNefrom Ph02 are included in the sample sets (5 – 13) as additional, if they were not extracted by theSSV in their sample; d<1 denotes that a set of PNehavedistances less than 1 kpc,1≤d<2, between 1 and 2 kpc,etc. † – 3 PNe (He1-5, HDW 6 and HaWe 5) are removed (see text). Table 2. The results of the Σ−D and L−D fits of PN samples defined using the different methods of distance determination. The parameters of the fit (β and α, respectively) and the correlation coefficient r are given for each sample. The number of PNe in the sample is given in the last column (N). No. Samplea β r α r N Σ−D Σ−D L−D L−D 1 Trigonometric Parallax 1.79±0.22 −0.90 −0.21 ±0.22 0.23 18 2 Trigonometric Parallax † 2.02±0.18 −0.95 0.02 ±0.18 −0.02 17 3 Geometrical (trig. parall. + expan.) 2.26±0.26 −0.86 0.26±0.26 −0.19 30 4 Geometrical (trig. parall. + expan.) † 2.44±0.22 −0.91 0.44±0.22 −0.36 29 5 Gravitational 1.42±0.37 −0.66 −0.58 ±0.37 0.33 22 6 Gravitational † 2.53±0.31 −0.89 0.53 ±0.31 −0.38 20 7 Spectroscopic 2.57±0.36 −0.91 0.57 ±0.36 −0.43 13 8 Extinction 1.69±0.38 −0.80 −0.31 ±0.38 0.25 13 9 USNO-PN 2.38±0.56 −0.90 0.38 ±0.56 −0.32 6 a For thefirst eight samples, PNeare extracted from the SSVand Ph02; for USNO-PN– see Table 3. Table 3. The data for the USNO PN sample. Name Trigonometric Parallax S5GHz Diameter [mas] [mJy] [′′] NGC7293 (036.1-57.1) 4.56±0.49 1292a 660b NGC6853 (060.8-03.6) 3.81±0.47 1325a 340a NGC6720 (063.1+13.9) 1.42±0.55 384a 60c A 21 (205.1+14.2) 1.85±0.51 157d 550e A 7 (215.5-30.8) 1.48±0.42 305a 760a A 24 (217.1+14.7) 1.92±0.34 36a 415d a b c d e References: Cahn et al. (1992); O’Dell (1998); George et al. (1974); Hua& Kwok (1999); Salter et al. (1984). D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae 7 -19 21 10 10 -20 20 10 10 1] -1-Hz sr 1-1 Hz] -2m -g s [W Hz [erGHz G L5 5 Σ -21 19 10 10 -22 18 10 10 0.1 1 0.1 1 D [pc] D [pc] Fig.1. (Left)Σ −D diagramat5GHzforsampleset9inTable2,having6USNO-PNewithaccuratetrigonometric ν distances determinations. (Right) Corresponding 5 GHz L −D diagram for the same set of PNe. ν slope may be different than given in Equations (26) and central star reaches a temperature of 40 000 K, the num- (27) (derived assuming the shell thickness ∆ ≈ const.). ber of ionizing photons becomes almost independent of If the filling factor f is constant, the slope β could be as temperature and radio luminosity no longer increases. As greatas3.Weconcludethatwhetherthesesimplyderived discussed by Van de Steene & Zijlstra (1995), this causes theoreticalrelations havethe correctformsor not,empir- a decrease in surface brightness with a gradualexpansion ical relations are under the influence of biases that make ofthenebula,whichthenbecomesopticallythinandthus their slopes different than they would be if one considers detectable. To keep our sample pure, we devised a set of the evolution of a single hypothetical PN. An empirical criteriawhichexcludedthreesuspectedthickPNeintheir relation can be useful for estimating distances only under rising phase ofevolution(as markedwith “†”in Table 1). strict conditions as discussed below. From Tables 1 and 2, we suggest that 24 slopes Van de Steene & Zijlstra (1995) and Phillips (2004) are most likely trivial Σ − D relations. Therefore, we showed that the evolution of PNe may not have a linear feel that the Malmquist bias is not so severe as in trend in log-log scales. When the radiation of the central the case of Galactic SNR samples, where the slopes star begins to ionize the nebula, the radio brightness will are significantly steeper. For example, the obtained increase rapidly. However, in this phase of PN evolution, slope for the sample of 132 bulge PNe collected by the observedincidence of radio nebulae is predicted to be Van de Steene & Zijlstra (1995)takesarelativelyshallow low.Therefore,our linear trend ofthe Σ−D formcanbe valueofβ ≈2.Thisisexpectedresultsinceall132PNeare justifiedandthequadraticformoftheR−T dependence located at the approximately same distance. Malmquist b givenbyPhillips (2004)isnotnecessaryforthestatistical bias does not exist in samples which contain objects at analysisgivenhere.Inaddition,differentdependenciescan the same distance. notbederivedfromthethermalbremsstrahlungradiation Both Schneider & Buckley (1996)and at a later time, formula (Equation 24). Bensby & Lundstr¨om (2001) attempted to improve the The subjects of our analysis are optically thin PNe. Bulge sample by removing foreground and background ThereasonforthelowincidenceofobservedradioPNdur- PNe. Bensby & Lundstr¨om (2001) defined stronger cri- ing the risingphase,discussedabove,is thatthe majority teria to extract a more pure Bulge sample (Set 14 in ofthesePNeareopticallythickat5GHz.Then,whenthe Table 1). We found a trivial Σ−D slope for this sam- 8 D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae -16 24 10 10 -17 23 10 10 -18 22 10 10 -19 21 10 10 1] -1-Hz sr 1-1 Hz] -2m 10-20 -g s 1020 [W Hz [erGHz G L5 5 Σ -21 19 10 10 -22 18 10 10 -23 17 10 10 -24 16 10 10 0.01 0.1 1 10 0.01 0.1 1 10 D [pc] D [pc] Fig.2. (Left) Σ −D diagram at 5 GHz for the set of 13 SSV PNe with distances < 1 kpc (from sample set 1 in ν Table 1). (Right) Corresponding L −D diagram at 5 GHz for the same set of PNe. ν ple. We also applied this stronger criteria directly to the By examination of sample sets 5 – 9 in Table 1 (and Van de Steene & Zijlstra (1995) sample. The change in Fig.3)wetrytodemonstratetheeffectofMalmquistbias. slope for this modified Van de Steene & Zijlstra sample Forset5withdistanceslessthan0.4kpc,thecorrespond- (Set 18 in Table 1) is negligible. This again supports the ingΣ−Dslopeisβ =1.47.Furtherdistancecomparisons idea that these sample sets are not strongly affected by (in the form of distance interval in kpc followed by the the Malmquist bias. Σ−D slope) include: 0.4 < d ≤ 0.6, 2.18; 0.6 < d ≤ 1.0, 2.37; 1.0 < d ≤ 2.0, 2.64; and > 2.0, 4.44. This implies that an increase in distance gives a Σ−D slope that is According to criteria defined in this paper, the corre- also greater. spondingslopesfortheBulgesamplesrepresenttheresult ofscatteringintheL−Dplane.Thisslope(β ≈2)wasob- Using both the Tukey-Kramer (T-K) and the Scheffe tainedforthe extragalacticsamplesofSNRs (exceptM82 (S) methods, we found the Σ− D slope for sample set sample) where Malmquist bias also does not have severe 9 in Table 1 (> 2.0, 4.44) is statistically different at a effect. Again, this is because all of the SNRs in each set confidence level of α = 0.05, if compared to the slopes areatapproximatelythesamedistance(seeUroˇsevi´c 2002 of sets 5, 6, and 7. However, since all of these samples and Uroˇsevi´c et al. 2005). consist of PNe with some uncertainty in their distances, D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae 9 any firmer conclusion about their Malmquist bias should 10-16 be avoided. Using the same methods, sets 1 (< 1, 2.61) SSV (d < 1) SSV (1 < d < 2) and 3 (≥ 2, 3.96) in Table 1 also represent a statistically SSV (d > 2) differentpopulationataconfidencelevelofα=0.05.The d < 1 kpc d > 2 kpc difference between their corresponding Σ − D slopes is statistically significant at a confidence level of α = 0.2 10-17 and can be explained by the Malmquist bias. Sample set 2 in Table 1 is omitted from our slope analysisbecause of its extremely high level of scattering in the L−D plane (see Table 1, rL−D =0.08). 10-18 InspectionofTable2revealsthatwhileallΣ−Dslopes arerelativelycloseto the trivialform,set4,6 and7area bit steeper. The increased slope of set 4 can be explained by perturbations induced by twelve very luminous and 10-19 highly Σ−D plane dispersed “expansion” PNe added to sample set 2. Since the gravitational and spectroscopic 1] -sr methodsaremostlyappliedfordistantPNe,slopesinsets 1 -z 6and7maybesteeperduetoaMalmquistbias.However, H alloftheL−DcorrelationcoefficientsinTable2arevery -2m 10-20 low and therefore poor for statistical investigations. W ThehighlevelofscatteringintheL−D plane,asseen [Hz G inFigs.1(right)and4,supportstheideathatmostofthe Σ5 slopesinTables1and2donothaveanyrealphysicalinter- 10-21 pretation.Itisthisluminosity-diameterscatteringartifact which produces the trivial Σ ∝D−2 form. Therefore, 24 of these trivial Σ−D slopes are not useful for the deter- mination of distances. Causes for this scattering include: imprecisely determined calibrator distances, mixtures of 10-22 different PNtypes in the same sample,limitations in sen- sitivity and resolution of radio surveys, source confusion and to a lesser extent, Malmquist bias. Although, our criteriaimply that the relationgivenin 10-23 Equation(30) could be useful for distance determination, we assigna low confidence that it represents a real evolu- tionary track. The Σ−D slope from this relation repre- sents the result of coupling of evolutionary and selection 10-24 effects. The corresponding sample set consist of only 13 0.01 0.1 1 10 PNe and they do not belong to the same morphological D [pc] type. Thus, this sample is not complete or homogeneous Fig.3. The Σ − D diagram at 5 GHz for all 39 SSV and distance estimates obtained by using Equation (30) ν PNe (sample set 4) with reliable individually determined should be viewed with appropriate caution. distances. The straight lines represent the least square fit PN radio surface brightness depends on a number of lines through the two distance limiting Galactic samples. factors that include: intrinsic nebula evolution, changes The difference in slope of corresponding Σ−D relations in density and temperature, changes in medium ioniza- are defined for PNe with distances < 1 kpc and ≥ 2 kpc tion andextremely largeluminosity variations in the cen- (represented with dash and tick line, respectively). These tral star during its evolution. In Section 2, we assumed couldbeexplainedasaresultofthevolumeselectioneffect that shell expansion was constant throughout each phase called Malmquist bias. of nebula evolution. On the contrary, PN shells are ac- celerated by the central star fast winds at the beginning of PN evolution, and only later are they driven at ap- proximatelyconstantspeedbythehotbubblescreatedby the crossing of the fast winds through the inner shocks. than 0.03 pc (see Fig. 4) and therefore in later phases of Shaw et al. (2001)discusses accelerationof PN shells, i.e. evolution where acceleration (or central star influence) is how the velocity of PN shells is related to the evolution negligible. This explanation is admittedly simplistic. For of the central star. They emphasized that the accelera- essential understanding of the theoretical and empirical tion of shells is nearly zero in larger (and probably older) Σ−D dependencies, the full explanation of the complex PNe having diameters > 0.03 pc. Almost all of the PNe co-evolution of a nebula and its central star would be from sample sets presented in Tables 1 and 2 are larger needed. 10 D.Uroˇsevi´c et al.: The Σ−D relation for planetary nebulae The most massive central stars will never fully ionize (NGC6720orRingNebula,seeGeorge et al. 1974 );sym- their environment and those PNe will never become op- metric (A21 or S274, see Salter et al. 1984 ); and S-type tically thin at 5 GHz since their central stars evolve so (A247).However,afurther searchofthe literaturereveals rapidly. Those observed PNe with smaller diameters and thatthe intrinsicmorphologiesforallfiveoftheseobjects surface brightnesses (optically thick PNe at 5 GHz) are is most likely to be bipolar, i.e. that the structures of the located significantly lower in the Σ−D or L−D planes emitting regions(or shells) rangefroma thick disk (Helix when compared to the average (thin) PNe. Their evolu- Nebula) to a triaxial ellipsoid (A24, and Ring Nebula) tion aresurely not defined by the Σ−D relations derived to a barrel-like shape (A21, and Dumbbell Nebula) here.Althoughthe majorityofthe thickPNewillbecome (see O’Dell 1998, Henry et al. 1999, Hua & Kwok 1999, thin,thequestionisWHEN?Thisdependsonanumberof O’Dell et al. 2007 and Meaburn et al. 2000). The sixth factors during their evolution. These effects significantly objectforwhichwefoundreliableradioobservations(A7) contribute to uncertainties in addition to the huge scat- appears to be a “classical” PN with a well defined spher- tering found in the L−D plane. ical shell (Xilouris et al. 1996). Other factors not included in our analysis may have We note that Kastner et al. (1996) and a large impact on the evolution of PN emissivity. These Zuckerman & Gatley (1988) detected a shock-excited include: compression of gas induced by the interaction of H2 emission from the Helix, Ring and Dumbbell fastwindwiththemediumonthemainshellinnerbound- nebulae, implying an ionization bounded case for ary, isothermal shocks associated with the passage of D- these objects. Additionally, no spectroscopic evi- ionizationfronts through neutral gas and progressiveion- dence for the existence of a high-velocity stellar wind izationofgreateramountsofsurroundingneutralgas(e.g. has been found for the previously mentioned ob- Shaw et al. 2006, and references therein). We do not an- jects, including A7 (Cerruti-Sola & Perinotto 1985, alyze ionization fronts in this paper; our derivations in Patriarchi& Perinotto 1991). Central stars of these ob- Section 2 are based only on the dynamics of correspond- jects are most likely well pass the hydrogen-shell burning ing shock waves.The surface brightnessevolutionmaybe phase and approaching the white dwarf cooling sequence. more closely connected to ionization than the dynamical A wide range of shell structures is predicted in the expansionoftheshell.Thismaybegoodenoughreasonto hydrodynamical models based on the ISW model and it challenge the quality of our derived theoretical relations is very likely that a large number of PNe classified by and anadditional explanationfor scattering in the L−D apparent morphology as spherical or elliptical, may in plane. Finally, we concede that some PNe could have dif- fact posses a bipolar structure. Even though morphologi- ferent initialconditions leading to independent evolution- calclassificationsbasedonopticalimaging favorelliptical ary paths. These paths could follow the same theoretical objects, Zuckerman & Aller (1986) predict that approxi- Σ−D curvebutwithvaryingintercepts,thereforeleading mately 50% of all PN in the Galaxy are actually bipolar. to the observed scatter. Nevertheless, we conclude that the USNO-PN sample is not representative of the majority of PNe because of its In order to define a valid PN sample for statistical “narrowness” in morphology (i.e. almost all appear to be investigation, one needs to separate a higher number of bipolar). nearby PNe with the same morphological characteristics. Wedefinesuchasample(USNO-PN)inthesubsectionbe- low(Section5.1).However,thequalityofradio-continuum 6. Summary observations today simply is not sensitive enough to pro- The main results of this paper may be summarized as vide a high number of similar PNe having reliably deter- follows: mined distances. (i) We have derived a new theoretical Σ−D relation for PNe in the form Σ ∝ D−1 by including the in- 5.1. The USNO-PN Sample teraction between the AGB star wind and the fast wind from the central star. This dependence is ob- We have examined the morphological properties of tained for both momentum conserving and energy PNe from the USNO-PN sample using the radio im- conserving phases of evolution. Only for the final ages of these objects found in the literature. We use stage,wedidderivethesametheoreticaldependence comparisons with radio morphological classes given by (Σ∝D−3) as in Paper I. Aaquist & Kwok (1996) which are based on the prolate (ii) We derivedempiricalΣ−D relationsfor 6updated ellipsoidal model (PES). The PES model describes diver- PN samples from the literature, 12 new sample sets gency in observed (empirical) morphologies through dif- extractedfromtheseand9additionalsetsasdefined ferent sky projections of relatively simple spherical shell inthispaper.Wediscusstheselectioneffectsthatin- modelswithbothradialandlatitudedensitygradientsand fluence these PN sample sets. Our results show that withdifferentionizationdepths.USNO-PNsampleobjects can be classified as follows: circular (NGC 7293 or Helix 7 For A24 we acquire a 1.4 GHz radio image from the Nebula, see Zijlstra et al. 1989); open elliptical (NGC NRAO VLA Sky Survey (NVSS) postage stamp server 6853 or Dumbbell Nebula, see Bignell 1983); elliptical (http://www.cv.nrao.edu/nvss/postage.shtml).