Deriving X-Ray Line Profiles for Massive-star Winds from Momentum-conserving Dynamical Working Surface Solutions

We present a general procedure for deriving a line-profile model for massive-star X-ray spectra that captures the dynamics of the wind more directly. The basis of the model is the analytic solution to the problem of variable jets in Herbig–Haro objects given by Cantó et al. In deriving our model, we generalize this jet solution to include flows with a prescribed nonzero acceleration for the context of radiatively driven winds. We provide example line profiles generated from our model for the case of sinusoidal velocity and mass-ejection variations. The example profiles show the expected shape of massive-star X-ray emission lines, as well as interesting but complicated trends with the model parameters. This establishes the possibility that observed X-rays could be a result of temporal variations seeded at the wind base, rather than purely generated intrinsically within the wind volume, and can be described via a quantitative language that connects with the physical attributes of those variations, consistently with the downstream momentum-conserving nature of radiatively cooled shocked radial flows.

1. INTRODUCTION X-ray generating dynamics in hypersonic massive-star winds have been studied using simple heuristic models, connecting the salient features of the dynamics to their consequences on observable line profile shapes and flux ratios (Ignace 2001;Owocki & Cohen 2001;Gunderson et al. 2022).The features of these past heuristic models connect with specific questions that those approaches are effective at answering.Sometimes these features included a local filling factor of hot gas (Owocki & Cohen 2001), while others focused on the spatial probability distribution of where the fast gas gets shocked (Gunderson et al. 2022).What has not been explored yet is whether a model can be derived such that the line shapes can be analyzed in a way that is dynamically consistent with the momentum conserving elements of the shocked gas sandwiched between forward and reverse shocks.In this paper, we derive a method of modelling the hypersonic dynamics as consequences of wind variations seeded at the lower boundary, in an effort to test if such a dynamically consistent model can succeed at generating the observed characteristics of the line profiles.
The "problem" with previous models is that they are heuristic by necessity.The dynamics of massive star winds are extremely complex to model in totality.Part of the complexity is the acceleration of the gas due to line-driving (Puls et al. 1996;Owocki 2004;Puls et al. 2008) that introduces instabilities such as the linedeshadowing instability (LDI; Owocki & Rybicki 1984).In addition, parameterizing the wind shocks requires deciding if the shocks originate very near the surface due to some variable boundary condition (e.g.VBC; Gunderson et al. 2022Gunderson et al. , 2024)), or if they require a stand-off distance as might be seen from the LDI (Owocki & Cohen 2001).
Simulations can address these issues in principle, but in practice the line driving requires the solving of radiation transfer equations while shocks require a fine grid in space and time, all in 3D.So any simulation requires compromises about what aspects of the physics to idealize (Puls et al. 2020;Moens et al. 2022).Even with the best compromises, extracting X-ray luminosities from the shocks is another difficulty in and of itself (e.g., Antokhin et al. 2004).
Given these difficulties, we take the approach of creating a model that is as simple and idealized as possible, yet still respects the basic momentum conserving nature of radiatively cooled hypersonic shocks.The goal of this model is to provide a natural language that allows for a more direct connection between wind dynamics and observables of the line profiles.
The foundation of this model is the analytic solution to the problem of Herbig-Haro (HH) jets given by Cantó et al. (2000).An HH jet is the collimated ejection of material from a young star that emits radiation from internal shocks, dubbed "working surfaces" (Reipurth & Bally 2001).For readers unacquainted with the working surface term, it was inducted into the study of HH jets by Mundt & Fried (1983) and Meaburn & Dyson (1987) to describe the shock front within the jet; or the "piston" if it is not a shock.
The Cantó et al. (2000) solution makes one very simplifying assumption: the jet is free-streaming, i.e., constant velocity.The wind of a hot star is of course not constant, and is often described by a β-velocity law where v ∞ is the wind's terminal velocity and β is a parameter controlling how fast the wind accelerates.So to use Cantó et al. (2000)'s model, we will need to generalize it to non-constant velocities, which is the majority of the modelling work done in this paper.
We will concern ourselves with non-constant velocities of the form v(r, τ ) = v 0 (τ )f (r), where v 0 (τ ) is variable in some chosen manner and f (r) is a unitless function describing the acceleration.By considering a separable velocity in space and time, we aim to capture the dynamics of the VBC theory from Gunderson et al. (2022).Specifically, wind shocks are seeded by variations at the wind base, and anything that happens spatially with radius traces back to the time-dependent variation that is seeded at the base.The spatial acceleration is thus similar for all gas parcels, subject only to an overall scale difference controlled by when they were emitted.
Building a line profile model for massive star winds from the solution to a variable jet sounds like an odd choice at first, but the variability of the wind provides an easy translation.To see how, let's fix our frame-ofreference to that of a radially directed cone extending from the base of the wind out to its termination.At the base of this cone, the star is assumed to be undergoing some kind of variability that will cause the wind to be variable in its speed and mass-loading.Inside the cone, the wind follows a prescribed acceleration law, until it encounters slower gas, at which point it joins a working surface (WS) piling up material between the two interacting flows.The gas is then assumed to radiatively cool, creating a narrow pancake-shaped shocked region of hot gas between the two shocks, one a forward shock at the front of the pancake of hot gas, and the other a reverse shock at its rear.The reverse shock decelerates the upstream fast wind to the speed inside the WS, and the forward shock accelerates the downstream slower wind to the speed of the overtaking WS.
The shock itself is the only high density gas needed to be considered, and is the source clumping within the wind.Note however the density of the gas does not need to be determined because efficient radiative cooling is assumed.The assumption of rapid radiative cooling of the hot gas also precludes the need to consider more complex wind geometries where high pressure in the shock region creates lateral blowouts from the otherwise predominantly radial flow.Whatever geometry is assumed, this picture is equivalent to a variable HH jet, so a massive star wind can be thought of as a collection of these "jets" spherically distributed around the star's surface.This is illustrated in Figure 1 where we give a (not-to-scale) schematic of an HH jet versus our "jets" in a massive star wind.
This paper is organized as follows.In § 2, we derive a solution for the dynamics and luminosity of a WS formed in a time-dependent velocity field in the same manner as Cantó et al. (2000).This is done to show that this generalization is possible and the algorithm for which later sections will use.In § 3, we give an example of a time-dependent solution by adding a time-dependent component to a case covered by Cantó et al. (2000).In § 4, we derive a solution for the dynamics and luminosity of a WS formed in a velocity field parameterized by position to connect our theory to the spatial shock distribution.In § 5 we detail how to turn the positiondependent WS luminoity into a line profile.In § 6, we give an example of a position-dependent case using the usual β = 1 velocity law, detailing both the WS luminosity and line profiles.Finally, we discuss our conclusions and future considerations in § 7. The most frequently used variables in our derivation are given in Table 1 for ease of reference.

A TIME-DEPENDENT WORKING SURFACE SOLUTION
To show that it is possible mathematically to generalize the solution by Cantó et al. (2000), and that it provides physically reasonable results, we will first derive the case of a time-dependent velocity.The existence of an accelerating jet can be debated since they are ob-Figure 1.Schematic of the difference between an HH jet from a young stellar object (left) and the "jets" from a massive star (right).In both objects, the faster material is the cyan line, the slower material is the purple line, and the WS is the thick red line.Note that neither drawing is to proper scale.The massive star also is representative of deep in the wind where expansion is negligible; a more detailed drawing would show the "jets" expanding as well.
served to flow with a constant velocity, but this exercise will be done to demonstrate the algorithm that will be followed for the more applicable case of a positiondependent velocity.
We will start with the kinematic statement of a nonconstant velocity where x is the position of the gas parcel, t is the time now, and τ is the time of ejection.As written, the velocity v(t ′ , τ ) is allowed to take any form; a more specific example will be done later, but in general one must define their velocity based on their system.As such, x here is meant as the general position.
The computation of the total mass in the WS does not change for accelerating flows, so where ṁ ≡ dm/dτ .What will change though, is the velocity of the WS, which now must have an explicit time dependence to it since the center of mass will be subject to the acceleration as well.This equation expresses the radial momen-tum conservation of the gas that enters the radiatively cooling WS, which is the heart of the model.The WS will form at This location will correspond to a time t c that must satisfy the relation where τ c is the ejection time of material that first enters the WS, satisfying the minimization condition If the velocity is periodic, then this minimization is taken over the first period.Equations ( 5) -( 7) are an implicit generalization based on the constant velocity equations derived by Raga et al. (1990).
Beyond the starting point, the position of the WS corresponds to the center of mass Since we assume a thin WS, the location of the unshocked wind parcels that are just about to enter the WS through its front and back must be the same x ws (t) = at all times.Thus defines the equations of motion for a WS formed from an accelerating flow.Just as prescribed in Cantó et al. (2000), one would first define the time-dependent velocity to then derive the velocity and position of the WS from Equations (3), (4), and (8).Then, assuming one wishes to have the system described in parametric form with either τ 1,2 , Equations ( 9) and ( 10) can be used to put the position and current time into other known quantities.
Now since we are interested in creating a line profile eventually, we need to propagate these through into the luminosity where ∆E = E 0 − E ws and are the energies before and after the shock occurs respectively.Differentiating these expressions and tracking the time dependencies results in a luminosity of In other words, the time-dependent version of the model Cantó et al. (2000) created amounts to simply adding an additional term to the luminosity of the WS that accounts for the acceleration of gas just entering the WS.
To check that we have the correct expressions, we can apply a constant velocity v(t, τ ) = v(τ ).This significantly simplifies the equations above as all integrals turn into multiplications of time differences τ 2 − τ 1 , but of special importance is the extra term in the luminosity.This term vanishes in the limit of time-independent velocities since the partial time derivative will be zero for constant velocities.So, in this limit we recover the correct expressions.To the extent that our assumptions prove useful, we have derived a time-dependent version of the HH jet solution from Cantó et al. (2000) and can apply it to a specific example.

EXAMPLE 1: DECAYING EXPONENTIAL AND SINUSOIDAL VARIABILITY
As an academic exercise to demonstrate the machinery derived above, consider the specific example of a multiplicative velocity with a sinusoidal velocity variability that accelerates with a decaying exponential and a constant mass ejection rate: These function were chosen for comparison against the first example of Cantó et al. (2000) as the exponential term is like a first order "correction." In this case, our choice of f (t, τ ) does not allow for analytic expressions for the location of the first WS, so we will not discuss that specific portion of this analysis.Moving immediately to following the steps laid out in Equations ( 3) -( 14), we can calculate all the respective quantities to describe the WS.The velocity and position, for example, are where 23) There are a few things to note about these equations.First, the velocity in Equation ( 19) is similar to what was found by Cantó et al. (2000) but with an extra term accounting for the assumed acceleration of the gas just entering the WS, which is not surprising given that we have chosen a velocity similar to those in the reference work.Secondly, while our velocity choice does not allow for a linear expression for the WS's position, Equation ( 20) can be written in condensed form as where which is of the same form as the WS position in the constant velocity case.
To fully solve this system, we need to choose one of our parameters as the parametric variable.We can start by eliminating time by using Equations ( 9) and ( 10), which results in where W (z) is the Lambert W -function, and for the sake of readability we define v 0i ≡ v 0 (τ i ).Due to the length of Equation ( 27), we will not be directly substituting it into any of the equations that follow, but there will be an inherent assumption that it has been inserted for the purposes of plotting.Given Equation ( 27), we can also reinsert it into, say, Equation ( 10) to get another statement of the WS position: meaning one has to solve for τ 1 , which then gives you the rest of the system.This would not yield solutions in closed form given our combination of sinusoidal and exponential functions, but it does follow the prescription outlined in the same way as Cantó et al. (2000).
We can provide an approximate form of the WS quantities if we change our parametric variables to and assume that δv w /v w < 1.In this limit, the average ejection time approaches τ ≈ π/ω, which means Applying these quantities to Equations ( 16), ( 17), ( 19), (28), and (29) simplifies their expressions to where The WS luminosity can also be written in an approximate form under this limit.If we write Equation ( 14) as L ws = L 0 + L accel , we can explore each of these quantities systematically.Starting with the former term, we find a similar form to that of Cantó et al. (2000), i.e., where is the luminosity of the constant velocity case.The latter term L accel was stated before to account for the acceleration on the gas just about to enter the WS.This manifests as the non-constant components in By inspection one can note that L accel ≪ L 0 for all ∆τ .So to first order, the luminosity of the WS can be taken as simply In other words, a WS formed from gas that experiences an exponentially decaying acceleration, in the limit of δv w /v w < 1, generates shocks observable to be similar to that of the constant velocity case.This luminosity is plotted in Figure 2. Focusing first on the top plot of Figure 2, we see that the exponential acceleration causes the WS to take longer to form, and as t 0 increases the emission is over a longer period of time.Of course the more that t 0 is increased the weaker the luminosity since the WS would have lower shock velocities.The time at which the emission starts does appear to only be weakly dependent on t 0 .Moving from the constant case to t 0 = 1 has the WS start at t ∼ 6 instead of t ∼ 5 (in the arbitrary time units used for illustration purposes), whereas going to t 0 = 3 has only moved the start to t ∼ 6.5.
The amplitude of the velocity variability δv w has a more outsized effect on the start time of emission over its whole range and for the luminosity's emission trend in general, as can be seen in middle plot of Figure 2.For both the accelerating (solid lines) and constant cases (dashed lines), the shape of the luminosity is essentially the same and all show the same approximate limiting behavior as time goes on.More interesting is what happens as δv w increases.At some point, there is a value for δv w such that both models will have the same luminosity but any larger than that the constant velocity case will emit more.
This behavior is also seen in the periodicity's frequency where for some values of ω we see the constant velocity case have more emission, such as ω = 2 in the bottom plot of Figure 2.This effect is interesting since it can highlight how important the role of acceleration in the strength of a shock is, depending on how weak/strong it is in comparison to the fluctuations of the injection.For example, when the oscillations are slower, the acceleration can cause a larger difference in the gas speeds when the shocks form.When the velocities are constant, the WS simply emits light for a longer timeframe but the peak emission does not change.
An important consideration to note is the peak of the emission.All three panels of Figure 2 show that the peak can increase quickly as parameters change.In some cases, though, this peak may be artificial, such as when t 0 → 0. A more realistic form of the chosen velocity would include a constant to the acceleration, say where 0 ≤ ϵ ≤ 1, to keep the emission from becoming artificially large.Such considerations are beyond the scope of this work, however, and left for possible future applications in studies of HH Jets.

A POSITION-DEPENDENT WORKING SURFACE SOLUTION
While the discussion of the time-dependent solution is interesting and worth further consideration, for the purposes of this paper it is only for the illustrative purpose of showing that is possible to find a solution to the problem of HH Jets with non-constant velocities.Our goal from the outset has been to use this methodology to create a line-profile model for an X-ray line generated within the wind of a massive star.Such a situation is one that requires an acceleration that is position-dependent.The derivation of the equations of motion for such a case are simplified by knowing the generalization procedure.
If the position x is now the dependent variable, then it is the time of flight t − τ of the parcel of gas that will be the independent variable.This can be written explicitly as While we rarely think of kinematics in this orientation, applying a constant velocity gives the expected result.Note that we are still using x as an arbitrary position of interest.
Many of the quantities we give below will be similar to our time-dependent solution, but the approach to the problem differs.The easiest example is how one finds the quantities x c and t c .To find these, we will look at the system formed from material ejected at τ c − dτ and τ c , which looks like This system has a solution when and τ c is still determined by the minimization condition d dτ where it is again taken within the first variability period.The WS velocity will not have any change since Equation (4) is only an integral over τ , but of course we now have an x-dependence Specific care should be taken to the context of the acceleration being applied to the gas in this case.If a multiplicative velocity v = v 0 (τ )f (x) is used, then equation Equation ( 53) will read similarly as v ws (x) = v 0ws (τ 1 , τ 2 )f (x), which says the acceleration is applying to the WS as well.This will be discussed further in § 6 in the context of our chosen velocity.
Since the WS is still assumed to be thin, its time of flight should correspond to the two ejecta as before The problem of a position-dependent velocity for an HH jet is thus the same as the time-dependent case.The equations given in this section give the equations of motion and can be put into parametric form with a chosen parameter.
The luminosity for position-dependent system can be easily derived if we apply a change of variable so that This significantly simplifies the problem of defining the luminosity because we can expect the same form to appear for the rate of energy loss.Thus the WS luminosity for a position-dependent velocity would be where we have noted that dx/dt = v ws since the position we are tracking is that of the WS.We have also kept the derivatives dτ i /x for reasons that will be clear in the examples that follow.
Applying the constant velocity check to these equations does recover the expect equations of motion, including the the acceleration term in the luminosity again going to zero.What we want to highlight is the similarity in the mathematical form of the two approaches to the WS solutions.While such similarities can be expected since we are investigating the same system, the connection between the two methods in the luminosity is remarkably only a change of variables Such a simple connection highlights that the emission of the WS is independent of the parameterization.So while the kinematics are non-standard, the results can be connected to the usual methods.

WORKING SURFACE LINE PROFILE
With the basic description of the expected dynamics of the shocks that produce the X-rays is defined, we can track the emission as the working surface traverses the wind and how that emission is attenuated by the cold, unshocked wind.Such attenuation is characterized at each radius r and emission directions, equally spread over µ = cos θ, by the optical depth D(µ, r), which is why we have generalized the Cantó et al. (2000) solution to position-dependent velocities.Note that the equations given previously we general to any position coordinate x, but the rest of the analysis that follows we will use the radial coordinate r as it is more appropriate for the derivation of the line profile.
Line profiles are a powerful probe of the shock signatures because the hypersonic wind Doppler shifts the line in such a way that the frequency dependence serves as a proxy for the spatial distribution.Hence it is the shape of the X-ray line profile that provides a test of our description of the shock physics.A profile can be generated by integrating the X-ray generation rate L ws (r) over radius and tracking how the emitted energy is concentrated into each dξ frequency bin.The methodology used in the defining of the VBC profile model in Gunderson et al. (2022) will be employed here.In this approach, the WS is assumed to be emitting isotropically in all directions µ = cos θ, and the emitting plasmas are assumed to be distributed spherically symmetrically around the wind.
The isotropic emission simplifies the integration for the line profile by removing the need to integrate around the azimuthal angles.This means that the line profile model will only be the radial integration where ξ = −µv/v ∞ is the line of sight velocity coordinate describing the Doppler shift of each photon and dµ/dξ is the local mapping of the solid angle to that Doppler shift space.It was noted in the introduction that our assumption in the use of a jet model requires some modifications of the input assumptions.One of those assumptions is whether the flow is expected to be collimated or spherically expanding.For the jet solution derived by Cantó et al. (2000), no expansion of any kind was included since jets stay collimated and do not accelerate.The "jets" in the stellar winds we are considering are expected to spherically expand as they flow out radial, and radially expand as they accelerate.The expansion affects the density, which has no independent significance in a purely radiatively cooled model, but will relate to when that assumption breaks down.
Determining when radiative cooling is not a good assumption requires comparing the radiative and adiabatic cooling times.We address that complicated issue with an efficiency factor ε eff (r).The choice in the functional form of ε eff is arbitrary but is expected to fall with radius.For example, a power law form is a natural choice.Here the index q describes the importance of adiabatic cooling.Alternatively, an exponential form ε eff = e −(r−R * )/ℓ0 (62) can provide similar weighting.In this case, ℓ 0 is describing the radial extent of the radiative cooling zone.
Note that these efficiencies are identical to the source descriptions from Owocki & Cohen (2001) and Gunderson et al. (2022Gunderson et al. ( , 2024) ) respectively.In the former, the power-law was used for a filling factor approach describing how much of the wind is heated through shocks at each radius.The latter used an exponential to describe the shocks through a mean-free path characterization.Both of these works choose these forms heuristically to match with the assumed spatial distribution of the shocked gas throughout the wind.In the present work, the adiabatic cooling cutoff will be found to be so important that it is appropriate to compare ε eff to these heuristic source terms.
We can do slightly better than these heuristic choices for the efficiency by directly comparing the cooling rates for radiative and adiabatic processes Here ρ is the density of the gas, mw is the average mass of the ions in the wind, Λ is the cooling function (Wang et al. 2014), k B is Boltzmann's constant, and T is the temperature.Note that for the adiabatic cooling rate (dE/dt) ad we assume a spherically expanding flow.
Both of the cooling rates require some knowledge of the shock strength to determine the temperature T , but if we note that ρ = Ṁ /4πr 2 v, we can estimate the radial dependence of the efficiency as where v ′ (r) = ∂v/∂r.All three efficiencies from Equations ( 61), (62), and (65) will be used in example profiles in § 6.
The optical depth can be derived from where p is the path-length of the photon's emission from the shock at radius r, r ′ is the photon's radial distance from the star, and is a fiducial optical depth describing the X-ray reabsorption (Cohen et al. 2010;Gunderson et al. 2022).Note that the definition of D * uses the average terminal velocity v ∞ , which describes the cool, unshocked wind that responsible for the absorption.
The lower bound of Equation ( 60) is This is the minimum radius r D capable of achieving a ξ Doppler shift or the radius r c where the WS forms.The latter term r c is derived from Equations ( 51) for the given choice of velocity v(r, τ ).The larger of the two is always taken because no X-rays are created below where the WS forms.The minimum Doppler shift radius is dependent on the direction cosine of the photon's emission µ, though more easily written in terms of ξ, and the velocity v(r) used to describe the flow.For the forward direction, the minimum occurs for µ = 1, which corresponds to The backwards direction, which corresponds to photons traversing around the star, must account for occulation.For this case, the minimum angle corresponds to µ = − r 2 − R 2 * /r, so one has to solve 6. EXAMPLE 2: OFFSET β-LAW AND SINUSOIDAL VARIABILITY In this example, we will consider velocity and massloading as it pertains to massive-star winds.Quantities will have appropriate units now, such as wind speeds in terminal velocity units v ∞ and distance in stellar radii R * .We will thus be measuring time in units of t fl = R * /v ∞ , i.e., flow time.
As mentioned before, the β-velocity law in Equation (1) is the ususal parameterization for the velocity field of massive star winds, so this example will use the specific case of β = 1 for the acceleration term.However, a β = 1 velocity introduces a logrithmic divergence in the integrals (e.g., Equation 48) involved.This divergence, physically speaking, is due to the β = 1 wind taking an infinite amount of time to step from R * to R * + dr as the wind launches.To account for this problem, we will instead use an offset-β-law where a = 0.995, or some similar value close to unity.
The differences in the speed of the fast and slow portions of the flow will be described by varying the terminal velocity The average value of this variability is the terminal velocity v ∞ as measured through P Cygni profiles (Puls et al. 1996).The variability amplitude δv ∞ will be restricted to 0 < δv ∞ < v ∞ , but in general is large enough such that given the frequency of the oscillation ω, the kinetic energy thermalizing in the WS is of order 1.We will pair this velocity with a non-constant mass ejection rate ṁ = ṁs + δ ṁs sin(ωτ ). ( The phase shift between v 0 and ṁ is chosen based on the density dependence of line-driving; higher density winds experience less acceleration than those of lower density (Gayley 1995).Thus it is actually the mass variability that causes the wind velocity to vary.Note that ṁ is not the mass-loss rate of the star Ṁ .It is only a small fraction describing the mass-flux along the ray.This velocity and mass variability will produce a WS at a radius of where is the position the WS when subject to only a constant velocity.Using this position in Equation ( 50), we can find the corresponding time of the WS formation Finally, the ejection time of the material that first forms the WS can be calculated from Equation ( 52), which gives Figure 3. Radius of WS formation as a function of the velocity variability amplitude δv∞ for an offset β = 1 flow.
It is interesting to note that Equations ( 76) and ( 77) are the same in both this example and those given by Cantó et al. (2000).This shows that for the case of an acceleration parameterized by position, the time that the WS is dependent on the shape of the variability, such as whether it is sinuisodal, saw-tooth, or some other periodic function.The acceleration will only change the location where it forms.
In Figure 3 we plot Equation ( 74) as a function of δv ∞ to illustrate the important feature of r c asymptoting to r c → R * as δv ∞ → v ∞ .Such behavior is expected because a stalled wind would never leave the surface and immediately shock with the next ejection.Additionally, while Equation ( 74) is analytic for δv ∞ > v ∞ , we will not consider such a case.Our system has a hard boundary at r = R * , so no wind can flow downward at the surface.This presents one limitation of our model since it cannot account for an ejection that stalls above the surface and falls back down.
In this example we will deviate significantly from the previous method of defining τ and ∆τ and instead show how the radius r can be used.This is achieved by differentiating Equations ( 54) and ( 55) with respect to r, forming the system of differential equations where The initial conditions of this system are τ 1 (r c ) = τ c − δτ and τ 2 (r c ) = τ c +δτ , where δτ is a small shift to avoid the singularity when τ 1 = τ 2 .Note that v0 = dv 0 /dτ .This system will have a unique solution as long as r > R * , but it has no analytic expression. 1 A different choice in v 0 and f may permit a solution that is an analytic expression but such choices would need to be determined on a case-by-case basis.
The velocity of the WS is significantly simplified in its computation since f (r) has no τ dependence, giving us where and τ 1,2 are the solutions to Equations ( 78) and ( 79).The velocity can be written more compactly as which implies the WS experiences acceleration just the same as the non-shocked gas.
We would not expect the line driving to apply to the shocked plasma since the ions will become highly ionized.A possible correction to this effect is to remove additional acceleration applied after the working surface forms However, after testing this method, we found that our differential equation system in Equations ( 78) and ( 79) is either highly numerically unstable or no longer has a solution.The exact nature of the numerical problem is not clear, so we elect to keep Equation (84) as the WS velocity for what follows and defer the removal of line driving to future considerations.Turning to the luminosity, inserting Equations ( 72), ( 71), ( 78), (79), and ( 84) into (57) gives us 1 Meaning solutions that can be written with known functions up to and including special functions. where is the variance of the velocity variability v 0 .This luminosity is plotted for various combinations of parameters in Figures 4.An interesting feature of the luminosity is the radial extent of the emission compared to its starting location.For all parameter combinations, there is only a small range of starting radii (illustrated in Figure 3) and the emission continues well past 10R * .Such a narrow range of starting location connects well with previous studies that inferred a small shock radii in the 1.5-2 R * range (Cohen et al. 2014;Owocki & Sundqvist 2018).
The overall extent of the emission, however, speaks to the slow nature of a β = 1 velocity law, which affects the timescale for kinetic energy to enter the working surface and be thermalized.The individual parameters also show interesting trends for the emission in Figure 4.The velocity variability amplitude δv ∞ (top panel) can significantly raise the total amount of the emission since more energy will be available for conversion and the shock speed will increase between the fastest and slowest portions of the flow.The latter reason is also why lowering the frequency ω (middle panel) also increases the total emission, though more subtly.The larger difference between the fastest and lowest flow is dependent on the the interplay between the acceleration and variability.A faster frequency results in the WS forming lower down where there is not as much time for acceleration to have occured, causing the fast and slow flows to be more similar in their velocity at the time of the shock.The last parameter δ ṁs , however, has very little effect on the profiles.It also can change the overall total emission, which is likely due to how the working surface velocity is inversely proportional to the total mass within it, but it is of smaller scale in how much it changes the overall emission.The changes these parameters have in the shape and amount of emission will propogate into the line profile formed by the emission, which we give below.To generate a line profile, we need three more things: an effiency ε eff , optical depth D(r, µ), and the minimum radius r m (ξ).Starting with the first of these, we could simply choose the exponential or power-law from Equations ( 62) and ( 61), but we also want to include the more physically-motivated choice of the cooling rate ratio from Equation (65).Using the offset β = 1 velocity, this efficiency is

Line Profile Plots
In this idealized wind, the radius where adiabatic cooling overtakes radiative cooling occurs at r ≈ 1.7R * .
The next thing we need for generating a line profile is the optical depth of the wind from Equation (66).For this example, we will use the ususal β = 1 velocity without the offset a.The luminosity curves in Figure 4 show that the vast majority of the emission is well beyond r = R * , so even though there is an infinite optical depth at the surface for this velocity, we will not be losing any substantial amount of emission.Given this, the optical depth for a regular β = 1 velocity is where This form of the optical depth was originally derived by Gunderson et al. (2022).Note that in our formulation, µ > 0 corresponds to directions of only increasing radius away from the star while µ < 0 is directions with initial decreasing radius.Finally, the minimum Doppler radius can be calculated from Equations ( 69) and (70).Starting with the forward direction, again, the minimum Doppler radius is and the backwards direction, when also taking into account occultation, is where These equations were originally derived in Gunderson et al. (2022) for a = 1.While the form is not exact to our velocity, such a small difference will not compromise the illustration of our model.
To analyze the derived line profile, we can first look at what the different efficiencies change in the profile shape.This is shown in Figure 5, where in the left panel we plot the profiles using Equations 62, (61), and (88).We have also included the case of ε eff = 1, i.e., the wind is equally efficient at cooling both radiatively and adiabatically at all radii, for comparison.This can be considered equivalent to assuming no expansion within our "jets", so it is also a check of the self-consistency in our assumptions.
The first thing to note is that for all non-constant efficiencies, the generated profile is Gaussian-like.Such a shape is necessary since the lines in observed massive star spectra are Gaussian-like, but it is also significant because the model was not designed purposely for this shape to occur.The Gussian-like shapes occur naturally, insofar as can be claimed based our minimal use of heuristics, from the kinematics considered.We can thus point to the exact source of the Gaussian-like shape: the increasing efficiency of adiabatic expansion as the cooling mechanism over radiative emission as the hot gas flows outward.The non-Gaussian profile generated from ε eff = 1 highlights this point because there is significant emission near the terminal velocity, i.e., coming from far far radii.
To explore this idea further, the right plot of Figure 5 shows how the shape of the profile changes as the characteristic radiative cooling efficiency length ℓ 0 from the exponential efficiency increases.For small ℓ 0 , the profiles are Gaussian-like, but if the region of radiative cooling is allowed to stay significant for longer, the profiles loses this shape.The balance between the two is not broad, however, because even at ℓ 0 = 2.5R * the profile becomes less Gaussian-like.Owocki & Cohen (2001) discussed how flat top profiles would occur for profiles generated by an optically thin, constant velocity wind, but in their derived model, these come from emission that starts at larger radii.Our model shows similar profiles would appear even from emission starting at deep radii.Thus the shape of a massive star's X-ray line profile appears to be caused by the combination of both optical depth and radiative cooling only being significant for a small region deep in the region.
In Figure 6, we show how varying the different parameters of the model affect the profile shape.Starting with the fiducial optical depth in the top left plot, increasing D * causes the profile to become more blue-shifted and asymmetric.This is expected behavior for increasing the total optical depth that the line would experience and is in line with how this parameter affects other model profiles (Owocki & Cohen 2001;Gunderson et al. 2022).The other three plots also describe properties of the WS, so we can connect their changes to those shown in the luminosity curves.The top right plot shows that as δ ṁs increases, the total flux is reduced with no changes in the overall shape.A similar property was shown in the bottom plot of Figure 4, implying that the total mass in the emitting hot gas will not change spectral diagnostics like width and centroid.
The velocity variability amplitude δv ∞ and frequency ω, on the other hand, have dramatic and similar effects on the profile shape, which the bottom plots of Figure 6 show respectively.When either parameter is small, the profiles take on a "shark fin"-like shape that has been observed in WR 6 (Ignace 2001;Huenemoerder et al. 2015).The luminosity curves for these parameter values (blue curves in the top and middle plots of Figure 4) show that the shark-fin shape occurs because the majority of the emission occurs well beyond the radius where adiabatic expansion becomes significant.For other plotted parameter combinations, the Gaussian-like shape is produced again.However, there does appears to be a limit to the profile shape as either parameter increases.The first time δv ∞ or ω is doubled, the profile shape changes dramatically.The second time produces a much less significant change of being slightly narrower.Further increases would produce even smaller narrowing that would become imperceptible in difference.Looking again at the luminosity curves, this would be due to the similar radial profile of the luminosity curves, particularly as increasing either parameter causes the curve to become more sharply peaked and start at lower radii.As the luminosity becomes more sharply peaked, it concentrated further in the core of the line profile but the WS can only squeeze so tightly under the constraints of our assumptions.
The importance of the luminosity's radial dependence to the line profile's shape is a surprising connection.Our model implies that to have a Gaussian-like profile from a β = 1 velocity wind, the embedded shocks must be occuring deep within the radiative region near the surface, such as Figure 3 shows for the WS formation radius.Previous studies have inferred such deep radii for the shock radii from He-like f ir line ratios (Cassinelli et al. 2001;Cohen et al. 2022).Wind profile models have also pointed to this as Cohen et al. (2014) found r ≈ 1.5 − 2 R * to hold across many lines and stars.Similarly, Gunderson et al. (2024) found that ζ Pup's X-rays can be well-fit with the VBC model, which assumes shocks can occur the moment after launching.So the produced Gaussians coming from a WS that formed near the surface is within expectations.The formation radius is not the only piece, though, as we had already discussed how the radial dependence in the radiative cooling efficiency is also necessary.This can also be put into context of previously models, specifically the models of Owocki & Cohen (2001) and Gunderson et al. (2022Gunderson et al. ( , 2024)).Both of these models can be described as having shell-like emission as the source of the X-rays, though they do so in different ways.In the former, a power-law filling factor r −q describes how much hot gas should be at any given radius after their onset radius.The latter assumes that the likelyhood of shock occurrence decays exponentially with radius e −(r−R * )/ℓ0 , using ℓ 0 to describe the mean-free path before a shock occurs.By assuming shell-like emission, these models have effectively hard-coded an assumed efficiency of radiative cooling into the model, which is why we chose to use these exact functional forms as a point of comparison.The difference now is that our model provides the insight that the shell-like emission is due to decreasing cooling efficiencies crossing with rising shock heating in our dynamical model.
The profiles discussed here are but one example of what can be generated from our modelling approach.In Appendix A, we provide general expressions for when β is left unspecified.

CONCLUSIONS
We derive a more dynamically consistent model for connecting observed massive-star wind X-ray line profile shapes with potential boundary variations that could be responsible for their generation.The goal was to test whether such boundary variations are generally consistent with, and can be constrained by, observations.The derived model is a modification of the HH jet WS solution derived by Cantó et al. (2000), for which we also generalized to non-constant velocities parameterized by either time or position.
For the specific case of sinusoidal velocity and massejection with a β = 1 acceleration, the resulting spatial heating distribution, described by L WS (r), is reminiscent of previous more heuristic models.The profiles generated by this distribution have appropriate Gaussianlike shape that is both blue-shifted and skewed due to the absorption.Additionally, certain combination of parameters, particularly slower variability frequencies, can generate profiles with non-Gaussian-like shapes reminiscent of the so-called shark-fin profiles.
The free parameters of the model, optical depth D * , variability frequency ω, velocity variability amplitude δv ∞ , and mass-ejection variability amplitude δ ṁs , provide a more direct connection between the dynamical origin of the wind shocks and the shape features of a line profile that contain information about the dynamics occurring up-to and into the shock.Moreover, they can inform us of the nature of the boundary variations that could possibly be the source of the wind clumping and shocks.This offers promise that within the context of models where X-rays are generated in hypersonic wind perturbations that are seeded by variations at the lower boundary, the observed line profile shapes and flux ratios can constrain the nature of the boundary variation.
Our model also reveals that the source of the Gaussian-like shape in the profiles is due to adiabatic expansion dominating as the cooling process at large radii.The core of a line's emission thus comes from a small radial shell around the star, not unlike the assumed geometry of emission from previous heuristic models.This shell like emission is achieved using the radiative cooling efficiency ε eff , which can be connected to a physically motivated cooling rate ratio.More detailed line analysis will want to include details left out in the cooling rate ratio we have used, though.
In general, the model needs to be implemented within spectral fitting software such as xspec (Arnaud 1996), ciao (Fruscione et al. 2006), or isis (Houck & Denicola 2000) and fit to a variety of single massive stars data.In particular, WR 6 is a candidate of interest to investigate since our model can naturally produce shark-fin-like profiles (Huenemoerder et al. 2015).If our model can fit both the Gaussian-like lines of OB and shark-fin-like lines of WR stars, it would provide evidence that the Xray production mechanism is the same in both spectral types, as posited by Gayley (2016).

Figure 5 .
Figure5.Peak-normalized line profiles generated from a WS defined by Equation (86) using different efficincies.In the left panel, the legend refers to what form of ε eff is used: no radial dependence (blue), a power-law from Equation (61; red), an exponential from Equation (62; green), and the rate ratio from Equation (88; orange).In the right panel, an exponential is used with different ℓ0 values to demonstrate the dependence the Guassian-like shape of the profiles have with the size of the region of efficient radiative cooling.Other parameters values used were: D * = 1, v∞ = 1, δv∞ = 0.2v∞, ṁs = 1, δ ṁs = 0.1 ṁs, ω = π/2 t −1 fl , and a = 0.995.

Figure 6 .
Figure 6.Peak-normalized line profiles generated from a WS defined by Equation (86) as parameters are varied.All generated profiles used Equation (88) for the efficiency.Unless stated otherwise in plot legends, the following parameter values were used: D * = 1, v∞ = 1, δv∞ = 0.2v∞, ṁs = 1, δ ṁs = 0.1 ṁs, ω = π/2 t −1 fl , and a = 0.995.Note that in the top right figure where δ ṁs is varied the profiles are peak-normalized to δ ṁs because no changes are made to the shape; only the total flux in the line is changed.