Observed epochal variations in X-ray lines from the O Supergiant $\zeta$ Puppis do not require substantial changes in the wind mass flux

We fit the high resolution \textit{Chandra} X-ray spectra of the O supergiant $\zeta$ Puppis using the variable boundary condition (VBC) line model to test the stability of its mass-loss rate between two epochs of observation: 2000 March and 2018 July -- 2019 August. At issue is whether the observed variations are induced by global changes in the cool (unshocked) wind itself or are isolated to the local pockets of hot gas (i.e., changes in the frequency and location of the shocks). Evidence in the literature favored the possibility of a 40 per cent increase in the mass flux of the entire stellar wind, based on X-ray reabsorption from a line-deshadowing-instability-inspired parameterization, whereas our fit parameters are consistent with a constant mass flux with a change in the velocity variations that determine the locations where shocks form. Our results suggest the shocks in the more recent data are formed at somewhat larger radii, mimicking the enhanced blueshifts and increased line fluxes interpreted in the previous analysis as being due to increases in both the X-ray generation and reabsorption from an overall stronger wind.


INTRODUCTION
At only 332 pc away (Howarth & van Leeuwen 2019),  Puppis (HD 66811) is one of the brightest O stars in the night sky.This has made it the the canonical O supergiant and one of the most well studied single massive stars across all wave bands.Its well studied nature has also made it an ideal calibration target for some of our most sensitive X-ray instruments (e.g., XMM-Newton's RGS and EPIC 1 ).As a massive star,  Pup's stellar wind is relatively well understood to be driven by radiative forces acting on hypersonically Doppler shifted UV line opacity from metal ions in its atmosphere, As such, it is expected to have a wind that, despite being clumpy (Martínez-Núñez et al. 2017), should on the whole be as steady as its luminosity.This is true even in the case of the star being a rapid rotator, as is the case for  Pup.The source of the high rotation rate, along with its runaway status, is thought to be evidence of a prior companion in  Pup's evolutionary history whose supernova resulted in  Pup's high speed ejection (Woermann et al. 2001).
While many of  Pup's stellar properties have been well constrained, a number of studies have also revealed that our picture ★ Contact e-mail: seang97@mit.edu 1 https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/uhb/routinecal.html is perhaps incomplete.Observations in both the optical (Ramiaramanantsoa et al. (2018) using BRITE) and X-ray (Nazé et al. (2018) using XMM-Newton and Swift and recounting earlier documented variabilities, and Nichols et al. (2021) using Chandra) regimes have highlighted apparent variabilities in the star.The main periodicity detected in the optical band has a period of 1.78 days (originally detected by Howarth & Stevens (2014)) and been argued to arise from bright spots on  Pup's surface and its rapid rotation (Ramiaramanantsoa et al. 2018;Nichols et al. 2021).There is also an apparent 2.56-day periodicity (Marchenko et al. 1998) that Howarth & van Leeuwen (2019) corroborated in the Hipparcos data, which interestingly does not show evidence of the aforementioned 1.78-day signal.
Adding another layer of complexity, the periodicities do not appear to be constant.Nazé et al. (2018) noted that the X-ray periodicity appear to change in phase and period across many XMM-Newton observations, at times correlating with the BRITE signal and other times not.The strength of the X-ray signal also appears to wax and wane: it was stronger during the 2007-2011 than before or after.The optical 1.78-day signal also appears to be varying as Howarth & van Leeuwen (2019) posited that its amplitude has grown with time to explain why Hipparcos does not show it.Older observations reported longer period variations than those today.These include periods of 5.21 days (Balona 1992) and 5.075 days (Moffat & Michaud 1981) using optical data, and 5.1 days (Howarth et al. 1995) in the UV.How these longer periods fit into the picture is less clear.They were originally proposed to be connected to the rotation of  Pup, but the rotation period is faster (<3 days; Howarth & van Leeuwen 2019) than these signal periods.Thus there appears to be multiple layers of periodicity which are themselves varying in frequency, amplitude, and phase relationship between various wavebands.
There is an important distinction to be made between X-ray variability and variability in other wavebands.As a massive star,  Pup generates X-rays not at its static photosphere but farther out in its hypersonic wind through embedded shocks (Lucy & White 1980;Feldmeier et al. 1997) as clumps experience differential acceleration.This X-ray generating mechanism is dependent on the number, size, speed, and many other factors of the wind clumps, so it has an inherent stochasticity.This randomness is associated with the flow time of the wind  flow =  * / ∞ , which is on the order of a couple of hours, a much shorter timescale than those discussed above.An example of a possible short-period signal is the 0.694-day period Berghoefer et al. (1996, reported as 1.44 cycles/day) found in simultaneous X-ray and H observations.X-ray observations are typically long enough compared with this timescale to average out any stochastic effects (see Table 2).The noted X-ray periodicities are more coherent signals that persist through the time averaging, so changes in the X-ray signal may correspond to global changes in the time-averaged wind.
In summary, many recent investigations of the relationship between optical and X-ray variations of this star on short time scales have revealed interesting links between the two wave bands.Hence it is timely to ask if there exists longer term changes, either in pockets of X-ray generating shocks, or more globally over the cool wind.Capitalizing on the decades of high resolution X-ray profile information that exists for  Puppis, our focus is exploring long term changes in the wind and its X-ray heating The first piece of evidence for changes in  Pup as an X-ray source is its brightening in the last few years as detected with both Chandra (Cohen et al. 2020) and XMM-Newton.The latter is shown in Figure 1, where the plotted flux values were measured by XMM-Newton's EPIC detector and reported in 4XMM-DR13 (Webb et al. 2020) Cohen et al. (2020) proposed that the increase in X-ray flux is due to the wind mass-flux increasing (more on this immediately below), but the exact origin is not entirely clear due to the the many factors which influence overall X-ray output from embedded wind shocks.More mass within the wind would mean more clumps and subsequently more shocks, but the increase in flux could also be caused by stronger shocks producing more light or a wind that is more porous to X-rays due to changes in clump size or density.
The second is changes in the wind as inferred from a line shape analysis, something that can only be measured using the higher spectral resolution of Chandra (though the time coverage is much less frequent).Constraints on the physical effects listed above can be extracted from the X-ray spectral line shapes via some appropriate parameterization of its features.Understanding what the increased flux and changes to the line shapes allude to requires an appropriate line model to parameterize quantities of interest.The inference of the physical changes to the wind using a line profile parameterization is a goal of this paper, for which the chosen parameterization will be discussed in § 2.
An important parameterization comes from Owocki & Cohen (2001), which Cohen et al. (2010Cohen et al. ( , 2020) ) used to calculate the massloss rate of  Pup from X-ray lines measured by Chandra.In this paper we will be referring to this model as simply the "Cohen model" to account for the number of papers that have made use of it (Owocki  2 for dates).Values were retrieved from 4XMM-DR13 (Webb et al. 2020).The reported errors in the catalogue are too small to visualize as vertical bars so are plotted in the sub-plot.& Cohen 2001;Cohen et al. 2010Cohen et al. , 2014Cohen et al. , 2020)).The most recent of these papers, Cohen et al. (2020) came to the remarkable conclusion that the mass-loss rate of  Pup has increased by 40 per cent in a roughly 20-year time frame, based on the changes to the line shapes (i.e., centroid and width) in Chandra High Resolution Transmission Grating System (HETGS) data from Cycles 1 and 19 (see Table 2 for observation dates).Such a large change in the mass-loss rate of a star brings to light the need to independently verify that observed profile changes are robustly interpretable as changes in the global wind, rather than in just its hot component, given that the hot component is a small fraction of the total wind (Cohen et al. 1997;Owocki et al. 2013;Gayley 2016) and hence might be more susceptible to larger variations.

The Challenge of Inferring Mass-loss Rates
The problem of determining mass-loss rates of massive stars is not a simple one.In Table 1, we list a number of calculated mass-loss rates for  Pup in the literature.The theoretical value expected from Vink et al. (2000) is unsurprisingly high as their theory does not account for the clumps that make up the wind.Puls et al. (2006) made detailed computations of the clumping-included mass-loss rates using H, IR, mm, and radio, finding values that bracket the predicted value from Vink et al. (2000).The wide difference between the values is a reflection of differing assumptions of where He recombination starts in the wind.These four wave bands provide precise values for the mass-loss rate, but the specific values determined are still heavily dependent on the assumed clumping factor of the wind.The X-ray values from Cohen et al. (2010) and Cohen et al. (2020) are free of this assumption since small scale clumping does not affect the reabsorption of X-rays by the cool wind.This allows measurements of X-ray line profiles to directly parameterize the mass-loss rate (see Section 2).Additionally, if the X-ray lines are fit individually in isolated regions, the (relatively flat) local continuum requires no correction for interstellar absorption, so there is no dependence on the distance to  Pup.Thus there is no need for distance corrections like those made by Howarth & van Leeuwen (2019) to the flux-dependent calculations done by Puls et al. (2006).However, the use of X-rays for inferring mass-loss rates is not without its own problems.Gunderson et al. (2022 determined that these estimates are dependent on the assumptions going into the wind profile model used to fit the data.

Expected Stability of the Global Wind Mass Flux
In classical theories of massive stars (Castor et al. 1975;Owocki 2004), the mass-loss rate is regarded as a stable quantity due to the nature of the line-driving of the wind.Therefore, assuming one's calculated value can be trusted, a large change in the mass-loss rate may be indcative of an incomplete understanding of the stability of these stars over longer periods.Variations on long timescales is not unheard of for stars; for example the Sun undergoes a 22 year cycle (Russell et al. 2016).The disks of Be stars also exhibit similar long timescales (Kee 2020).
To know if Cohen et al. (2020)'s reported changes are part of a long-period variations, we need much more data to see more than one period to determine if it is indeed even periodic, or stochastic change over a longer time scale.As a first step, however, we need to know if the secular changes that have been highlighted in our above discussion correspond to real changes in the global wind properties.
Thus is the goal of this paper.We will investigate the X-ray line shapes of  Pup for the Cycle 1 and 19 Chandra observations to investigate the large change in mass-loss rate reported by Cohen et al. (2020).To do so, we will use the VBC model from G2022.While the specific values of mass-loss rates derived through X-ray lines are dependent on the model used, any large change in a global parameter such as the mass-loss rate should be relativity model independent.
To explore this, we will be comparing many of our modelling results with those of Cohen et al. (2010) and Cohen et al. (2020), so we will hereafter refer to these works as C2010 and C2020 respectively.
This paper is organized into the following sections.In § 2, we give a short summary of the VBC model while including some new insights.In § 3, we provide details on the data reduction and modelling.Finally, in § 4 and 5, we give our results and conclusions.

VARIABLE BOUNDARY CONDITION MODEL
For a more detailed derivation of the VBC model, readers are referred to G2022.Here we give a summary of the model and its parameters, which starts with the wind's velocity.As with most massive-star wind models, the velocity field is assumed to follow the usual -velocity law We choose  = 1 to allow for analytic results in subsequent equations.For example, with  = 1 the optical depth takes the form where  = cos  and The variable is a fiducial optical depth describing the amount of absorption a photon produced at  =  * experiences in a constant velocity wind ( = 0) with the same mass-loss rate as an accelerating wind.For the case of  = 1, it is instead the optical depth at a radius of / * = e/(e − 1) ≈ 1.58 for a photon on the same central ray. 2 For information on how  * can be used to calculate a mass-loss rate, see C2010.Additionally, details on the role of complex numbers in the derivations of Equations ( 2)-( 5) are given in Appendix A.
The main distinguishing feature of the VBC model is its direct parameterization of the heating of the gas.This is in contrast to the use of a hot-gas filling factor, a number which implicitly conflates the effects of the heating and cooling rates.This is achieved through tracking pockets of gas that are assumed to follow a slightly faster velocity field with terminal velocities of  ∞,  than the ambient, slower pockets.Note that this velocity is not used in the definition of  * because it is the ambient, unshocked gas described by  ∞ that absorbs the photons generated in shocks.The faster clumps are assumed to have a differential probability d d of shocking with a slower clump.The above equation contains the final model parameter ℓ 0 which parameterizes the mean free path between shocks.Thus for a given line can be described by how much optical depth  * the photons traveled through, a characteristic formation radius  * + ℓ 0 , and the terminal speed  ∞,  of the faster clump that caused the shocked.The line shape described by these parameters is modelled using the relative luminosity per frequency bin where  ≡ −()/ ∞,  is the frequency shift from line center in terminal speed units and is the local mapping from solid angle to frequency space.The lower bound of the integral is the minimum radius capable of producing enough Doppler shift to reach the  in question.It is a complicated function dependent on whether the emission is toward the forward ( > 0,  < 0) or backward ( < 0,  > 0) hemisphere of the star.
Readers are encouraged to see G2022 for details on this minimum radius.In Appendix B, we give details on approximating Equation ( 8) with a Gauss-Laguerre quadrature.

Terminal Velocity Parameter
In G2022, a constant, frozen fast-gas terminal velocity of  ∞,  = 2500 km s −1 for each line was assumed.This was done to simplify the model fitting.We use the same process in this work but provide a more mathematical backing for this choice.
Let us assume that the fast gas terminal velocity is a simple scaling constant with respect to the slow gas,  ∞,  =  ∞ .Then, for a given pocket of fast gas that shocks within a mean-free path, its velocity is Thus the true parameter to investigate is ℓ 0 .This is difficult to interpret and compare with existing work however, as other line profile models only parameterize the slow gas terminal velocity.It is for this reason that we choose to fix the fast-gas terminal velocity, i.e. assume a fixed , and vary ℓ 0 to better compare against spatial parameters in other models.

Source Function Description
In this section we will discuss how the VBC model compares with previously published models when put into the standard language.In general, a line profile function should have the form where (, ) is the emissivity function of the light produced in the volume  of the wind.An example of this  function is the emissivity   used in the Cohen model.If we change the integration of Equation ( 8) to that of volume, it follows that our model's equivalent emissivity function (when converted back to radial form) is Example curves for this emissivity function are plotted for various ℓ 0 values in Figure 2.For comparison, in the bottom half of that figure we have plotted the emissivity function of the Cohen model.The full emissivity function is given in Owocki & Cohen (2001), but for our purposes, for a photon emitted at line center, it scales as where  is an X-ray emission filling factor.The authors of this model assume, based on the work of Ignace ( 2001   of .The velocity function () which appears in  = /4() 2 is calculated using a -law with the value  = 1. Figure 2 highlights an important difference between these model parameterizations: the extent of the X-ray emission.Using a constant or slowly varying filling factor causes the emission to stem from a rather narrow radius interval where the cool wind emission measure is high.But since the VBC approach parameterizes the shock heating and assumes efficient radiative cooling, the cool wind emission measure is of no consequence, and the X-ray generation can be more extended if the shock distribution is more extended.Ultimately, only observations and more detailed simulations can clarify which approach better characterizes the spatial distribution of the heating.Therefore our current intent is to understand how robust the conclu- sions inferred from line profile shapes are when working within these two schemes.

DATA AND MODELLING
The observations of  Pup included all archival HETGS data available in the Chandra Archive.The list of observations with their Observation ID (Obs ID) are given in Table 2.Each observation used the Chandra HETGS instrument which simultaneously provides data from two grating arrays: the medium energy grating (MEG) and high energy grating (HEG).These gratings have resolutions of 0.023 Å and 0.012 Å respectively (Canizares et al. 2005).
Each observation was retrieved from the Chandra archive and reprocessed using the standard pipeline in ciao version 4.13 (Fruscione et al. 2006).This process produced two first-order spectra for both the HEG and MEG arrays, corresponding to the positive and negative diffraction orders.Each of the positive and negative orders were subsequently co-added to produce a single first order data set for their respective grating array and Obs ID.For the Cycle 19 datasets, we additionally summed the 42 datasets (21 for the HEG and MEG respectively) to produce a single total observation for each of individual grating array.Thus in the end we had 1 HEG and 1 MEG spectra for Cycle 1 and 1 HEG and 1 MEG spectra for Cycle 19.These final datasets were then rebinned by a constant factor of 3 for model fitting.This corresponds to bins sizes of 0.015 Å for the MEG and 0.0075 Å for the HEG.
When comparing data taken by Chandra separated by almost two decades, we needed to account for changes in the detector properties during that time.The most important change is a result of contamination build up on Chandra's ACIS-S detector used for HETGS observations (Marshall et al. 2004;O'Dell et al. 2017).This contamination decreases the detector's effective area, particular in the long-wavelength region of HETGS spectra, reducing overall count rates.However, this change in the detector should not have an important affect on the analysis carried out here.The contamination is calibrated in the standard data products (i.e., date-dependent response files), and we have useful diagnostics out to 17 Å as indicated in Figure 3.In addition, it is worthwhile to note that the Cycle 19 observations have a much greater exposure time, so even with the greater contamination effects there are many more total counts in most of the X-ray lines compared to the Cycle 1 observation.
Using xspec version 12.12 (Arnaud 1996), through the python wrap around pyxspec, we fit the same list of lines as in G2022 in isolated regions.The total source model used consisted of the VBC model3 discussed above plus a constant continuum component.For the line features Si xiii He ,4 He-like Lyman  is shortened to Si xiii He .Mg xi He , Fe xvii at 15.01 Å, and Fe xvii at 16.78 Å, a total of three line model components were used to model the three lines that occur in close spacing at these wavelengths.The  * and ℓ 0 were the same for each line of these multi-component fits, being tied to the line with the shortest wavelength of the feature's fitted region.The shortest wavelength emission line was chosen for being the strongest single line in these features.Each line component in a feature had a freely fit normalization .We assumed Poisson statistics and used a Cash (maximum likelihood) statistic (Cash 1979).

Model Fitting Behavior
The fitting algorithm used for our analysis was the Markov chain Monte Carlo (MCMC) function provided in pyxspec.For each fit, the resulting chain was analyzed for all steps after a convergence was achieved, giving parameter distributions like those in Figure C1 in Appendix C for the Mg xii H  line.This line is within an isolated region with little continuum, so it provides a well constrained example of the model behavior, There are two interesting trends that can be seen in the parameter pairs.First is the inverse relationship between  * and ℓ 0 .This is similar to the discussion in G2022 when comparing observableparameter pair contours.A more interesting effect is between ℓ 0 and the normalization, which show a slight positive correlation.This is likely due to the factor of 1/ℓ 0 that applies to the entire integral in Equation ( 8) through the relationship defined in Equation ( 7).Since such a constant factor would normally be considered a part of the normalization, it is not surprising that these would be linked in this way.

RESULTS
The results of our modelling are summarized in Table 3 where we report the best-fit values of free parameters in our model.Note that for line features fit with multiple components we report the total flux in the feature.These best-fit values were extracted from the MCMC chains after they stabilized along with their 68 per cent confidence intervals.The actual model fits are plotted in Appendix D. Figures 3  and 4 show these parameter results visually.

Determination of Mass-loss Rate
Focusing first on Figure 3, we can draw two important conclusions.First, it appears that mass-loss rates calculated from X-ray line profiles is independent of the assumed heating model.The best-fit values for our Cycle 1 analysis (blue dots) and C2010 (green up-triangles) show an over-all consistency when looking at 68 per cent confidence intervals.This is further bolstered by the values from our Cycle 19 analysis (orange squares) as the consistency is maintained.For all three fit groups, we would find the same mass-loss rate, since  * is a simple rescaling of  based on Equation ( 6), save for deviations accounted for by errors in a trend-line fitting process.
Given this overall consistency, the second conclusion we draw from Figure 3 is that  Pup's mass-loss rate has not changed.If such a change occurred, particularly one suggested to be so large, our independent analysis on the same datasets would have found a similar result.However, we do confirm that there is a real change in the lines.The flux in all but one of our fitted line features has increased, just as was found by C2020.Combined with the XMM-Newton flux in Figure 1, there is high confidence in more X-rays being produced and/or escaping the wind.

Issues with Epochal Analysis
The result that mass-loss rate determinations from X-ray lines is model independent is in contrast to G2022 that concluded the opposite by analyzing only Cycle 19 data.The choice in only looking at Cycle 19 data may have been the source of the erroneous conclusion as G2022 was comparing against the results of C2020, which is the only set to break from the noted consistency.Why this is the case is not immediately clear.For lines below 10 Å, all four groups are in agreement except for Si xiii He , though only marginally so.For the longer wavelengths, there are large differences in the values found for the Cohen model.It is possible that the limited fitting region of the Ne x H  and Fe xvii at 15 Å lines used by C2020 is responsible for some of the inflated values of  * .The fits presented here contrast this by fitting the entire red-half of the profile of Ne x H  and all three lines in the Fe xvii at 15 Å feature (Figures D3 and D6).This is simply speculation as to what differences may exist between the two modelling approaches that would result in the different  * values found.Whatever the source is, there may be anomalous behavior in the C2020 fitting analysis that is also responsible for their claim of a 40 per cent increase in  Pup's mass-loss rate.
From their values in Table 3, we argue that this difference stems from the ℓ 0 parameter, which is also the source of the line shape differences that were noted by C2020.Put another way: the distribution of the embedded wind shocks has changed.Figure 4 more clearly shows that the mean-free path between shocks in  Pup's wind has increased by a significant amount.Note that in this figure we plot  * + ℓ 0 so as to compare directly with the turn-on radius  0 of the Cohen model.What this means is that the shocks are occurring farther out in the wind, allowing for more X-ray emission to escape, thereby raising the amount of flux that we see.
This presents two competing interpretations of the same data: either the mass-loss rate has changed or the shock distribution.In comparing these, it is worth considering the implications each has over the entire wind.Given C2020's increased mass-loss rate, the entire wind is required to be changed to account for the differences we see in the line profile shapes.Without evidence of  Pup having large scale variabilities, i.e. those that persist after time-averaging during an observation, that alter the stellar parameters, such a large increase in the mass-loss rate would imply either that  Pup has undiscovered properties or our understanding of massive stars, particularly Otypes, is missing an important piece.Given the importance of  Pup as a prototype, guide star, and calibration target, it has many further implications in every other wave band.The conclusion from the VBC model analysis that only the shock distribution has changed, and not the global mass-loss rate, has less drastic implications involving a much smaller fraction of the wind.

Shock-Temperature Spatial Distribution
The primary motivation for using the VBC analysis is to understand how the gas is heated, not to determine the cool wind mass loss rate.So now we turn to the issue of using the VBC parameterization to interpret the shock distribution, and contrast that with the picture from Cohen's parametrization.For example, part of the motivation of using a turn-on radius  0 in the Cohen model is the idea that the line deshadowing instability requires a certain stand-off distance in order to have room to take effect, whereas the VBC picture holds that fast gas is destined to overtake slow gas over the length scale ℓ 0 .In the latter picture, it is natural that more significant boundary variations should yield more violent encounters between fast and slow wind, all at lower ℓ 0 , whereas in the former picture, there could be more of a tendency for  0 to be independent of shock strength, or even rise with shock strength (Driessen et al. 2021).
G2022 visually identified a trend in ℓ 0 vs. the temperature of maximum emissivity  max .That trend can be demonstrated quantitatively using the Pearson correlation coefficients for the plotted points in Figure 4.The values of the coefficients  ℓ 0  max are given in Table 4.Note that this correlation coefficient is invariant under linear transformations, so the calculated values are the same for both ℓ 0 and  * + ℓ 0 .Additionally, we calculated  ℓ 0  max using log( max ), but for readability of the coefficient subscripts, we omit the logarithm.We used a random sampling method of the ℓ 0 parameter distributions to calculate the reported values.The central values are the most frequent of the resulting  ℓ 0  max distribution while the uncertainties are then taken as the 68 per cent levels.
We are interested in investigating trends related to the shock strength, for which the temperature of maximum emissivity  max is a useful proxy.From the values in Table 4, a negative correlation does exist between the ℓ 0 and the shock strength for both Cycle 1 and Cycle 19.The trend is not exactly linear, which would correspond to  ℓ 0  max = −1, but it does imply an inverse trend nonetheless.
Of note is how the Cycle 1 results were handled.Two of the lines fit with our model, Si xiv H  and Ne ix He , show significantly more uncertainty in their best fit ℓ 0 compared to the rest of the lines.When these points are included, the correlation is consistent with 0, meaning no linear relationship exists.However, the rest of the Cycle 1 results show a tighter trend with log( max ), so we also calculated the correlation coefficients with these two lines removed.The restricted coefficient is constrained to only negative values and is statistically Notesa Predicted line center wavelength as reported by AtomDB (Foster et al. 2012;Smith et al. 2001).
b Temperature of maximum emissivity as reported by AtomDB (Foster et al. 2012;Smith et al. 2001).c This line is not confidently detected in the Cycle 1 data.d The reported flux for these line features is the total flux in the three components used.Note that we have applied constant offsets to the different result groups to make the error bars more visible.The plotted wavelengths correspond to  p for the blue dots,  p + 0.1 Å for the orange squares,  p + 0.2 Å for the green up triangles, and  p + 0.3 Å for the red down triangles.Data points that do not have a corresponding group of four are due to that line feature not being fit in some of the works considered here.consistent with the computed value with all lines considered.If the uncertainties in the outlier lines can be better constrained, the true value is likely to be within the overlap of these two values.These coefficients confirm the prediction that under a VBC initiation for a massive star's wind, the stronger shocks will occur at lower radii.There are interesting ramifications for how the shock heating distribution appears to have changed with time.While  0 parameter in the Cohen model is broadly unchanged, the ℓ 0 parameter of the VBC model has increased significantly.Such a difference between the temporal changes in the shock distributions is likely explained from the fact that in the VBC model it is only the hot gas that has changed.This requires a more significant adjustment of the heating distribution instead of the entire wind.The latter of which is required by the Cohen model based on the wind mass flux increase (C2020).Thus we can conclude that the hot gas has moved to farther radii in recent years while the cool, unshocked gas has not changed, though some changes can not be ruled out.

CONCLUSIONS
We present new and better constrained results for the VBC line model fit to both epochs of X-ray line emission data on  Pup as observed by Chandra.Based on the best-fitting values of the  * parameter, we find that  Pup's mass-loss rate has not undergone any large-scale changes within the near two decade between of the two observations.Additionally, the consistency between our Cycle 1 and 19  * values and those of C2010 implies that mass-loss rates inferred using X-ray line emission are independent of model parameterizing and thus robust as a technique.
It is concerning then that a 40 per cent change in the mass-loss rate shows up in one parameterization and not in the other.However, rather than reflecting an overall incompatibility between the two models, the difference is possibly from the specifics of the line analysis.For example, there appears to be a dependence on whether an entire line feature, such as the three resolved lines of the Fe xvii at 15 Å, is fit versus a single line of the feature.This is based on our modelling accounting for the three lines in this feature for both Cycle 1 and 19 and finding a smaller  * value than C2010 and C2020, the former still being within error, which fit only the Fe xvii line at 15.014 Å.It remains unclear what systematic uncertainties are present in the mass flux calculations due to the complexity of the modelling choices and ad-hoc nature of the heating distributions used in the two models discussed.Further work should be done on two fronts to confront these problems.First, the source of uncertainties in modelling should be investigated as the potential source for the discrepancies noted in this work.Secondly, models with more self-consistent shock physics should be developed to better approach this problem.
If indeed  Pup's mass-loss rate has not changed, we must then explain why there are observed increases in the line fluxes that are not due to changes in Chandra's detectors.Our spatial parameter ℓ 0 provides such an explanation: the embedded wind shocks are occurring at farther radii now than before.The further out a shock occurs, the more photons can escape since they travel through less material.The increase in ℓ 0 can be traced back to the surface where the lag between the boundary variations is getting longer, giving the slower gas a longer head start before the fast gas is launched.So if the shocks are seeded by boundary fluctuations, this would suggest the timescale of the fluctuations has gotten longer.That any such change is possible on a 20-year timescale is already interesting, and a challenge to interior models, but some sort of stellar cycle seems potentially relevant.Future monitoring of  Pup in the X-ray band is needed to test the apparent long-term variations.LGL|/((LQ + LGL)/2) Relative difference between the general purpose quadrature  Q and Gauss-Laguerre quadrature  GL integration methods available in the VBC model codes.
relative difference between using scipy's general purpose quadrature algorithm  Q and the Gauss-Laguerre algorithm  GL for  = 20 in Figure B1.At the greatest, the difference between these two methods are of order 1.7 per cent, far below any uncertainties from statistical sources.Both of these methods are included in the available files for implementing the VBC model on personal distributions of pyxspec.

APPENDIX C: VBC PARAMETER DISTRIBUTION EXAMPLE APPENDIX D: FIT PLOTS
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Year-to-year X-ray flux of  Pup measured by XMM-Newton's EPIC detector.Vertical lines denote the Cycle 1 Chandra observation and the start of the Cycle 19 campaign (see Table2for dates).Values were retrieved from 4XMM-DR13(Webb et al. 2020).The reported errors in the catalogue are too small to visualize as vertical bars so are plotted in the sub-plot.
), that this filling factor is a power law  () ∝  − .Calculations of  Cohen as a function of radius are shown in the bottom plot of Figure2for a range of values

Figure 2 .
Figure 2. Plot of the line emissivities for the VBC model (top panel) and Cohen model (bottom panel; Owocki & Cohen 2001), both starting at a radius of 1.5 * .Note that for the VBC model, the equivalent of the emissivity function Equation (12) is describing emission over a distance ℓ 0 , not a local emissivity.Both  VBC and  Cohen are normalized at  = 1.5 * .

Figure 3 .
Figure3.Plot of best-fit  * values for the four analyses discussed.Error bars represent 68 per cent uncertainty.Note that we have applied constant offsets to the different result groups to make the error bars more visible.The plotted wavelengths correspond to  p for the blue dots,  p + 0.1 Å for the orange squares,  p + 0.2 Å for the green up triangles, and  p + 0.3 Å for the red down triangles.Data points that do not have a corresponding group of four are due to that line feature not being fit in some of the works considered here.

Figure C1 .Figure D1 .Figure D2 .Figure D4 .Figure D5 .
Figure C1.Example parameter-pair confidence contours for the VBC model fit to Mg xii H  with  ∞,  frozen.Dashed vertical lines in the marginal distributions are the 68 per cent levels for that parameter.Parameter-pair contours are representative of model behavior for most lines fit.

Table 1 .
Archival values of  Pup's mass-loss rate
a CXO Proposal Cycle 1 observation.b CXO Proposal Cycle 19 observation.

Table 3 .
VBC model fitting results for Cycles 1 and 19.

Table 4 .
Plot of best-fit spatial parameters values for the four analyses discussed against the temperatures of maximum emissivity taken from the APED.The VBC model's parameter is plotted as  * + ℓ 0 in units of stellar radii.Error bars represent 68 per cent uncertainty.Points with large errors are not plotted in entirety for readability and are represented by error bars with arrows.Pearson correlation coefficients.
a Correlation coefficient when outlier values from Si xiv H  and Ne ix He  are included.b Correlation coefficient without outlier values from Si xiv H  and Ne ix He .