1. Introduction
This paper forms part of a program addressing relatively neglected close binary systems in the Southern Hemisphere (Idaczyk et al. Reference Idaczyk, Blackford, Budding and Butland2013). A recent example was that of Erdem et al. (Reference Erdem2022) dealing with the system V410 Pup, which has some similarities in its physical properties to NO Pup. Close binary systems provide a recognised source of fundamental data on stellar parameters, notably their masses and radii. The advent of space based facilities, such as the Transiting Exoplanet Survey Satellite (TESS: Ricker et al. Reference Ricker, Oschmann, Clampin, Fazio and MacEwen2014 and Reference Ricker2015) with its wealth of high-precision photometry, together with available high-resolution spectrometry, have transformed our knowledge of stellar astrophysics.
This knowledge becomes interestingly augmented in the case of eclipsing binaries such as NO Pup (HD 71487, HIP 41361, CoD –38
$^{\circ}$
4462, HR 3327), whose orbital motion includes a steady additional rotation of its main elliptical form. This is associated with the departure from simple Keplerian motion mainly due to the perturbation from the sphericity of the component stars. Data on the difference between observed and calculated times of minimum light from these eclipsing binaries, ‘O – C’s, allow tests of the physics underlying this apsidal motion (Wolf et al. Reference Wolf, Zejda and de Villiers2008). Substantiation of this point can be found in the reviews of Giménez (Reference Giménez, Kondo, Sistero and Polidan1992), Claret & Giménez (Reference Claret and Giménez1993), Tohline (Reference Tohline2002), Horch (Reference Horch, Oswalt and Barstow2013), and others. The short-period (
$\sim$
1.2569 d) of NO Pup’s orbit makes for a relatively fast apsidal motion. In fact, that such a short period system should have retained a noticeable eccentricity seems prima facie surprising (cf. the case of
$\zeta$
TrA; Skuljan et al. Reference Skuljan, Ramm and Hearnshaw2004). As well as the classical two-body discussion of apsidal motion (Sterne Reference Sterne1939), the case of NO Pup may present effects associated with the close binary’s wide orbit companions.
The early-type close pair forms the ‘A’ component of the multiple star WDS J08263-3904Footnote
a
(Figure 1). The ‘B’ component, discovered by John Herschel in 1835, is now at a separation of about 8.1 arcsec and position angle 124 deg. The AB double star, at a distance of 172 pc, has apparently closed in by a few arcsecs over the last
$\sim$
200 yr, although without a significant change of position angle. The B component is itself an astrometric binary (Bab), with a relatively short period of 103 yr. This pair has a combined V magnitude of about 7, and the source would have been included in conventional photometry of NO Pup A. A third component (C) apparently has optical closeness only (Veramendi & González Reference Veramendi and González2014), while companion D was identified by Tokovinin (Reference Tokovinin, Chalabaev, Shatsky and Beuzit1999), with a separation of 5.4 arcsec and position angle 265
$^{\circ}$
. The whole group lies to the south of a small asterism, about half a degree south of the Galactic plane.
Figure 1. Schematic of the main four stars of the NO Pup system.
The close pair’s combined brightness is given as V
$\approx 6.49$
(SIMBAD), with
${B-V} \approx -0.06$
, which is slightly bluer than the reported combination of Main Sequence spectral types B8V+A7V would suggest (for the Aa and Ab stars). The visible companion Bab is reported with V mag 7.04 and
${B-V} = 0.07$
(SIMBAD: Wenger et al. Reference Wenger2000).
NO Pup (Aab) was discovered to be variable by Grønbech (Reference Grønbech1976) during uvby photometry of bright southern stars in 1972 (Jørgensen Reference Jørgensen1972). He derived the first light elements, revealing the relatively short orbital period. A diaphragm of 30 arcsec was used in the photometric measurements of Grønbech (Reference Grønbech1976), so the light contribution of the visible companion was always included in the flux measurements. This four-colour photometry was analysed by Giménez et al. (Reference Giménez, Clausen and Jensen1986), who gave the following linear ephemerides for the primary minimum:
\begin{align} {\mathrm{Min}\ I} = \ & {\mathrm{HJD}} \, 2441752.6576 \pm 0.0012 \nonumber \\ & + (1.2569 \pm 0.0027 ) E\end{align}
and for the secondary:
\begin{align} {\mathrm{Min}\ II} = \ &{\mathrm{HJD}} \, 2441753.3780 \pm 0.0004 \nonumber \\ &+ ( 1.2569 \pm 0.0004 ) E\end{align}
Giménez et al. (Reference Giménez, Clausen and Jensen1986) improved our knowledge of the orbit and its behaviour, announcing the short period of apsidal motion
$U = 37.2 \pm 0.2$
yr.
2. Photometry
2.1. WinFitter fits to TESS data
We downloaded the TESS data in sectors 34, 35, 61, and 62 with short cadence (120-s sampling) from the Mikulski Archive for Space Telescopes (MAST). Our adopted procedure for extraction of TESS data has been spelled out by Blackford (2025). While noting the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) (Jenkins et al. Reference Jenkins, Chiozzi and Guzman2016), we concentrated LC analysis on the Simple Aperture Photometry (SAP) fluxes, since the PDCSAP detrending produces artificial side-effects associated with the search for planetary transits.
LC modelling for close binary stars often refers to the numerical integration procedure of Wilson & Devinney (Reference Wilson and Devinney1971) (WD), which represents the distorted component surfaces as equipotentials, according to the classical point-mass formulation attributed to Roche (Reference Roche1873), recalled in Ch. 3 of Kopal (Reference Kopal1959). Both the WD and WiNFitter (WF) methods converge to the same approximation for the surface perturbation when the internal structural constants
$k_j$
are neglected, implying disregard of the effects of tides on tides. The relevant formula, Eqn 1–11 in Ch 2, or Eqn 2-6 in Ch. 3 of Kopal (Reference Kopal1959), is:
\begin{equation} \frac{ \Delta^{\prime} r }{r_0} = q \sum^{4}_{j = 2} r_0^{j+1} (1 + 2k_j) P_j(\lambda) + n r_o^3(1 - \nu^2) \,\,\, , \end{equation}
where r is the local stellar radius expressed as a fraction of the orbital separation of the components with mean value
$r_0$
.
$ P_j(\lambda)$
is the Legendre polynomials and
$n = (1+q)/2$
, where q is the mass ratio.
$\lambda$
and
$\nu$
are the direction cosines for an arbitrary point on the star’s surface with respect to the line of centres (
$\lambda$
) and the spin axis (
$\nu$
). In any case, the difference between the stellar distortions in WD and WF become small at greater relative separations of the two stars (
$r_0 \rightarrow 0.$
). The direction cosine of the angle between the radius vector
$\hat{r}$
and the line of centres is here
$\lambda$
, and
$\nu$
is the direction cosine of the angle between
$\hat{r}$
and the rotation (‘spin’) axis. The coefficients
$k_j$
(in WF) can be taken from suitable stellar models, for example, Inlek et al. (Reference Inlek, Budding and Demircan2017). They are set to zero in WD.
Sector 34 gathered observations of the photometric flux during the period 2021 Jan 14 to 2021 Feb 08; Sector 35 2021 Feb 09 to 2021 Mar 06; Sector 61 2023 Jan 18 to 2023 Feb 12; and Sector 62 2023 Feb 12 to 2023 Mar 10.Footnote
b
Rhodes (Reference Rhodes2023) provides background and a user manual for WinFitter, Banks & Budding (Reference Banks and Budding1990) discussed its optimisation methods that build on chapter 11 of Bevington (Reference Bevington1969). WinFitter numerically inverts the Hessian of the
$\chi^2$
variate in the vicinity of its minimum to derive estimates for resulting parameter uncertainties, that include effects of inter-correlations between these parameters.
Later, we binned the 17 385 individual observations of Sector 34 by phase to produce 1 021 representative points. Similarly, 14 156 observations in Sector 35 were binned to 1 010 points, 17 771 to 1 045 in Sector 61, and 18 025 to 1 045 in Sector 62. Such binning tends to remove photometric structure with periods not synchronised with that of the orbit (see Section 2.3). However, certain systematic, but unmodelled effects remain in the data, working against the idea of one clear and unequivocal LC ‘solution’. Parameter estimates from the best fitting model are given in Table 1, while Figure 2 displays the model fits to each sector’s data.
Table 1. Parameter values for WF models to TESS Sectors 34, 35, 61, and 62, folding each sector’s data by the ephemeris given by Veramendi & González (Reference Veramendi and González2014). The mass ratio q was adopted as 0.47 (Veramendi & González Reference Veramendi and González2014). The linear limb-darkening coefficient for the primary star was set at 0.29, and for the secondary 0.39. The latter parameters depend on the assigned effective temperatures and wavelength. These were set as
$T_1 = 12\,000$
K;
$T_2 = 7\,700$
K;
$\lambda_{\textrm{eff}} = 0.835\,\unicode{x03BC}$
m. The fractional luminosities
$L_i$
are the mean relative fluxes from each star, normalised so that their sum is unity. The radii
$r_i$
are mean radii of the two stars in the close binary system divided by the semi-major axis of the relative orbit. Angles are given in degrees.
$M_0$
is the mean anomaly at phase zero. The phase bin size is then close to 0.4 deg. See Figure 2 for plots of the model fits to the four data sets.
Figure 2. Binned TESS data are plotted for sectors 34, 35, 61, and 62 with optimal WinFitter models. Fluxes for Sector 35 are offset by
$-0.2$
from their actual values, Sector 61 by a further
$-0.2$
, and Sector 62 by an additional
$-0.2$
. The model fluxes are presented as red continuous curves. The TESS data for each sector have been folded by the orbital period and then binned to 3 600 points. Model parameters are given in Table 1.
The main optimal parameters derived from the WF analysis of the TESS light curves are as follows (the angular parameters i,
$M_0$
and
$\omega$
in degrees) :
$L_1 = 0.595\pm 0.013$
,
$L_2 = 0.128\pm 0.008$
,
$L_3 = 0.277\pm 0.019$
,
$r_1 = 0.255\pm 0.001$
,
$r_2 = 0.181\pm 0.004$
,
$i = 78.8 \pm 0.4$
,
$e = 0.136\pm 0.012$
,
$M_0 = 340.5\pm 6.6$
,
$\omega = 119.8 \pm8.5$
.
2.2. WinFitter fits to BVR photometry
Multicolour photometry of NO Puppis was carried out over 6 nights in January and February 2019 from the Congarinni Observatory, NSW, Australia (
$152^{\circ} 52'$
E,
$30^{\circ} 44'$
S, Alt. 20 m). Images were captured with an ATIK
$^{\textrm{TM}}$
One 6.0 CCD camera equipped with Johnson–Cousins BVR filters attached to an 80 mm f6 refractor stopped down to 50 mm aperture. MaxIm DL
$^{\textrm{TM}}$
software was used for image handling, calibration, and aperture photometry.
HD 71932 was the main comparison star. Its magnitude and colours were measured as
$V = 8.274$
,
${B - V}$
= 0.486, and
${V} - R$
= 0.283. The magnitudes and colours of NO Pup just before and after primary eclipse were then determined as
$V = 6.086(6)$
,
$B - V = -0.019(10)$
,
$V - R = -0.002(9)$
. These measures, coupled with the determinations of the relative fluxes of the components from the LC analysis, allow the magnitudes and colours of the components to be derived, and thence the surface temperatures checked.
Table 3 presents the best fit parameter estimates for the three light curves, while Figure 3 shows the model LC fits to the observations. Table 2 presents the magnitudes of the components in the standard BVR system. From this we see that the primary star (NO Pup Aa) has
$B-V$
of
$-0.11$
, in agreement with the assigned B8V spectral type (Section 1). Ab, with
$B-V = 0.30$
, corresponds to an F0 main sequence star (Table 9.2 in Budding & Demircan Reference Budding and Demircan2022). We could expect NO Pup B to be characterised by a mid-A dwarf spectral type. The WDS catalogue gives a slightly brighter combination V magnitude for NO Pup B as 7.23. The separate magnitudes given in the same catalogue would correspond to A5V and A6V spectral types, with a total mass of about 3.8 M
$_{\odot}$
.
Table 2. BVR magnitudes of NO Pup Aab and B.
Table 3. Parameter values for WinFitter models to the BVR photometry. The parameter symbols carry the same meaning as in Table 1. See Figure 3 for plots of the model fits to the three data sets. Angles are in degrees. The eccentricity (
$e = 0.127$
), adopted after checking the results of numerous optimisation estimates, has been used in these fittings (see Section 5).
Table 4. Final parameter values for WD+MC model to the BVR and TESS light curves. r (volume) is the radius of a sphere having the same volume as the tidally distorted star.
$l_3$
is the third light contribution to the total light at phase 0.25.
Figure 3. WinFitter model lightcurves for the ground-based BVR photometry. The V and R light curves are offset by
$-0.1$
and −0.2, respectively, in normalised flux for display purposes. Optimal parameter values are listed in Table 3.
2.3. WD+MC fits to TESS data
We separately downloaded TESS data in sectors 34, 35, 61 and 62 in short cadence (120-s sampling rate) from the MAST,Footnote c as recalled in Subsection 2.1. We preferred to work with the SAP fluxes, as mentioned above. Data with a quality flag of zero were selected and used approximately 17 000, 14 000, 16 500, and 17 000 points to define the LCs for sectors 34, 35, 61, and 62, respectively.
In the light curves of NO Pup, the light level of maximum I remain higher than that of maximum II. The maximum I light level was then selected for reference in all sectors of the TESS data and additional detrending was applied with a low-order polynomial fit in order to normalise the LCs. As an example, the maximum light levels in two consecutive orbital light curves selected from Sector 35 are shown in Figure 4. Apart from the asymmetry in the maximum light levels, NO Pup shows low amplitude pulsations. We have performed a frequency analysis of these oscillations that will be discussed in Section 7.
Figure 4. Maximum light levels in the light curve of NO Pup along two consecutive orbits from TESS Sector 35 data. Asymmetry between the maximum light levels is evident, and it is also seen that NO Pup shows pulsations with very low amplitude.
In the SAP data of all sectors the CROWDSAP parameter, i.e. the ratio of target to total flux in the photometric aperture, is given as 0.62. This indicates that approximately 38% of the observed flux in the SAP aperture does not come from the target star. This is most probably due to other background starlight collected in the rather large pixels of the TESS detector. The relative contribution of third light (
$l_3$
) was therefore taken into account in the LC analysis of NO Pup A (see Table 4).
We used the numerical integration method of Wilson & Devinney (Reference Wilson and Devinney1971) (WD) which models the light curve of a given binary star by taking into account proximity effects, regarding the surfaces of the components as Roche equipotentials. The original program, wd, has been combined with a Monte Carlo (MC) optimisation procedure, as discussed in Zola et al. (Reference Zola2004). Representative values and uncertainties of the adjusted parameters were derived in this way.
WD requires a preliminary value of the primary’s effective temperature, which was deduced as follows: Veramendi & González (Reference Veramendi and González2014) assigned the spectral type of the system as B5V + B9.5V from their spectral analysis;. from that, they adopted that
$T_{1}$
= 13 000 K. In our later Section 4, we set
$T_{1}$
= 13 500 and 13 700 K (Table 10) by applying Kurucz atmospheric modelling to the disentangled spectrum of the primary. However, the colour index of
$B-V = -0.11$
mag from the photometric analysis in Section 2.2 yields the spectral type as B8V for the primary star. Furthermore, we estimated
$T_{1}$
= 12 000 K from the measured equivalent widths of He I 6678 lines in the UCMJO spectra (see Section 3.4). As a result, the values trialled for
$T_{1}$
range from 12 000 to 13 700 K.
To determine which
$T_{1}$
value is the most suitable, we followed the following procedure: TESS Sector 62 light curves were fitted with
$T_1$
values between 11 000 and 14 000 K in steps of 500 K. We then calculated the photometric parallax from each LC solution. As a result, we adopted the value
$T_{1}$
= 13 000 K in our LC analysis. This gives the photometric parallax value (5.747 mas, see Section 5) closest to the Gaia DR3 parallax. Gaia DR3 gives the trigonometric parallax for NO Pup as 5.799 mas (Gaia Collaboration et al. Reference Collaboration2023). If we add to this the 0.015 mas to account for the zero-point offset suggested by Lindegren et al. (Reference Lindegren2021), the trigonometric parallax of NO Pup becomes 5.814 mas. The effective temperature of the secondary (
$T_2$
) was adjustable in the range of 5 000–9 000 K.
The input range of the orbital inclination was set to
$70^\circ \lt i \lt 90^\circ$
, considering the WinFitter LC fittings in Section 2.1. For changes in
$T_0$
and P, the input range for phase shift (
$\Delta \phi$
; which allows the WD code to adjust for a zero-point error in the ephemeris used to compute the phases; the unit is the orbital period) was set to
$-0.01 \lt \Delta \phi \lt0.01$
.
Based on the RV results in Tables 7 and 9 in Section 3, the mass ratio (q) was fixed at 0.473. Consequently, the input range of the surface potentials (
$\Omega_1$
,
$\Omega_2$
) was set to 3.5–5.0. Taking into account the RV and
$O-C$
analyses, the input range for eccentricity (e) was set to 0.09–0.15. The orbital cycle number of NO Pup A corresponding to the start and end times of the TESS data was calculated using the linear ephemeris given in Table 6, and the argument of periastron (
$\omega$
) corresponding to these cycle numbers were calculated from Equation (6). Accordingly, the input range for
$\omega$
was entered as
$100^\circ$
to
$140^\circ$
.
Table 5. Magnitudes and colours of NO Pup Aab and B from WD+MC results. Errors are on the order of 0.02 mag.
Table 6. Best-fit estimates for the apsidal motion elements of NO Pup. See Figure 7 for plots of the model fits to the times of minima. The parameter estimates from Wolf et al. (Reference Wolf, Zejda and de Villiers2008) are given for easy reference.
Table 7. Best-fit WD modelling for the RV curves measured from selected HARPS spectra of NO Pup A.
A range from 0.3 to 0.8 was set for the fractional luminosity of the primary component (
$L_1$
). The input range for the third contribution of light to the total light of the system (
$l_3$
) was set to 0.30–0.50, based on the CROWDSAP parameter discussed above with the recognition that NO Pup is a multiple star.
A quadratic limb-darkening law was assumed; the coefficients were taken from Claret (Reference Claret2017) according to the effective temperatures and the filter used. The bolometric gravity-darkening exponents and albedoes for both components were taken as 1.0, assuming that the components have radiative atmospheres, in accordance with the regular procedure of the WD program.
The values of the adopted WD + MC model for all the sectors’ data are listed in Table 4. A comparison of the sectors’ LCs with the WD + MC fits is provided in Figure 5.
Figure 5. TESS light curves with the WD model fitting. Residuals to the LC model are plotted in the lower figure. The fluxes for sectors 35, 61, and 62 and their residuals are shifted downward to enhance visibility.
2.4. WD+MC fits to BVR photometry
We also used the wd+mc program for simultaneous fittings of our ground-based BVR data. The input ranges of the adjustable parameters (e,
$\omega$
, i,
$T_2$
,
$\Omega_1$
,
$\Omega_2$
,
$L_1$
,
$L_2$
, and
$L_3$
) were entered into the program as with the TESS data (Section 2.3). For the periastron longitude parameter,
$\omega$
, 80 to 110-degree limits were entered, according to the observation epochs, calculated from Equation (6).
The orbital eccentricity, e, was found to drop below the reasonable limit of 0.10 when allowed to be set by the optimiser program. We have associated this with the effects of data irregularities, especially around the secondary minimum. Subsequently e was fixed at
$e = 0.13$
.
The adopted parameters of the WD+MC model for the the BVR LCs are presented in Table 4. However, when these results were compared with the TESS results given in Table 4, it transpired that allowing freely adjustable third light contributions (
$l_3$
) could change their nominal values by a factor of up to 2. In an alternative approach to the BVR fittings, the geometric parameters (e, i,
$\Omega_1$
and
$\Omega_2$
) were fixed at those found in the TESS LC analyses and just
$T_2$
,
$\omega$
, phase shift,
$L_1$
,
$L_2$
, and
$l_3$
were left free.
We have thus referred to the normal approach as producing Model I, and the alternative, with the geometric elements taken as constant, as Model II. The results of Model I and II fittings are given together in Table 4. The third light contribution (
$l_3$
) in Model II is consistent with the
$l_3$
values obtained in other solutions. A comparison of Model II with the WD+MC and BVR LCs is shown in Figure 6.
Figure 6. BVR light curves with the wd+mc model fitting. Residuals to the LC model are plotted in the lower figure.
We assumed that the third light came from the B component of NO Pup and calculated the relative light contributions of each component to the total light of the system (as
$l_1+l_2+l_3$
) in each band at 0.25 phase from the Model II results in Table 4. These
$l_1$
,
$l_2$
, and
$l_3$
values, and their corresponding magnitudes and colours, are given in Table 5. In this calculation, the magnitude
$V = 6.086$
and colour indices
$B-V = -0.019$
and
$V-R = -0.002$
given in Section 2.2 for NO Pup were used. The dereddened colour indices were obtained from the colour excess
$E(B-V) = 0.020$
derived from the SED analysis in Section 5. Spectral types correspond to the
$(B-V)_0$
colour index according to the calibration in Budding & Demircan (Reference Budding and Demircan2022) and in Eker et al. (Reference Eker2018).
When comparing the WinFitter results (Table 2) with the WD+MC results (Table 5) for the magnitudes and colours of the components, there are small discrepancies of a few percent within the error limits in the magnitudes and colours of the Aa and Ab components, while the slightly larger discrepancy in those of the B component, and thus its spectral type shifts slightly to the early A spectral type in the WD+MC estimation. However, since component B is probably a binary star (see Section 1), it does not seem possible to make a definitive estimate for this component.
2.5. Times of minima
We followed the methodology of Zasche et al. (Reference Zasche, Liakos, Niarchos, Wolf, Manimanis and Gazeas2009) to analyse times of minima (ToM)s, extending the analysis of Wolf et al. (Reference Wolf, Zejda and de Villiers2008), who analysed data up to HJD 2452284.261 covering 25 minima. Our analysis is based on 53 ToMs extending from HJD 2441351.7094 to 2460013.952957 with the inclusion of estimates from TESS data (20 ToMs), our BVR photometry (2 ToMs), and 6 ToMs from Kreiner (Reference Kreiner2004).
To find the parameters related to the apsidal motion using the ToMs of NO Pup A, the methods given by Giménez & Garcia-Pelayo (Reference Giménez and Garcia-Pelayo1983) and Giménez & Bastero (Reference Giménez and Bastero1995) were used. Accordingly, each observed ToM of the eclipsing binary star with apsidal motion can be represented by the equation:
\begin{equation} T = T_0 + EP_s + \Delta \tau, \end{equation}
where
$P_s$
is the sidereal period and is related to the anomalistic period
$P_a$
and apsidal motion rate
$\dot{\omega}$
by:
\begin{equation}P_s = P_a \left(1 - \frac{\dot{\omega}}{360}\right),\end{equation}
where
$P_s$
and
$P_a$
is in days,
$\dot{\omega}$
is in degrees per cycle, and 360 is the number of degrees in one cycle.
$\Delta \tau$
in Equation (4) is the term which shows that the primary and secondary ToMs shift periodically in antiphase about the linear ephemeris and as a function of the eccentricity e, argument of periastron
$\omega$
and orbital inclination i of the relative orbit of the binary star (see Giménez & Bastero Reference Giménez and Bastero1995, their Eqs (15–21)).
We performed our O-C analysis under the assumption of apsidal motion, taking the initial values for
$T_0$
and
$P_s$
from Kreiner (Reference Kreiner2004) and the input values for other apsidal motion parameters from Wolf et al. (Reference Wolf, Zejda and de Villiers2008). Our results from the best-fitting model are given in Table 6, while Figure 7 plots the model fit against the ToMs, and also the rescaled residuals. Our findings are in good agreement with those of Wolf et al. (Reference Wolf, Zejda and de Villiers2008), with the apsidal period U also in accord with that of Giménez et al. (Reference Giménez, Clausen and Jensen1986). However, although there are more data in our O-C analysis, the error estimates in the parameters of Wolf et al. (Reference Wolf, Zejda and de Villiers2008) are generally smaller. A possible reason for this, apart from the different accuracies of the data samples used, is that different methods have been followed in the O-C analysis and therefore in the error estimation. We utilised a numerical method based on the incomplete gamma function for error estimation (see Zasche Reference Zasche2008, and his references), while Wolf et al. (Reference Wolf, Zejda and de Villiers2008) use the least squares method to derive formal uncertainty values.
Figure 7. In upper panel we plot the optimal model (shown as a black curve) for the primary minima (black filled circles), along with the optimal model fit to the secondary ToMs (blue curve and unfilled circles). Units are days. Lower panel shows the residuals from the optimal models.
Keeping in mind the variation of the argument of periastron (
$\omega$
) determined from the O-C analysis, as an alternative approach, we plotted the
$\omega$
values found from the solutions of the RV and LCs of the system according to the cycle number in Figure 8. The
$\omega$
values in the lower left of the figure are taken from the analysis of the HARPS RV curves (Table 7), the RV1 solution from the UCMJO spectra and that from the REOSC spectra (Table 9), and the
$\omega$
values in the upper right of the figure are from the BVR and TESS LC fittings (Table 4). In this way, we obtained the linear time dependence of
$\omega$
shown in Figure 8:
\begin{equation} \omega = \omega_0 + \dot{\omega} E,\end{equation}
where E is the cycle number of the mid-times of the related observations. These were calculated from the linear ephemeris given in Table 6. The plot of the
$\omega$
values against cycle number in Figure 8 reveals that
$\dot{\omega}$
=
$0.0340 \pm0.0003$
deg/cycle and
$\omega_0$
=
$-373 \pm4$
deg as the linear best fit parameters. This value of
$\dot{\omega}$
, which was confirmed by the O-C analysis (see Table 6), results in a value of
$U = 36.4 \pm 0.3 $
yrs.
Figure 8. Variation of the argument of periastron
$\omega$
of the eccentric binary NO Pup A.
3 Spectrometry
3.1. HARPS spectra
The European Southern Observatory (ESO) Science Archive Facility (SAF)Footnote
d
contains 43 HARPS (High-Accuracy Radial-velocity Planet Searcher) spectra of NO Pup, taken between December 09, 1996, and June 21, 2015.Footnote
e
Of these spectra, we selected 34 from the nights of April 2–7, 2009. We used the HARPS cross-dispersed échelle spectrograph at the 3.6-m telescope at La Silla Observatory (Mayor et al. Reference Mayor2003). We chose to operate HARPS in the EGGS mode, which has a larger fibre entrance on the sky than the standard HAM mode (1.4 versus 1.0 arcsec), resulting in a higher throughput (by a factor of approximately 1.75), a lower resolving power (80 000 as against 115 000) and a lower but still excellent RV precision (3 m s
$^{-1}$
rather than 1 m s
$^{-1}$
).
Each spectrum consists of 72 orders incident on two CCDs that cover the range 3 780–6 900 Å with a gap at 5 304–5 337 Å between the CCDs. The data were reduced using the standard HARPS pipeline, with the pipeline products retrieved from the ESO SAF.
3.2. FEROS spectra
Upon reviewing source material we found that the ESO SAFFootnote f maintains several datasets from the Fiber-fed Extended Range Optical Spectrograph (FEROS) on the 2.2-m telescope at La Silla (Elkin, Kurtz, & Nitschelm Reference Elkin, Kurtz and Nitschelm2012). This échelle spectrograph has a resolving power of 48 000 and each exposure covers the full optical range (3 600–9 200 Å). The majority of these data were obtained under ESO Program ID 088.D-0080(B) (PI: Hełminiak). However, since the HARPS spectra have a higher spectral resolution and signal-to-noise (S/N) ratio (on average 170), we preferred to use only the HARPS spectra in the present study. We discuss these data next.
3.3. HARPS spectral analysis
Although all 43 HARPS spectra available at ESO SAF were taken into account for RV measurements and atmosphere modelling (see Section 4), the consequences of the rapid apsidal motion caused us to select only RVs from the 34 HARPS spectra during a 5-day interval for the spectroscopic orbit model.
For the RV measurements, we employed the cross-correlation technique using the IRAF FXCOR task (Tody Reference Tody and Crawford1986). As template material, we generated synthetic spectra based on the known spectral types of the binary components (Sections 1 and 2). For each binary component, two different theoretical templates were generated based on their spectral types and the spectral type–effective temperature relationship (Cox Reference Cox2000), using ATLAS9 model atmospheres (Kurucz Reference Kurucz1993) and the synthe code (Kurucz & Avrett Reference Kurucz and Avrett1981) as the two components have significantly different spectral types. These synthetic templates were then broadened according to the spectral type–projected rotational velocity relationship (Cox Reference Cox2000). RV variations were determined and are listed in Table A1. These HARPS RVs were phased according to the linear ephemeris given in Table 6. Additional spectral analysis to estimate atmospheric parameters is provided in Section 4.
However, due to the orbital eccentricity tending to drop below the assigned limit of 0.10, presumably due to insufficient optimal parameter resolution, e was fixed to the weighted average of the values in the previous RV and LC solutions, that is, 0.127. The optimal values are given in Table 7, and a comparison of the model with the observations is shown in Figure 9. The amplitudes (
$K_1$
and
$K_2$
) were calculated accordingly. Similar results were found using the Winfitter program.
Figure 9. RV curves measured from selected HARPS spectra of NO Pup A with the WD model fitting. Residuals to the RV models are plotted in the bottom figure. RVs of the primary and the secondary components are marked as filled and hollow symbols, respectively.
3.4. UCMJO spectra
Relevant spectroscopic observations include those made with the HERCULES spectrograph (Hearnshaw et al. Reference Hearnshaw, Barnes, Frost, Kershaw, Graham, Nankivell, Ikeuchi, Hearnshaw and Hanawa2003), together with the 1m McLellan telescope at the University of Canterbury Mt John Observatory (UCMJO).
Observations were recorded with a 4k
$\, \times \,$
4k Spectral Instruments (SITe) camera (Skuljan et al. Reference Skuljan, Ramm and Hearnshaw2004). Wavelength and relative flux calibration was performed using the latest version of the software package hrsp (Skuljan Reference Skuljan, Kurtz and Pollard2004, Reference Skuljan2021) that outputs measurable data in fits (Wells, Greisen, & Harten Reference Wells, Greisen and Harten1981) formatted files. Typical exposures lasted for
$\sim$
700 s. Further information on the spectroscopic arrangements at UCMJO were given by Bakiş et al. (Reference Bakiş2024).
24 exposures of NO Pup were made during the interval December 3–13, 2009, although, unfortunately, a number of them were affected by technical difficulties. However, 10 spectra were selected from these spectral images for analysis. The individual lines detected in the orders in the observed spectra are listed in Table A2. Since the light contribution of the secondary component is relatively low – on the order of
$\sim$
10 per cent (see Section 2), mostly only the primary features could be distinguished.
Apart from H
$\alpha$
and H
$\beta$
, the best-defined line is probably the primary He I
$\lambda$
6678, where there is no significant contribution from the secondary since its effective temperature is
$\sim$
7 500 K (see Table 4). The He I
$\lambda$
6678 lines in these 10 spectra were studied, and the RVs of the primary determined, using the program Prof (latest version, Erdem et al. Reference Erdem2022).
If the resolution is sufficiently high, Prof models spectral line profiles with a parameter set that determines the RV of the centre of light, as well as the rotation rate of the source and the turbulence scale in the surrounding plasma. The application of Prof to the He I
$\lambda$
6678 features resulted in the values of RV, rotation parameter (r) and equivalent width (EW) listed in Table 8, and an example of such a profile fitting is displayed in Figure 10.
Table 8. Values of RV, rotation parameter (r) and equivalent width (EW) of the primary component of NO Pup A derived from the He I lines in the UCMJO spectra.
Table 9. Best-fit modelling for the RV curves of Veramendi & González (Reference Veramendi and González2014) of NO Pup.
Table 10. The results of atmospheric parameter analyses for the primary and secondary components of No PUP using the KOREL and FDBINARY disentangled spectra. The metallicity values given in the table are the [M/H] and [Fe/H] for the KOREL and FDBINARY analyses, respectively.
*represents the fixed parameters.
Figure 10. Convolved rotation Gaussian fitting to the He I
$\lambda$
6678 line profile in UCMJO spectrum of NO Pup.
The mean value of the rotation parameter (r) for the data sets in Table 8 yields a projected mean equatorial rotation speed of
$82 \pm2$
km s
$^{-1}$
. If we assume that the inclination of the primary star’s rotation axis is equal to that of the system’s orbit (
$i_{rot}=i_{orb}$
) and we neglect effects related to the low orbital eccentricity, the primary exhibits a synchronised rotation with the mean orbital revolution (
$P_{rot}=P_{orb}$
), we find the synchronous projected rotation speed of the primary as
$86 \pm5$
km s
$^{-1}$
using the formula
$v_{rot}\sin i_{rot}$
=
$(2\pi R_{1}\sin i_{rot})/P_{rot}$
and the values of i and
$R_1$
in Table 11. In this way, the projected rotational velocity of the primary component, derived from the He I
$\lambda$
6678 line profile fitting, and the projected synchronous rotational velocity computed from Table 11 agree within the uncertainty limits, supporting that NO Pup Aa rotates synchronously. The possibility of pseudo-synchronous rotation in the NO Pup A system is considered in Section 8.
Table 11. Absolute parameters of the eclipsing binary NO Pup A.
The fitting function in Prof has two main components: uniform rotation and a Gaussian turbulence broadening. Here Prof applied to the He I line in the UCMJO spectra of NO Pup A gives the turbulence parameter of the order of a few km s
$^{-1}$
for the surface of the primary star.
Prof also estimates the equivalent width of a line by numerical integration. The results are given in Table 8. The mean value of the measured equivalent widths (EWs) for the data sets in Table 8 is
$0.07 \pm0.01$
Å. The relative noise in the He I lines in the observed spectra makes for a fairly uncertain mean value of the EWs. However, in comparison with the calibration data of Leone & Lanzafame (Reference Leone and Lanzafame1998), the effective temperature of the primary was estimated to be
$12\,000 \pm1\,000$
K, which corresponds to a spectral type B8/9.
The UCMJO observations, adopted as of suitable quality, provide coverage for the first half of the full radial velocity (RV) cycle, consistent with the spectroscopic data of Veramendi & González (Reference Veramendi and González2014) (see Figure 11). Their observations were made with the 2.15 m telescope and the REOSC échelle spectrograph at the Complejo Astronómico El Leoncito (CASLEO) during 10 allocations between 2008 and 2013. RVs were determined by the cross-correlation technique, with spectral disentangling of double-lined systems.
The spectroscopic results, as a whole, have confirmed the type classifications for NO Pup and obtained good orbital coverage, with twenty individual points having signal/noise ratios of greater than 100. Veramendi & González (Reference Veramendi and González2014) went on to provide absolute parameters of the system, making use of Grønbech (Reference Grønbech1976)’s uvby photometry and the well-known Russell paradigm.
Figure 11. RV curves of NO Pup A with the WD model fitting. Black circles denote REOSC RV data of Veramendi & González (Reference Veramendi and González2014), whereas orange circles denote the RVs of the primary component derived from He I lines in the MJUCO spectra. Residuals to the RV models are plotted in the bottom figure. RVs of the primary and the secondary components are marked as filled and hollow symbols, respectively.
We remodelled the RV data set of Veramendi & González (Reference Veramendi and González2014) using the program suites Winfitter (Rhodes Reference Rhodes2023) and WD and presented our results in Table 9. Although the main parameters of the optimal RV curve fittings are essentially similar to those of Veramendi & González (Reference Veramendi and González2014), the values of the longitude of the periastron (
$\omega$
) appear different. Noting the decline in precision of the determination of
$\omega$
at low e, the estimates of Winfitter and WD are tolerably in agreement; however, the
$\omega$
value of Veramendi & González (Reference Veramendi and González2014) is quite closer to zero. This may reflect a different reference epoch for the data of Veramendi & González (Reference Veramendi and González2014), given the fast rate of apsidal motion (
$\sim$
10 deg y
$^{ -1}$
).
4. Determination of the atmospheric parameters
The atmospheric parameters of a star, in particular the effective temperature (
$T_{\textrm{eff}}$
) and surface gravity (
$\log g$
), are crucial for understanding the nature of the source. This is especially true for binary stars, where accurately estimating the
$T_{\textrm{eff}}$
values of both components is essential for performing reliable binary modelling analysis.
There are two methods that can be applied in advance to determine the atmospheric parameters of binary systems. The first involves creating composite synthetic spectra and comparing them with the observed binary spectra. The second method is spectral disentangling, which allows us to separate the spectra of both binary components. For both methods, the flux ratios of the binary components must be known. However, for spectral disentangling, reliable results require spectra spread over the orbital phase. Since our data have such a distribution, we preferred applying the spectral disentangling. This is explained more fully in the next section.
4.1. Disentangling components’ spectra and atmosphere modelling
We have selected two spectral regions in the HARPS data to disentangle the component spectra. Both regions are centred around the H
$\beta$
line, as shown in Figure 12. This region was chosen due to the presence of a relatively high number of metal lines for both components. This enables the extraction of their individual spectra and facilitates the modelling of the spectroscopic orbit. The codes korel (Hadrava Reference Hadrava2004) and fdbinary (Ilijic et al. Reference Ilijic, Hensberge, Pavlovski, Freyhammer, Hilditch, Hensberge and Pavlovski2004) were used to cross-check the results.
Figure 12. Spectral regions around H
$\beta$
in the HARPS data used for disentangling.
Spectra obtained during eclipses, except at mid-eclipse, were excluded from the analysis, because non-Keplerian effects, such as the Rossiter–McLaughlin effect, are not accounted for in the codes. In the korel analysis, each spectral region was re-binned to 2 048 bins, resulting in a resolution of approximately 7 km s
$^{-1}$
per pixel. This re-binning smoothed the data without significantly compromising resolution. FPor the fdbinary analysis, the spectra were used directly, with no additional processing applied before the disentangling. During the analysis, some orbital parameters, such as P and
$K_{1,2}$
were held constant. Additionally, the flux ratio of the binary components was examined.
The disentangled spectra determined from both korel and fdbinary programs were modelled using synthetic spectra computed from Kurucz’ atlas9 model atmosphere grids (Kurucz Reference Kurucz1993). In the atmospheric determination with the korel disentangled spectra, the grid parameters and intervals for the primary and secondary components were chosen as given below. During this analysis, the microturbulence parameter,
$\xi$
, was set at 2 km/s for both components.
-
• Primary:
-
- Effective temperature (
$T_\mathrm{eff}$
): 11 000–14 000 K, in steps of 100 K -
- Surface gravity (log g): 4.3–4.4 cgs, in steps of 0.05 cgs
-
- Metallicity ([M/H]):
$-0.5$
to
$+0.5$
, in steps of 0.5 -
- Projected rotational velocity (
$v \sin i$
): 60–100 km s
$^{-1}$
, in steps of 10 km s
$^{-1}$
-
-
• Secondary:
-
-
$T_\mathrm{eff}$
: 7 000–9 000 K, in steps of 100 K -
- Surface gravity (log g): 4.3–4.4 cgs, in steps of 0.05 cgs
-
- Metallicity ([M/H]):
$-0.5$
to
$+0.5$
, in steps of 0.5 -
-
$v \sin i$
: 60–100 km s
$^{-1}$
, in steps of 10 km s
$^{-1}$
-
with fdbinary. the input parameters used in the spectral analysis of the separated spectra were selected a follows:
-
• Primary:
-
- Effective temperature (
$T_\mathrm{eff}$
): 10 000–14 000 K, in steps of 100 K -
- Surface gravity (log g): 3.8–4.4 cgs, in steps of 0.1 cgs
-
- Metallicity ([Fe/H]) :
$-0.5$
to
$+0.5$
, in steps of 0.1 -
-
$v \sin i$
: 30–150 km s
$^{-1}$
, in steps of 1 km s
$^{-1}$
-
-
• Secondary:
-
-
$T_\mathrm{eff}$
: 6 800–8 500 K, in steps of 100 K -
- Surface gravity (log g): 3.8–4.4 cgs, in steps of 0.1 cgs
-
- Metallicity ([Fe/H]):
$-0.5$
to
$+0.5$
, in steps of 0.1 -
-
$v \sin i$
: 30–150 km s
$^{-1}$
, in steps of 1 km s
$^{-1}$
.
-
In this analysis, the spectrum synthesis method (Niemczura & Polubek Reference Niemczura, Polubek, Fletcher and Thompson2006) was used to estimate the atmospheric parameters (
$T_\mathrm{eff}$
, log g,
$\xi$
) and iron (Fe) abundances of both close binary components by taking into account the Kurucz line list.Footnote
g
Based on the Saha–Boltzmann equation, the atmospheric parameters were estimated during the analysis in the way presented by Kahraman Aliçavuş et al. (Reference Kahraman Aliçavuş2016).
The resulting atmospheric parameters, after examination of the KOREL and FDBINARY disentangled spectra, are given in Table 10. As can be seen in the table, the results of both analyses agree within the adopted error limits.
Our results are consistent with those of Veramendi & González (Reference Veramendi and González2014) within the error bars; however, the
$v \sin i$
of the secondary component shows a slight discrepancy (58.7
$\pm$
0.6 km s
$^{-1}$
from Veramendi & González Reference Veramendi and González2014). Given that the secondary component is fainter and the spectral data used by Veramendi & González (Reference Veramendi and González2014) have lower resolving power and signal-to-noise ratio, we consider this difference acceptable. The disentangled spectra of both components, along with the best-fitting synthetic spectra, are displayed in Figures 13 and 14. The derived effective temperatures in Table 10 appear slightly greater than those given in Table 4, but the differences are within reasonable uncertainty estimates of each other.
Figure 13. Disentangled metal lines and best-fitting synthetic spectrum for the primary (top) and secondary (bottom) components, respectively.
Figure 14. Disentangled H
$_\beta$
line and best-fitting synthetic spectrum for the primary (top) and secondary (bottom) components, respectively.
5. Absolute parameters
Applying the well-known rearrangement of Kepler’s third law as:
\begin{equation} (M_1 + M_2) \sin^3 i = CP (1 - e^2)^{3/2} (K_1 + K_2) ^3 ,\end{equation}
where the constant
$C = 1.03615 \times 10^{-7} $
, period P is in days, and the RV amplitudes
$K_1$
and
$K_2$
are in km s
$^{-1}$
. Adopting the weighted average values of eccentricity and orbital inclination from Section 2 as
$e = 0.127$
and
$i = 81.33^\circ$
and taking the RV amplitudes
$K_1$
and
$K_2$
from Table 7, we find the component masses of NO Pup A to be
$M_1 = 3.58 \pm 0.11$
and
$M_2 = 1.68 \pm 0.09$
in solar units. These masses are slightly larger than those found using the RV data of Veramendi & González (Reference Veramendi and González2014) in Table 9, perhaps due to the larger values of
$K_1$
and
$K_2 $
from the HARPS analysis.
With the total mass of the close binary at
$5.26 \pm 0.20$
M
$_{\odot}$
and the period 0.0034411 yr, Kepler’s third law yields the semi-major axis of the orbit as
$0.03956 \pm 0.00023$
AU, or
$8.51 \pm 0.05$
solar radii. From the results given in Section 2 we find
$R_1 = 2.17 \pm 0.03$
and
$R_2 = 1.51 \pm 0.06$
.
The high accuracy of this determination allows us to constrain the temperatures used in the photometric parallax (see Section 2.3 and the following paragraphs). The mean effective temperatures, taking into account the results given in Sections 2–4, are, for the primary,
$\sim \! 13\,300 \pm 500$
; and for the secondary
$\sim \! 7\,400 \pm 500$
K.
The surface gravitational accelerations (
$g_1$
,
$g_2$
) are directly related to solar values through
$g/g_{\odot}= (M/M_{\odot}) / (R/R_{\odot})^2$
. The bolometric magnitudes (
$M_{bol}$
) and luminosities (L) of the component stars are calculated using Pogson’s formula and the absolute radii and effective temperatures listed in Table 11. We thus write
$M_{bol}=M_{bol,\odot}+10\log T_{\odot} - 10\log T - 5\log (R/R_{\odot})$
, and
$L/L_{\odot}=10^{0.4(M_{bol,\odot}-M_{bol})}$
. The solar values, adopted by IAU 2015 Resolutions B2 and B3, were used in our calculations.
Bolometric corrections in the V band for the components (
$BC_{1,2}$
) were taken from Flower (Reference Flower1996), while those in the TESS band were taken from Eker & Bakiş (Reference Eker and Bakiş2023), according to their effective temperatures used in the conversion from bolometric magnitudes to V and TESS-band absolute magnitudes (
$M_{V,1,2}$
and
$M_{TESS,1,2}$
). The absolute magnitude of the eclipsing binary NO Pup A is also computed from following equation:
\begin{equation} M_{band, \text{system}}=M_{band, 2} -2.5 \log \left(1+10^{-0.4\left(M_{band, 1}-M_{band, 2}\right)}\right).\end{equation}
Finally, the distance to the system is calculated from the distance modulus (
$d = 10^{m_{band} - M_{band} + 5 - A_{band}}$
). Here
$m_{band}$
is the apparent magnitude and
$A_{band}$
is the interstellar extinction in the given band.
The
$A_{band}$
extinction is estimated from the SED analysis following the method described in Bakiş & Eker (Reference Bakiş and Eker2022), later refined by Eker & Bakiş (Reference Eker and Bakiş2023) for the TESS pass-band. The best-fitting SED model determines the reddening as
$E(B-V)=0.02\pm 0.01$
mag, corresponding to
$A_V$
=0.062
$\pm$
0.031 mag and
$A_{TESS}$
=0.041
$\pm$
0.019 mag. Figure 15 presents the SED data alongside synthetic spectra computed using the system parameters, demonstrating a strong agreement between the model and observations. The absolute parameters thus obtained are given in Table 11 with their errors.
Figure 15. SED data (black dots) and the combined synthetic spectra of the components, which are calculated using the absolute parameters of the components and the distance of the system given in Table 11.
On the other hand, the following equations given by Budding & Demircan (Reference Budding and Demircan2007) can also be used for the photometric parallax:
\begin{equation}\log \Pi = 7.454 - \log R - 0.2 \text{V} - 2 F'_{\text{V}} \\\end{equation}
where
$F'_{\text{V}}$
(flux scale) is equal to 0.25
$\times$
the logarithm of the surface flux in the V band and is specified by:
\begin{equation}F'_{\text{V}} = \log T_{\rm eff} + 0.1 BC\end{equation}
Here, if we use the equations given above for the NO Pup Aa component, that is,
$V = 6.630$
from Table 5, R = 2.17 R
$_{\odot}$
and
$A_V = 0.062$
from Table 11, and
$T_{\rm eff} = 13\,000$
K and
$BC = -0.89$
, we find the photometric parallax of this component as
$\Pi = 5.747$
mas (i.e. its distance
$d = 174$
pc).
Figure 16. Location of the components of NO Pup A in the H-R diagram. The Geneva evolutionary tracks for 3.58 M
$_{\odot}$
(blue line) and 1.68 M
$_{\odot}$
(red line), corresponding to the primary and secondary stars, are plotted for
$Z_{\odot}=0.020$
. The Geneva isochrone of 20 Myr for Z = 0.020 is also indicated by the dashed black curve. Filled and open circle symbols represent primary and secondary components, respectively. Vertical and horizontal lines show error bars of the measured quantities.
One way to assess the accuracy of the absolute parameters in Table 11 is to compare them with those in the Gaia DR3 catalogue. An absolute magnitude of
$M_{G}=0.225$
mag for the system is calculated from the distance modulus, using the G-band apparent magnitude
$G=6.525$
mag, distance
$d=172$
pc and extinction
$A_{G}=0.122$
mag from the Gaia Archive.Footnote
h
This absolute Gaia magnitude,
$M_{G}=0.225$
mag, is then converted to the bolometric magnitude
$M_{bol}=-0.595$
mag from the formula
$M_{bol}$
=
$M_{G}$
+
$BC_{G}$
, using the bolometric correction
$BC_{G}=-0.820$
mag for
$T=13\,300$
K, derived from Eker & Bakiş (Reference Eker and Bakiş2023) (their Eqn. 4).
The bolometric magnitude and distance values derived from the BVR & TESS-band LCs + HARPS RVs solutions for NO Pup A (
$-0.617 \pm0.426$
mag and
$171 \pm20$
pc) agree with those computed from the Gaia Archive (
$-0.595\pm0.026$
mag and
$172\pm1$
pc) within the error limits. Here, the difficulties in calibrating bolometric corrections, especially for hot stars with
$T\gt12\,000$
K, should be kept in mind.
We used the Geneva evolution models (Yusof et al. Reference Yusof2022) to study the evolutionary status of the eclipsing binary NO Pup A. According to the mass values of the components listed in Table 11, the Geneva evolutionary tracks and isochrone were created using the interpolation interface on the Geneva group website.Footnote i
The H-R diagram shown in Figure 16 was used to estimate the ages of the components. In this diagram, the positions of the two components are plotted on the evolutionary tracks for their measured masses. As a result, the isochrone of log(age) = 7.30 (i.e. 20 My) with supersolar metallicity (Z=0.020) matches the positions of both components within the error limits. We may note that the atmosphere models in Section 4.1 produced metallicities close to solar, while according to the evolution models at 20 Myr, the metallicities of the components would be somewhat larger than that. This result (age and metallicity) was also confirmed using the Padova evolution models (e.g. Nguyen et al. Reference Nguyen2022).
6. Astrometry
Astrometric data for the visual binary WDS J08263-39044 (NO Pup Bab) were taken from the Washington Double Star catalogue (WDS, Mason et al., Reference Mason, Wycoff, Hartkopf, Douglass and Worley2024). We made use of the Markov Chain Monte Carlo (MCMC) optimisation technique as described by Ersteniuk et al. (Reference Ersteniuk, Banks, Budding and Rhodes2024) to fit the apparent orbit. The best fit parameter estimates together with formal uncertainties are given in Table 12, while Figure 17 shows the model orbit based on these parameter values plotted against the fitted data.
Table 12. Best-fit estimates for the orbital parameters of NO Pup B. WDS refers to the orbit parameters listed in the Washington Double Star catalogue and MCMC refers to the best-fitting parameters estimate derived by this study using Hamiltonian Markov Chain Monte Carlo. Angles are in degrees, semi-major axis a in arc-seconds, period in years, and the epoch (time of periastron passage) is in fractional Besselian year.
Figure 17. Model fit to WDS astrometric data for NO Pup B. The red curve plots the model orbit (see Table 11 for the listed parameter values), the black dots show the observational data, and the short blue lines connect the fitted data points with their expected positions along the model orbit. The primary is indicated by the star symbol at the origin. East is to the right, and down is northwards. Labels give the observation dates.
Comparison between the parameters published in the WDS from Josties & Mason (Reference Josties and Mason2019) and those derived from this study are in reasonable agreement, particularly for the period (P in years), semi-major axis (a in arcsecs), and eccentricity (e).
$\omega$
is the argument of periastron.
$\Omega$
gives the position angle of the ascending node relative to the north direction in the sky plane. The Eulerian angles (
$\omega$
, inclination i, and
$\Omega$
) are less well constrained by the MCMC fit.
The WDS catalogue lists the primary component of NO Pup B to be about 1.42 V mag fainter than that of the B8 primary of NO Pup A. The secondary component of NO Pup B is 0.20 V mag fainter than the primary. Using the data in Table 9.2 of Budding & Demircan (Reference Budding and Demircan2022) we can then surmise that component B likely comprises two Main Sequence dwarfs of spectral types A5 and A6 respectively. Their masses, from that same Table 9.2, are around 2 and 1.8 M
$_{\odot}$
. If we correspondingly adopt the mass sum of NO Pup B to be 3.8 M
$_{\odot}$
we require the inclination to be close to 115 deg to allow Kepler’s third law to be satisfied with the orbit’s semi-major axis being near to 34.4 AU. Although this inclination is appreciably different from the MCMC estimate in Table 12, the difference is comparable to that between the WDS and MCMC results. Further monitoring of the pair of stars in NO Pup B is clearly needed for derivation of more reliable orbital parameters (see also the related discussion on model fitting and the use of historic data in Tokovinin Reference Tokovinin2024).
7. Pulsation analysis
The out-of-eclipse ranges of the TESS LCs of NO Pup show systematic variations that we attributed to pulsations. We investigated these using the 120-s cadence data. With close binary variations having a dominant effect, we first separated the binary LC from the full flux. This was achieved by fitting the harmonics of the orbital frequency of NO Pup A using a method similar to the study of Kahraman Aliçavuş et al. (Reference Kahraman Aliçavuş, Çoban, Çelik, Dogan, Ekinci and Aliçavus2023).
The residual light curve was then analysed using the Period04 program, which derives pulsational frequencies based on a discrete Fourier transform algorithm. To estimate the significant frequencies, a 4.5
$\sigma$
significance limit was applied, as outlined in the study of Baran & Koen (Reference Baran and Koen2021). As a result, we obtained pulsation frequencies between 1.26 and 36.95 d
$^{-1}$
. The list of derived frequencies is provided in Table 13, and the amplitude spectrum is shown in Figure 18.
Table 13. Results of the frequency analysis.
Figure 18. Power spectrum of NO Pup. The vertical line represents the 4.5
$\sigma$
level.
Taking into account the estimated frequencies and
$T_{\textrm{eff}}$
values of the binary components, as well as the visual companion, we conclude that there are probably two types of oscillating star present in the source. One is of
$\delta$
Scuti type, and the other exhibits pulsations characteristic of a slowly pulsating B (SPB) star.
$\delta$
Scuti stars, which range from A to F spectral type, show oscillations in the frequency range of approximately 5–80 d
$^{-1}$
(Chang et al. Reference Chang, Protopapas, Kim and Byun2013). In contrast, SPB stars are hotter objects, typically of spectral types B3–B9. Their oscillations occur with frequencies ranging from roughly 0.3 to 1.3 d
$^{-1}$
(Aerts et al. Reference Aerts, Christensen-Dalsgaard and Kurtz2010). These two types of pulsating stars have their own instability regions in the HR diagram, on evolving into which, oscillatory behaviour is expected from current theory.
We have plotted the close binary components Aa and Ab on the instability strips of both
$\delta$
Scuti and SPB stars, as shown in Figure 19. The diagram was plotted with the same Z value as Figure 16 assuming single-star evolution. In Figure 19, the component Aa lies within the SPB instability strip, while the cooler component Ab is located within that of
$\delta$
Scuti. On this basis, the hotter component should be an SPB and the secondary a
$\delta$
Scuti star. However, we must take into account also the B component, which has mid-A spectral type components (Veramendi & González Reference Veramendi and González2014).
Figure 19. Position of the primary (Aa blue dot) and secondary (Ab red dot) binary components on the instability strips of
$\delta$
Scuti (below solid red lines Murphy et al. Reference Murphy, Hey, Van Reeth and Bedding2019) and SPB (above solid blue line Pamyatnykh Reference Pamyatnykh1999) stars. The theoretical evolutionary tracks (faint dotted lines) were taken from the MESA Isochrones and Stellar Tracks (MIST; Dotter Reference Dotter2016) and were generated using the same input parameters as in Figures 16.
Since this companion binary is located 8" away from the binary companion and the TESS pixel size is larger than this separation, we are at present unable to distinguish whether the
$\delta$
Scuti -type oscillations originate from this visual companion. Additionally, there are hybrid
$\delta$
Scuti
$-$
$\gamma$
Doradus stars (Uytterhoeven et al. Reference Uytterhoeven2011) that exhibit both low and high frequency pulsational effects. The component Ab could also be a
$\delta$
Scuti
$-$
$\gamma$
Doradus star. In such systems, there is generally a gap between the low-frequency
$\gamma$
Doradus-type oscillations and the higher-frequency
$\delta$
Scuti ones (Grigahcène et al. Reference Grigahcène2010). In our data, however, no such gap appears in the amplitude spectrum, and we may exclude this possibility.
To summarise, we infer that the hotter component of the A system is likely to be an SPB object. The
$\delta$
Scuti -type oscillations could originate either from the cooler component of the close binary or the B component of NO Pup.
As a check, we examined the position of the
$\delta$
Scuti pulsator on the proposed relationship between orbital (
$P_o$
) and pulsation periods (
$P_{puls}$
) (Kahraman Aliçavuş et al. Reference Kahraman Aliçavuş, Soydugan, Smalley and Kubát2017). This relationship was derived from data on eclipsing binaries with
$\delta$
Scuti pulsating components. Due to binary effects on oscillation, a correlation between pulsation and orbital periods is observed. This relationship could be used to determine whether the
$\delta$
Scuti pulsator is in the A or B binaries.
The position of the
$\delta$
Scuti star on the log(
$P_o$
)-log(
$P_{puls}$
) relationship is shown in Figure 20. As seen in the figure, the position of the
$\delta$
Scuti pulsator follows the general trend of the relationship. We conclude that the secondary in the eclipsing binary system, that is, Ab, is probably the
$\delta$
Scuti variable.
Figure 20. The log
$P_o$
- log
$P_puls$
relationship and the position of the
$\delta$
Scuti star (triangle) on it. The dots represent the known
$\delta$
Scuti stars in eclipsing binaries, taken from Kahraman Aliçavuş et al. (Reference Kahraman Aliçavuş, Soydugan, Smalley and Kubát2017).
8. Discussion and concluding remarks
NO Pup belongs to a small, young association (Tokovinin Reference Tokovinin, Chalabaev, Shatsky and Beuzit1999), with an estimated age of
$\sim \! 20$
Myr (see Section 5). The Aab component has separation (0.03956 AU), that is, three orders of magnitude smaller than its astrometric binary B component. In its present configuration, therefore, the dynamical behaviour of A is essentially governed by rotational and tidal distortions between the Aab components alone. Given the short period (
$P \approx 1.256$
d) and close separation (
$R/a \approx 0.2$
), the relatively high orbital eccentricity (
$e \approx 0.13$
) appears anomalous. Tidal friction should tend to circularise the orbit in the available time. On the other hand, very close encounters (of the order of the Aab separation) with other stars in the system afford a possible way to increase orbital eccentricity – for example, through the Kozai–Lidov mechanism (Naoz et al. Reference Naoz, Farr, Lithwick, Rasio and Teyssandier2013). Therefore, it is necessary to determine the timescales over which tidal friction and close 3-body encounters with Aab occur, and whether they may have a substantial effect within the system’s age.
The components of NO Pup B, from the information given in the WDS Catalogue, are Main Sequence stars with V-magnitudes corresponding to spectral types A5V (
$ M \approx 2.0\,{\rm M}_\odot$
) and A6V (
$M \approx 1.8\,{\rm M}_\odot$
). From standard stellar modelling they should possess convective cores and radiative envelopes. Within these envelopes, internal (
$g-$
mode) gravity waves can propagate and efficiently transport angular momentum from the convective-core boundary to the outer layers (Zahn Reference Zahn1977; Goldreich & Nicholson Reference Goldreich and Nicholson1989; Aerts, Christensen-Dalsgaard, & Kurtz Reference Aerts, Christensen-Dalsgaard and Kurtz2010). Such waves are partly dissipated through thermal damping in the thin, non-adiabatic surface region (Khaliullin & Khaliullina Reference Khaliullin2010; Townsend, Goldstein, & Zweibel Reference Townsend, Goldstein and Zweibel2018). Tidal forcing (Zahn Reference Zahn1975; Zahn Reference Zahn1977) operates mostly in these outer layers with radiative damping, creating spin-up/down in the visible surface, which may thus be a poor indicator of the internal rotational angular momentum.
Orbital circularisation by tidal friction depends on the age of the system. As a measure of its efficacy, one may specify an upper limiting value of the ratio
$a/R$
of orbit size to stellar radius, below which circularisation can be achieved within an interval of one-quarter of the star’s Main Sequence lifetime (Zahn Reference Zahn1977, Table 2). For NO Pup Aab, using the Table 11 masses, one finds this upper limiting value
$(a/R)_{lim}$
to be
$\sim \! 4.0$
. Since for the Aa component of NO Pup,
$(a/R) \sim \! 4.2$
, the case for effective tidal circularisation of the orbit appears marginal, especially given our lack of knowledge of the rotational structure of the stars. Zahn’s analysis indicates that the upper limiting value of
$a/R$
for synchronisation would be
$\sim \! 5.6$
, suggesting that tidal friction may have de-spun the components into synchronicity. The spectroscopic measurements set out in Tables 7, 9, and 10 in Section 3, including values of
$v\sin i$
, enables assessment of synchronicity of the Aab stars to the available accuracy of measurement. If the pseudo-synchronisation formula proposed by Hut (Reference Hut1981) (in his Equation 42) is used for the components of NO Pup A, rotation velocities 95 and 66 km s
$^{-1}$
for the components Aa and Ab, respectively, are obtained. Comparison of these estimated velocity values with the measured rotation velocities (82 and 64 km s
$^{-1}$
for the components Aa and Ab; see Sections 3.4 and 4), would suggest that, in this eccentric binary, the rotation and orbital motion are pseudo-synchronised for Ab but not for Aa, as noted by Veramendi & González (Reference Veramendi and González2014). However, Aa appears to have a rotation period of 1.32 d – close to the orbital period, as noted in Section 3.4 – and thus nearly synchronised to the mean orbital motion. Confirmation of these suggestions would be possible with the availability of more extensive and more accurate data.
As noted in Section 7, the Ab component lies in the instability region of the Hertzsprung–Russell diagram for the
$\delta$
Scuti variables (Aerts et al. Reference Aerts, Christensen-Dalsgaard and Kurtz2010, Figure 1.12). These are likely to experience p-mode oscillations, whereas in the higher-mass Aa component, one would expect g-modes excited at the convective core/radiative envelope boundary. In addition, the high eccentricity of the orbit means that a large number of harmonics of the orbital frequency could be excited, thereby providing possibilities for resonance with one or more of the p- and g- stellar oscillation modes, that in turn may be split by stellar rotation. The increase in amplitude associated with these resonances could provide additional tidal damping. However, harmonics of the orbital frequency
$(f_{orb} = 0.7956\,{\rm d}^{-1})$
become separated out with the general binary effects in the pre-whitening process.
A possible alternative explanation for the observed orbital eccentricity is by a three-body interaction. Since the orbital size of the Aab eclipsing binary system is much smaller than that of the Bab astrometric binary, we may regard Aab as a single point mass. The resulting hierarchical 3-body system may be subject to the classical Kozai–Lidov mechanism, mentioned above. Consider the Aab–Bab system, in which Aab is approximated by a single point mass The Aab component (mass =
$M_{Aab} = 5.26\,{\rm M}_\odot$
) is in a wide orbit about the centre-of-mass of the Bab (astrometric) system (masses
$M_{Ba}, M_{Bb}$
). The Kozai–Lidov cycle time
$t_{KL}$
for such a system, to the quadrupole level of approximation (Naoz Reference Naoz2016, Equation (27)), is:
\begin{align*}{t_{KL}} \approx \frac{{16}}{{30\pi }}\frac{{{M_{tot}}}}{{{M_{Aab}}}}\frac{{P_{wide}^2}}{{{P_{close}}}}{(1 - e_{wide}^2)^{3/2}}\end{align*}
where ‘close’ and ‘wide’ refer to the (Ba,Bb) astrometric binary and the wide orbit of Aab about Bab’s centre-of-mass, respectively.
$M_{tot} = M_{Ba}+M_{Bb}+M_{Aab}$
is the total mass.
Neglecting the value of
$e^2_{wide}$
, and from Table 11 using
$M_{Aab} = 3.58+1.68 = 5.26\,{\rm M}_\odot$
, then with the values given above for the masses
$M_{Ba} = 2\,{\rm M}_\odot, M_{Bb} = 1.8\,{\rm M}_\odot$
, so
$M_{Ba}+M_{Bb} = 3.8\,{\rm M}_\odot$
; the distance d being 171 pc, we find
$a_{close} = d\times 0.179$
arcsec
$= 30.609$
AU. Now the B component has orbital period
$P_{close} = 101.3$
yr, while for the wide orbit, assuming the apparent (projected) angular separation of 8.1 arcsec corresponds to the semi-major axis, we derive
$a_{wide} = 1\,385$
AU. Kepler’s Third Law then yields
$P_{wide} \sim 17\,100$
yr, on using
$M_{tot} = 9.06\,{\rm M}_\odot$
. For the NO Pup system, we find
$t_{KL} \approx 8.5\times 10^5$
yr – or at least an order of magnitude smaller than the estimated age of the system.Footnote
j
This would clearly imply the possibility of the Kozai–Lidow effect operating in the NO Pup system, to produce apparently anomalous parameters.
In addition to classical Kozai–Lidow processes in this quadruple system’s history, since NO Pup is a member of an association (Tokovinin Reference Tokovinin, Chalabaev, Shatsky and Beuzit1999), it may previously have comprised several more components that have since escaped via close binary-binary or 3-body interactions (Aarseth & Mardling Reference Aarseth, Mardling, Podsiadlowski, Rappaport, King, D’Antona and Burderi2001; Mikkola Reference Mikkola2010). Of interest is the short-term dynamical behaviour of the system.
To gain insights, we may replace Aab by a point mass, as before and use the available observational data on NO Pup. Initial conditions for a 3-body model system were then chosen to be
$a_{wide} = 1\,385$
AU; the orbit was assumed to be circular:
$e_{wide} = 0$
, and edge-on (
$i_{wide} = 90 \deg$
); and masses
$(M_{Ba},M_{Bb},M_{Aab}) = (2,1.8,5.26)\,{\rm M}_\odot$
. This system was integrated over a time interval
$T = 1\times 10^7$
yr, implementing an Aarseth–Zare (AZ) regularised scheme (Aarseth & Zare Reference Aarseth and Zare1974) with Kustaanheimo–Stiefel regularisation (Kustaanheimo & Stiefel Reference Kustaanheimo and Stiefel1965) to remove singularities in the equations of motion due to possible 2-body collisions. The AZ scheme was adapted so that, after each set of 2 000 time-steps of integration, the Keplerian orbital parameters corresponding to the current smallest two of the three interparticle separations were computed. The time transformation t(s) between ‘regularised time’ s and physical time t was chosen to be the Lagrangian, which gives the most accurate integration scheme and has the advantage that the function t(s) is known explicitly (Alexander Reference Alexander1986). Constancy of total energy and angular momentum was maintained to one part in
$10^{14}$
and
$10^{12}$
, respectively. In addition, integrations of the Kozai–Lidov equations, taking into account quadrupole and octupole level perturbations (Naoz Reference Naoz2016) were performed. The results of AZ integrations for
$e_{close}, I_{close}$
are shown in Figure 21.
Figure 21. Top: Eccentricity of astrometric binary system (Ba,Bb) over a time interval of 10 Myr, from integrating the reduced 3-body system (Aab,Ba,Bb) using the Aarseth–Zare regularisation scheme. MCMC data from Table 11 was used. The Kozai cycle is clearly evident. Bottom: Inclination between the close orbital and invariable planes.
Numerical integrations of the AZ scheme over 10 Myr revealed the existence of Kozai cycling with respect to the B components, with a cycle time of 4.5 Myr. During the span of the integrations, for the astrometric binary (B components),
$e_{close}$
twice reached a maximum of
$e_{close} = 0.902$
, corresponding to a minimum separation between the B components of 3.0 AU (
$\approx 645\,{\rm R}_\odot$
). The perturbation approach (quadrupole or octupole) may become suspect at such close encounters. However, the distances between Aab and both of the B components remain large (
$\approx$
1 400 au) throughout the Kozai cycle, so that tidal effects between Aab and Bab are negligible. By the same token, apsidal motion in the Aab eclipsing system will be governed exclusively by tidal and rotational distortions of its components.
With regard to the B components, the quadrupole formula above for
$t_{KL}$
yields a value
$\sim 5$
times smaller than was found with the much more accurate AZ numerical integrations. In the immediate vicinity of closest approach, where
$e_{close}, I_{close}$
were rapidly changing, there was poor agreement between the perturbative and AZ results; however, away from these events, there was fairly good agreement.
In this paper, we have marshalled evidence from the available data streams – photometric, spectrometric, and astrometric – so as to quantify parameters as accurately as possible. This allows well-defined and comprehensive analyses. To this end we have combined various curve-fitting techniques to determine optimal values and uncertainties. Corresponding results are presented in the preceding tables and diagrams. In turn, this permits more detailed theoretical discussion.
Acknowledgement
This paper includes data collected by the TESS mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. This research has made use of the Washington Double Star (WDS) Catalogue maintained at the U.S. Naval Observatory. We thank Dr. Rachel Matson for extracting data from the WDS for us, and also the University of Queensland for collaboration software. FKA thanks the Scientific and Technological Research Council (TUBITAK) project 120F330 for supporting the study. The authors thank Prof. G. Handler for his valuable comments on pulsational analysis. This work uses the VizieR catalogue access tool, CDS, Strasbourg, France; the SIMBAD database, operated at CDS, Strasbourg, France. We thank the UCMJO time allocation committee for observing time with hercules. We recently learned of the death in December 2024 of Sverre Aarseth. He was a pioneer in N-body (both small and large N) numerical integrations and inspired many people in this field, including one of the authors (MA). He will be greatly missed.
Data availability statement
The TESS data are available from the MAST data archive (https://archive.stsci.edu/). The astrometric data may be obtained from United States Naval Observatory on request. The other data supporting this study are available upon reasonable request to the corresponding author.
Funding statement
This research was partially supported by a grant from the Scientific and Technological Research Council of Türkiye (TUBITAK project 120F330).
Competing interests
None.
Appendix A
Table A1. Radial velocity values of NO Pup derived from the HARPS spectra.
Table A2. Identified spectral lines for NO Pup based on comparison with the ILLSS Catalogue (Coluzzi Reference Coluzzi1993; Coluzzi Reference Coluzzi1999). The lines are confidently detected mainly for the primary (Aa).
Open access
