Modelling of the B-type binaries CW Cep and U Oph
Key Words.:stars: keyword 1 – stars: keyword 2 – stars: keyword 3
Context:Intermediate-Mass stars are often overlooked as they are not supernova progenitors but still host convective cores and complex atmospheres which require computationally expensive treatment. Due to this, there is a general lack of such stars modelled by state of the art stellar structure and evolution codes.
Aims:This paper aims to use high-quality spectroscopy to update the dynamically obtained stellar parameters and produce a new evolutionary assessment of the bright B0.5+B0.5 and B5V+B5V binary systems CW Cep and U Oph.
Methods:We use new spectroscopy obtained with the Hermes spectrograph to revisit the photometric binary solution of the two systems. The updated mass ratio and effective temperatures are incorporated to obtain new dynamical masses for the primary and secondary. With these, we perform isochrone-cloud based evolutionary modelling to investigate the core properties of these stars.
Results:We report the first abundances for CW Cep and U Oph as well as report an updated dynamical solution for both systems. We find that we cannot uniquely constrain the amount of core boundary mixing in any of the stars we consider. Instead, we report their core masses and compare our results to previous studies.
Conclusions:We find that the per-cent level precision on fundamental stellar quantities are accompanied with core mass estimates to between %. We find that differences in analysis techniques can lead to substantially different evolutionary modelling results, which calls for the compilation of a homogeneously analysed sample to draw inference on internal physical processes.
The model independent estimates of the absolute dimensions of and distances to stars provided by eclipsing binary systems serve as a fundamental calibrator in modern astrophysics. In the best cases, such systems offer dynamical mass and radius estimates better than one per-cent (torres2010). Such precise measurements combined with the powerful constrains of co-evolution and identical initial chemical composition have allowed the thorough investigation of the importance of rotation in stellar evolutionary theory (brott2011a; brott2011b; ekstrom2012; deMink2013; schneider2014; ekstrom2018), the calibration of pre through post main-sequence evolution (torres2013; higl2017; beck2018b; kirkbyKent2018), the critical investigation of magnetic fields in stars (takata2012; grunhut2013; torres2014a; pablo2015; kochukhov2018; wade2019), the calibration of distances (guinan1998; ribas2000c; ribas2005; hensberge2000; bonanos2006; pietrzynski2013; gallenne2016; suchomska2019), the investigation of abundances and rotational velocities (pavlovski2005; pavlovski2009a; pavlovski2009b; pavlovski2018; simondiaz2017), as well as the calibration of asteroseismic modelling (decat2000; decat2004; aerts2004; schmid2015; schmid2016; beck2018a; beck2018b; johnston2019). Additionally, the advent of such precise measurements has led to the uncovering of the reported systematic discrepancy between masses obtained via dynamics or empirical spectral relations and fitting theoretically calculated evolutionary tracks (herrero1992; ribas2000b; tkachenko2014). This discrepancy has served as the centrepiece of intense debate over the importance of convective core boundary mixing in stellar evolution theory (ribas2000b; torres2010; torres2014b; tkachenko2014; stancliffe2015; claret2018; constantino2018; johnston2019).
In general, both element and angular momentum transport processes throughout a star are poorly calibrated (aerts2019). It is a well known short-coming of most 1-D theoretical descriptions of convection that convective boundaries are not well described (hirschi2014). Proposed as a means to remedy this short-coming, the inclusion of convective core overshooting as a way to increase near-core mixing in evolutionary models is now highly debated, with several competing studies claiming that models with and without overshooting can reproduce observed binaries across different mass ranges and evolutionary stages (andersen1990; schroder1997; pols1997; claret2007; stancliffe2015; claret2016; claret2017; higl2017; constantino2018). Convective core overshooting is a phenomenon theoretically predicted in intermediate- to high-mass stars with a convective core where the inertia of a convectively accelerated mass element propels said mass element beyond the convective boundary described by the Schwarzchild stability criterion into the stably stratified radiative region (zahn1977; roxburgh1978; zahn1991; maeder2009). The mathematical form of the implementation into stellar evolutionary codes is not universally agreed upon, with different forms having been successfully used to describe both binary (ribas2000b; guinan2000; tkachenko2014; claret2016; claret2017) and asteroseismic observations (briquet2007; moravveji2015; moravveji2016; vanReeth2016; johnston2019). To date, two such descriptions have been implemented in 1-D stellar evolution codes: i) convective penetration where the temperature gradient in the overshoot region is the adiabatic one, , and ii) diffusive overshooting where the temperature gradient is the radiative one, . This difference results in a fully chemically and thermally mixed extended region in the case of penetration, effectively meaning the core is extended and thus more massive. In the case of diffusive overshooting, the extended region is only partially chemically mixed, and any increase in core mass is due to the transport of chemicals into the convective core via this chemical mixing. In either case, the convective core will thus have more hydrogen available to burn (or He in the He-core burning phase), thus extending the main-sequence (MS) lifetime of the star, having a pronounced effect on the morphology of evolutionary tracks. Alternatively, some studies have used near-core rotational mixing to enhance the core mass, effectively producing the same situation where more rotational mixing corresponds to a more massive core. We point out that in 1D diffusive codes, the implementations of overshooting and rotational mixing are seemingly different but both are able to function as a proxy for the total amount of near-core mixing, whatever the physical cause, and that both prescriptions contain un-calibrated parameters. We adopt the approach of using overshooting as a general proxy for the total amount of near-core mixing, whatever its physical cause (rotation, convection, magnetism, waves). The mass-discrepancy reported between either spectroscopic (herrero1992) or dynamical masses (guinan2000; ribas2000b; claret2007; tkachenko2014) and evolutionary masses has traditionally been resolved by increasing the amount of overshooting in a stellar model. This increase in overshooting effectively increases the core mass at a given age, mimicking a more massive star.
It was theoretically outlined that the extent of an overshooting region would be limited by the total energy (mass) of the core (roxburgh1992), and hence the mass of the star. This theoretical prediction has been investigated by numerous studies, some claiming no significant mass dependence (schroder1997; pols1997; stancliffe2015; constantino2018), while others claim a statistically significant mass dependence (claret2007; claret2016; claret2017; claret2018; claret2019). Yet another body of work suggests caution at the ability to constrain overshooting from classical observable quantities given the sensitivity of the data and methodologies (valle2016; higl2017; valle2017; valle2018; johnston2019; constantino2018). On the theoretical side, valle2016 studied the ability for models to uniquely describe a set of observables, revealing an inability to uniquely constrain overshooting. This result was supported by the findings of valle2018 and constantino2018 who show that traditional observed quantities do not provide enough discriminating power to uniquely constrain overshooting, with constantino2018 being unable to reproduce the mass-dependence of overshooting reported by claret2016. Furthermore, johnston2019 showed that even with the inclusion of asteroseismic information, the extent of overshooting, stellar mass, and age cannot be uniquely constrained when properly accounting for correlated nature of stellar model parameters. Instead, johnston2019 suggest that the mass and radius of the convective core should be reported and considered in place of the overshooting.
In this paper, we follow the paradigm of johnston2019 to investigate the ability of well detached double-lined eclipsing binaries (EBs) to probe the core mass. Additionally, we investigate the comparatively sparsely sampled mass range of 4-6 M, connecting the mass ranges intensively covered by claret2016; claret2017 and pols1997; higl2017. We revisit the intermediate- to high-mass double-lined EBs CW Cep and U Oph with new spectroscopy and radial velocities to obtain updated mass and radii estimates. In Section 2, we will provide an overview of both systems, including past modelling efforts. In Sections 3 and 4 we discuss the new spectroscopy, the newly determined orbital elements from spectral disentangling, and the determination of spectroscopic quantities from the disentangled spectra, respectively. Section 5 details the modelling procedure and results for both systems with the mass ratio fixed as derived in the previous section. Section LABEL:section:evol_models covers our evolutionary modelling procedure. Sections LABEL:section:discussion and LABEL:section:conclusions discuss the newly determined mass and radii estimates for each system, the modelling results, and places them in the context of the larger modelling efforts of the community. Following the results of constantino2018 and johnston2019, we report and discuss the estimated core mass and overshooting from our modelling procedure.
2 Literature overview on CW Cep and U Oph
2.1 CW Cephei
The detached double-lined EB CW Cephei (HD 218066, mag) is an intensively studied system. The component masses have reported values ranging from and (popper1974; popper1980; clausen_gimenez_1991; han2002), placing this system at the lower end of the high-mass sequence. This spread in masses results in an uncertainty of % compared to the median value (solution by han2002). The quality of the photometric light-curve solution, in particular the determination of the masses and radii, has been restricted by uncertainty in the mass and light ratio, respectively. This problem has been extensively discussed by clausen_gimenez_1991,who found that the spread in ratio of radii is also accompanied by a significant spread in the sum of the radii. Subsequent analysis of their own new photometry by han2002 and erdem2004 did not settle issue as they used a different methodology from clausen_gimenez_1991. Namely, han2002 did not prefer the photographically determined light ratio over that returned from the lightcurve modelling, and erdem2004 allowed for the possibility of a-synchronous rotation in the components, which alters the light ratios derived from photometric modelling. Comparing the radii derived by different previous analyses (a complete set of the references are given in han2002), a spread of % is found.
Apsidal motion was detected in CW Cep by nha1975 with improvements to the apsidal period made by han2002, erdem2004, and wolf2006. The last authors settled the apsidal period to yr, with an eccentricity of . The relatively short apsidal period, coupled with the brightness of the system have made it an ideal target for dynamical and evolutionary studies. Currently, the nature of the mechanism that drives the apsidal motion is not well understood. New space-based, high-precision, high duty-cycle observations from the NASA TESS mission (ricker2015) promise to provide hitherto unseen constraints on the apsidal motion observed in this system.
Due to a distinct lack of constraints on the metallicity of CW Cep, the unique determination of evolutionary models for CW Cep has proven difficult (clausen_gimenez_1991). To date, several age estimates for CW Cep exist, with clausen_gimenez_1991 reporting an age of , placing both components in the first half of the main-sequence. In their fitting work, ribas2000a derived a much younger system with . However, ribas2000a varied the metallicity and helium content, which introduces a near perfect degeneracy with age and mass in evolutionary modelling. Thus, their solution with , , and a high helium contents albeit with a large uncertainty, , is entirely consistent with that of clausen_gimenez_1991 given this degeneracy. Recently, CW Cep has been modelled by schneider2014, who used a Bayesian modelling framework wrapped around Bonn evolutionary tracks to derive an age of , with best fit for an initial rotational velocity km s for both components. It should be noted that schneider2014 used a less massive solution in their modelling than ribas2000a by for the primary and for the secondary and fixed the metallicity of their tracks to solar. Furthermore, blaauw1959 identified CW Cep to be a member of the Cep OB3, one of the smaller associations in the Orion arm. blaauw1961 also indicated that this association is composed of two subgroups. CW Cep is located in the older subgroup for which clausen_gimenez_1991 found an average age of about 10 Myr in perfect agreement with the age they obtained for CW Cep. However, in a comprehensive study of a new homogeneous photometry and membership jordi1996 obtained ages of 5.5 and 7.5 Myr for the two subgroups, in disagreement with the ages derived by both clausen_gimenez_1991 and schneider2014.
CW Cep is also characterised as an intrinsically variable polarized object (elias2008). Both CW Cep and another early B+B binary system AH Cep, were observed with the Chandra X-ray Telescope in search for evidence of a wind-wind collision (ignace2017). Although CW Cep and AH Cep are comprised of stars with similar properties (c.f. pavlovski2018), X-rays were only detected for AH Cep, despite it being nearly a factor 2 further away than CW Cep. The authors could not disentangle, however, whether the X-rays detected from AH Cep were caused by colliding winds, or perhaps from magnetic activity originating in one of the other components of the quadruple system of AH Cep (ignace2017).
2.2 U Ophiuchi
U Oph (HD 156247, mag) is a detached double-lined EB comprised of two B5V components. Much like in the case of CW Cep, the dynamical solution of U Oph suffers from uncertainties in the light and mass ratios from spectral analysis. The reported masses for U Oph span from and , whereas the reported radii span from and (holmgren1991; vaz2007; wolf2006; budding2009). This represents an uncertainty of and in and and an uncertainty of and in and when compared to the most recent solution by budding2009. Additionally, a wide range of effective temperatures has been reported for both components, with differences up to 3 000 K (clements1979; eaton1973; holmgren1991; andersen1990; budding2009).
A majority of the past lightcurve solutions are based on either the unfiltered photoelectric measurements of huffer1951, the OAO-2 spacecraft photometery of eaton1973, or both. However, the work of vaz2007 and budding2009 relies on new photometric and spectroscopic data. While all of these analyses used different modelling methodologies, codes, and assumptions, they produce derived quantities within a rather small range, as discussed above, and with high precision, which is promising. U Oph displays a very rapid apsidal motion with a period of yr attributed to a distant third body (koch1977; kaemper1986; wolf2002). The apsidal motion has been studied thorougly with several proposed apsidal periods, some as large as 55 yr (frieboes1973; panchatsaram1981; wolf2002; vaz2007). Several recent studies have tried to constrain the nature of the tertiary component, reporting an orbital period of and (kaemper1986; wolf2002; vaz2007; budding2009).
Largely due to uncertainties in its metallicity, there have been several discrepant ages reported for U Oph. holmgren1991 first reported an age of Myr for U Oph when compared to evolutionary tracks without overshooting, and an age Myr when compared to evolutionary tracks with overshooting. Later, vaz2007 compare their solution to evolutionary tracks of different metallicities, considering the apsidal constant as an additional constraint in their modelling and find the best agreement with isochrones for Myr, Myr, and Myr calculated at , and , respectively. budding2009 perform their own evolutionary analysis, again with different codes and solutions compared to the previous evolutionary modelling attempts, and arrive at an average age estimate of Myr between the two components for tracks calculated at . The authors also note that a younger solution is found at Myr from tracks calculated at . Most recently, schneider2014 model U Oph with the Bonnsai code, assuming rotational mixing in their models (calculated at Z=0.014) and find an average age of Myr for the system. budding2009 provide a comprehensive and detailed discussion of U Oph, to which we refer the reader for additional information.
3 Orbital elements from new high-resolution spectroscopy
For both CW Cep and U Oph, we obtained a new series of high-resolution échelle spectra using the High Efficiency and high Resolution Mercator Échelle Spectrograph hermes on the 1.2 m Mercator telescope at the Observatorio del Roque de los Muchachos, La Palma, Canary Islands, Spain. The hermes spectrograph covers the entire optical and NIR wavelength range (3700 - 9100 Å) with a spectral resolution of (raskin2011). CW Cep was observed a total of 18 times over 13 nights. Three observations were taken in January 2015 with the remaining 15 taken in August 2016. The argument of periastron progressed between these two subsets, and less than one degree within either subset. U Oph, was observed 11 times over 10 nights from April to August 2016, during which time the argument of periastron progressed . The resulting spectra have an average S/N of 110 in a range 51-144 and 145 in a range 117-163 for CW Cep and U Oph, respectively.
The basic reduction of the spectra was performed with the hermes pipeline software package. This pipeline delivers merged, un-normalised spectra. Therefore, before disentangling the spectra, we performed normalisation via spline function.
|Param.||Unit||CW Cep||U Oph|
Spectral disentangling (hereafter spd) models the Doppler shift of spectral lines from a time-series of double-lined stellar spectra to determine the spectroscopic orbital elements as well as simultaneously reconstruct the individual spectra of the components (simon1994). Since the orbital elements are directly optimized in spd the determination of radial velocities for each individual exposure is side-stepped. This removes the dependence on template spectra as are commonly used in the cross-correlation function (CCF) radial velocity (RV) determination method, which is often a source of systematic error due to mismatches between the spectral type of the star and that of the template (hensberge2007). Moreover, the resulting disentangled spectra of each component have an increased signal-to-noise compared to single-shot spectra, since disentangling acts as co-addition of the input spectra (c.f. pavlovski2010). To perform spd, we employ the FDBinary code (ilijic2004), which performs spd in Fourier space in order to efficiently solve the large and over-determined system of linear equations represented by the data, through applying discrete Fourier transforms to the spectra (hadrava1995).
FDBinary calculates the RV curve for each component using the standard set of orbital elements: period , time of periastron passage , eccentricity , the argument of periastron , and the semi-amplitudes of the RVs variations for the components , and . FDBinary simultaneously optimises all orbital parameters across the entire set of spectra utilising the simplex algorithm. Although the Balmer lines dominate the optical spectra of hot stars, these lines are broad and usually cover a majority of a single échelle order, thus, any imperfections in the order-merging and normalisation procedure would propagate into the optimisation and affect both the orbital elements and the resulting disentangled component spectra. Therefore helium and metal lines are more suitable for our purposes. The resulting optimised orbital parameters for CW Cep and U Oph are listed in Table 1.
The orbital parameters of CW Cep and U Oph have been derived from fitting RVs in numerous previous studies. For CW Cep, strickland1992 determined km s, and km s () by fitting RVs extracted from three days of IUE spectra using a CCF method. However, since they had a relatively small number of spectra (21), the authors chose to fix the eccentricity to following clausen_gimenez_1991. Alternatively, popper1991 obtained km s and km s () by fitting RVs obtained via the CCF method from digitised plate spectra of CW Cep obtained with the Lick Observatory 3m telescope. The value for derived by popper1991 is substantially larger than that obtained by popper1974 who used the very same data, but employed the oscilloscopic method to determine RVs as opposed to the CCF method that was used by popper1991. By comparison, our results for CW Cep, listed in Table 1, place our estimates within of the solution presented by strickland1992 and within of popper1991.
popper1991 also re-fit the orbital parameters of U Oph on RVs determined via CCF from Lick Observatory plate spectra, reporting km s, and km s (). Additionally, holmgren1991 reported km s and km s () from 31 RV measurements extracted via CCF from 31 Reticon spectra obtained at the DAO 1.2m telescope. Later, vaz2007 obtained slightly different estimates of km s and km s () from 34 plate spectra obtained by the ESO 1.5m telescope. Until this work, the only results based on échelle spectra were presented by budding2009 who found km s and km s () from 30 RV measurements determined via CCF from spectra obtained with the hercules spectrograph attached to the 1m Canterbury University McLellan Telescope located at Mt. John University Observatory in New Zealand. Our results are in rough agreement with the literature values, but highlight the increased precision provided by spd which inherently minimises uncertainties presented by line-blending and template-mismatches that plague CCF techniques.
|[K]||[dex]||[km s]||[km s]|
|CW Cep A||28 300460||4.079||2.00.5||105.22.1|
|CW Cep B||27 550420||4.102||1.50.5||96.21.9|
|U Oph A||16 580180||4.073||2.0||1106|
|U Oph B||15 250100||4.131||2.0||1086|
|CW Cep A||8.30 0.07||7.79 0.08||8.71 0.07||-0.51 0.11||-0.92 0.11||7.55 0.08||7.49 0.06|
|CW Cep B||8.24 0.07||7.70 0.08||8.70 0.06||-0.54 0.11||-1.00 0.10||7.53 0.09||7.45 0.07|
|OB binaries||8.26 0.05||7.70 0.04||8.71 0.04||-0.56 0.06||-1.01 0.06||7.59 0.08||7.57 0.10|
|B single stars||8.33 0.04||7.79 0.04||8.76 0.05||-0.54 0.06||-0.97 0.06||7.56 0.05||7.50 0.05|
4 Atmospheric parameters from disentangled spectra
4.1 CW Cep
CW Cep consists of two early-B spectral type stars with K (popper1974; popper1980; clausen_gimenez_1991; han2002). These temperature estimates place both components in the temperature range where the strength of lines starts to grow, thus allowing us to obtain precise effective temperature estimates through fine tuning the helium ionisation balance. As such, we apply the same methodology as described in pavlovski2018, which we briefly summarise here.
As an observed spectrum of a binary is a composite of spectra of the two components, the disentangled spectra are equal to the intrinsic components’ spectra multiplied by the respective light factors, i.e. the components’ fractional light contribution to the total light of a binary system, such that their co-addition reaches unity in the continuum. Generally, the fractional light contribution of each component can be determined either in the light curve analysis, or extracted from disentangled spectra. In the case of partially eclipsing binary systems where the components have similar radii, the light ratios are degenerate with the radii ratio and inclination. Therefore, it is advantageous to use the light ratio derived from disentangled spectra in the lightcurve modelling. We follow an iterative approach, where we first vary both the light factors and surface gravities, and then impose the light factors derived from spectroscopy as priors in our lightcurve modelling. To obtain atmospheric parameters, an optimised fit to the disentangled spectra of each component, which are re-normalised by their light-factor, is performed over a grid of pre-calculated non-local thermodynamic equilibrium (NLTE) models using the starfit code (tamajo2011; kolbas2014). These theoretical NLTE spectra were calculated using Atlas9 model atmospheres and the NLTE spectral synthesis suite detail/surface (giddings1981; butler1984). The synthetic spectra grid used in the optimisation contains models with , and dex, and solar metallicity [M/H] = 0. However, we are able to fix the for each component according to the values listed in Table LABEL:tab:derived_pars, since high precision, independent estimates for the surface gravities were derived from the light curve modelling. Fixing the surface gravity effectively lifts the degeneracy between the effective temperature and surface gravity, and enables us to use the Balmer lines as constraining information in our fit found by the helium ionisation balance. By fixing the surface gravity and micro-turbulence per component, we reduce the optimisation to eight free parameters: the effective temperature per component, projected rotational velocity per component, a relative Doppler shift between disentangled spectra, and laboratory reference frame, as well as the light-factors of the disentangled components. The optimisation across this parameter space is performed via a genetic algorithm modelled after that of the PIKAIA subroutine by charbonneau1995, with the errors calculated via Markov Chain Monte Carlo technique as described by ivezic2014 and implemented by kolbas2014. The optimisation was carried out over the spectral segment from 4000-4700 Å, and includes the Balmer lines H and H, in addition to helium lines from both ionisation stages. Other spectral lines were masked. Due to the strong interstellar absorption band which effects the red wing of the line, we were unable to use this spectral segment which covers the filter. However, since the effective temperatures of CW Cep A & B are similar, the wavelength dependence of the light-ratio is very small. The final analysis with fixed surface gravities and variable light ratios returned K, and K with light-factors of and , for the primary and secondary, respectively. We note that these light factors are the same as those determined from the initial iteration, within errors The full optimised parameters are listed in Table 2. The best fit for the and lines for both components is shown in Fig. 1.
The reported values for the effective temperature of the primary of CW Cep have a broad range of almost 3 000 K, from K in clausen_gimenez_1991, to a lower extreme in terrell1991, and with intermediate values K in han2002 (terrell1991 and han2002 fix and do not report formal uncertainties for these values). It should be noted, however, that terrell1991 adopt their value for the primary effective temperature from a spectral type classification of B0.5, and clausen_gimenez_1991 determine a mean value from different color calibrated photometric relations. In these studies, the effective temperature of the secondary was then determined from the light curve solution. The reported spread in secondary effective temperature is only 1 300 K, with the hottest solution being only 600 K (clausen_gimenez_1991) cooler than the primary and the coolest solution being 1 900 K cooler than the primary (terrell1991). If we compare the values of our spectroscopically determined effective temperatures for the components of CW Cep, and the difference of their optimal values, 460 K and K, to the various estimates in previous analyses, we find the closest agreement with the estimates of clausen_gimenez_1991. Comparatively, we are able to reduce the uncertainties considerably due to our methodology combining the spectral disentangling, ionisation balancing, and fixing the surface gravity.
Following our atmospheric analysis, we determine a detailed photospheric composition for both stars as well. We calculate ATLAS9 model atmospheres for the atmospheric parameters derived above, from which theoretical spectra are calculated with the detail/surface suite. Details on the model atoms used can be found in pavlovski2018. The abundances are then varied and optimised against the disentangled spectra, from which we report abundances for carbon, nitrogen, oxygen, magnesium, and silicon, as listed in Table 3. Additionally, we are able to derive the microturbulence velocity from the oxygen-lines and the condition of nul-dependence of the oxygen abundance on equivalent width. The derived values for CW Cep A and B are listed in Table 2. For comparison, the ’present-day cosmic standard’ abundance pattern for sharp-lined early-B type stars of nieva2012 is provided in the bottom row of Table 3, with which we find general agreement. We also note that the abundances of CW Cep are in close agreement with the abundance pattern and ratios derived for OB binaries by pavlovski2018 as listed in the third row of Table 3.
Since iron lines are not visible in early-B type stars, the iron abundance can not be directly measured and used as a proxy for stellar metallicity. Instead, lyubimkov2005 determined the magnesium abundance from the Mg ii line in a sample of 52 un-evolved early- to mid-B type stars and used this as a proxy for stellar metallicity. lyubimkov2005 determined the mean abundance to be in close agreement with the solar magnesium abundance, as determined in asplund+2009. Exploiting the Mg abundance as a proxy for metallicity, lyubimkov2005 find that the metallicty of young MS B-type stars in the solar neighbourhood and the Sun are the same. Following this work, we find that our reported magnesium abundance suggests that CW Cep has solar metallicity.
Additionally, we note that we observe H to be in emission in the new spectra assembled for this work. Fig. 2 displays spectra at roughly quarter phases as labelled, all of which show clear double-peaked emission with central absorption. The corresponding velocity difference between the blue () and red () peaks remains constant at km s through the orbital phase. Similarly, we find that the intensity ration between the two peaks remains roughly stable at throughout the orbit as well. For comparison, in Fig. 2 we show also synthetic H profile for 0.25 phase. Although H emission is typical for Be stars or mass-transfer binaries we cannot reliably attribute the emission to a given component. Moreover, as there is no clear evidence of variability of the emission with the orbital phase, we postulate that the emission originates from a circumbinary envelope or the nebula of the Cep OB3 association in which CW Cep is located. An extensive H i nebula in which the Cep OB3 association is embeded is well documented (e.g. simonson1976).
4.2 U Oph
U Oph consists of two main-sequence components of spectral type (mid-)B. Given that our NLTE grid discussed in the previous section is limited to stars hotter than 15 000 K, and that the use of the LTE formalism is overall justified for un-evolved stars in this temperature range, we employ the Grid Search in Stellar Parameters (GSSP, tkachenko2015) code for the analysis of the disentangled spectra of the U Oph’s stellar components. The GSSP algorithm is based on a grid search in basic atmospheric parameters (, , , , and [M/H]) and, if necessary, individual atmospheric abundances, and utilizes a merit function and statistics to judge the goodness-of-fit between the grid of synthetic spectra and the observed spectrum and to compute 1 confidence intervals. Synthetic spectra are computed by means of the SynthV radiative transfer code (tsymbal1996) based on the pre-computed grid of LLmodels atmosphere models (shulyak2004). Both atmosphere models and synthetic spectra can be computed for arbitrary chemical compositions, where one, several, or all chemical elements’ abundances can be set, also with an option of a vertical stratification in the stellar atmosphere. Similarly, the effect of the microturbulent velocity can be taken into account, if necessary assuming its vertical stratification.
GSSP is a multi-function software for spectrum analysis that is able to deal with spectra of single stars (GSSP_single module), and those of spectroscopic double-lined binaries, either with their observed composite spectra (GSSP_composite module) or with the disentangled spectra of individual stellar components (GSSP_single or GSSP_binary module). In the former of the two binary cases (GSSP_composite module), a composite spectrum of a binary is fitted with a grid of composite synthetic spectra that are built from all possible combinations of grid points for the primary and secondary star. Individual radial velocities can also be optimised along with all the aforementioned atmospheric parameters of the two stars, where individual flux contributions are taken into account by means of the stellar radii ratio factor. In the latter case, the distinction is made whether the spectra are analysed as those of a single star with a certain light dilution factor (GSSP_single module, so-called unconstrained fitting where the light dilution factor is assumed to be independent of wavelength) or they are fitted simultaneously by optimising radii ratio to account for individual light contributions (GSSP_binary module, so-called constrained fitting with wavelength dependence of individual light contributions taken into account). A simultaneous fit of the two disentangled spectra is essential when a binary consists of two stars which are significantly different from each other in terms of their atmospheric properties. In this instance, their relative light contributions will strongly depend on wavelength. In the instance where the two stars have similar atmospheric parameters, independent fitting of the disentangled spectra is justified, while still enforcing that the two (wavelength-independent) light factors ultimately add-up to unity (see tkachenko2015, for detailed discussion).
As with CW Cep, the atmospheric parameters of U Oph A & B are similar enough that we fit the disentangled spectra individually. Again, we use an iterative approach where the light-ratios are first determined from the disentangled spectra, then used as priors in the light curve solution. The photometric surface gravities are then fixed and the light-ratios are re-optimized along with the other atmospheric parameters from the disentangled spectra. We found the light factors to be 0.5750.007 and 0.4250.008. The final solution is presented in Table 2, while the quality of the fit is demonstrated in Fig. 3.
5 Revised Photometric Models
Both CW Cep and U Oph have been studied extensively in the literature for several decades, with a heavy focus on the rapid apsidal motion displayed by the systems (holmgren1991; clausen_gimenez_1991; han2002; wolf2002; erdem2004; wolf2006; vaz2007; budding2009). This study aims to use the updated mass ratio, semi-major axis and effective temperatures of the primary and secondary obtained in Sections 3 and 4 to determine updated dynamical masses, radii, and surface gravities from photometric modelling with PHOEBE (prsa2005; prsa2011).
5.1 Photometric Data
For CW Cep we revisit the photometry initially analysed by clausen_gimenez_1991. These data consists of 21 nights of observations spanning 3 years in the Stromgren photometric system, totalling 1396 measurements in the filters, and 1318 in the filter. Both HD 218342 and HD 217035 served as photometric comparison stars, from which the final differential magnitudes were obtained. Extinction corrections were applied to the data as were determined by nightly coefficients determined across the listed comparison stars and other standard objects (gimenez1990). According to clausen_gimenez_1991, the observations were constant to 0.004 mag in all filters, which we adopt as the uncertainty on each point.
We also revisit archival data for U Oph, initially analysed by vaz2007. These data consist of 25 nights of observations spanning 1992-1994 in the Stromgren photometric system, totalling 645 measurements, however, due to a trend in the data, we do not use the -band lightcurve. The data were taken with the 0.5m ESO SAT telescope in La Sille, Chile. HR 6367, HR 6353, and SAO 122251 were all used as comparison stars, from which the final differential magnitudes were obtained. As with CW Cep, extinction corrections were calculated each night from the comparison stars used. For more information on the comparison targets and observations, we refer to vaz2007. Finally, vaz2007 report a standard deviation of 0.0037 mag in the filters, which we adopt as the uncertainty on each point.
5.2 Photometric modelling methodology
Both CW Cep and U Oph are well detached systems, exhibiting mild out of eclipse variability and slow apsidal motion, which for the purposes of our modelling is effectively mitigated by phase-binning the data. Our photometric modelling uses the PHOEBE binary modelling code, which is a modern extension of the original WD code but also incorporates new physics such as dynamic effects, the light travel time effect, and the reflection effect (prsa2005; prsa2011). Given that all components considered are expected to have radiative envelopes, we fix the gravity darkening exponent to unity for all components (vonZeipel1924). To obtain statistically robust estimates for the fit parameters, we wrap PHOEBE into a Bayesian Markov Chain Monte Carlo (MCMC) framework using the emcee affine-invariant ensemble sampler MCMC code (foremanMackey2013), which has already been successfully applied by schmid2015; hambleton2016; pablo2017; johnston2017; kochukhov2018.
MCMC procedures numerically evaluate Bayes’ Theorem, given by:
to estimate the posterior probability, , of some varied parameters given the data .We can see above that is proportional to the product of the likelihood function and the prior probability of the parameter vector . We write the likelihood function as:
where, is each individual model point and are the individual uncertainties associated with the data. We have written the log-likelihood function above as this is what is used in practice. To make efficient use of the information obtained via the spectroscopic analysis, we apply Gaussian priors on the light factors (per-cent contribution per component) and estimates per component, as well as the projected binary separation , the mass ratio , the effective temperature of the secondary , and the eccentricity of the orbit. However, since PHOEBE does not directly sample all of these, we calculate the separately for each component and for every considered. By including the spectroscopic light factors and simultaneously fitting all filters, we arrive at a more robust solution than if we were to fit them all individually and mitigate any degeneracies between the temperatures, light-factors, and potentials of each component (clausen_gimenez_1991). Furthermore, inclusion of priors on for each component helps constrain the spin paramters and , which are otherwise largely unconstrained.
We draw parameter estimates and uncertainties as the median and () Highest Posterior Density (HPD) intervals of the marginalised posterior distribution for each sampled parameter. As both systems undergo apsidal motion, we bin each lightcurve such that each phase bin covers 0.0033 phase units, which covers the entire periastron advance in a single binned point for either system. Although PHOEBE accepts and directly, we sample and in our MCMC analysis and solve for and afterwards. To aid in the discussion and provide additional constraints, we also report the relative radii in the bottom panel of Table 5.4.
5.3 PHOEBE Model: CW Cep
To propagate our newly derived spectroscopic and orbital information into updated dynamical masses and radii, we fix the effective temperature of the primary () to the value listed in Table 2. As mentioned above, we apply Gaussian priors on the mass ratio, the eccentricity, the projected binary separation, the per component, and the light-factor per component in the - and -band lightcurves, since these correspond to the spectral range for which we derived the light factors. The light factor for the - and -bands are given a uniform prior. For each sampled , we interpolate limb-darkening coefficients for the square-root law from the provided PHOEBE girds. Finally, given the radiative envelope of hot stars such as CW Cep A & B, we fix the albedo to unity in both components.
CW Cep is known to suffer from third light which scales the apparent eclipse depths across each filter. Accounting for this scaling is non-trivial as there is a degeneracy between inclination and third light levels. However, this degeneracy is crucial to account for when determining the derived masses. We use a uniform prior on the third light contributions per filter. Additionally, we sample the reference date (), period (), the inclination, the total binary separation, the secondary effective temperature, as well as potentials and synchronicity parameters per component ( and , respectively), giving all uniform priors. All sampled values, and their type of prior, are noted in Table 5.4.
The analysis of clausen_gimenez_1991 states that the argument of periastron changes from to over the course of the photometric campaign, which corresponds to the secondary minima shifting phase units. To mitigate this change in periastron, we phase bin our data to 300 points, with each bin covering 0.0033 phase. Thus, the argument of periastron that we sample does not correspond to the value provided in the literature of the zero-point, but rather to the mean periastron during the photometric campaign.
For a consistency check, we compare the luminosities derived from the binary modelling with the luminosity derived from the Gaia parallax for CW Cep: mas (luri2018; lindegren2018). We take following the reported value of from clausen_gimenez_1991 and calculated as the average correction between CW Cep A and B (reed1998). The summed luminosity derived from our binary model yields , while the GAIA derived luminosity yields . Given the large ( uncertainty on the Gaia parallax, we also check the luminosity derived from the Hipparcos parallax ( mas; vanLeeuwen2007) which yields . We find that all of these agree within the uncertainties.
5.4 PHOEBE Model: U Oph
As with CW Cep, we fix the effective temperature of the primary to the value listed in Table 2 and impose Gaussian priors on , , , per component and the light factors per component in the - and -band lightcurves. Although we can safely ignore the small eccentricity and set it to zero to perform spd, we cannot ignore the eccentricity in the lightcurve. As such, we apply a Gaussian prior according to the values taken from vaz2007. Limb-darkening coefficients are interpolated from PHOEBE tables at every model evaluation. The albedos of both components are fixed to unity as both stars are expected to have radiative envelopes.
Since U Oph is also known to suffer from third light, we take the same approach as with CW Cep, using uniform priors for the third light per filter and uniform priors in all other parameters listed in Table 5.4. To mitigate the effects of the apsidal advance, we phase bin into 300 bins, which effectively covers the apparent change in superior / inferior conjunction. Again, this means that the argument of periastron reported is an average over the photometric campaign when the data was collected. The best model according to the median estimates listed in Table 5.4 is shown in Fig. 5. Derived parameters for both CW Cep and U Oph are reported in Table LABEL:tab:derived_pars alongside other solutions from the literature.
As with CW Cep, we compare the total luminosity obtain from binary modelling with the luminosities derived from GAIA ( mas; luri2018; lindegren2018) and Hipparcos ( mas vanLeeuwen2007), assuming as taken from vaz2007 . The luminosity we calculate as does not agree with the GAIA derived luminosity as within , but does agree with the Hipparcos derived luminosity as within .