1. Introduction
Winds are a common and dynamically important feature of many astrophysical systems, including OB-type and protostellar winds (Castor et al. Reference Castor, Abbott and Klein1975a; Weaver et al. Reference Weaver, McCray, Castor, Shapiro and Moore1977; Vink et al. Reference Vink, de Koter and Lamers2001), Wolf-Rayet (WR) stars (Avedisova Reference Avedisova1972; Crowther Reference Crowther2007), white dwarfs (Cordova & Mason Reference Cordova and Mason1982), stellar clusters (Chevalier & Clegg Reference Chevalier and Clegg1985), active galactic nuclei (King & Pounds Reference King and Pounds2015; Murray et al. Reference Murray, Chiang, Grossman and Voit1995), and starburst galaxies (e.g. Thompson et al. Reference Thompson, Quataert, Zhang and Weinberg2016). Related hydrodynamic structures may also arise in some tidal disruption events (TDEs) (Rees Reference Rees1988; Gezari Reference Gezari2021), where the tidal destruction of a star by a supermassive black hole leads to rapid fallback of stellar debris and the formation of an often super-Eddington accretion flow (Evans & Kochanek Reference Evans and Kochanek1989; Wu et al. Reference Wu, Coughlin and Nixon2018). That accretion flow can in some cases launch a relativistic jet or disk wind (e.g. Strubbe & Quataert Reference Strubbe and Quataert2011). While only a small number of TDEs show clear relativistic jets or spectroscopically confirmed outflows, many now exhibit late-time radio emission consistent with an outflow interacting with the circumnuclear medium (Giannios & Metzger Reference Giannios and Metzger2011; Horesh et al. Reference Horesh, Cenko and Arcavi2021; Cendes et al. Reference Cendes2024; Matsumoto & Piran Reference Matsumoto and Piran2024). The shocks driven by that interaction can produce a distinct observable flare, separate from the thermal optical, ultraviolet, and X-ray emission associated more directly with the accretion episode itself. Similar considerations may apply to at least some luminous fast blue optical transients (LFBOTs), for which central-engine activity has been proposed in a number of events and whose radio/X-ray observations likely imply fast ejecta interacting with dense circumstellar material (Margutti et al. Reference Margutti2019; Ho et al. Reference Ho2019, Reference Ho2020; Coppejans et al. Reference Coppejans2020; Yao et al. Reference Yao2022; Perley et al. Reference Perley2026).
In each case, the wind inflates a hot, pressurised ‘bubble’ that is bounded by both a forward and reverse shock, where the former sweeps up and heats the ambient material and the latter (also called a ‘wind termination shock’) decelerates and pressurises the otherwise-cold wind material (Castor et al. Reference Castor, McCray and Weaver1975b). These wind-driven bubbles play a central role in regulating feedback in the interstellar medium (ISM) (Weaver et al. Reference Weaver, McCray, Castor, Shapiro and Moore1977; Kim et al. Reference Kim, Ostriker and Kim2013; Fierlinger et al. Reference Fierlinger, Burkert, Ntormousi, Fierlinger, Schartmann, Ballone, Krause and Diehl2016; Fielding et al. Reference Fielding, Ostriker, Bryan and Jermyn2020), and accurately characterising their hydrodynamic evolution is therefore essential for understanding problems ranging from stellar feedback to large-scale galactic outflows.
When the wind ‘turns on’ and initially impacts its surroundings, radiative losses are dynamically insignificant and the system evolves adiabatically. The shocked regions then expand and cool while the dynamical timescales increase and radiative losses become more significant. As outlined in Castor et al. (Reference Castor, McCray and Weaver1975b) and Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977), radiative cooling is initially only important for the forward-shocked medium, which compresses into a thin and dense shell while the reverse-shocked interior remains hot and continues to drive the expansion through
$p-dV$
work. As the system continues to evolve and the interior pressure declines due to continued expansion and cooling, radiative losses become important for the reverse-shocked wind as well, and the entire system ultimately approaches a momentum-conserving ‘snowplow’ phase (Ostriker & McKee Reference Ostriker and McKee1988).
In the original analyses of Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977) and Avedisova (Reference Avedisova1972), motivated primarily by stellar feedback into a uniform ISM, the early adiabatic phase was found to be comparatively short-lived (
$\sim2\times10^3\,$
yr) and therefore of limited physical significance. This phase corresponds to the classic energy-conserving stage, in which the shocked wind has relaxed to an approximately uniform-pressure interior, the reverse shock lies well interior to the contact discontinuity, the shocked ejecta density is well below that of the shocked ambient medium, and the flow is therefore well approximated as isobaric. In such environments, the dense swept-up shell cools efficiently, rapidly driving the flow into the subsequent radiative phases. However, for decreasing ambient medium density profiles (
$\rho \propto r^{-n}$
where
$n\gt 0$
), reduced post-shock densities are expected to increase the cooling time relative to the expansion time, allowing the adiabatic phase to persist substantially longer. Such scenarios could arise, for instance, when a newly-launched, faster stellar wind from a later evolutionary phase overtakes and sweeps up the slower, denser outflow from an earlier phase (e.g. Garcia-Segura et al. Reference Garcia-Segura, Langer and Mac Low1996), when a galactic outflow propagates through a gravitationally- or hydrostatically-stratified halo or circumgalactic medium (CGM) (Mac Low & McCray Reference Mac Low and McCray1988), or when an outflow from a tidally disrupted star interacts with the circumnuclear medium.
Furthermore, Coughlin (Reference Coughlin2024) recently identified additional self-similar solutions that characterise the flow at even earlier times, before the onset of the aforementioned energy-conserving phase. At this stage, the freely expanding wind still interacts more directly with the ambient medium: the ratio f of the wind density to the ambient density is high, the reverse shock is not well separated from the contact discontinuity, and the shocked ejecta is not well-approximated as isobaric. As with the later adiabatic stage, steeper ambient density gradients prolong this early regime by delaying the attainment of this smaller overdensity ratio (e.g. see the figures and analysis in the following sections).
Taken together, these considerations suggest that wind-driven systems evolve through multiple distinct similarity regimes whose durations and relative dynamical importance depend sensitively on the ambient medium density profile and initial density contrast. In lower-energy systems, this initial density contrast may be vanishingly small, so that the shocked wind rapidly becomes underdense relative to the shocked ambient medium, and the interior flow reaches the approximately isobaric structure of the energy-conserving phase with little or no extended interaction regime. By contrast, in higher-energy outflows, such as those associated with WRs, TDEs, LFBOTs, etc., the initial overdensity ratio can be extremely high, allowing for a more prolonged interaction-dominated phase as discussed in detail in Section 4.2.
While the interaction solutions and the later energy-conserving solutions are each well understood in their respective limits, the transition between these regimes (and the accuracy with which they describe the full hydrodynamic evolution) has not been systematically explored, particularly for stratified ambient media with
$\rho \propto r^{-n}$
. In this work, we present a unified study of the adiabatic evolution of wind-driven bubbles expanding into power-law density profiles, comparing analytic and semi-analytic scalings to one-dimensional hydrodynamic simulations. Because we are interested in applications to high-energy systems, we largely consider
$f \gg 1$
. We quantify the ‘relaxation timescales’ toward each similarity state and assess the accuracy of the solutions during transitions between self-similar regimes in different ambient density profiles.
This paper is structured as follows. In Section 2, we derive scalings for relevant timescales and shock positions in both similarity regimes. We use Section 3.1 to describe the novel shock-capturing numerical scheme we implement to resolve the narrow shocked region of interest, and we state our simulation parameters and initial conditions. We compare our hydrodynamic simulations to the relevant self-similar solutions in their respective limits in Section 3.2, and interpret our results in the broader context of astrophysical winds in Section 4. We conclude with a summary of our findings in Section 5.
2. Basic considerations
We consider a time-steady, spherically symmetric wind launched at radius
$R_{w,0}$
with constant mechanical luminosity
$L_w=\frac{1}{2}\dot{M}v_w^2$
, where
$v_w$
is the wind velocity and
$\dot{M}$
is the wind mass flux. The wind density is then
$\rho_w(r) = \rho_{w,0}(r/R_{w,0})^{-2}$
, where
$\rho_{w,0} = \dot{M}/(4\pi v_w R_{w,0}^2)$
. This wind expands into an ambient medium with density
$\rho_a(r) = \rho_{a,0}(r/R_{w,0})^{-n}$
, where
$\rho_{a,0}$
is a constant. The collision of the wind and the ambient medium produces a double-shock structure: a reverse shock (or wind termination shock) at radius
$R_{{r}}(t)$
, a forward shock at
$R_{{s}}(t)$
, and a contact discontinuity at
$R_{{c}}(t)$
that separates the shocked wind from the shocked ambient gas. At the earliest times, when the wind density can greatly exceed that of the ambient medium, the dynamics are governed by direct ram-pressure interaction and rapidly approach the interaction-dominated, or ‘coasting,’ self-similar regime described by Coughlin (Reference Coughlin2024). Once the ambient and ejecta densities become comparable, it is expected that the system transitions to an energy-conserving expansion described by the first stage of the solutions presented in Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977), Ostriker & McKee (Reference Ostriker and McKee1988), Koo & McKee (Reference Koo and McKee1992).
2.1. Interaction phase
At sufficiently early times, the freely expanding wind can be much denser than the ambient medium,
$\rho_w \gg \rho_a$
, and the shocked wind speed is comparable to the free-wind speed,
$v_{w}$
. The forward shock propagates faster than the wind speed, such that the relative thickness of the forward-shocked shell,
\begin{equation} \Delta_{{s}} \equiv \frac{R_{{s}}}{R_{{c}}}-1, \end{equation}
is nearly constant. Conversely, the reverse shock gradually separates from the contact discontinuity, so that the relative thickness of the reverse-shocked shell,
\begin{equation} \Delta_{{r}} \equiv 1 - \frac{R_{{r}}}{R_{{c}}},\end{equation}
is fundamentally time-dependent and evolves in a way that is set primarily by the smallness parameter
\begin{equation} \delta(R_{{c}}) \equiv \left( \frac{\rho_{a,0}}{\rho_{w,0}} \right)^{1/2}\left[ \frac{R_{{c}}(t)}{R_{w,0}} \right]^{(2-n)/2}.\end{equation}
In the regime of interest,
$\delta \ll 1$
, and the relative shell thickness
$\Delta_{r}$
remains small. In this regime,
$\Delta_{{r}}$
can also be expressed as
\begin{equation} \Delta_{{r}} = \kappa_{{r}}\delta(R_c),\end{equation}
where
$\kappa_{{r}}$
is a dimensionless constant set by the self-similar solution.
Despite the asymmetry between the forward- and reverse-shocked regions, the shocked fluid variables within both the forward- and reverse-shocked shells admit approximately self-similar forms. Specifically, in the forward-shocked shell,
\begin{equation}\begin{split} v_{{s}}(r,t) &= V_{{c}}f_{{s}}(\eta)\,,\\ \rho_{{s}}(r,t) &= \rho_{a,0}\left( \frac{R_{{c}}}{R_{w,0}} \right)^{-n}g_{{s}}(\eta)\,,\\ p_{{s}}(r,t) &= \rho_{a,0} V_{{c}}^2\left( \frac{R_{{c}}}{R_{w,0}} \right)^{-n}h_{{s}}(\eta), \end{split}\end{equation}
and in the reverse-shocked shell,
\begin{equation}\begin{split} v_{{r}}(r,t) &= V_{{c}}\{1+\Delta_{{r}}f_{{r}}(\eta)\}\,,\\ \rho_{{r}}(r,t) &= \rho_{w,0}\left( \frac{R_{{c}}}{R_{w,0}} \right)^{-2}g_{{r}}(\eta)\,,\\ p_{{r}}(r,t) &= \rho_{w,0}\left( \frac{R_{{c}}}{R_{w,0}} \right)^{-2}V_{{c}}^2\Delta_{{r}}^2h_{{r}}(\eta),\end{split}\end{equation}
Here
$V_{{c}} = dR_{{c}}/dt$
, and
$\eta = (r - R_{{c}}) / (R_{{c}} - R_{{r}})$
in the reverse-shocked shell, while
$\eta = (r - R_{{c}}) / (R_{{s}} - R_{{c}})$
in the forward-shocked shell. Substituting these expressions into the Euler equations reduces the flow to a coupled system of ordinary differential equations for the dimensionless functions
$f_{{s,r}}$
,
$g_{{s,r}}$
, and
$h_{{s,r}}$
, which determine the structure of the shocked gas in each region. We make use of additional corrections for the forward-shocked variables to account for the effects of finite f, outlined in Coughlin (Reference Coughlin2024).
In these solutions, the position of the contact is given explicitly by
\begin{equation} R_{{c}}(t) = R_{w}[1-\kappa_{{c}}\delta(R_{{c}})],\end{equation}
where
$R_w = R_{w,0} + v_w t$
, and
$\kappa_{{c}}$
is a dimensionless constant fixed by the similarity solution.
2.2. Transition to energy-conserving phase
The interaction solution remains valid so long as the reverse-shocked region remains thin relative to the forward-shocked shell. A natural estimate for the end of this phase is obtained by requiring that the reverse-shocked shell no longer remain thinner than the forward-shocked shell, i.e.
$\Delta_{{r}} \sim \Delta_{{s}}$
. Letting
$t_{\mathrm{ss,end}}$
represent the time at which this occurs, then
\begin{equation} \kappa_{{r}}\delta(t_{\mathrm{ss,end}}) = \Delta_{{s}}.\end{equation}
Using the expression for
$\delta(t)$
outlined earlier, this becomes
\begin{equation} \frac{R_{w}(t_{\mathrm{ss,end}})}{R_{w,0}} = \left[ \frac{\Delta_{{s}}}{\kappa_{{r}}}\left( \frac{\rho_{w,0}}{\rho_{a,0}} \right)^{1/2} \right]^{2/(2-n)}.\end{equation}
Assuming
$R_{w}(t_{\mathrm{ss,end}}) \sim v_wt_{\mathrm{ss,end}}$
, one finds
\begin{equation} t_{\mathrm{ss,end}} = \frac{R_{w,0}}{v_w}\left( \frac{\Delta_{{s}}}{\kappa_{{r}}} \right)^{2/(2-n)} \left( \frac{\rho_{w,0}}{\rho_{a,0}} \right)^{1/(2-n)}. \end{equation}
For
$n \lt 2$
, this yields a finite transition time at which the interaction solution ceases to apply. As
$n \rightarrow 2$
, the exponent diverges, indicating that the interaction regime no longer transitions to a distinct later energy-conserving phase.
Once the reverse-shocked shell width becomes comparable to the width of the forward-shocked shell, or equivalently, once the density of the swept ambient medium becomes comparable to that of the shocked wind, the flow approaches the pressure-driven, energy-conserving regime first analysed by Avedisova (Reference Avedisova1972), Falle (Reference Falle1975), Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977), in which the shocked wind remains hot and approximately isobaric, and the expansion is driven primarily by the conversion of wind kinetic energy into thermal energy at a radius deeply interior to the contact. For a power-law ambient medium, the corresponding generalisation of the forward-shock evolution was given by Ostriker & McKee (Reference Ostriker and McKee1988), while the reverse-shock scaling was discussed in Koo & McKee (Reference Koo and McKee1992).
The temporal scaling of the forward shock (and contact discontinuity) in this regime can be derived by assuming that (1) the rate at which energy is injected into the forward-shocked shell scales in proportion to the (assumed time-independent) energy injection rate from the wind, which should be upheld when
$R_{r} \ll R_{c}$
(or when the sound speed of the shocked ejecta is large in comparison to the speed of the reverse shock itself; see Weaver et al. Reference Weaver, McCray, Castor, Shapiro and Moore1977), (2) the pressure at the contact discontinuity scales with the post-forward-shock pressure, and (3) the radius of the contact discontinuity scales in proportion to the radius of the forward shock. Taken together, these assumptions predict
\begin{equation}\begin{split} \dot{E} &= 4\pi R_{c}^2V_{c} p_{c} \propto R_{s}^2 V_{s} R_{s}^{-n}V_{s}^2 = \mathrm{const.} \\ &\Rightarrow V_{s} \propto R_{s}^{(n-2)/3},\end{split}\end{equation}
such that the temporal power-law solutions for the radius of the forward shock and the pressure of the shocked wind (equal to the pressure at the contact discontinuity) are given by
\begin{equation} R_{{s}}(t) \propto t^{3/(5-n)}\,, \quad p_{c}(t) \propto t^{-(4+n)/(5-n)}. \end{equation}
Within the shocked ambient shell, the flow admits the usual self-similar structure, which can be defined in terms of the standard self-similar variable
$\xi = r/R_{s}(t)$
that treats the forward shock radius as the fundamental radial scale. On the other hand, because we know that the solution must terminate in a contact discontinuity, it is more natural to introduce the same self-similar variable that is used for the interaction phase, i.e. the energy-conserving solutions for the forward shock can be written in precisely the same way as Equation (5). Because this simply corresponds to a change of variables, the functions are the same as those derived in Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977), now written in terms of
$\eta$
and rescaled by the appropriate factors of
$V_{c}/V_{s}$
.
The structure of the shocked wind can be determined under the isobaric assumption: with the temporal scaling of the pressure given in Equation (12), the gas-energy equation becomes
\begin{equation} \frac{1}{r^2}\frac{\partial}{\partial r}\left[r^2 v_{{r}}\right] = -\frac{1}{\gamma}\frac{1}{p_{c}}\frac{\partial p_{c}}{\partial t} = \frac{n+4}{3\gamma}\frac{V_{c}}{R_{c}},\end{equation}
which is integrated alongside the boundary condition
$v_{{r}}(R_{c}) = V_{c}$
to give
\begin{equation} v_{{r}} = V_{c}\left(\frac{n+4}{9\gamma}\frac{r}{R_{c}}+\left(1-\frac{n+4}{9\gamma}\right)\left(\frac{r}{R_{c}}\right)^{-2}\right). \end{equation}
This shows that the velocity structure of the shocked wind is self-similar, and depends only on
$r/R_{c}$
. Since the reverse shock has already been assumed to be at small radii relative to the contact discontinuity, this also shows
\begin{equation} v(R_{r}) \simeq \left(1-\frac{n+4}{9\gamma}\right)V_{c}\left(\frac{R_{r}}{R_{c}}\right)^{-2} \simeq \frac{\gamma-1}{\gamma+1}v_{w},\end{equation}
where the final equality follows from the jump conditions at the reverse shock and the assumption that the reverse shock speed is much less than the wind speed. The radius of the reverse shock therefore satisfies
\begin{equation} R_{r} = R_{c}\sqrt{\frac{\gamma+1}{\gamma-1}\left(1-\frac{n+4}{9\gamma}\right)}\sqrt{\frac{V_{c}}{v_{w}}} \propto R_{c}^{\frac{n+4}{6}} \propto t^{\frac{n+4}{2\left(5-n\right)}}, \end{equation}
where the last proportionality follows from Equation (12) (recall that the radius of the contact discontinuity is assumed to scale in proportion with the radius of the forward shock).
Finally, we can determine the density profile of the shocked wind by noting that the boundary condition at the reverse shock is
$\rho(R_{r}) \propto R_{r}^{-2} \propto R_{c}^{-(n+4)/3}$
. Since the velocity profile (14) is self-similar and
$R_{r}/R_{c} \simeq 0$
in this regime, the density permits a self-similar solution of the form
\begin{equation} \rho_{{r}} = \frac{\gamma+1}{\gamma-1}\rho_{w,{r}}\left(\frac{R_{c}}{R_{c, 0}}\right)^{-\frac{n+4}{3}}g\left(\frac{r}{R_{c}}\right),\end{equation}
where
$R_{c, 0}$
is the radius of the contact discontinuity at some initial time,
$\rho_{w,{r}}$
is the wind density at the reverse shock when
$R_{c}$
coincides with that scale radius, and
$g(0) = 1$
to satisfy the boundary condition at the reverse shock (again assuming
$R_{r}/R_{c} \simeq 0$
). Then using the velocity profile in the continuity equation shows that
\begin{equation} g\left(\frac{r}{R_{c}}\right) = \left(1-\frac{r^3}{R_{c}^3}\right)^{-\frac{\left(n+4\right)\left(\gamma-1\right)}{9\gamma-n-4}},\end{equation}
and hence the density is
\begin{equation}\begin{split} \rho_{{r}} &= \frac{\gamma+1}{\gamma-1}\rho_{w,{r}}\left(\frac{R_{c}(t)}{R_{c, 0}}\right)^{-\frac{n+4}{3}}\left(1-\frac{r^3}{R_{c}(t)^3}\right)^{-\frac{\left(n+4\right)\left(\gamma-1\right)} {9\gamma-n-4}} \\ &\propto t^{-\frac{n+4}{5-n}}\left(1-\frac{r^3}{R_{c}(t)^3}\right)^{-\frac{\left(n+4\right)\left(\gamma-1\right)} {9\gamma-n-4}}.\end{split}\end{equation}
These solutions were presented in Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977) for the specific case of
$n=0$
and
$\gamma=5/3$
; while the velocity profile was later generalised in Koo & McKee (Reference Koo and McKee1992), we have not encountered these generalisations for
$\rho$
elsewhere. Note that they break down for
$n \ge 2$
, as they are based on the assumption that
$R_{r} \ll R_{c}$
, which – from Equation (16) – will only be upheld when
$n \lt 2$
. They also illustrate that, in the limit that the radius of the reverse shock is much less than the contact discontinuity, the inner shocked-wind flow is also self-similar in terms of the same variable as the interaction phase for
$R_{r} \lt r \lt R_{c}$
. Specifically, because
$R_{r} = 0$
is just a limiting case, the preceding solutions can be written
\begin{equation} p_{{r}} = \rho_{a,0} V_{c}^2\left(\frac{R_{c}}{R_{c, 0}}\right)^{-n}h_{s}(0),\end{equation}
\begin{align} v_{{r}} &= V_{c}\Bigg(\frac{n+4}{9\gamma}\left(1+\left(1-\frac{R_{r}}{R_{c}}\right)\eta\right) \nonumber\\ & \quad +\left(1-\frac{n+4}{9\gamma}\right)\left(1+\left(1-\frac{R_{r}}{R_{c}}\right)^{-2}\right)\Bigg)\end{align}
\begin{align} \rho_{{r}} &= \frac{\gamma+1}{\gamma-1}\rho_{w,{r}}\left(\frac{R_{c}}{R_{c, 0}}\right)^{-\frac{n+4}{3}} \nonumber\\ & \quad \times\left(1-\left(1+\left(1-\frac{R_{r}}{R_{c}}\right)\eta\right)^3\right)^{-\frac{\left(n+4\right)\left(\gamma-1\right)}{9\gamma-n-4}}\!,\end{align}
where
$h_{s}(0)$
is the self-similar pressure at the contact discontinuity as determined from the forward shock solution, and
$\eta$
is the same variable as introduced in Section 2.1 for the reverse-shock shell. This demonstrates that the Weaver solution is simply the one that approximates
$R_{r} = 0$
, and this representation can be used to directly compare to the interaction-phase solutions given any (finite) value of
$R_{r}/R_{c}$
.
As an aside, the equivalent parametrisation (i.e. substitute the standard self-similar variable
$\xi$
for
$\eta$
) can be used for the Sedov–Taylor blastwave (Sedov Reference Sedov1946; Taylor Reference Taylor1950) in scenarios where the contact discontinuity approaches the origin (relatively small values of n for a given
$\gamma$
; see Goodman Reference Goodman1990). Analogously to the present case where the reverse shock reaching the origin constitutes the “true” self-similar solution that is reached at asymptotically late times, the exact Sedov–Taylor solution holds when the contact discontinuity reaches the origin. Nevertheless, and just as the self-similar velocity profile of the shocked wind can be used to constrain the propagation rate of the reverse shock, i.e. Equation (16), the velocity profile of the Sedov–Taylor blastwave can be used to predict the propagation rate of the contact discontinuity: in the limit that
$r\rightarrow 0$
, the Sedov–Taylor velocity profile satisfies
\begin{equation} v(r,t) = \frac{2}{\gamma+1}V_{s}(t)\frac{r}{R_{s}(t)}.\end{equation}
Because the fluid velocity must be continuous across the contact discontinuity, it follows that its velocity and time dependent position are
\begin{equation} V_{c} = \frac{2}{\gamma+1}\frac{V_{s}}{R_{s}}R_{c} \,\,\, \Rightarrow \,\,\, R_{c} \propto R_{s}^{\frac{2}{\gamma+1}} \propto t^{\frac{4}{\left(5-n\right)\left(\gamma+1\right)}}.\end{equation}
Thus, and even though this is a distinct power-law dependence from the forward shock, the time-dependent evolution of the contact discontinuity can be gleaned from the approximate self-similar behaviour of the shocked ambient fluid.
3. Hydrodynamic simulations
3.1. Numerical methodology
Although this work concerns the properties of non-relativistic winds (
$u_w \lesssim 0.01$
), our numerical simulations employ a one dimensional finite-volume special relativistic hydrodynamics (SRHD) code originally developed for modelling gamma-ray burst afterglows, as its moving-mesh framework proves particularly well-suited to this problem. Specifically, the domain boundaries are designed to track discontinuities (shocks and contact discontinuities), which ensures consistently high resolution in the narrow shocked regions without wasting grid cells on the unshocked wind and ambient medium, both of which are analytically prescribed. The details of this code, including the moving mesh algorithm, are presented in Appendix A.
All simulations involve the collision of a cold, spherical, constant-velocity wind (
$u_w=0.01$
) with a cold, power-law ambient medium, initialised at radius
$R_{w,0}=1$
and time
$t_{w,0}=1$
. It is useful to express t in dynamical times
$t_{\mathrm{dyn}}=R_{w,0}/u_w$
, noting that for our parameter choices,
$t_{\mathrm{dyn}}=100\,t_{w,0}$
. The ambient medium and unshocked wind densities are defined according to Section 2.
Comparisons of our simulated fluid profiles for
$f=10^5$
to those of the interaction solutions and of the two-shock Riemann problem (RP) solutions at two dynamical times (right) and at
$t \ll t_{\mathrm{dyn}}$
(left). For
$t \ll t_{\mathrm{dyn}}$
, only minimal radial structure has formed in each of the shocked shells; as such, the system is well-approximated by the solution to the two-shock RP. By
$t\sim t_{\mathrm{dyn}}$
, gradients have steepened in every fluid variable within each shocked shell, and the simulation has fully relaxed into the interaction regime.
The discontinuities that bound the subdomains are initialised by solving a two-shock Riemann problem (Zou et al. Reference Zou, Wu and Dai2005) over a small interval of time
$\delta t$
, where
$\delta t / t_{w,0} = 10^{-6}$
. The domain edges then consist of the reverse shock at
$R_{{r}}$
, the contact discontinuity at
$R_{{c}}$
, and the forward shock at
$R_{{s}}$
, with
$R_{{r}} \lt r \lt R_{{c}}$
comprising the reverse-shocked subdomain, and
$R_{{c}} \lt r \lt R_{{s}}$
comprising the forward-shocked subdomain. At the inner boundary, inflow conditions are imposed using the wind primitive state
$(\rho_w, u_w, p_w)$
evaluated at
$R_{{r}}$
. At the outer boundary, the ambient medium state
$(\rho_a, u_a, p_a)$
is prescribed at
$R_{{s}}$
.
We performed 25 simulations with
$\gamma=5/3$
,
$n \in [0, 0.5, 1, 1.5, 2]$
and
$f \in [10^1, 10^2, 10^3, 10^4, 10^5]$
. The two shocked subdomains were resolved with 800 zones each, which was sufficient to resolve the overall profiles; see Appendix B for a discussion of convergence. Each simulation was evolved to
$t/t_{\mathrm{dyn}}=10^{18}$
.
3.2 Results
3.2.1 Interaction phase
In the limit
$f \gg 1$
, the early-time evolution should be described by the interaction phase solutions presented in Coughlin (Reference Coughlin2024), following a transient that originates from the initial conditions – i.e. the state first has to ‘relax’ from the two-shock Riemann problem (RP) to the self-similar solution. Given the initial rate at which the CD-RS shell expands, the estimate for when the self-similar state is reached is
$t_{\mathrm{dyn}}$
(see the discussion in Section 4.1 of Coughlin Reference Coughlin2024). Figure 1 shows convergence to the interaction solutions as the system evolves from
$t \ll t_{\mathrm{dyn}}$
to
$t \sim t_{\mathrm{dyn}}$
, also depicting the simultaneous divergence from the two-shock RP solution. This transition from the initial two-shock RP to the interaction regime occurs many orders of magnitude earlier than the transition from the interaction to energy-conserving stage, which occurs approximately when the wind density declines to that of the ambient density (see Section 2.2).
Figure 2 shows the convergence of our simulations and the interaction solutions for different f at
$t\sim t_{\mathrm{dyn}}$
in the case of
$n=0$
: as the overdensity parameter is increased from
$f=10^2$
to
$10^5$
, the simulated profiles approach the interaction solutions throughout the shocked structure. For all other diagrams in this work, we use
$f=10^5$
to ensure sufficient convergence to the self-similar solutions in the interaction regime.
This agreement extends to all other n tested, as shown in Figure 3. For each value of n, the snapshot across all three fluid variables (
$\rho$
, v, and p) was chosen at the time for which the relative
$L_1$
error in the forward-shocked velocity profile, measured with respect to the interaction self-similar solution, is minimised. Across the full range of ambient density profiles considered here, the interaction solutions effectively capture the structure of the simulated profiles in both shocked regions given
$f \gg 1$
.
Comparison of the density (top), velocity (middle), and pressure (bottom) profiles from our simulations to those predicted by the interaction solutions for
$n=0$
. The representative time at which each snapshot was taken for any given n corresponds to the minimum relative
$L_1$
error of the forward-shocked shell pressure between the simulation output and the interaction-phase self-similar solutions. Since these solutions assume
$f \gg 1$
, results converge as f increases.
Radial fluid profiles from our simulations and the interaction solutions for
$f=10^5$
and
$n\in [0, \frac{1}{2}, 1, \frac{3}{2}, 2]$
. The representative time at which each snapshot was taken for any given n corresponds to the minimum relative
$L_1$
error of the forward-shocked shell velocity between the simulation output and the interaction-phase self-similar solutions. Agreement is robust across all tested ambient medium density profiles.
Relative
$L_1$
error between our hydrodynamics simulations and self-similar solutions in both interaction and energy-conserving regimes for the radial pressure profiles in the forward-shocked shell. The initial overdensity parameter f was set to
$10^5$
for each case, and
$n \in [0,\frac{1}{2},1,\frac{3}{2},2]$
. The change in error is shown vs. time (top) and instantaneous overdensity
$f(R_{{c}})\equiv f\cdot (R_{w,0}/R_{{c}})^{2-n}$
(bottom). We measure
$t_{\mathrm{ss,end}}$
to be the time at which these error curves intersect. This transition from the interaction to energy-conserving stage takes place at various times for different n, but at comparable
$f(R_{{c}}) \sim 1$
. For
$n=2$
, no such transition occurs since the system maintains a perpetual similarity state with f constant throughout its evolution.
3.2.2 Transition to energy-conserving phase
The shocked regions only remain in the interaction regime so long as the reverse-shocked shell remains thin compared to the forward-shocked shell. Once these shell widths become comparable, so too do the wind and ambient densities, and the flow is expected to depart from this interaction-dominated state and resemble the later, energy-conserving regime. To quantify this transition in the simulations, we compare the radial pressure profile in the forward-shocked shell to both sets of self-similar solutions and compute the relative
$L_1$
error as a function of time.
The top panel of Figure 4 shows the resulting errors vs. time for
$n \in [0, \frac{1}{2}, 2, \frac{3}{2}, 2]$
. At the earliest times, this error is high since we do not initialise the simulations in the interaction self-similar state; as such, the system relaxes into the interaction regime over the course of
$\sim$
a few dynamical times. Once this relaxation has occurred, the interaction solutions provide highly accurate descriptions of the fluid behaviour, until the error again begins to rise as the system gradually departs from these solutions and moves toward the energy-conserving ones. We define
$t'_{\mathrm{ss,end}}$
– the measured analogue to the derived
$t_{\mathrm{ss,end}}$
in Section 2.1 – to be the time at which the two error curves intersect, i.e. the epoch at which the interaction and energy-conserving descriptions are equally (in)accurate. For
$n=2$
, no such intersection occurs: the interaction solution remains the more accurate description throughout the evolution, consistent with the expectation that this case does not transition to a distinct later energy-conserving similarity state (see the discussion toward the end of Section 2). For all n, the error does not reach 0 for two reasons: (1) the numerical simulations are not infinitely accurate, exhibiting (e.g.) discrepancies at the CD due to the finite nature of the numerical method (see Appendix B), and (2) the interaction solutions used here do not include any of the higher-order terms discussed in Coughlin (Reference Coughlin2024), i.e. the self-similar solutions are only the leading-order terms in an infinite series in
$\delta$
.
The bottom panel of Figure 4 shows the L1 error vs. instantaneous
$f(R_c) = f \cdot (R_{w,0}/R_c)^{2-n}$
for
$n \in [0, \frac{1}{2}, 2, \frac{3}{2}, 2]$
. The transition from the interaction to energy-conserving stage occurs around
$f(R_c) \sim 1$
for all
$n \lt 2$
, demonstrating that this change takes place when the densities of the wind and ambient medium are comparable. This is an expected result since
$f(R_{c}) = 1/\delta^2(R_{{c}}) = (\kappa_{{r}} / \Delta_{{r}})^2$
, and because
$\Delta_{{s}} \sim \Delta_{{r}}$
at the end of the interaction regime,
$f(R_{c}) \sim (\kappa_{{r}} / \Delta_{{s}})^2 \sim 1$
.
The measured dependence of
$t'_{\mathrm{ss,end}}$
on the overdensity parameter is shown in Figure 5. For each
$n \lt 2$
, the measured transition time increases with both f and n, following the scaling
$t_{\mathrm{ss,end}} \propto f^{1/(2-n)}$
predicted by Equation (10). The duration of the interaction phase therefore grows rapidly with increasing ambient density slope, and as
$n \rightarrow 2$
, the exponent in the predicted scaling relation diverges, while the measured behaviour approaches a single persistent similarity state. Thus, the interaction regime can remain dynamically important for a substantial fraction of the pre-radiative evolution, especially in steep ambient density profiles.
The duration
$t'_{\mathrm{ss,end}}$
of the interaction phase as a function of initial overdensity parameter f. The duration is proportional to
$f^{1/(2-n)}$
as predicted by Equation (10), and is therefore extended for steeper ambient medium density profiles. For
$n=2$
, the system is characterised by a single perpetual self-similar state, corresponding to
$t_{\mathrm{ss,end}} \rightarrow \infty$
.
Furthermore, depending on the criterion for an acceptable threshold of error, Figure 4 shows that for
$n \lt 2$
, the window of time in which the system is not well-modelled by either set of self-similar solutions is relatively broad. Even adopting a conservative criterion in which a model is regarded as ‘valid’ once its relative error lies within an order of magnitude of its global minimum, the system can remain in this intermediate phase until times several orders of magnitude larger than the characteristic transition timescale itself. This indicates that the transition between the interaction and energy-conserving regimes is itself dynamically extended, and it may not be a suitable approximation to model the system with an abrupt change from one asymptotic solution to the other, though extending the interaction solutions to higher-order terms may shorten this transition window (as discussed further in Section 5). Additionally, we interpolate between these two regimes and provide analytic time-dependent functions for the forward shock, contact discontinuity, and reverse shock positions in Appendix C.
As the shocked regions evolve away from the interaction phase and into the energy-conserving phase, the fluid profiles and shock structure change accordingly. Figure 6 shows the trajectories of the reverse shock, contact discontinuity, and forward shock for each value of n, together with the corresponding interaction and energy-conserving solutions. At early times, all three discontinuities are well described by the interaction-stage solutions. For
$n \lt 2$
, the simulated shock structure subsequently departs from this early-time behaviour at
$t \sim t_{\mathrm{ss,end}}$
and approaches the later energy-conserving scalings, adhering specifically to the power-law shock temporal scalings in Equations (12) and (16). For
$n=2$
, the discontinuity positions remain consistent with a single similarity state, and differ in normalisation from the energy-conserving predictions since the system never leaves the interaction regime.
The positions of the reverse shock, contact discontinuity, and forward shock as a function of time for different values of n from our hydrodynamic simulations, the interaction solutions, and the energy-conserving solutions. The dotted vertical lines represent the measured timescales characterising the transitions from the interaction to energy-conserving regime. For
$n=2$
, these are not distinct regimes, and the system exists in a single perpetual similarity state.
The evolution of the fluid profiles is shown explicitly in Figure 7 for the representative case
$n=0$
,
$f=10^5$
, at epochs before, during, and after the transition. At early times (
$t \ll t_{\mathrm{ss,end}}$
), the simulated profiles are well described by the interaction solution across both shocked regions, and at late times (
$t \gg t_{\mathrm{ss,end}}$
), the same quantities are instead well reproduced by the energy-conserving solution. Near
$t \sim t_{\mathrm{ss,end}}$
, however, neither asymptotic description provides a comparably good match to the full shocked structure.
Comparisons of our simulated fluid profiles to those of the interaction and energy-conserving solutions across various epochs in the system’s evolution. The initial overdensity ratio
$f=10^5$
for this example, and the ambient medium density profile is uniform (
$n=0$
). The snapshots here correspond to instantaneous overdensities of
$f(R_{{c}})=9.227\times 10^2$
(left),
$f(R_{{c}})=1.109$
(middle), and
$f(R_{{c}})=4.843\times 10^{-5}$
(right), again reinforcing the criterion that the transition period is characterised by
$f(R_{{c}}) \sim 1$
. During the interaction (
$t \ll t_{\mathrm{ss,end}}$
) and energy-conserving (
$t \gg t_{\mathrm{ss,end}}$
) phases, the radial profiles of all fluid variables are well-modelled by one set of self-similar solutions, while the transition phase at
$t \sim t_{\mathrm{ss,end}}$
is not well-modelled by either set.
4. Discussion
4.1. Duration of the energy-conserving phase
The energy-conserving phase lasts as long as radiative losses remain dynamically negligible – e.g. as long as the forward-shocked ambient shell cannot cool on a timescale comparable to the expansion time. Beyond this point, the shell rapidly compresses into a thin, dense layer, while the reverse-shocked wind interior remains hot, marking the onset of the first radiative stage described in Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977).
An equivalent way to estimate this transition is through the comparison of the rate at which the contact discontinuity does work on the swept-up shell (i.e. the rate at which energy is added to the shocked ambient gas via the contact discontinuity) to the rate at which the shell radiates energy. The former rate is, from energy conservation and the isobaric structure of the shocked wind,
\begin{equation} \dot{E}_{{c}} = \frac{9(\gamma-1)}{9\gamma-n-4}L_w,\end{equation}
while the radiative luminosity
$L_{\mathrm{rad}}$
of the forward-shocked shell may be estimated as
\begin{equation} L_{\mathrm{rad}} \sim n_{{s}}^2\Lambda(T_{{s}})\mathscr{V}_{{s}},\end{equation}
where
$n_{{s}}$
and
$T_{{s}}$
are the characteristic post-shock number density and temperature,
$\mathscr{V}_{{s}}$
is the shell volume, and
$\Lambda(T_{{s}})$
is the cooling function evaluated at the post-shock temperature. Here we take the ion, electron, and total number densities to be equal to within factors of order unity. Using the strong shock compression ratio
$n_2/n_1 = (\gamma+1)/(\gamma-1)$
, together with
$n_{{s}} \sim \rho_{{s}}/m_p$
, gives
\begin{equation} L_{\mathrm{rad}} \sim \frac{4}{3}\pi R_{w,0}^3\left(1-\xi_{c}^3\right)n_{a,0}^2\Lambda(T_{{s}}) \left( \frac{\gamma+1}{\gamma-1} \right)^2 \left( \frac{R_{{s}}}{R_{w,0}} \right)^{3-2n},\end{equation}
where
$\xi_{c} = R_{c}/R_{s}$
is the ratio of the contact discontinuity radius to the forward shock radius, which is constant during the adiabatic (Weaver/energy-conserving) phase. The transition to the radiative stage is then estimated by requiring
$L_{\mathrm{rad}} \sim \dot{E}_c$
, which yields a shock position of
\begin{equation} \frac{R_{{s}}}{R_{w,0}} \sim \left[ \frac{27(\gamma-1)^3L_w}{4\pi n_{a,0}^2R_{w,0}^3(1 - \xi_{{c}}^3)\Lambda(T_{{s}})(9\gamma-n-4)(\gamma+1)^2} \right]^{1/(3-2n)}. \end{equation}
This expression gives the characteristic shock radius at which the energy-conserving phase ceases to apply. For typical cooling curves, the dominant dependence enters through the drop in post-shock density with increasing radius, so the transition radius scales approximately as
$\mathrm{const.}^{1/(3-2n)}$
, with
$\mathrm{const.} \gt 1$
. The transition is therefore pushed to larger radii for steeper ambient density profiles. Since the interaction phase is also lengthened for steeper ambient density profiles, the cumulative effect is that steeper ambient density profiles result in a longer total non-radiative stage, delaying the onset of the second and third stages outlined in Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977).
Furthermore, a critical slope occurs at
\begin{equation} n_{\mathrm{crit}} = \frac{3}{2},\end{equation}
above which the flow is not expected to ever become dynamically dominated by radiative losses. This is particularly relevant for winds embedded in Bondi-like accretion flows, as such systems are still expected to transition to the energy-conserving regime since
$n \lt 2$
, while steeper systems where
$n \geq 2$
are expected to remain interaction-dominated in perpetuity.
4.2. Observational implications
For a sufficiently overdense outflow embedded in an ambient medium with a relatively steep density gradient, our results demonstrate that the interaction phase can occupy a substantial fraction of the total non-radiative lifetime of the outflow, and in some cases may be the hydrodynamic regime most relevant to the observed emission. A useful estimate for the duration of this regime is given by Equation (10); since
$\Delta_{{s}} / \kappa_{{r}}$
is of order unity, this equation can be re-expressed as
\begin{equation} t_{\mathrm{ss,end}} \sim \left[ \frac{v_w^{(n-3)}\dot{M}}{4\pi\rho_{a,0} R_{w,0}^n} \right]^{1/(2-n)},\end{equation}
where
$R_{w,0}$
again represents the characteristic radius at which the freely expanding wind becomes established and begins driving the double-shock structure. This transition occurs for
$n \lt 2$
if the wind remains active indefinitely.
On the other hand, realistic winds are launched over a finite time interval that depends on the driving mechanism, resulting in ejecta shells of finite thickness. The dynamics of the shocked regions will then deviate from this framework once the reverse shock reaches the interior edge of the shell, which occurs at a time
\begin{equation} t_{\Delta} \sim (t_{\mathrm{dyn}} + t_{w,e}) f^{1/(4-n)},\end{equation}
where
$t_{w,e}$
is the duration of the wind. This could be particularly important in scenarios where
$t_{w,e}$
is much shorter than relevant observing windows.
In the following subsections, we estimate the duration
$t_{\mathrm{ss,end}}$
of the interaction phase for various astrophysical phenomena. For TDEs and LFBOTs specifically, depending on how long the wind is sustained, the reverse-shock crossing time may be shorter than
$t_{\mathrm{ss,end}}$
, though perhaps still long enough to be observationally relevant.
4.2.1. Massive-star winds
For an O-type stellar wind, Falle (Reference Falle1975) and Weaver et al. (Reference Weaver, McCray, Castor, Shapiro and Moore1977) estimated the duration of the entire non-radiative phase to be
$\sim 2\times 10^3\,$
yr, and we obtain a comparable number using similar parameters
$v_w=2\,000\,\mathrm{km}\,{s}^{-1}$
,
$\dot{M}=10^{-6}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$
, and
$n_{a,0}=1\,\mathrm{cm}^{-3}$
, and
$\Lambda = 8\times 10^{-22}$
in Equation (28), yielding
$t_{\mathrm{rad}}\sim 2112\,$
yr. Adopting a ‘wind-launch’ radius of
$R_{w,0}\sim 10^{12}\,\mathrm{cm}$
, consistent with standard calibrations for Galactic O stars (Martins et al. Reference Martins, Schaerer and Hillier2005), we estimate an inner stellar wind density of
$\rho_{w,0} \sim 10^{-14}\,\mathrm{g}\,\mathrm{cm}^{-3}$
, which corresponds to
$f \sim 10^{10}$
. Since
$f \gg 1$
, the early wind evolution should be interaction-dominated, and we can estimate the duration of the interaction-specific component of this regime using the same parameters,
\begin{align} t_{\mathrm{ss,end,O}} &\sim 19.4\,\mathrm{yr} \, \times \, \left( \frac{v_w}{{2\,000\,\mathrm{km}\,{\mathrm{s}^{-1}}}} \right)^{-3/2} \nonumber\\& \quad \times \left( \frac{\dot{M}}{10^{-6}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}} \right)^{1/2}\times \left( \frac{n_{a,0}}{1\,\mathrm{cm}^{-3}} \right)^{-1/2}. \end{align}
Taking
$R_{\mathrm{ss,end}} \equiv v_wt_{\mathrm{ss,end}}$
, this corresponds to
$R_{\mathrm{ss,end}}=1.73\times 10^{17}\,$
cm, equivalent to the transition radius predicted in Falle (Reference Falle1975), which was obtained by solving for the radius at which the unshocked wind density is equivalent to that of the unshocked ambient medium. This is consistent with our numerical result that the transition occurs when the instantaneous overdensity has declined to only
$f(R_{{c}}) \sim$
a few (Figure 4), i.e. when the wind is no longer overwhelmingly denser than the ambient medium.
4.2.2. White dwarf winds
Another application of this framework is the white-dwarf-wind-driven remnant Pa 30 (Gvaramadze et al. Reference Gvaramadze, Gräfener, Langer, Maryeva, Kniazev, Moskvitin and Spiridonova2019; Oskinova et al. Reference Oskinova, Gvaramadze, Gräfener, Langer and Todt2020), the likely remnant of SN 1181 (Ritter et al. Reference Ritter, Parker, Lykou, Zijlstra, Guerrero and Le Dû2021). The central white dwarf appears to drive a very fast wind,
$v_w \geq 10^4\,\mathrm{km}\,{s}^{-1}$
, and Coughlin et al. (Reference Coughlin, Salvesen and Pasham2026) argue that this wind interaction is central to the remnant’s unusual radial filamentary morphology. Although we only consider spherically symmetric winds here, and therefore cannot describe the Rayleigh–Taylor growth invoked in that work, their interpretation is nevertheless embedded in the same interaction-phase framework considered here. Filament growth shuts off once the wind and ambient densities become comparable near the contact discontinuity, marking the end of the strongly overdense interaction regime. Coughlin et al. (Reference Coughlin, Salvesen and Pasham2026) use this same framework to infer a characteristic duration of
$\sim1-10\,\mathrm{yr}$
for the filament-growth phase in Pa 30, placing that system naturally within the broader class of interaction-dominated outflows discussed here.
4.2.3. Galactic winds
As highlighted in earlier sections, these interaction timescales increase for steeper n. Large-scale galactic outflows interact with the circumgalactic medium (CGM), adhering approximately to a power-law density profile of index
$n=1-2$
(Li et al. Reference Li, Bregman, Wang, Crain, Anderson and Zhang2017; Bregman et al. Reference Bregman, Anderson, Miller, Hodges-Kluck, Dai, Li, Li and Qu2018; Singh et al. Reference Singh, Majumdar, Nath and Silk2018). For such an outflow, we take
$R_{w,0}$
to be the gas disk scale height
$\sim$
$100\,\mathrm{pc}$
, i.e. the characteristic radius at which clustered supernova-driven bubbles break out of the disc and begin venting into the halo as a galactic wind (e.g. Chevalier & Clegg Reference Chevalier and Clegg1985; Fielding et al. Reference Fielding, Quataert and Martizzi2018). Adopting a wind velocity
$v_w = 100\,\mathrm{km}\,{s}^{-1}$
and mass flux
$\dot{M}=10^{-1}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$
(e.g. Fox et al. Reference Fox, Richter, Ashley, Heckman, Lehner, Werk, Bordoloi and Peeples2019), as well as a local density normalisation
$n_{a,0} = 10^{-2}\,\mathrm{cm}^{-3}$
, we estimate an initial wind overdensity
$f \sim$
tens. This is a somewhat modest overdensity, implying a relatively short-lived interaction phase – however, the steepness of the ambient density profile counteracts this.
To estimate the duration, we adopt a fiducial index of
$n=3/2$
, giving
\begin{align} t_{\mathrm{ss,end,gal}} \sim 9.75\times 10^{8}\,\mathrm{yr} \, \times \, \left( \frac{R_{w,0}}{100\,\mathrm{pc}} \right)^{-3} \\\times \left( \frac{v_w}{{100\,\mathrm{km}\,{\mathrm{s}^{-1}}}} \right)^{-3}\\\times \left( \frac{\dot{M}}{10^{-1}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}} \right)^{2}\times\left( \frac{n_{a,0}}{10^{-2}\,\mathrm{cm}^{-3}} \right)^{-2},\end{align}
though this result is highly sensitive to the normalisation of the CGM density profile, which is less constrained than
$v_w$
and
$\dot{M}$
. Nevertheless, this estimate seems to indicate that the interaction phase is highly relevant for galactic winds.
4.2.4. TDE outflows
An
$n=3/2$
scenario is also typical of the expected environments of some TDE outflows, given the natural stratification of accreting material around the central black hole (Bondi Reference Bondi1952) (though radio modelling often favours even steeper ambient profiles; e.g. Alexander et al. Reference Alexander, Berger, Guillochon, Zauderer and Williams2016; Cendes et al. Reference Cendes, Alexander, Berger, Eftekhari, Williams and Chornock2021). We use
$R_{w,0} = 10^{13}\,\mathrm{cm}$
, consistent with the simplified scenario presented in Strubbe & Quataert (Reference Strubbe and Quataert2009) where the outflow launch radius is
$\sim$
twice the pericenter distance of the stellar orbit around the black hole. Extrapolating typical densities of
$n_a \sim 100\,\mathrm{cm}^{-3}$
at
$\sim$
$0.01\,\mathrm{pc}$
(Alexander et al. Reference Alexander, Berger, Guillochon, Zauderer and Williams2016) back to
$10^{13}\,\mathrm{cm}$
for
$n=3/2$
yields
$n_{a,0} \sim 10^7\,\mathrm{cm}^{-3}$
. For total disruptions, peak fallback accretion rates range from a few
$\mathrm{M}_{\odot}$
to tens of
$\mathrm{M}_{\odot}$
per year (Bandopadhyay et al. Reference Bandopadhyay2024; Fancher et al. Reference Fancher, Bandopadhyay, Coughlin and Nixon2025), but not all of this material will be ejected in the outflow and the peak fallback rate is (clearly) greater than the average fallback rate. Partial disruptions can also yield substantially smaller fallback rates (e.g. Guillochon & Ramirez-Ruiz Reference Guillochon and Ramirez-Ruiz2013; Miles et al. Reference Miles, Coughlin and Nixon2020; Bandopadhyay et al. Reference Bandopadhyay, Coughlin and Nixon2026). Taking
$\dot{M} = 10^{-1}\,\mathrm{M}_{\odot}\,\mathrm{yr}$
and
$v_{w} = 10^4\,\mathrm{km}\,{s}^{-1}$
, we estimate an initial overdensity
$f \sim 10^5$
, indicating these outflows are heavily interaction-dominated.
For the duration of the interaction phase, we then estimate
\begin{align} t_{\mathrm{ss,end,TDE}} \sim 2.86 & \times 10^7\,\mathrm{yr} \, \times \, \left( \frac{R_{w,0}}{10^{13}\,\mathrm{cm}} \right)^{-3} \left( \frac{v_w}{{10^4\,\mathrm{km}\,{\mathrm{s}^{-1}}}} \right)^{-3} \nonumber\\& \times \left( \frac{\dot{M}}{10^{-1}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}} \right)^{2}\left( \frac{n_{a,0}}{10^7\,\mathrm{cm}^{-3}} \right)^{-2}.\end{align}
Of course, the accretion-stratified structure of the ambient medium would not extend far enough to support such a long interaction phase, so in reality these systems would be limited by the radial extent of the surrounding gas profile. Furthermore, the outflow may shut off much earlier than
$t_{\mathrm{ss,end}}$
. If the wind is tied to the super-Eddington phase of the TDE accretion flow, then the duration of such a wind can be reasonably estimated to be
$\sim$
months (Strubbe & Quataert Reference Strubbe and Quataert2011; Wu et al. Reference Wu, Coughlin and Nixon2018). Conservatively adopting a duration of 1 month,
$t_{w,e} \gg t_{\mathrm{dyn}}$
based on the same parameters used to estimate
$t_{\mathrm{ss,end}}$
. Thus, the time at which the reverse shock is expected to fully cross the finite-width ejecta shell is
\begin{equation} t_{\Delta} \sim 8.24\,\mathrm{yr}\,\times \left( \frac{t_{w,e}}{2.6\times 10^6\,{s}} \right) \times \left( \frac{f}{10^5} \right)^{2/5}.\end{equation}
While this timescale is much shorter than
$t_{\mathrm{ss,end}}$
, radio observations of TDE outflows are carried out within
$\sim$
months to years after the initial disruption. As such, these systems are expected to remain interaction-dominated throughout much of their observationally relevant evolution.
4.2.5 LFBOTs
The analysis here may also be relevant to LFBOTs, whose radio and X-ray properties in at least some well-studied cases suggest fast outflows interacting with dense circumstellar material (e.g. Margutti et al. Reference Margutti2019; Ho et al. Reference Ho2019). Recent modelling of the LFBOT AT2024wpp suggests that the circumstellar density structure approximately adheres to a power law where n is as high as
$\sim$
3 (Nayana et al. Reference Nayana2025; Perley et al. Reference Perley2026). Such steep power-laws are consistent with modelling assumptions for other similar events such as AT2018cow (Ho et al. Reference Ho2019), where it is often assumed that the emission is coming from an outflow propagating through a wind-like (
$n \sim 2$
) circumstellar medium (Guarini et al. Reference Guarini, Tamborra and Margutti2022; Margutti et al. Reference Margutti2019). These systems are therefore especially compelling from the standpoint of the present work, because they occupy precisely the regime in which the interaction phase is most prolonged.
Indeed, for
$n \geq 2$
the flow never relaxes to an energy-conserving similarity state; instead, the interaction-dominated structure persists indefinitely. In these scenarios, the ambient density declines more rapidly with radius than that of the freely expanding wind, and hence the reverse shock never propagates to small radii relative to the contact discontinuity. Therefore, if the circumstellar media for these events truly are wind-like or steeper, then these systems could remain in the interaction regime perpetually. Direct application of the interaction solutions to emission modelling in such transients is therefore a particularly natural next step, though outside the scope of this work. As discussed in prior sections, the relevance of these interaction solutions could potentially be limited by the duration of the outflow. However, the properties of these systems are not very well-constrained, so it is difficult to assess the importance of this.
5. Summary and conclusions
We presented a unified study of the early adiabatic evolution of wind-driven bubbles expanding into power-law ambient media, focusing on the relationship between the interaction-dominated early phase and the later energy-conserving phase, each of which has a self-similar solution. Using analytic scalings together with one-dimensional shock-capturing hydrodynamic simulations, we quantified both the timescales over which the flow relaxes toward these asymptotic states and the degree to which the corresponding solutions reproduce the full shocked structure. Our principal conclusions are:
-
1. In the limit of highly overdense ejecta, the interaction solutions provide an excellent description of the shocked flow at early times. For all density profiles considered here,
$0 \leq n \leq 2$
, the numerical solutions rapidly relax toward the interaction-phase profiles when
$f \gg 1$
, where f is the initial ratio of the wind density to the ambient density. -
2. For
$n \leq 2$
, both the analytic estimate derived in Section 2 and the numerical results in Section 3.2 show that the duration of the interaction regime increases strongly with both f and the ambient density power-law slope, following
$t_{\mathrm{ss,end}} \propto f^{1/(2-n)}$
. Physically, the interaction solutions cease to apply once the shocked-wind shell (i.e. the fluid bounded by the reverse shock and the contact discontinuity) is no longer thinner than the forward-shocked ambient shell, which in practice occurs when the instantaneous overdensity has fallen to
$f(R_{{c}}) \sim$
a few. At later times, the flow approaches the energy-conserving phase, with the shock positions converging toward the expected power-law scalings. -
3. The case
$n=2$
is qualitatively distinct. In this limit, the predicted transition timescale diverges, and our simulations confirm that the system does not relax to a distinct energy-conserving similarity state. Rather, the interaction-dominated structure persists throughout the evolution, which is consistent with the fact that the ratio of the reverse shock radius to the contact discontinuity radius does not approach zero at late times. More generally, our results imply that steep ambient density gradients lead to prolonged interaction-dominated evolution and delay – or even entirely prevent – the onset of the later energy-conserving regime. -
4. The transition between these regimes is itself dynamically extended. Near
$t \sim t_{\mathrm{ss,end}}$
, neither the interaction nor the energy-conserving solution provides an accurate representation of the full shocked structure, especially in the velocity and pressure profiles. Even under a conservative error criterion, the interval over which the system is not particularly well described by either asymptotic solution can extend for several orders of magnitude in time beyond the nominal transition timescale itself. This indicates that it may generally be insufficient to model the evolution as an abrupt hand-off from one self-similar regime to the other, and we developed a set of closed-form, analytic equations for the discontinuity positions as functions of time in Appendix C. To adequately account for the evolution of the full radial profiles, there is room for future work here in exploring the higher-order corrections to the interaction solutions. -
5. These results have important implications for engine-driven transients. Although the interaction phase is shorter than the total non-radiative lifetime for classical O-star wind bubbles, it can still persist for decades. Galactic winds in stratified halos, TDE outflows, and at least some LFBOTs may remain interaction-dominated throughout much of their observed lifetimes.
There are several natural extensions of this work. On the theoretical side, it would be useful to incorporate higher-order corrections to the interaction solutions and to develop a more quantitative description of the broad intermediate regime between the two asymptotic states. On the astrophysical side, coupling these hydrodynamic solutions to explicit radiation post-processing would allow more direct comparisons with radio, X-ray, and other observational diagnostics in systems such as TDEs and LFBOTs. Such extensions would help clarify when the interaction phase is merely an early hydrodynamic transient and when it is, in practice, the dominant regime governing the observable evolution of the source.
Acknowledgements
We thank the anonymous referee for their useful and constructive comments. B.A. and E.R.C. acknowledge support from NASA through the Astrophysics Theory Program, grant 80NSSC24K0897, and through Chandra Award Number 25700383 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.
Appendix A. 1D SRHD code
We solve the one-dimensional special relativistic hydrodynamics equations in spherical coordinates using a moving-mesh finite volume method. The equations are written in conservative form,
\begin{equation} \frac{\partial \mathbf{U}}{\partial t} + \frac{1}{r^2} \frac{\partial (r^2 \mathbf{F})}{\partial r} = \mathbf{S},\end{equation}
where the source term
$\mathbf{S}=(0, 2p/r, 0)^T$
accounts for geometric factors arising from spherical coordinates, following Zhang & MacFadyen (Reference Zhang and MacFadyen2006), and
$\mathbf{U} = (D, S, E)^T$
contains the lab-frame mass density
$D = \Gamma\rho$
, momentum density
$S = \Gamma^2 \rho h u$
, and energy density
$E = \Gamma^2 \rho h - p - D$
. Here,
$\rho$
is the comoving mass density,
$u = \Gamma\beta$
is the radial four-velocity, p is the gas pressure,
$\Gamma = \sqrt{1 + u^2}$
is the Lorentz factor, and
$h = 1 + \hat{\gamma} p / [(\hat{\gamma}-1)\rho]$
is the specific enthalpy for an ideal gas with adiabatic index
$\hat{\gamma}$
. The corresponding flux vector is
\begin{equation} \mathbf{F} = \begin{pmatrix} Dv \\ Sv + p \\ Ev + pv \end{pmatrix},\end{equation}
where
$v = u/\Gamma$
is the three-velocity.
The spatial discretisation uses piecewise linear reconstruction to achieve second-order accuracy in smooth regions. Slopes are limited using the minmod limiter,
\begin{equation} \sigma = \frac{1}{4}|\mathrm{sgn}(a) + \mathrm{sgn}(b)| \cdot (\mathrm{sgn}(a) + \mathrm{sgn}(c)) \cdot \min(|a|, |b|, |c|),\end{equation}
where
$a = \theta(q_i - q_{i-1})$
,
$b = (q_{i+1} - q_{i-1})/2$
,
$c = \theta(q_{i+1} - q_i)$
, and
$\theta = 1.5$
controls the limiter’s compressiveness. This prevents spurious oscillations near discontinuities while maintaining accuracy in smooth flow regions.
Numerical fluxes at cell interfaces are computed using the HLLE approximate Riemann solver. Given left and right states at an interface, the relativistic wave speeds are estimated as
\begin{equation} a^{\pm}_{L,R} = \frac{v_{L,R} \pm c_{s}}{1 \pm v_{L,R} c_{s}},\end{equation}
where
$c_s = \sqrt{\hat{\gamma} p / \rho h}$
is the relativistic sound speed. The bounding wave speeds
$a^- = \min(a^-_L, a^-_R)$
and
$a^+ = \max(a^+_L, a^+_R)$
define the extent of the Riemann fan, and for a face moving with speed
$v_F$
, the ALE-HLLE numerical flux is given by
\begin{equation} \hat{\mathbf{F}} = \begin{cases} \mathbf{F}_L - v_F\mathbf{U}_L\,, & v_F \leq a^-\,,\\ \mathbf{F}^* - v_F\mathbf{U}^*\,, & a^- \lt v_F \lt a^+\,,\\ \mathbf{F}_R - v_F\mathbf{U}_R\,, & v_F \geq a_+\,, \end{cases}\end{equation}
with
\begin{equation} \mathbf{F}^* = \frac{a^+ \mathbf{F}_L - a^- \mathbf{F}_R + a^+ a^- (\mathbf{U}_R - \mathbf{U}_L)}{a^+ - a^-}\end{equation}
and
\begin{equation} \mathbf{U}^* = \frac{a^+\mathbf{U}_R - a^-\mathbf{U}_L - (\mathbf{F}_R - \mathbf{F}_L)}{a^+ - a^-}.\end{equation}
Recovery of primitive variables
$(\rho, u, p)$
from the conserved quantities (D, S, E) requires solving a nonlinear equation for the pressure. We employ Newton-Raphson iteration with a convergence tolerance of
$10^{-12}(D + E)$
.
A key feature of our method is the moving-mesh formulation, wherein the computational domain is bounded by tracked discontinuities rather than fixed spatial boundaries. The global domain is divided into two subdomains, each consisting of N uniformly-spaced grid cells. The outer edge of the outer subdomain coincides with the forward shock, while the inner edge of the inner subdomain tracks the reverse shock, and the boundary separating the two subdomains follows the contact discontinuity separating the shocked ejecta from the shocked ambient medium. Shock positions are evolved using speeds determined from the exact Rankine-Hugoniot jump conditions, given the upstream (unshocked) and downstream (shocked) states. The upstream state is supplied as a boundary condition representing either the unshocked ambient medium (for the forward shock) or the freely-expanding ejecta (for the reverse shock). At the contact discontinuity, the interface moves at the common fluid velocity shared by the adjacent zones.
Time integration is performed using a third-order strong stability-preserving Runge–Kutta (SSP-RK3) scheme:
\begin{align} \mathbf{U}^{(1)} &= \mathbf{U}^n + \Delta t \, \mathscr{L}\,(\mathbf{U}^n), \end{align}
\begin{align} \mathbf{U}^{(2)} &= \tfrac{3}{4} \mathbf{U}^n + \tfrac{1}{4} \mathbf{U}^{(1)} + \tfrac{1}{4} \Delta t \, \mathscr{L}\,(\mathbf{U}^{(1)}), \end{align}
\begin{align} \mathbf{U}^{n+1} &= \tfrac{1}{3} \mathbf{U}^n + \tfrac{2}{3} \mathbf{U}^{(2)} + \tfrac{2}{3} \Delta t \, \mathscr{L}\,(\mathbf{U}^{(2)}),\end{align}
where
$\mathscr{L}\,$
denotes the spatial discretisation operator. The timestep is determined by a CFL condition,
$\Delta t = C_{\mathrm{cfl}} \, \Delta r_{\mathrm{min}} / \lambda_{\mathrm{max}}$
, where
$\lambda_{\mathrm{max}}$
is the maximum signal speed and
$C_{\mathrm{cfl}} = 0.3$
.
Appendix B. Numerical convergence of simulations
The shock-capturing method outlined in Appendix A allows for a high degree of numerical accuracy with relatively few zones. All fluid profiles are generally very well converged even for resolutions as low as
$N=200$
zones per subdomain, though the behaviour near the contact is extremely challenging to resolve. Because this study is primarily concerned with broader agreement with existing self-similar solutions and associated transitional timescales, as well as shock dynamics specifically, we adopt a baseline resolution of
$N=800$
while acknowledging that the density cusps adjacent to the contact are slightly under-resolved in both the interaction and energy-conserving regimes. The under-resolved features near the contact can be seen in the insets in Figure B1; outside of this narrow region, however, the simulated profiles are broadly numerically converged even for lower N.
A numerical convergence test for the profiles of the fluid mass density, velocity, and pressure vs. self-similar coordinate
$\eta$
at six different resolutions. This simulation was initialised with a wind velocity
$v_w=0.01$
and an overdensity
$f=10^3$
, impacting a uniform (
$n=0$
) ambient medium. These snapshots were taken at
$t/t_{\mathrm{dyn}}=10$
(left) and
$t/t_{\mathrm{dyn}}=10^6$
(right) to explore convergence in both the interaction and energy-conserving regimes. Broadly, the profiles are extremely well-converged with resolutions even as low as
$N=200$
zones per subdomain, though the immediate vicinity of the contact (within
$\sim$
10 zones) remains under-resolved even for large N.
Appendix C. Analytic functions for discontinuity positions
The time-dependent positions of all three relevant discontinuities – the forward shock, contact discontinuity, and reverse shock – have analytic expressions in both the interaction and energy-conserving limits. We can smoothly connect these expressions to produce a set of formulae for these positions that can be explicitly evaluated at any time t, given the necessary parameters. We adopt the smooth interpolation function
\begin{equation} S(t) = \frac{t^{2p}}{t^{2p} + qt_{\mathrm{ss,end}}^{2p}},\end{equation}
where
$t_{\mathrm{ss,end}}$
is the transition timescale in Equation (10),
$q=1.11-0.64n$
, and
$p=1.16-0.53n$
(these latter two values were obtained via fitting). Using this function to conjoin the analytic expressions in each limiting regime, we obtain
\begin{align} R_{s}(t) &= [1-S(t)](R_0 + V_{\mathrm{ej}}t)(1+\Delta_{s})\nonumber\\ &\quad \times \left[ 1 + \left( \frac{\kappa_{s}}{1 + \Delta_{s}} - \kappa_{c} \right)\frac{1}{\sqrt{f}} \left( \frac{R_0 + V_{\mathrm{ej}}t}{R_0} \right)^{(2-n)/2} \right]\nonumber\\ &\quad + S(t)K\big( 2\pi R_0^{2-n}V_{\mathrm{ej}}^3 \,f \big)^{1/(5-n)}t^{3/(5-n)} \end{align}
\begin{align} R_{c}(t) &= [1-S(t)](R_0 + V_{\mathrm{ej}}t)\nonumber\\ & \quad \times \left[ 1 - \frac{\kappa_{c}}{\sqrt{f}} \left( \frac{R_0 + V_{\mathrm{ej}}t}{R_0} \right)^{(2-n)/2} \right]\nonumber\\ & \quad + S(t)\xi_{c}K\big( 2\pi R_0^{2-n}V_{\mathrm{ej}}^3 \,f \big)^{1/(5-n)}t^{3/(5-n)} \end{align}
\begin{align} R_{r}(t) &= [1-S(t)](R_0 + V_{\mathrm{ej}}t)\nonumber\\ &\quad\times \left[ 1 - \frac{(\kappa_{c} + \kappa_{r})}{\sqrt{f}} \left( \frac{R_0 + V_{\mathrm{ej}}t}{R_0} \right)^{(2-n)/2} \right]\nonumber\\ &\quad + S(t)\sqrt{\frac{\gamma+1}{\gamma-1}\left( 1 - \frac{n+4}{9\gamma} \right)\frac{3}{(5-n)V_{\mathrm{ej}}}}\nonumber\\ &\quad \times(\xi_{c}K)^{3/2}\big( 2\pi R_0^{2-n}V_{\mathrm{ej}}^3 \,f \big)^{3/(10-2n)}t^{(n+4)/(10-2n)}. \end{align}
Here,
$f_0$
,
$R_0$
, and
$V_{\mathrm{ej}}$
are initial conditions. All other parameters are extracted from either the self-similar solutions presented in Coughlin (Reference Coughlin2024) or by numerically integrating to obtain the contact position in the energy-conserving limit. Values for these parameters are tabulated in Table C1 for several different values of n and a fixed value of
$\gamma=5/3$
. Note that
$n=2$
is not included here since the system will remain in the interaction-dominated state indefinitely, and would therefore be described solely by the solutions presented in Coughlin (Reference Coughlin2024). To further minimise the error in the transition window for the full radial profiles for the density, velocity, and pressure, higher-order terms can be included in these same interaction solutions.
A direct comparison between these equations and the simulation output is shown in Figures C1 and C2 for
$f=10^3$
and
$n=0,1.5$
.
Self-similarity parameters (
$\gamma=5/3$
).
Table C1. Long description
A table with 10 rows and 6 columns comparing different parameters for various values of n and a fixed value of K. The columns are labeled n, K, xi_c, Delta_s times 10, kappa_s times 10, kappa_c times 10, and kappa_f times 10. The table presents the following data: Row 1: n 0.00, K 0.8837, xi_c 0.8584, Delta_s times 10 0.9572, kappa_s times 10 0.3350, kappa_c times 10 3.9526, kappa_f times 10 0.4852. Row 2: n 0.25, K 0.8531, xi_c 0.8557, Delta_s times 10 1.0199, kappa_s times 10 0.3405, kappa_c times 10 4.2295, kappa_f times 10 0.5263. Row 3: n 0.50, K 0.8184, xi_c 0.8531, Delta_s times 10 1.0914, kappa_s times 10 0.3431, kappa_c times 10 4.5495, kappa_f times 10 0.5751. Row 4: n 0.75, K 0.7790, xi_c 0.8505, Delta_s times 10 1.1740, kappa_s times 10 0.3412, kappa_c times 10 4.9240, kappa_f times 10 0.6341. Row 5: n 1.00, K 0.7340, xi_c 0.8480, Delta_s times 10 1.2703, kappa_s times 10 0.3314, kappa_c times 10 5.3691, kappa_f times 10 0.7069. Row 6: n 1.25, K 0.6826, xi_c 0.8457, Delta_s times 10 1.3843, kappa_s times 10 0.3084, kappa_c times 10 5.9078, kappa_f times 10 0.7990. Row 7: n 1.50, K 0.6238, xi_c 0.8435, Delta_s times 10 1.5214, kappa_s times 10 0.2621, kappa_c times 10 6.5756, kappa_f times 10 0.9186. Row 8: n 1.75, K 0.5566, xi_c 0.8416, Delta_s times 10 1.6898, kappa_s times 10 0.1731, kappa_c times 10 7.4298, kappa_f times 10 1.0845.
f=105
t≪tdyn
t≪tdyn
t∼tdyn
n=0
L1
f≫1
f=105
n∈[0,12,1,32,2]
L1
L1
105
n∈[0,12,1,32,2]
f(Rc)≡f⋅(Rw,0/Rc)2−n
tss,end
f(Rc)∼1
n=2
tss,end′
f1/(2−n)
n=2
tss,end→∞
n=2
f=105
n=0
f(Rc)=9.227×102
f(Rc)=1.109
f(Rc)=4.843×10−5
f(Rc)∼1
t≪tss,end
t≫tss,end
t∼tss,end
η
vw=0.01
f=103
n=0
t/tdyn=10
t/tdyn=106
N=200
∼
γ=5/3
n=0
f=103
n=1.5
f=103
Open access
