{'fraction_non_alphanumeric': 0.03067987700717458, 'fraction_numerical': 0.012606764605398018, 'mean_word_length': 4.586068702290076, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time Dependent Radiation Driven Winds

We study temporal variability of radiation driven winds using one dimensional, time dependent simulations and an extension of the classic theory of line driven winds developed by Castor Abbott and Klein. We drive the wind with a sinusoidally varying radiation field and find that after a relaxation time, determined by the propagation time for waves to move out of the acceleration zone of the wind, the solution settles into a periodic state. Winds driven at frequencies much higher than the dynamical frequency behave like stationary winds with time averaged radiation flux whereas winds driven at much lower frequencies oscillate between the high and low flux stationary states. Most interestingly, we find a resonance frequency near the dynamical frequency which results in velocity being enhanced or suppressed by a factor comparable to the amplitude of the flux variation. Whether the velocity is enhanced or suppressed depends on the relative phase between the radiation and the dynamical variables. These results suggest that a time-varying radiation source can induce density and velocity perturbations in the acceleration zones of line driven winds.


INTRODUCTION
Line driving is a viable mechanism for launching winds in massive stars, cataclysmic variables (CVs) and active galactic nuclei (AGN). Simulations have shown that in all but the simplest, spherically symmetric geometries (Castor, Abbott & Klein 1975, hereafter CAK), line driven winds exhibit complex, time-dependent behaviour. Perturbations at the base of the wind may grow due to the line deshadowing instability (LDI), leading to a clumpy, non-stationary flow (Owocki, Castor & Rybicki 1988, hereafter OCR88;Sundqvist, Owocki & Puls 2018 hereafter SOP18). In the case of a disc wind, differences in geometry between the flat, equatorial source of photons and matter and the spherically symmetric gravitational potential lead to the growth of density structures at the base of the wind in both 2D axisymmetric simulations (Proga, Stone & Drew 1998;Proga, Stone & Drew 1999) and in 3D simulations (Dyda & Proga 2017 hereafter DP17;Dyda & Proga 2018, hereafter DP18). Line driving is also sensitive to the ionization state of the gas, which is dependent on the aformentioned density features in the flow, which can act to shield the outermost gas from ionizing radiation (Proga, Stone & Kallman 2004;Matthews et al. 2016).
In the case of massive stars, one can correlate changes sdyda@physics.unlv.edu in apparent magnitude, which probe stellar luminosity, to changes in emission lines, which probe the outflows. This has been used to show increases in the luminosity of ζ Puppis lead to increases of the mass flux of its stellar wind (see Ramiaramanantsoa et al 2018 and refrences therein).
In the case of AGN, studies of BAL QSOs have shown that approximately 5% of systems exhibit disappearing BALs on time-scales of months to years (see De Cicco et al. 2018 and refrences therein). Though such variations are typically explained via changes in covering fraction or ionization state of the absorbers, the dynamics of the absorbers or clouds may be highly sensitive to intrinsic source variability. By comparing the time-variability of spectral lines relative to the continuum using reverberation mapping techniques (for a review see for example Peterson 2001), one can model the distribution and velocity of gas in the system. A key question is how to distinguish between changes in flux due to intrinsic source variations and due to variations in the column density of interposing gas between source and observer. By self consistently modelling the source variation and its effects on the outflow, it may be possible to break the degeneracy of observed flux variations due to intrinsic source variability and due to interposing gas clouds.
To develop a physical intuition for the effects of temporal source variability on winds, we explore time-dependent models of CAK type 1D spherical line driven winds. This approach is motivated by work by Waters and Proga (2016) that showed that a time varying radiation source can accelerate AGN clouds more efficiently than a stationary source. Other groups have shown that properly treating the line driving may require a more accurate description of the radiation transfer to capture the LDI (SOP18), multiple resonance surfaces in the wind (Springmann 1994) and correctly inferring line opacities from ionization states of the gas (Higginbottom et al. 2014). As a first step to understanding time variability, we use the simplest model for radiation transfer and work in the Sobolev approximation, but vary the driving flux. This will allow us in future work to study time varying models in 2D, and eventually 3D, and compare to our previous work on disc winds, DP17 and DP18, using the Sobolev approximation.
The CAK solution is a particularly simple case to study because winds are stationary and fully characterized by their physical parameters at the critical point, where the monotonicaly increasing velocity satisfies the relation, dv/dr|c = vc/rc, where v and r are the velocity and radial position, respectively and the c subscript denotes the critical point. Thus the solution is characterized by a single physical scale rc and dynamical time scale τc = rc/vc.
For time-dependent winds we expect different solutions depending on the ratio between the dynamical time and the period of the source TS. Driving the wind on long time scales relative to the dynamical time, TS τc, the solution oscillates between the high and low flux stationary states. Driving the wind on short time scales, TS τc, the solution behaves like the mean flux stationary solution. When time scales are comparable TS ≈ τc, velocity perturbations with amplitude comparable to variation in driving radiation are induced at the base of the wind. Surprisingly, these velocity perturbations are greater (smaller) than the velocity of stationary solutions, when the phase between the driving radiation and the velocity perturbations are negative (positive).
In Section 2, we describe our simulation setup and our model for intrinsic source variation. In Section 3, we describe our results for a fiducial run and the dependence of winds on driving source frequency. In Section 4 we discuss some implications for observations of massive stars, CVs and AGN. We conclude in Section 5 where we discuss some of the limitations of our approach and possible improvements in future work.

SIMULATION SETUP
Consider a point source of radiation, surrounded by a spherically symmetric, isothermal shell of gas that is optically thin to the continuum radiation. The basic equations for single fluid hydrodynamics are where ρ is the fluid density, v the velocity, P a diagonal tensor with components P the gas pressure, Φ = −GM/r is the gravitational potential of the central object and E = 1/2ρ|v| 2 + E is the energy where E = P/(γ − 1) is the internal energy. The radiation force is F rad and is described below. We take an equation of state P = ρ γ where γ = 1.01. The isothermal sound speed is a 2 = P/ρ and the adiabatic sound speed is c 2 s = γa 2 . We compute the temperature from the internal energy via T = (γ −1)Eµmp/ρk b where µ = 0.6 is the mean molecular weight, k b is the Boltzmann constant and mp is the proton mass.
We model the radiation force as a sum of contributions from electron scattering and line driving where the force due to electron scattering is with Γ * = L * σe/4πcGM the stellar Eddington parameter, L * the stellar luminosity and σe the Thompson cross section. The radiation force due to lines is where the force multiplier uses the CAK parametrization with the optical depth paramter and line driving parameters k = 0.2 and α = 0.6, consistent with a gas thermal velocity v th = 4.2 × 10 5 cm/s. To be able to discuss applications of our results to various objects (from O stars and CVs to AGN), we perform our simulations in dimensionless units and later convert to scalings for the objects of interest. The central object has a mass M = M * and the simulation region extends from r * < r < 20 r * . This dynamical range allows the wind to fully accelerate and reach its terminal velocity. The velocity scale is then v * = GM * /r * and the time scale t * = r 3 * /GM * . We use a logarithmically spaced grid of Nr = 128 points and a scale factor ar = 1.08 that defines the grid spacing recursively via drn+1 = ardrn. We choose a hydrodynamic escape parameter HEP = GM * /r * c 2 s = 10 3 , which ensures that thermal driving is negligible (Stone & Proga 2009;. At the inner boundary, we impose outflow boundary conditions on v and E while keeping the density fixed at ρ * = 10 −10 g/cm 3 in the first active zone. This density ensures that we are sufficiently resolving the atmosphere at the base of the wind, but also providing sufficient mass to launch a wind given our choice of Eddington parameter. At the outer boundary, we impose outflow boundary conditions on ρ, v and E.
We explore cases where the Eddington parameter varies as a function of time, oscillating between a high luminosity state Γ+ = 0.12 and a low luminosity state Γ− = 0.08. We consider a model where where the amplitude A = (Γ+ − Γ−)/2Γ0 and t0 is chosen so the initial wind reaches a steady state. We explore cases where 1.16 × 10 −2 t * TS 1.16 × 10 2 t * . We use the 0, + and -subscripts to refer to the steady state winds with fiducial, high and low luminosity respectively. The dynamics of CAK line driven winds are determined by the mass flux at the critical point (denoted by the subscript c) where dv/dr |c = vc/rc (see for example Lamers and Cassineli 1999). Our fiducial solution with Γ * = Γ0 allows us to determine the characteristic time-scale of the flow, τc = rc/vc = 1.57 r * /0.92 v * = 1.71 t * . The characteristic angular frequency is then ωc = 2π/τc = 3.67 t −1 * .

RESULTS
We first determined steady state solutions for Eddington parameters Γ−, Γ0 and Γ+ by running simulations for t ≈ 200 t * . The fiducial run steady state was used as the initial condition for the time varying model (see eq. 7). We find that after a transient period set by the dynamical time, the time-varying solution settles into a quasi-stationary state where the dynamical variables oscillate sinusoidally in time about the fiducial solution where ∆X and φX are the amplitude and phase shift of the density, velocity and mass flux. Crucially, both the amplitude and phase shift are functions of the driving frequency, which we further describe below.
In Fig 1, we plot the Eddington parameter Γ * , the density ρ, outflow velocity v and mass fluxṀ as functions of time at r = 2r * (light blue), r = 5r * (blue) and r = 20r * (dark blue) for the typical case TS = 1.16t * = 0.68τc. For reference, we plot the high (upper dashed line) and low (lower dashed line) flux solutions. We note the critical point of the fiducial solution, rc = 1.57 r * , so these points are all exterior to the critical point. We plot times near t0 around which the source begins to oscillate.
The transient period lasts approximately for a time τc, during which the solution becomes periodic, and dynamical variables begin to oscillate with angular frequency of the source, ω = ωS. Oscillations begin at the innermost radii r = 2 r * (light blue) and propagate outwards to r = 5 r * (blue) and to the outermost radius r = 20 r * (dark blue). The phase difference between peaks in the Eddington parameter and in the dynamical variables increases as we go out in radius. In this case, the amplitude of density oscillations are bounded by the high and low flux stationary solutions, as is the mass flux. On the other hand, the velocity oscillations are not bounded by the high and low flux solutions. At r = 2 r * , the amplitude of these velocity oscillations is a factor of a few times the velocity difference in the high and low flux solution. This difference decreases as we move out in radii, with the oscillation amplitude damping with radius.
The magnitude of the oscillations ∆ρ, ∆v, ∆Ṁ as well as their phase shifts φρ, φv and φṀ depend non-trivially on ωS. In the upper panels of Fig. 2, we plot the maximum and minimum values for the density ρ (triangles), velocity v (circles) and mass fluxṀ (diamonds), normalized to the steady state values ρ0, v0 andṀ0 respectively for r = 2 r * (light blue) and r = 20 r * (dark blue). The dashed black lines indicate the high and low luminosity steady state solutions. In the lower panels we plot ∆X/(X+ − X−), the difference between the maximum and minimum amplitudes for each dynamical variable, normalized to the difference between the high and low luminosity solutions.
For small driving frequency ωS ωc, the wind oscillates between the low and high luminosity steady state solutions. This is most apparent in the density and mass flux profiles. At large driving frequencies ωS ωc, the wind resembles the fiducial steady state wind. Near the critical frequency wS ∼ ωc, we see a smooth transition between these two regimes. Most interestingly, the velocity amplitudes exhibit a resonance effect whereby the difference between the maxima and minima grows beyond the velocity differences between the high/low luminosity states. To understand why velocity is enhanced, consider a stationary wind solution, that begins to receive less (more) driving flux. A parcel of gas just beyond the critical point in such a flow will now be more (less) dense than the same parcel of gas in the stationary wind with reduced (enhanced) driving flux. Therefore, it will feel a smaller (larger) acceleration, compared to the reduced flux stationary wind. The mass flux has already been set at the critical point, so this reduced (enhanced) acceleration leads to a slower (faster) wind. Therefore, the wind will be slower (faster) than for a steady accelerating flux. Because we are sinusoidally varying the flux rather than monotonically changing it, we can induce oscillations in the outflow velocity. These fluctuations tend to damp out as you move further out in the wind. The effects described above are most pronounced in the innermost parts of the flow (r = 2 r * ) and smaller as we move out in radius to (r = 20 r * ). At the outermost radius, the source must be varied on timescales roughly 5 times larger to generate a commensurate change in the ∆X/(X+ − X−) profile.
For driving frequency ω 1/2 ωc, the velocity fluctuations are greater than the fiducial solution, whereas for ω 1/2 ωc fluctuations are less than the fiducial solution. This enhancement can be understood by considering the phase difference between the Eddington paramerter and the dynamical variables. In Fig. 3, we plot the phase shift for each variable as a function of driving frequency at r = 2 r * . At small driving frequency, density perturbations are in phase (φρ = 0) and velocity perturbations are φv = π/2 out of phase. This is analogous to the case of a simple harmonic oscillator, where the force/acceleration are π/2 out of phase with the velocity. As the driving frequency approaches ω = ωc, the velocity phase difference φv ≈ 0 before becoming negative at larger driving frequency.
Consider what happens if density is enhanced (suppressed) at the critical point in the stationary flow. Because mass flux is solely determined by the Eddington parameter at the critical point, increasing (decreasing) density at the critical point will decrease (increase) velocity. When ωS 1/2ωc, both density and velocity lag the driving radiation i.e. φρ, φv < 0. This means that the positive acceleration phase occurs when the flow is denser, and hence the negative acceleration phase occurs when the flow is less dense. Hence the velocity minimum will be enhanced. When ωS 1/2ωc, the density lags the driving radiation φρ < 0 but the velocity leads the radiation φv > 0. The positive accelertion phase therefore occurs when the flow is less dense and the velocity maximum is enhanced. At ωS ≈ 1/2 ωc, the velocity is nearly in phase with the radiation, φv = 0 and there is no significant velocity enhancement. As with the amplitude of the fluctuations, the turnover in the phase angle plot (Fig. 3) is shifted to smaller frequencies by a factor of 5 at the outermost radius. . Phase shift φ between peaks in the Eddignton parameter Γ(t) and peaks in the density ρ (red circles), velocity v (green triangles) and mass fluxṀ (blue diamonds) as a function of the driving angular frequency ω S at r = 2 r * .

DISCUSSION
We investigated the effect of varying the amplitude of the luminosity variation with A = 0.01 and 0.9. For A = 0.01 our results are qualitatively unchanged. Dynamical variables oscillate sinusoidally at late times and we find changes in the velocity as high as ∆v/(v+ − v−) ≈ 10. The change in the velocity is approximately given by the amplitude of the luminosity oscillations, ∆v ≈ ±0.5% v0 ∼ A/2 v0, as with our fiducial case A = 0.2. In the case A = 0.9 results are qualitatively different. The density and mass flux vary sinusoidally, but the velocity profile behaves like a sawtooth wave. The velocity rise-time for TS = 1.16 t * is ∆trise ≈ 1 TS whereas ∆t fall ≈ 5.3 TS. Ulmschneider (1970Ulmschneider ( , 1971) has shown that saw-tooth waves can produce shocks, which will heat the gas. Our isothermal treatment may therefore be insufficient for these larger amplitude cases. The solution reaches a quasi-periodic state, with fluctuations occasionally producing spikes in the density or velocity profile. This case is of less physical interest, as the amplitude is large and thus we do not expect the solution to vary as smoothly as for cases with small amplitude. We do not see a clear trend in boosts to the velocity and only mention this case to confirm that at large amplitude the behaviour of our fiducial case breaks down. Time variability is a ubiquitous feature of outflowing systems. An important question is determining whether this variability is due to changes in the central object or the result of changes in the wind. Fortunately, the mechanism responsible for driving the flow couples the central object to the wind, allowing for the possible lifting of this degeneracy. In the context of AGN, studies of BAL QSOs have shown that ∼ 5% of systems exhibit disappearing X-ray absorption features (De Cicco et al. 2018). BAL variability is seen simultaneously, at multiple velocities and hence at multiple radii in the flow. This suggests the need for a global mechanism for producing BAL variability, and not one that only affects individual clouds. One idea is that variability in the X-ray flux produces global changes in the ionization state of the clouds, which in turn changes the acceleration of clouds due to line driving. We explore this scenario in a limited way in the present work, since the CAK model for line driving depends only on the Eddington fraction from the source. More sophisticated model for line driving may determine the ionization state from the X-ray flux, and the line driving from the UV flux. One may then expect correlations betweeen variability in the X-ray and UV features and therefore require multi-wavelength observation campaigns. The other possible scenario involves clouds in the innermost part of the flow shielding gas in the outermost parts as they orbit at the base of the wind. Simulations in DP17 and DP18 have shown that clumps naturally form in line driven disc winds due to symmetry breaking of the velocity field. They found clumps were primarily restricted to the base of the wind and therefore disappearing absorption troughs should appear on the orbital time-scale, over which clumps orbit into the line of sight. Intrinsic source variability may provide an additional mechanism for clump generation, which may change the size, physical distribution and formation time-scale of clumps. This may be important in the shielding gas scenario.
Another application in the context of AGN is for reverberation mapping studies. The key variable in the reverberation mapping formalism is the transfer function, which describes the emission line response of the gas as seen by a distant observer due to a delta-function pulse of continuum radiation from the source. The simplest models assume a stationary transfer function i.e. any time variability in the line emission is solely due to variability in the continuum. A more complete model would account for a time dependent distribution of gas, accounted for in a time dependent response function. Such non-stationary transfer functions require knowledge of the radial velocity of the outflow. Since in line driven winds the outflow velocity depends on the radiation source, this could lead to enhanced coupling between the continuum and emission line response. Such coupling may be weaker in the case of thermal winds for instance and absent in the case of magnetic driving.
Time variability of outflows is is also important in the context of CVs, such as BZ Cam which exhibits a fast and rapidly varying wind (Ringwald and Naylor 1998). CIV absorption line equivalent widths vary in time on ∼ 10Å, corresponding to velocity dispersions of ∆v ∼ 2000 km/s, comparable to the outflow velocity (Prinja et al. 2000). DP18 found that velocity dispersion of clumps is ∼ 5% of the outflow velocity. The non-stationarity of line driven disc winds may therefore be insufficient to explain highly variable systems such as BZ Cam. As pointed out by Prinja et al (2000), for compact systems such as CVs the variability time-scale is comparable to the flow time. Therefore variability may be due to the structure of the wind changing, due to changes in the driving luminosity for example, and not because clumps in the wind have evolved. This is different from the model considered in this work, where variability is on time-scales shorter than the flow time and induces perturbations in the wind rather than altering the entire wind.
Other CV systems such as IX Vel and V3885 Sgr have much more stationary winds (Hartley et al 2002). On timescales shorter than the flow time, they exhibit variability of v ∼ 90 km/s and v ∼ 130 km/s respectively, roughly 10% of the outflow velocity. However, they show little evidence of any density structure. CV systems thus exhibit different types of temporal variability, which suggests the need to explore different mechanims to break stationarity of wind solutions.
Another possible application of this model is to spectral variations in O stars. Ramiaramanantsoa et al. (2018) used the BRITE-Constellation's high-precision, time-dependent photometry to study the early O-type supergiant ζ Puppis. After subtracting out a 1.78 d periodic signal due to the stellar rotation, they found stochastic variation of ∼ 20 mmag. They showed that these variations in magnitude, which probes the stellar photosphere, are positively correlated with variations in the He II 4686Å emission line, which probes clumps in the wind.
A variation of 20 mmag corresponds to fractional changes in intensity A = ∆I/I ≈ 0.02. These variations were found to be coherent on time-scales of hours. With a stellar radii r * ≈ 18 R and wind terminal velocity v∞ = 2300 km/s, we estimate the dynamical time τc ≈ r * /v∞ = 1.5 h. This is comparable to the time scales of variability explored in our model. Previous study of the He II emission line of this system by Eversberg, Lepine & Moffat (1998) had found that velocity perturbations are largest within R 2R * and tend to decrease at larger radii, consistent with our findings that perturbations are largest near the base of the wind. We stress that direct comparison of our model to observations is overly optimistic, since previous work has shown the need to properly account for radiative transfer in massive stars. However, it serves as further motivation to study models with intrinsic source variability.
Additional observational campaigns of O and Wolf-Rayet stars are required to further investigate the relationship between intrinsic stellar variability and variability of accelerating clumps. Given the faintness of such stars and the need for high S/N observations, this is well suited for upcoming Thirty Meter Telescope observations.

CONCLUSION
Time variability of emission and absorption features is a fairly ubiquitous feature of outflows from O stars, CVs and AGN. Numerical modeling has shown that when outflows are radiation driven, variability can be induced by instabilities (such as the LDI), geometric effets (as in disc wind models) or by intrinsic variability of the source (as in this work). As a first step, in understanding the effects of intrinsic source variability, we have investigated how a timevarying source alters the otherwise stationary CAK solution.
Our model makes several simplifying assumptions and is only a first step in understanding time variability of line driven outflows. Within the Sobolev approximation, it may be important to have a more complete treatment of the radiation field. The simplest approach may be to account for optical depth effects, whereby continuum radiation is absorbed as radiation propagates outwards through the wind. At large radii, accounting for optical depth of the over/under densities should be small since they will average out over many oscillations. At small radii, the effect may be more important, as flux will depend on the sign of the density perturbation. Beyond treating optical depth effects, one could perform full radiation transfer and allow for the density perturbations in the wind to couple directly to the radiation field. This may introduce additional time-scales into the problem, beyond the source frequency, and generate additional perturbations in the dynamical variables. Such a treatment may be important for line driven disc winds, relevant to CVs and AGN, where clouds form at the base of the wind.
Beyond how the radiation field is treated, we must also consider the parametrization of the line driving force itself. The CAK model for line driving is the simplest possible parameterization of the line force. In more sophisticated models of line driving (see for example OCR88) the UV flux determines the line driving but the X-ray flux determines the ionization state, and hence the force multiplier. Using photoionization codes, state of the art line lists may be used to compute the force multiplier from the ionization state and temperature of the gas, rather than relying on approximate powerlaw fits of the force multiplier in terms of the optical depth parameter as in CAK (Dannen, Proga et al. in prep).
Our model assumes the Sobolev approximation is valid. As such, our non-stationary solution is fundamentally different from models where time-dependence is triggered by the LDI. Owocki and Rybicki (1984) showed that for perturbations generated by the LDI, φv = π/2 and φρ = −π/2, with additional corrections due to thermal effects of order 1/HEP. These points do not appear in our phase diagram, which suggests that simulations which resolved the LDI would behave differently and exhibit more complex behaviour.
This model finds that density and velocity perturbations are induced by a time varying radiation source in a non-rotating, spherically symmetric flow. In spherical geometries rotational effects have been shown to increase mass flux and outflow velocities (Friend & Abbott 1986;Pauldrach et al. 1986) and recent work by Araya et al. (2018) have studied time-dependent transitions between slow and fast rotating solutions. In more complex geometries, such as line driven disc winds, complex structure form due to symmetry breaking between the spherically symmetric gravitational potential and the axisymmetry of the matter reservoir and radiation field. DP17 and DP18 characterized the formation of 3D non-axisymmetries at the base of such line driven disc winds, so called clouds, assuming the wind is optically thin to the continuum and isothermal. A time varying radiation field, due to a non-stationary rate of accretion for instance, may lead to an enhancement in the rate of cloud formation. Furthermore, simulations by Waters and Proga (2015) have demonstrated that clouds are sensitive to thermal effects. Clouds may grow due to the thermal instability, and their opacity, which controls the strength of the line driving, depends sensitively on their temperature. Therefore understanding the dynamics at the base of the wind requires proper treatment of the thermodynamics and radiation transfer of the clouds. Ultimately the correct treatment of the clouds at the wind base is necessary to resolve the question of whether variability of lines in the wind is a consequence of variability in the ionization state of the gas or variability in the distribution of clouds along the line of sight, the so called covering factor.
Additional driving mechanisms may also be operating in these systems, which may also affect variability.  showed that thermal driving is sensitive to the relative intensity of UV and X-ray photons, as with line driving. Therefore we expect thermal and line driving effects to be sensitive to the particular SED in the system. Further, magnetic driven effects may be important, particularly in AGN systems.