A Consistent Spectral Model of WR 136 and its Associated Bubble NGC 6888
We analyse whether a stellar atmosphere model computed with the code CMFGEN provides an optimal description of the stellar observations of WR 136 and simultaneously reproduces the nebular observations of NGC 6888, such as the ionization degree, which is modelled with the pyCloudy code. All the observational material available (far and near UV and optical spectra) were used to constrain such models. We found that even when the stellar luminosity and the mass-loss rate were well constrained, the stellar temperature at = 20, can be in a range between 70 000 and 110 000 K. When using the nebula as an additional restriction we found that the stellar models with 70 000 K represent the best solution for both, the star and the nebula. Results from the photoionization model show that if we consider a chemically homogeneous nebula, the observed N/O ratios found in different nebular zones can be reproduced, therefore it is not necessary to assume a chemical inhomogeneous nebula. Our work shows the importance of calculating coherent models including stellar and nebular constraints. This allowed us to determine, in a consistent way, all the physical parameters of both the star and its associated nebula. The chemical abundances derived are 12 + log(N/H) = 9.95, 12 + log(C/H) = 7.84 and 12 + log(O/H) = 8.76 for the star and 12 + log(N/H) = 8.40, 12 + log(C/H) = 8.86 and 12 + log(O/H) = 8.20. Thus the star and the nebula are largely N- and C- enriched and O-depleted.
keywords:circumstellar matter – stars: Wolf-Rayet stars.
NGC 6888 is an emission nebula associated with the Wolf-Rayet star WR 136 with spectral type WN6(h), and it is a typical example of a ring nebula resulting from the continuous interaction between a strong stellar wind and the interstellar medium. Projected on the sky, this nebula seems to be a filamentary and clumpy gaseous ellipse emitting mainly in H and optical [N II] lines, with an axis-size of arcmin (Chu, Treffers, & Kwitter, 1983), and it is surrounded by a homogeneous sphere which emits mainly in [O III] 5007 Å (Gruendl et al., 2000; Moore, Hester & Scowen, 2000). The combination of its filamentary morphology and its nitrogen enhanced abundance suggests that the ejections, coming from a previous evolutionary stage of the star WR 136, have contributed in an important way to the formation of the ring nebula that we see today (Esteban & Vílchez, 1992, hereafter EV92; Fernández-Martín et al., 2012, hereafter FM12; Mesa-Delgado et al., 2014, hereafter M-D14).
Both the nebula NGC 6888 and its central star WR 136 have been already studied in several wavelength ranges from X-rays to radio emission. For example, the stellar parameters, such as the stellar temperature, mass-loss rate and terminal velocity, have been calculated by Crowther & Smith (1996) and Hamann, Grafener & Liermann (2006, hereafter HGL06) by using stellar atmosphere codes and comparing the results with observed ultraviolet, optical and infrared stellar spectra; while the wind velocity law and clumpiness were studied by Ignace, Quigley & Cassinelli (2003). Skinner et al. (2010) showed that WR 136 emits in X-ray energies. In particular, these authors found soft emission ( keV) coming from outside the wind and hard emission ( keV) coming from the gas inside the wind.
The complete hot bubble in the interior of NGC 6888 was observed with the Space Telescope Suzaku (Zhekov & Park, 2011), confirming the existence of nitrogen enriched X-ray emitting plasma which consists of two thermal components, one emiting soft X-rays and the other one, hard X-rays. Based on this, the authors proposed a new X-ray emission mechanism based on the thermal interaction between heavy particles, given that the known mechanisms are not capable to explain simultaneously these two X-ray emitting components. However, as these authors explain, it is necessary to have better observational restrictions to confirm this idea.
On the other hand, complete hydrodynamical simulations performed by Toalá & Arthur (2011), attempting to reproduce the X-ray observations, indicate that the nebular morphology plays a key role in the presence of X-rays in ring nebulae. These authors conclude that the X-ray emission is higher in bubbles with shell-like morphology than the one present in clumpy filamentary bubbles. However, as the authors indicate, the stellar evolutionary models are not able to reproduce simultaneously the observed stellar parameters and the spectral signatures of the observed X-ray luminosity.
The physical properties in different zones of NGC 6888, such as electronic temperature and density, expansion velocity and chemical composition, have been investigated by several authors since Kwitter (1981) to recently Stock & Barlow (2014). EV92 made a chemodynamical study and they found that NGC 6888 has a chemical composition revealing the material processed in the interior of a star with an initial mass between 25 M and 40 M. These results were confirmed by M-D14 who detected for the first time the C II 4267 Å recombination line, allowing to estimate the C abundance. In this way they presented a complete study of the trace of CNO cycle in the nebula.
FM12 made a very complete study of NGC 6888 through integral field spectroscopy observations performing 1D and 2D analyses. Their results show indeed that NGC 6888 has a complex structure with shells and filaments moving with different velocities. In particular they proposed an “onion-like” morphology composed mainly by three shells. Concerning the plasma properties, they found that the electron density ranges from less than 100 to 360 cm and the electron temperature presents variations from 7 700 to 10 200 K. In addition, they reported important variations in the nitrogen abundance across the nebula, the N/O ratio being much higher in the inner region. Recently Stock & Barlow (2014) derived N/O values from the infrared [N III] 57-m and [O III]-88 m lines and their results would indicate that the edge region of NGC 6888, over the major axis, has lower N abundance with respect to the inner region near the star.
All the previous works focus only on the modelling of the star or the nebula but not both simultaneously. Further, when a single study of these systems is made (star or nebula) it is difficult to assign a value to the physical parameters only based on generic data. This is why we focus in a semiempirical method based in the self-consistent modelling of the star and its associated nebula.
It is logical to proceed in this way because the nebula associated to the star is naturally an additional restriction to the stellar atmosphere model, because the nebula can tell us about the ionizing ultraviolet photons coming from the star and then reprocessed and reemitted with low energy through recombination and collisional processes. In this way we can realize how well a stellar model can produce the ionization stage of the nebula and for this purpose we build a photoionization model of NGC 6888 based on spectroscopic observations where the modelled stellar spectral energy distribution (SED) is used as the ionization source.
This methodology will give us the advantage to obtain more constrained and coherent physical values for all the components of the entire system than those resulting from a single study of each component as it was made until now. The results from this work will help to better constraint other kind of studies as models of stellar evolution, astrophysical bubble formation and, in particular, the presence of X-rays in the bubbles since this hot gaseous component imposes an additional and important restriction to understand the whole system, given that is there where the stellar wind energy is stored.
This paper is structured as follows: In Section 2, observations and data reduction process are described. Sect. 3 contains the information of the line flux measurements and the error estimates. In Sect. 4 we describe the stellar and nebular models as well as our methodology to constraint the free parameters. A discussion about our best self-consistent model is presented in Sect. 5. Finally, in Sect. 6, we present our conclusions.
|Observation ID||A. R||Dec.||Date||Exp. Time||Slit dimensions|
|NGC 6888A||20:12:36.5||38:26:56.5||1st sept. 2011||41800||413.3|
|NGC 6888B||20:12:03.7||38:25:07.8||30 Aug. 2011||41800||413.3|
|NGC 6888C||20:11:51.2||38:21:45.8||30 Aug. 2011||41800||413.3|
|WR 136||20:12:06.5||38:21:17.7||28 Aug. 1997||3300||9.313.3|
|HR 7596||19:54:44.7||00:06:25.1||1st sept. 2011||3420||9.313.3|
2 Observations and data reduction
The main analysis in this work is based on the optical observations carried out the night 1997 August 28, for the star WR 136; and the 2011 August 31 and September 1, for the nebula and the standard star HR 7596. We used the REOSC echelle spectrograph (Levine & Chakrabarty, 1994) attached to the 2.1-m telescope at the Observatorio Astronómico Nacional (OAN) in San Pedro Mártir, B. C. México. With this spectrograph we obtained high-resolution spectra () with a wavelength coverage from 3700 Å to 7300 Å.
Three different positions of the nebula were observed using a slit with size 4 13.3 arcsec oriented east-west. For each position we obtained four spectra with an exposure time of 1800 s each. The slit positions are shown in Figure 1. Together with slit positions of other authors that will be described in the following sections.
To perform the spectroscopy of WR 136 and the standard star the slit size was 9.3 13.3 arcsec. For each star three spectra were obtained, with a exposure time of 420 s in the case of the WR star and 300 s for the standard star. The log of observations is shown in Table 1.
After the bias correction was made, the corresponding different exposures were combined and then background and cosmic rays were subtracted using the filter/cosmic task of the astronomical software MIDAS111http://www.eso.org/sci/data-processing/software/esomidas. Then, once the wavelength calibration was performed (using a ThAr lamp), the orders of the echelle images were extracted.
To correct the extracted orders by the blaze function and to make the flux calibration in the case of the nebular spectra, we modelled the star HR 7596. First the individual orders were normalized to a curve which is theoretically the convolution between the blaze function and the sensitivity function. In this way, the normalized orders were merged to obtain a complete normalized stellar spectrum to be modelled with the code SPECTRUM (Gray, 1999). This code assumes a LTE atmosphere with a plane-parallel geometry and is adapted for modelling standard stars with a spectral type between B and M. The spectrum can be calculated in a normalized way or the model can give the emergent flux per wavelength unit. First we obtained the normalized modelled spectrum to directly compare it with the observed spectrum and once this one was reproduced, the not normalized spectrum was calculated.
The model spectrum of the standard star was then reddened considering the interstellar and atmospheric extinction as well as the radiation dilution. In this way, we constructed the synthetic standard star spectrum. Then, the extracted orders of the standard star observation can be divided by this spectrum to obtain the single order function where the information about the instrumental sensitivity and flux calibration is contained. To correct for the CCD sensitivity and to flux-calibrate the nebular spectra, these orders were divided by their corresponding single order function. Finally the corrected extracted orders were analysed as individual spectra. These resulted spectra have a total exposure time of two hours.
In the case of the orders from the WR 136 frames, we proceeded with the rest of the reduction similarly as we made with the standard star and given that observations of the standard star and WR 136 were not made at the same date, the resulting optical spectrum was just normalized to the continuum.
In addition to the optical spectra of WR 136, we analyzed a spectrum in the UV region between 900 Å and 1200 Å obtained from the FUSE satellite observational program C097 where some interesting and useful lines, such as S V 1063 and P V 1118 lie. To extend the wavelength range to 3300 Å we retrieved several reduced IUE spectra222see https://archive.stsci.edu/iue/file_formats.html for more information. To improve the signal to noise ratio we combined the individual IUE short-wavelength high-dispersion spectra SWP57618, SWP57622, SWP57625, SWP57626, SWP57630, SWP57635, SWP57638, SWP57645, SWP57647 and SWP57663; and the individual IUE long wavelength high dispersion spectra LWP32398, LWP32407, LWP32409, LWP32413, LWP32416, LWP32418, LWP32423, LWP32440, LWP32441 and LWP32446.
3 Flux measurement and error estimates
In Figure 2 we show the resulting H and [N II] 6548, 6584 Å lines present in each spectrum of NGC 6888. In this figure it is apparent that the spectra show line emissions coming from three components at different velocities. This is more obvious for the zone B in whose spectrum we can clearly distinguish these components located at around 45, 15 and 65 km s. According to the wavelength calibration performed, the radial velocities have an accuracy of 12 km s.
In the spectrum of position C only two components are present with radial velocity of 57 and 12 km s. Probably the redshifted one is also present but very weak. In the case of position A, corresponding to the nebular edge, it is only evident the more intense blueshifted component at around 36 km s. We can see that the line profiles are not gaussian at all, they exhibit a wing on the red side where surely the others components, not resolved, are contained.
The radial velocities of the blueshifted and redshifted components agree with the values in the literature. For example, in a position near the star, EV92 detected two components at 64 and 78 km s and M-D14 reported two blueshifted components at 60 and 25 km s, all of them associated with the emission of NGC 6888. Recently, performing a 2D study in a region close to our position A, FM12 found that the emission of the [N II] 6584 Å line comes from plasma with radial velocities between 40 to 30 km s, in agreement with our value at position A.
The kinematical component at 15 km s is surely associated with the extended H II region around NGC 6888 and it has been also identified by EV92 and M-D14 with a velocity of 18 km s and 12 km s, respectively.
In general, we found radial velocities which follow the trend previously reported, i.e, the highest values appear near the star and they decrease as moving away from the centre of the nebula.
For the spectra at positions B and C we measured the line fluxes by fitting a gaussian profile to each kinematical component, but for the spectrum at position A, all the emission was integrated for each identified line. The errors assigned to our measurements correspond to a two sigma value over the mean noise measured at each side of the lines detected.
The fluxes were then dereddened with a logarithmic extinction coefficient c determined from the H/ H ratio using Cardelli, Clayton & Mathis (1989) extinction law. The uncertainties in c are the result of the flux-error propagation.
The c value at position A (0.70) agrees with the mean value of 0.76 obtained by FM12 in a region close to this observation. In the case of position B the extinction values are 0.36, 0.38 and 0.77 for the blueshifted, foreground and redshifted components, respectively. The first and second values are similar to the value 0.34 obtained by M-D14 for a region near our position B.
The high c value of the redshifted component seems to indicate that the emission coming from the back shell of NGC 6888 is more reddened surely by the presence of inner extinction. This is supported by the extinction value for WR 136, mag, reported by van der Hucht (2001). Using and , as derived from the extinction law mentioned above, we find c for the star. Given that WR 136 is embedded in NGC 6888, inner extinction could be contributing to redden its radiation.
Finally, observations of position C do not seem to be reddened. The only observations to compare with a neighbour zone are those performed by EV92. They found , which is the lowest value reported in the works cited above, indicating this region of NGC 6888 does not present high extinction.
In Table 2 the whole set of identified lines with their respective dereddened flux-integrated values are shown. In this table we present one set of values for position A, three (B, B and B) for position B, and two (C and C) for position C, where the subscripts “b”, “c” and “r” correspond to the blueshifted, central and redshifted components, respectively.
|Model||333 and are defined at||log||log||X||Y|
|This work star1||1.50||5.45||-4.63||1550||110||55||2||0.1||0.016||0.95|
|This work star2||2.10||5.40||-4.95||1550||90||48.30||1||0.175||0.068||0.91|
|This work star3||3.38||5.40||-4.95||1550||70||45.91||1||0.25||0.12||0.86|
|This work star4||3.38||5.40||-4.95||1550||70||42.26||2||0.25||0.12||0.86|
|Hamann et al.(2006)||3.34||5.40||-4.95||1600||70.8||…||1||0.25||0.12||0.86|
|Crowther & Smith (1996)||5.00||5.30||-3.91||1750||55.7||28|
4 The models
4.1 The stellar model
To model the observed stellar spectra of WR 136, in the UV and optical ranges, the CMFGEN code (Hillier & Miller, 1998) was used. CMFGEN (CoMoving Frame GENeral) is a non-LTE line-blanketed atmospheric code designed for spectral analysis of stars with stellar winds. It solves the radiative transfer and statistical equilibrium equations considering spherical symmetry. This code is adequate to model the atmosphere of massive stars with an extended wind that is, in general, larger than the hydrostatic radius of the star.
To start-up CMFGEN code it is required to set the stellar hydrostatic structure through input parameters such as the gravity and stellar radius or, in an equivalent way, to set a velocity law that approaches in depth to the hydrostatic structure. In this study we chose the later option and used the velocity law given by:
as proposed by Hillier et al. (2003). In expression (1) is the velocity at the hydrostatic stellar radius , is the photospheric velocity which controls the match between the hydrostatic atmosphere and the wind, is the terminal velocity, is the scale height of photosphere, in units of , which specifies the density structure at low velocities and is the so-called velocity law exponent. The values used for all our models for , and were 1 km s, 100 km s and 0.002, respectively, as it is suggested for WR stars444See the CMFGEN Documentation and WR models at kookaburra.phyast.pitt.edu/hillier/web/CMFGEN.htm. The values used for the models were 1 or 2, which are usual values for this kind of stars (Ignace, Quigley & Cassinelli, 2003).
Then, once the code calculates the hydrostatic values, it solves the radiative transfer equations in the wind zone and the result is the emergent flux (as emitted at 1 kpc from the observer) as a function of the wavelength. To fit the UV spectra, we used for the modelled spectrum a distance of 1.45 kpc, based on the HIPPARCOS parallax (van Leeuwen, 2007), and a final extinction value of 0.52 (or c = 0.75) which agrees with the value by van der Hucht (2001) shown in §3. In this way we fit the stellar continuum and the absorption feature around 2200 Å. For the optical spectrum, the model was analysed using the continuum normalized spectrum.
For the stellar models computed in this study, the input values used to constrain the model atmosphere for WR 136 were the mass-loss rate , the terminal wind velocity , the stellar luminosity , the hydrostatic radius of the star (at optical depth equal to 20) and the chemical composition.
According to the stellar wind theory, it is expected that the recombination line luminosity is proportional to the mass-loss rate through the relation , therefore the H recombination line was used to adjust the value of (Lamers & Cassinelli, 1999).
On the other hand, when we consider an extended stellar atmosphere, it is not possible to assign a surface temperature to the star as the temperature varies with the radius, therefore the temperature cannot be a global input value in the model. However, it is possible to define a temperature value for = 20 and such a parameter is set in the model through the luminosity and the stellar radius with the common expression , where is the Stefan-Boltzmann constant. Therefore, for constructing the model, the and values were changed to fit the line intensities corresponding to consecutive ionization degrees of the same element, such as the He I and He II lines.
While constructing a model and when the He lines temperature indicators (He II 4686 and He I 5876) pointed that we were very close to the best fit, we found that the predicted N III 4640 line was weaker than the observed line, while the predicted N IV 4058 was more intense than the observation, implying the temperature should be smaller (see Fig. 3). It has been stated in the literature that the ion N III is very dependent on the detail of its atom model (Hillier et al., 2003). We tried to correct this issue by using a more complex N III ion model, i.e., more levels were considered when decollapsing a superlevel. However, the problem was not solved although the fit of these nitrogen lines improved slightly.
The value of the terminal velocity was constrained by fitting the blue edge of the saturated C IV 1550 P-Cygni profile. The obtained value (1550 km s) fits well other P-Cygni profiles such as S IV 1063,1073 and P V 1118,1128 in the far UV spectrum of the star (see Fig. 3) and it is in excellent agreement with values reported in previous works. Prinja, Barlow & Howarth (1990) used the common diagnostic line, the saturated UV C IV resonance doublet at 1548, 1550 Å and obtained km s; Eenens & Williams (1994) derived km s by using the He I 10830 Å line width; HGL06 measured km s by fitting their model to their observed spectrum; Crowther & Smith (1996) determined a km s; finally, through an IR study, Ignace et al. (2001) used the Ca IV 3.207 m forbidden line width, reporting km s. As our fit was done measuring the velocity at the edge of the satured component of the blue wing, instead of measuring it where the wing reaches the continuum, this renders a terminal velocity km s lower than when measured at the continuum. This excess in velocity can be explained invoking a turbulent velocity or by the clumping factor of the wind.
The chemical composition in the models was initially set at typical WN composition (van der Hucht, Cassinelli & Williams, 1986) and it was modified to generate some of the models shown here. The final H and He abundances (X, Y by mass-fraction) employed for all the stellar models are shown in Table 3, while for the other elements such as C, N, O, S and Fe, the abundances derived from CMFGEN are shown in Table 6 to compare then with the resulting nebular abundances from the modelled nebula.
As the modelled spectral features do not respond linearly to the input parameters, once a value is found it is not necessarly fixed but only constrained while changing the other input values. Therefore we calculated a grid of models in order to find the more appropriated models which reproduce most of the observed spectral signatures.
To properly fit the luminosity through the stellar continuum without changing the stellar temperature, it was necessary to use the transformed radius
(where ) because it has been shown that two models with the same and values have the same spectral behaviour (Schmutz, Hamann & Wessolowski, 1989).
From the grid of models calculated by varying the input parameters we obtained star1 which fits very well most of the stellar spectral characteristics. This model is presented in Table 3. In addition we also show in this table models from the literature, as HGL06 and Crowther & Smith (1996) ones. The mass-loss rate value of star1 is slightly higher but of the same order as HGL06 value and it is lower than the value published by Crowther & Smith (1996). We are confident in the mass-loss rate value of our models because as it can be seen in Figure 3 the recombination lines, like the Balmer’s series in the optical spectrum, are very well fitted.
In Figure 3 we can see that the He II 1640 and 5411 lines are overestimated in star1 only by a factor . It is important to say that it has been stated in the literature that the He lines in an observed spectrum cannot be fitted simultaneously (Hillier et al., 2003). From Figure 4 we can also see that star1 reproduces in general very well the observed spectra from the UV to the optical wavelength range.
The main difference between star1 and the model reported by HGL06 is the stellar radius which in our case has a low value of 1.5 (compared to the value of 3.34 by HGL06) leading to a higher stellar temperature of 110 kK. Additionally the chemical composition, represented by the values of X and Y, is different. These differences have important implications when modelling the nebula around the star (see §4.2.2) and are indicating that it could exist more than one solution for a stellar model of a WR star as both, star1 and HGL06 model, reproduce similarly well the stellar spectrum (see Fig. 68 by HGL06). This apparent degeneracy has already been remarked in a study of WR stars modelling by Hillier (1991). In this work the origin of this degeneracy is attributed to our misunderstanding of the dynamics in the inner wind regions.
To investigate this degeneracy we searched for other CMFGEN models that could fit the stellar observations. The results are plotted in Figure 4 where the general behaviour of the models is compared with the observed spectra.
In a simple manner we used interpolated values between the model star1 and the HGL06 model for the radius, the clumping factor and the hydrogen and helium fractions by mass, X and Y. The result is model star2 which has an effective temperature of 90 000 K (its characteristics are presented too in Table 3). Model star2 reproduces the stellar spectrum similarly well as star1. It is presented in Fig. 4, in green.
We computed two additional models. The first is star3 which basically uses HGL06 model parameters but it was calculated with the CMFGEN code. Model star3 reproduces very well the results of HGL06 model, which shows that CMFGEN and the Potsdam WR model atmosphere code are equivalent. The last one is model star4 which is almost the same as star3 but calculated with the value equal to 2, leading to a lower value because the wind region becomes more extended. Models star3 and star4 reproduce fairly well the observed stellar continuum and many spectral features (lines yellow and blue in Fig. 4) although the emission lines indicative of the stellar temperature are worse fitted than with star1(see Fig. 3). All the stellar model characteristics are presented in Table 3.
From this, it is then clear that there is not a unique stellar model reproducing the stellar spectrum. As said above, we are in presence of a degeneracy occurring when modelling WR stellar atmospheres and evidently more constraints are required to decide for the best model. Additional constraints can be found by analyzing the ionized nebula.
|Parameter||Position A||Position B||Position C||EV92 Position||Position B2 (FM12)||M-D14|
|T[O III](K)||…||…||…||9500:||(6855639)555Temperature derived from Pérez-Montero & Contini (2009) formula.||9050810|
|n[S II] (cm)||240||180||215||371||108||480:|
|n[O II] (cm)||190||…||…||301||…||31030|
|N||Galavis, Mendoza & Zeippen (1997)||Tayal (2011)|
|Wiese, Fuhr & Deters (1996)|
|O||Zeippen (1982)||Pradhan et al. (2006)|
|Wiese, Fuhr & Deters (1996)||Tayal (2007)|
|O||Wiese, Fuhr & Deters (1996)||Aggarwal & Keenan (1999)|
|Storey & Zeippen (2000)|
|S||Podobedova, Kelleher & Wise (2009)||Tayal & Zatsarinny (2010)|
|S||Podobedova, Kelleher & Wise (2009)||Galavis, Mendoza & Zeippen (1995)|
|Ar||Kaufman & Sugar (1986)||Galavis, Mendoza & Zeippen (1995)|
4.2 The nebular model
4.2.1 Abundances and plasma analysis
The starting point to make a photoionization model for NGC 6888 is through the spectral analysis of the observations taken for different zones on this nebula.
As we said before, the fluxes obtained from the spectrum at position A contains not only the nebular emission but also the foreground emission. This foreground contribution is known from the spectra at positions B and C, where the flux values corresponding to the 15 km s component are similar. This enforces the assumption that these components have their origin in the interstellar medium surrounding the nebula and therefore, we can assume that these foreground fluxes are similar on the whole nebula NGC 6888. By comparing the total fluxes from spectrum A with these values it is found that the foreground contribution is only around 1 per cent for all the identified lines, except for the faint [O III] 4959,5007 lines where this flux contribution reaches 5 per cent. We can therefore assume that the integrated fluxes of the spectrum at this position represent well the nebular behaviour.
From the data in Table 2 we determined the electronic temperature and density from the available diagnostic ratios. In the case of positions B and C only the intense blueshifted components were used. In all the cases, only derived from the [N II] line flux ratio could be determined because the [O III] 4363 line (crutial for determining [O III]) was not detected in the spectrum. The value comes from the usual [S II] diagnostic ratio because the lines coming from this ion are present in all the spectra. The temperature and density as well as the ionic abundances, derived by using these physical conditions and the line intensities from Table 2, are shown in Table 4. All these quantities as well as the errors were determined with the software PyNeb 1.0.1b9 (Luridiana, Morisset & Shaw, 2015) using the atomic data shown in Table 5. The total abundances were only derived for position A where [O II] 3726,3729 were detected. The ionization correction factors (ICF) by Kingsburgh & Barlow (1994) were used with this purpose. The results are listed in Table 4.
In addition to our results, in Table 4 we are including the results obtained by EV92 for their blueshifted component, and the data by FM12 corresponding to the inner zone observed by them, i.e., that labeled in their paper as B2. Both zones (EV92, and B2 FM12) are clearly not associated with the rim region of NGC 6888 as it is shown in our Fig. 1 where their positions are marked. Also the results by M-D14, for their position in the rim, have been listed.
From these data we can make a direct comparison of the chemistry in different zones of NGC 6888, as obtained by different authors.
In Table 4 it is found that our values of [N II] electronic temperature and density for positions A, B, and C agree with the values derived by EV92, FM12 and M-D14. According to our observations, the plasma in position B has the highest electronic temperature (although with a large uncertainty); similar values were reported by FM12 for the gas located farther than the edge of NGC 6888 and although their regions and our B position are not spatially close an explanation for the high temperature in position B could be that the main nebular emission comes from a low-density component of NGC 6888 if we assume this ring nebula is made of two gaseous components (see §4.2.2).
The abundances we obtain for position A in the rim are in good agreement with the previously published abundances (except for position B2 of FM12, see below). The fact that the values of O/H and N/H are compatible within the error bars for spatially different regions of the nebula indicates that the hypothesis of a chemically homogeneous nebulae need to be explored.
The abundances found for regions slightly farther than the edge of NGC 6888 (positions B1, E1, E2 and MB3 of FM12) are not taken into account because, according to FM12, there is no nitrogen enrichment there as this component was probably ejected during an early evolutionary stage of the star WR 136 as proposed by García-Segura, Langer & Mac Low (1996).
The region B2 of FM12 presents an O/H abundance much larger than the values of the other positions. Note that FM12 did not obtain directly the electron temperature for the [O III] zone in their B2 region, but rather defined it using the relation (2) from Pérez-Montero & Contini (2009) which was derived for other kind of objects as galaxies and H II regions, and given that NGC 6888 is an extented nebula which has a complex structure, probably made of two gaseous components, the changes in temperature we detect across the nebula cannot be assumed to follow the same behavior as found for global observations, since in addition we must deal with the ionization structure.
Therefore, their very low [O III] temperature value (not found by EV92 who directly determined it in a region close to B2) explains the high value they consequently obtain for O/H ionic abundance and their high total O/H. If we use the T[O III] of EV92 the O abundance is only 8.01 and the total abundance of oxygen is 8.29 which is in good agreement with the rest of the values shown in Table 4.
In addition FM12 obtain a low log N/O ratio of 0.33 for this zone, which is also questionable, as this is contrary to the chemical enrichment predicted (log N/O 0.20) by the stellar evolution models (discussed below), where N-enrichment should be detected also near the star. This is actually the case for the also inner position observed by EV92 who obtained log N/O = 0.29. Note that the N/O abundance ratio is defined as equal to the N/O ratio (ICF=1), while we will see in §5.3 and Fig. 8 that this is not the case for observations obtained through small slits in the inner zone of this extended nebula.
The values by FM12 could be indicating a chemically inhomogeneous nebula, as these authors suggested. However given the above discussion, it is possible instead to assume that the main structure of NGC 6888 has a homogeneous chemical abundance. The existence of such homogeneous nebulae are supported by the stellar evolution models for WR stars. For example, from Figure 16 of Toalá & Arthur (2011) we can see that along the evolution of a WR star, the atmospheric log N/O ratio reaches a high value of 0.20 in 10 000 years and remains almost constant for around 300 000 years, in the case of a non-rotational star with initial mass of 40 .
In a similar way M-D14 followed the time evolution of the N/O ratio for NGC 6888, concluding that the material of the rim was ejected during the Red Super Giant phase highly nitrogen-enriched and the subsequent WR phase did not provide additional enrichment (see their Figure 5).
Therefore, given that the dynamical age of NGC 6888 has been estimated to be of about 20,000 - 40,000 yr (van Marle, Langer & García-Segura, 2005), the idea of a chemically homogenous nebula for NGC 6888 is consistent with the evolutionary models of a massive star becoming a WR star.
In the following we explore this possibility, by means of a photoionization model.
4.2.2 Photoionization model of NGC 6888
To reproduce the optical spectra of NGC 6888 we computed a photoionization model using the pyCloudy package (Morisset 2013) based on the 1D photoionization code Cloudy (Ferland 2008). It allows us to generate 2D images for each emission line from a pseudo-3D model and to apply a mask on them, reproducing the effect of the finite aperture through which the observations were performed.
For our nebular model we used the version C10.0 of Cloudy, last described by Ferland et al. (2013). We assumed, as it is apparent from the optical images (Gruendl et al., 2000; Moore, Hester & Scowen, 2000), that the NGC 6888 projected morphology is the result of two gaseous components, as shown in Figure 5: one bright, granular region emitting mainly low ionization lines (e.g., [N II]), and the other one, more diffuse, mainly emitting the higher ionized emission lines (e.g., [O III]).
The denser (brighter) component has an ellipsoidal shape with an inner minor semi-axis of 80 arcmin, a minor to major semi-axis ratio of 0.33 with a H density of 400 cm; the other component has a lower density of 1 cm, and a spherical geometry with a central cavity with a radius of 20 arcmin. In addition we set a constant density profile for each component, with a filling factor = 1.0 for the low-density component and = 0.1 for the high-density component.
As we aim to build a consistent model of the system we used as ionizing sources the stellar atmosphere models of the star WR 136, as calculated with the CMFGEN code and described in the previous section. We therefore ensure to have well constrained stellar parameters, the ionization parameter and the ionization degree of the nebula through the stellar luminosity and temperature.
In a simple way, we calculate a chemically homogeneous model for NGC 6888 and we verify if observations can be reproduced.
The chemical abundances were first set to the values derived for position A in §4.2.1 (Table 4) which were completed with the abundances of C, Ne and Fe published by M-D14. Then the abundances were slightly modified until the photoionized model reproduces the observed line intensities.The final values which describe better the chemical composition in NGC 6888 are given in Table 6.
The coordinates for the slit positions, from where the model integrated fluxes were obtained, were set in angular units with respect to the central star adopting a distance of 1.45 kpc for NGC 6888 (van Leeuwen, 2007). In the optical images the star does not seem to be in the geometrical center of the nebula, but as we are only modelling the northwest part of it, the geometry of the model is an ellipsoid with the geometrical center on the central star position and whose edge matches the observed region.
The final intensities of the emission lines are obtained by summing up the contributions of each component, assuming that the high-density region weight on the total intensity is . Thus:
where and are the line intensity values from the low-density and high-density components, respectively. Once the abundances were fixed, only this weight value was changed, from one position to another, to fit the observed fluxes.
|w (see text)||7.6 10||1.7 10||3.8 10|
The comparison between the model predictions and the observed dereddened line fluxes is presented in Table 7 and Table 8 corresponding to the photionization model results by using star1 and star3 SED’s as ionizing source, respectively. The q parameter, in these tables, shows the goodness of our model and is defined following Morisset & Georgiev (2009) as:
where is the value returned by the model, is the observed flux and is the accepted tolerance in dex of this observable. We define the tolerance to be
where is the measured error. When a line flux is well reproduced the q value is between 1.0 and 1.0.
The results in Table 7 show that, under the assumptions made for this model and using star1 as ionizing source, the observed line fluxes are not well reproduced. In particular the lines from highly ionized elements (O, Ar) are extremely overpredicted. The best case is for the values in position A where only half of the whole line fluxes (the low excitation ones) detected in this region are well reproduced. The model results are worse for the other positions.
The fact that the ionization of the nebula is strongly overpredicted could be due to: 1) a too high ionization parameter 666 where is the ionizing photon rate, the hydrogen number density and the distance to the ionizing source., 2) a too high stellar effective temperature. Notice that this model was obtained using only the high density component (). Adding any contribution of the low density component would made the situation even worse.
As , and are well constrained through the HIPPARCOS distance to the object, its observed angular size, the observed stellar luminosity and the nebular density coming from the plasma diagnostic, then it follows that is also well constrained. The only remaining way to reproduce the observed low value of the [O III]/[O II] line ratio is to consider a star with a temperature lower than the value derived for model star1, i.e., we must not deal with the number of ionization photons but rather with their energy.
The fit improves very much when we use the model star3 which has 70 000 K. The ionization of the dense gas is even lower than observed, leading to the full use of our hypothesis of the 2-components model. The predicted emission lines from the corresponding photoionization model are compared with observations in Table 8. It can be seen that almost all the fluxes are well reproduced, indicating that the star WR 136 emits the UV ionizing photons corresponding to a star with this low effective temperature.
5.1 Effective temperature and mass-loss rate of WR 136
As it can be seen in Figure 4 and it was discussed in §4.1, at least three stellar models, (star1, star2 and star3) with different and therefore a different fit similarly well the stellar observations, in particular the stellar continuum and several emission lines. The differences among these models appear mainly in the He lines, where the He II/He I line ratio is well reproduced when considering an effective temperature of 110 000 K but it is understimated when the effective temperature is reduced, i.e., at a lower stellar temperature the stellar He I lines are stronger and the stellar He II lines become weaker than observed. As we discussed before, it seems that there is not a unique solution for the stellar model.
In Figure 6 we show a comparison between the UV region of the stellar models from Table 3. We also overplot the UV behaviour of two blackbodies with temperatures of 70 000 and 46 000 K which correspond to the temperatures and of model star3.
We can see that the key difference between the stellar models comes from the far UV region of the spectra, whose UV photons are absorbed by the nebula (between 1.0 and 7.0 Ryd). This difference is well appreciated when we compare the models with the blackbodies. From this figure it is clear that a blackbody with the same temperature than the modelled spectra, 70 000 or even 46 000 K, shows an excess of UV photons which would lead to a much more ionized nebula. It is also evident the high absorption of the UV photons emitted by the star, due to the presence of the stellar wind.
Owing to there are not stellar observations in the far UV to discriminate which stellar model is reproducing the real spectral energy distribution of the star WR 136, it is necessary to add an external restriction. In this case, such additional restriction is the associated nebula NGC 6888 that is highly sensitive to this UV ionizing radiation.
In addition, the combination of the temperature and mass-loss effects in the stellar wind, considered for this model, can reproduce the line profiles coming from recombination process. Therefore, if we are close to a good value for , it is the same for the mass-loss rate, . This is important because represents a good constrain for the hydrodynamical simulations.
As in this study we have found a well constrained value for (log = 4.95), which is similar to the values found by other authors (Table 3), we can rule out the assumption made recently by Zhekov & Park (2011) where they conclude that it is necessary a 3 - 4 times lower than observed for this type of stars to reproduce the X-ray flux of the emitting hot bubble.
5.2 Is NGC 6888 mainly photoionized?
The presence of a very fast stellar wind coming from the star WR 136 could imply that part of the ionization of the gas is due to shocks. Diagrams made by Veilleux & Osterbrock (1987) and Baldwin, Phillips & Terlevich (1981) are usually employed to determine if a line emitting region is mainly ionized by OB stars (like H II regions) or by other source (AGN, LINERS, shocks). See for example Allen et al. (2008) for the case of shocks.
In Figure 7 we show an [O III] 5007/H vs. [N II] 6584/H diagram. The solid line shows the separation between global photoionization models for Star Forming regions (NSF galaxies) and AGNs as defined by Kauffmann et al. (2003). The dashed line represents the same separation, but shifted on the log ([N II] 6584/H) axis by 0.96 dex, which is the difference between the N/O adopted by Kauffmann et al. (2003) and the N/O of NGC 6888 deduced in this work. The three black diamonds correspond to observations for positions A, B and C and are obviously located to the right side of the solid line, but well inside the photoionized zone with N/O ratio similar to NGC 6888 value.
The effect of a limited aperture size, which is not including the whole nebula is illustrated by the position of the color points in our Fig. 7. They show the emission line ratios of every spaxel of our 2D images obtained from our models. The colors correspond to different values of (the weight of the dense component), the blue symbols are obtained for the case of considering only the low-density component, and the red ones for the case of considering only the high-density component, the green symbols being an intermediate situation. As we can see, the points spread in an ample region of the diagram. The black square symbol is the position in the diagram of the global model, as if it were fully included in the slit.
Therefore if just a small portion of a whole nebula is observed, one must be careful in using diagnostic diagrams to decide if the emission of a nebula is caused by the presence of shocks, because as we demonstrated in Fig. 7 only photoionization and an enhanced N abundance are sufficient to obtain points located out of the regions represented by global photoionized models (solid line).
5.3 Oxygen and nitrogen abundances in NGC 6888
To deeper investigate our assumption of chemical homogeneity of NGC 6888, in Fig. 7 we added the observations by EV92, M-D14. It is found that they lie inside or close to our model points. Also we added the observations by FM12 (marked as black dots and stars) resulting from their 1D analysis of different regions of NGC 6888. We can see that their points show a large spread in this diagram. To explain the behaviour of their observations and to derive values of the N and O abundances, FM12 performed a grid of photoionization models using the code Cloudy varying the ionization parameter and the N/O ratio. As they used global models for each observation without taking into account that they were observing only a part of the nebula they concluded that there are important differences in the N/O ratio across NGC 6888.
However, our figure shows that this is not necessarily the case. In particular, let us focus on the data of FM12 marked as black dots, which correspond to their regions inside NGC 6888 with large nitrogen enrichment. We can see that all these points, but the point labeled X2 in their paper, are close to our model points and are following their trend. The rest of the points (stars) show a large dispersion compared to our model points and it is not expected for the model to reproduce their positions because they belong to regions located out from the edge of NGC 6888 with low (solar) nitrogen abundance.
Therefore, Figure 7 illustrates the fact that all the observed points within the nebula can be reproduced by a rather simple chemically homogeneous model.
To continue the discussion about the chemical homogeneity in NGC 6888, we analyse the behaviour of the ICF used to determine the N/O ratio. Let us remember that it is generally assumed that N/O is equal to N/O.
From our photoionization model, we are able to determine in each region of the nebula the ICF needed for the N/O ratio to obtain the N/O value. We show in Figure 8 a map of this ICF for three different values of w: 10, 10 and 10. It is evident that for a given w, this ICF varies over the nebula showing a difference as large as a factor of two. In adittion, we can see that the value of this ICF changes as the contribution of the high-density component is modified.
For example, in the case when w = 10, the ICF for position A is around 1.0. In general, we can see in Figure 8 that for positions closer to the star the ICF is around 1.5 or larger. Therefore, in the zones near to the star an ICF greater than one is needed to estimate the N/O ratio. Hence, this indicates that N/O could be understimated when it is simply assumed that N/O is equal to N/O. In this way, to calculate the nitrogen abundance for position B2 of FM12 and for EV92 a better value for this ICF is 1.5. Once this ICF is used for these positions the N/O ratio increases relative to the values in Table 4, to 0.47 for EV92 and 0.15 for B2 FM12. Therefore these values are still discrepant. In particular, for the position B2 the N/O ratio reaches a value of 0.15 lower than the others for the rim, which is contrary to the predictions of stellar evolution models.
The differences seen of the ICF in different zones of the map are due to the combination of two effects: the ionization structure of the nebula and the relative contribution to the fluxes of the high density clumpy region.
A similar problem with the ICFs is affecting the results of Stock & Barlow (2014) who used Herschel observations of the [O III] 88 m and [N III] 57 m lines, emitted in two positions of NGC 6888, to determine the N/O abundance ratio under the assumption that N/O is equal to N/O. For NGC 6888 this is a very rough approximation because only a small fraction of the nebula was included in their slits and Stock & Barlow (2014) did not take into account the nebular ionization structure, which indicates that the ion N is only dominant over N in the zones near to the star.
The previous discussion illustrates very well that detailed modelling (even the simplest as we do here: 2 almost roundish homogeneous components) is a very powerful way to identify effects that, if not considered, can easily drive the observer to wrong conclusions.
We computed several atmosphere models using the code CMFGEN to reproduce the stellar observations of WR 136. The characteristics of the star such as , and agree with the values found by other authors, i.e., HGL06 and Crowther & Smith (1996). We also derived the stellar abundances of He, C , N, O, Si, P, S and Fe (listed in Table 6) which agree with HGL06’s model.
We find that models with very different effective temperature, from 70 000 K to 110 000 K, can reproduce fairly well the stellar continuum and several emission lines, although the He lines, indicative of temperature, are worse reproduced by the cooler models. Therefore there is a degeneracy in atmosphere models, than apparently, as pointed out by Hillier (1991), is due to our misunderstanding of the structure of the wind velocity field.
To decide which stellar model produces the correct UV amount of ionizing photons coming from WR 136, we computed a photoionization model for the associated ring nebula NGC 6888, using the SED produced by the stellar models. This led us to conclude that WR 136 follows better the spectral behaviour of a star with a of 70 kK.
The photoionization model was built assuming chemical homogeneity and including two components: a low-density region enclosing a denser region.
When these two gaseous components are considered and assuming only photoionization, we were capable to reproduce the observed ionization structure in both shells, contrary to the presence of shocks as the probable mechanism producing the observed ionization structure of NGC 6888. We cannot immediatly rule out the importance of shock effects, but given that some diagnostic diagram such as the [O III] 5007/H vs. [N II] 6584/H indicate emission of H II type, and that we can reproduce several emission features of the nebula by only assuming photoionization, we established that the photoionization heating is more important than the presence of shocks. This latter idea was already mentionned by EV92 combining the diagnostic diagrams of Sabbadin, Minello & Bianchini (1977) with their observations in a zone close to the central star and shown in Figure 1.
The nebular abundances derived from our photoionization model are similar to the values derived from our position A and the values reported by EV92 and M-D14.
It is worth to notice that both the stellar atmosphere of WR 136 and the nebular plasma present a large N- and C-enhancement and an O depletion, as it can be seen in Table 6. This contributes to the idea that the observed gas was ejected by the central star through the mass-loss driven by the stellar wind.
Considering only photoionization and a chemically homogeneous nebula, we were able to explain the pretended differences found in the N/O ratio (derived from N/O or N/O in different zones of the nebula) by FM12 and Stock & Barlow (2014), who proposed a chemically inhomogeneous nebula. These pretended differences are due to effects of the ionization structure within the nebula. With our simple model we show that the ICFs for global nebulae should not be use for observations of small zones observed in extended nebulae.
From our results we have learned that there are several stellar models with different set of values that could reproduce the observables, like the models star1 and star3 presented in this study. Therefore external constraints are crucial to determine which atmosphere models represent the real stellar properties. The nebulae around stars can help to deal with the degeneracy mentioned above.
Finally, in this work we show the advantages of modelling simultaneously the stellar and nebular spectra of a Wolf-Rayet star plus nebula system. We verify whether a stellar atmosphere model of the star WR 136 provides a good description of the stellar observations and is able to satisfactorily reproduce the nebular characteristics such as the ionization degree. This allow us to better constraint all the physical parameters of both the star and the nebula, including their physical conditions and their chemical abundances.
This work received partial finantial support from UNAM PAPIIT IN109614 and CONACyT CB-2010/153985. The FUSE and IUE data presented in this paper were obtained from the Multimission Archive at the Space Telescope Science Institute (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. We thank G. Koenigsberger, A. Peimbert and J. García-Rojas for their helpful comments.
- Aggarwal & Keenan (1999) Aggarwal K. M., Keenan F. P. 1999, ApJS, 123, 311
- Allen et al. (2008) Allen Mark G., Groves Brent A., Dopita Michael A., Sutherland Ralph S., Kewley Lisa J., 2008, ApJS, 178, 20
- Baldwin, Phillips & Terlevich (1981) Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
- Cardelli, Clayton & Mathis (1989) Cardelli Jason A., Clayton Geoffrey C., Mathis John S., 1989, ApJ, 345,245
- Chu, Treffers, & Kwitter (1983) Chu Y.-H., Treffers R. R., Kwitter K. B. 1983, ApJS, 53, 937
- Crowther (2007) Crowther P. A. 2007, ARA&A, 45, 177
- Crowther & Smith (1996) Crowther P. A., Smith L. J. 1996, A&A, 305, 541
- Eenens & Williams (1994) Eenens P.R.J., Williams P. M., 1994, MNRAS, 269, 108
- Esteban & Vílchez (1992, hereafter EV92) Esteban C., Vílchez J. M. 1992, ApJ, 390, 536, (EV92)
- Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
- Ferland et al. (2013) Ferland G. J., Porter R. L., van Hoof P. A. M., Williams R. J. R., Abel N. P., Lykins M. L., Shaw G., Henney W. J., Stancil P. C. 2013, RMAA, 49, 13
- Fernández-Martín et al. (2012, hereafter FM12) Fernández-Martín A., Martín-Gordon D., Vílchez J. M., Pérez Montero E., Riera A., Sánchez S. F., 2012, A&A, 541, 119, (FM12)
- Galavis, Mendoza & Zeippen (1995) Galavis M.E., Mendoza C., Zeippen C. J. 1995, A&AS 111, 347
- Galavis, Mendoza & Zeippen (1997) Galavis M. E., Mendoza C., Zeippen C. J. 1997, A&AS, 123, 159
- García-Segura, Langer & Mac Low (1996) García-Segura G., Langer N., Mac Low M.-M. 1996, A&A, 316, 133
- Gray (1999) Gray R. O. 1999, ascl.soft10002G
- Gruendl et al. (2000) Gruendl R. A., Chu Y.-H., Dunne B. C., Points S. D. 2000, AJ, 120, 2670
- Hamann, Grafener & Liermann (2006, hereafter HGL06) Hamann W.-R., Grafener G., Liermann A., 2006, A&A, 457, 1015, (HGL06)
- Hillier (1991) Hillier D. J., 1991, IAUS, 143, 59
- Hillier & Miller (1998) Hillier D. J., Miller D. L. 1998, ApJ, 496, 407
- Hillier et al. (2003) Hillier D. J., Lanz T., Heap S. R., Hubeny I., Smith L. J., Evans C. J., Lennon D. J., Bouret J. C. 2003, ApJ, 588, 1039
- Ignace et al. (2001) Ignace R., Cassinelli J. P., Quigley M., Babler B. 2001, ApJ, 558, 771
- Ignace, Quigley & Cassinelli (2003) Ignace R., Quigley M. F., Cassinelli J. P., 2003, ApJ, 596, 538
- Kauffmann et al. (2003) Kauffmann G., Heckman T. M., Tremonti C., Brinchmann J., Charlot S., White S. D. M., Ridgway S. E., Brinkmann J., Fukugita M., Hall P. B., Zeljko I., Gordon T. R., Donald P. S. 2003, MNRAS, 346, 1055
- Kaufman & Sugar (1986) Kaufman V., Sugar J. 1986, JPCRD, 15, 321
- Kingsburgh & Barlow (1994) Kingsburgh R. L., Barlow M. J. 1994, MNRAS, 271, 25
- Kwitter (1981) Kwitter K. B., 1981, ApJ, 245, 154
- Lamers & Cassinelli (1999) Lamers H. J. G. L. M., Cassinelli J. P. Introduction to stellar winds / Cambridge U Press, 1999
- Levine & Chakrabarty (1994) Levine S., Chakrabarty D. 1994, Tevhnical Report MU-94-04, Instituto de Astronomï¿½a, Universidad Nacional Autónoma de México.
- Luridiana, Morisset & Shaw (2015) Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573,42
- Mendoza (1983) Mendoza C. 1983, I.A.U. Symp. No. 103, p. 143
- Mesa-Delgado et al. (2014, hereafter M-D14) Mesa-Delgado A., Esteban, C., García-Rojas J. Reyes-Pérez J. Morisset C., Bresolin F., 2014, ApJ, 785, 100, (M-D14)
- Moore, Hester & Scowen (2000) Moore Brian D., Hester J. Jeff, Scowen Paul A., 2000, AJ, 119, 2991
- Morisset (2013) Morisset C. 2013. sites.google.com/site/cloudy3d
- Morisset & Georgiev (2009) Morisset C., Georgiev L., 2009, A&A, 507, 1517
- Pérez-Montero & Contini (2009) Pérez-Montero E., Contini T. 2009, MNRAS, 398, 949
- Podobedova, Kelleher & Wise (2009) Podobedova L. I., Kelleher D. E., Wiese W. L. 2009, JPCRD, 38,171
- Porter et al. (2013) Porter R. L., Ferland G. J., Storey P. J., Detisch M. J. 2013, MNRAS 433, 89
- Pradhan et al. (2006) Pradhan A. K., Montenegro M., Nahar S. N., Eissner W. 2006, MNRAS, 366, 6
- Prinja, Barlow & Howarth (1990) Prinja Raman K., Barlow M. J., Howarth Ian D. 1990, ApJ, 361,607
- Sabbadin, Minello & Bianchini (1977) Sabbadin F., Minello S., Bianchini A. 1977, A&A, 60, 147
- Schmutz, Hamann & Wessolowski (1989) Schmutz W., Hamann W.-R., Wessolowski U., 1989, A&A, 210,236
- Skinner et al. (2010) Skinner Stephen L., Zhekov Svetozar A., Güdel Manuel, Schmutz Werner, Sokal Kimberly R., 2010, AJ, 139, 825
- Stock & Barlow (2014) Stock D. J., Barlow M. J., 2014, MNRAS, 441, 3065
- Storey & Zeippen (2000) Storey P.J., Zeippen C.J., MNRAS, 312, 813, 2000
- Tayal (2007) Tayal S. S., 2007, ApJS, 171, 331
- Tayal (2011) Tayal S. S. 2011, ApJS, 195, 12
- Tayal & Zatsarinny (2010) Tayal S. S., Zatsarinny, O. 2010, ApJS, 188, 32
- Toalá & Arthur (2011) Toalá J. A., Arthur S. J. 2011. ApJ, 737, 100
- van der Hucht (2001) van der Hucht K.A., 2001, New Astron. Rev., 45, 135
- van der Hucht, Cassinelli & Williams (1986) van der Hucht K.A., Cassinelli J.P., Williams P.M., 1986, A&A, 168, 111
- van Leeuwen (2007) van Leeuwen F., 2007, A&A, 474, 653.
- van Marle, Langer & García-Segura (2005) van Marle A. J., Langer N., & García-Segura G. 2005, A&A, 444, 837
- Veilleux & Osterbrock (1987) Veilleux S., Osterbrock D., E., 1987, ApJS, 63, 295
- Wiese, Fuhr & Deters (1996) Wiese W. L., Fuhr J. R., Deters T. M., JPCRD Monograph(7) pp. 1-522 (1996)
- Zeippen (1982) Zeippen C. J. 1982, MNRAS, 198, 111
- Zhekov & Park (2011) Zhekov S. A., Park S. 2011, ApJ, 728, 135