# Galactic winds and bubbles from nuclear starburst rings

Dustin D. Nguyen Center for Cosmology and Astro-Particle Physics (CCAPP)
The Ohio State University, Columbus, Ohio, 43210, USA  Department of Physics,
The Ohio State University, Columbus, Ohio, 43210, USA  Todd A. Thompson Center
for Cosmology and Astro-Particle Physics (CCAPP)
The Ohio State University, Columbus, Ohio, 43210, USA  Department of Physics,
The Ohio State University, Columbus, Ohio, 43210, USA  Department of
Astronomy, The Ohio State University, Columbus, Ohio, 43210, USA

###### Abstract

Galactic outflows from local starburst galaxies typically exhibit a layered
geometry, with cool $10^{4}\,$K flow sheathing a hotter $10^{7}\,$K,
cylindrically-collimated, X-ray emitting plasma. Here, we argue that winds
driven by energy-injection in a ring-like geometry can produce this
distinctive large-scale multi-phase morphology. The ring configuration is
motivated by the observation that massive young star clusters are often
distributed in a ring at the host galaxy’s inner Lindblad resonance, where
larger-scale spiral arm structure terminates. We present parameterized three-
dimensional radiative hydrodynamical simulations that follow the emergence and
dynamics of energy-driven hot winds from starburst rings. We show that the
flow shocks on itself within the inner ring hole, maintaining high $10^{7}$ K
temperatures, whilst flows that emerge from the wind-driving ring unobstructed
can undergo rapid bulk cooling down to $10^{4}\,$K, producing a fast hot bi-
conical outflow enclosed by a sheath of cooler nearly co-moving material. The
hot flow is collimated along the ring axis, even in the absence of pressure
confinement from a galactic disk or magnetic fields. In the early stages of
expansion, the emerging wind forms a bubble-like shape reminiscent of the
Milky Way’s eROSITA and Fermi bubbles. The bubble is preceded by a fast
transient flow along the minor axis that can reach velocities usually
associated with AGN-driven winds ($\gtrsim 3000-10000\,\mathrm{km\,s^{-1}}$),
depending on the density of the surrounding medium. We discuss the physics of
the ring configuration, the conditions for radiative bulk cooling, and the
implications for future X-ray observations.

Active galaxies (17); Starburst galaxies (1570); Galaxy evolution (929);
Hydrodynamic simulations (767); Stellar feedback (1602); Astronomical
simulations (1857); Galactic winds (572)

††software: numpy , scipy , Cholla, pyatomdb , yt ††thanks: 2d and 3d movies
of the simulations can be found online at dustindnguyen.com.

##  1 Introduction

Galactic winds are important to galaxy evolution. They drive gas into the
circumgalactic medium, help truncate star formation, and set the mass-
metallicity relation (see reviews by Veilleux et al., 2005; Zhang et al.,
2018). Energy-driven galactic winds, powered by massive star supernovae (SNe)
and stellar winds, may generate the hot ($T\approx 10^{6}-10^{7}\,\mathrm{K}$)
gas seen in diffuse X-ray emission from local starbursts (Strickland et al.,
2004; Strickland & Heckman, 2007). While the hot gas kinematics have not yet
been directly measured, current models predict velocities of $\sim
500-2000\,\mathrm{km\,s^{-1}}$, depending on the mass-loading of the hot phase
(Chevalier & Clegg, 1985; Strickland & Heckman, 2009; Nguyen & Thompson,
2021). In the M82 superwind, other local starbursts, and high-redshift rapidly
star-forming galaxies, the cooler neutral and ionized gas is observed to be
moving at velocities of order $\simeq 500-1000\,\mathrm{km\,s^{-1}}$ (Shopbell
& Bland-Hawthorn, 1998; Shapiro et al., 2009). A motivating fact for this
paper is that in local systems where the different thermal components can be
spatially resolved, the hot and cool outflows exhibit a layered geometry, with
the more cylindrical hot gas surrounded by a cooler gas “sheath” of confining
medium (e.g., as seen in M82, NGC 253, NGC 3079, and others; Westmoquette et
al. 2009; Leroy et al. 2015; Hodges-Kluck et al. 2020).

A critical open issue in galactic wind physics is the acceleration mechanism
for the cool component. One well-known picture is that the cool phase is
composed of clouds outside of the star-forming core that are ram pressure
accelerated by the emergent hot supersonic wind. Simulations of this
interaction often lead to the destruction of clouds by hydrodynamical
instabilities (Cooper et al., 2009; Scannapieco & Brüggen, 2015; Schneider et
al., 2018). In some regimes of cloud parameters, the wind-cloud interaction
leads to cloud growth and acceleration (Gronke & Oh, 2018), although this
process has not yet been shown to operate in global galactic wind simulations
(Tanner et al., 2016; Schneider et al., 2020). Another possible solution is
that the cool phase condenses directly from the hot phase via radiative
instability (Wang, 1995; Silich et al., 2004; Thompson et al., 2016). The
condition for rapid cooling is set primarily by the physical size of the wind-
driving region and the hot wind mass-loading rate.

Here, we argue that the characteristic layered outflow geometry seen in nearby
starbursts may be produced by galactic winds driven from a ring-like geometry.
Rings of massive star clusters form as a result of gas inflow to the inner
Lindblad resonance (Athanassoula, 1992; Binney & Tremaine, 2008). These rings
often dominate the total star formation rate within their host galaxies
(Mazzuca et al., 2008). For example, the circum-nuclear ring of the archetypal
nearby starburst galaxy M82 is well studied. Nakai et al. (1987) discovered a
double-lobed distribution of CO, which is interpreted as a molecular torus.
Hughes et al. (1994) studied thermal reradiated emission from dust in sub mm
wavelengths and found that the intensity profile can be modeled as an edge-on
nuclear starburst torus. Weiß et al. (2001) provided a multi-transition
analysis of the molecular gas and showed the regions of most violent star
formation are confined to the molecular lobes, and are arranged in a toroidal
topology around the nucleus. Star-forming rings are also observed in the
nearby starburst NGC 253 (Arnaboldi et al., 1995; Leroy et al., 2018), other
local galaxies (see Buta & Crocker, 1993; Böker et al., 2008; Comerón et al.,
2010; Leroy et al., 2021), and the Milky Way galaxy (see Henshaw et al.,
2022).

3d ring simulations have focused on the dynamical processes associated with
stellar populations and bars, and are typically carried out with parameters
characteristic of the Galaxy’s Central Molecular Zone (Armillotta et al.,
2019; Tress et al., 2020). Moon et al. (2021) found that individual SN shocks
can merge together to create large hot outflows resembling galactic winds.
These works motivate a focused study of the general properties, dynamics, and
morphology of energy-driven winds from star-forming rings. Here, we posit the
ring structure and focus on the evolution of large-scale outflows that can be
produced from parameterized energy and mass injection within the ring volume.
As nuclear rings are observed to substantially vary in size, ranging from
diameters of tens of pc to several kpc (Buta & Crocker, 1993; Comerón et al.,
2010), we consider several different ring sizes and mass-loading rates. In
addition, we consider two ring geometries: the first is an ideal azimuthally
uniform ring, which is used to explore the basic structures that can be
obtained from this configuration. The second is a series of uniform spheres
distributed in a ring-like geometry that is meant to mimic individual super
star clusters formed at the inner Lindblad resonance (Elmegreen, 1994). For
both sets of simulations, we find that a ring-like distribution of sources can
produce an overall wind morphology akin to that seen in observations, where
the hot and cool phases are layered, and where the hot flow is confined to a
more cylindrical structure.

Figure 1: Schematic of the nuclear starburst ring model. In this model,
supernovae inject energy and mass into a ring geometry, leading to a steady
state flow that has an initial temperature within the ring of $\sim
10^{7}\,$K, which accelerates to high velocity as it leaves the ring volume.
For high mass-loading, the wind emanating from the outer-ring can become
radiative and quickly cool to $10^{4}\,$K at the cooling radius,
$r_{\mathrm{cool}}$ (see Section 3). Within the ring hole, along the minor
axis, the flow shocks on itself, thermalizing its kinetic energy and
maintaining high temperature. As the hot flow emerges from the inner ring
region, it is collimated by the surrounding cooler flow.

##  2 Model and Methods

In Fig. 1 we illustrate the model problem. We assume SNe deposit energy and
mass at total rates $\dot{E}_{\mathrm{T}}$ and $\dot{M}_{\mathrm{T}}$,
uniformly throughout a symmetric ring with volume $V$ (see Appendix A). The
controlling parameters of the problem are the source terms, the energy and
mass-loading rates, the ring width $R_{A}$, and the ring radius $R_{B}$. These
quantities are shown in Table 1 and are chosen for Simulations ‘A’ and
‘$\texttt{A}_{\star}$’ such that the unobstructed flow emanating from the
outer ring undergoes bulk cooling.

We use the GPU accelerated radiative hydrodynamics code Cholla (Schneider &
Robertson, 2015) to run a suite of simulations to study the properties of the
starburst ring model. The simulations and parameters are summarized in Table
1. The simulation box is $10\times 10\times 15\,$kpc for the
$\hat{x}\times\hat{y}\times\hat{z}$ directions, respectively. There are 2 high
resolution simulations, ‘A’ and ‘$\texttt{A}_{\star}$’, with $768\times
768\times 1152$ cells in each spatial dimension, giving a cell resolution of
$\Delta x\simeq 13$ pc. These two serve as the fiducial simulations for the
study. The uniform ring in Simulation ‘A’ (see Section 3.1) has a ring width
of $R_{A}=25\,$pc. Simulation ‘$\texttt{A}_{\star}$’ (see Section 3.2) has
injection regions modeled as 8 identical spherical star clusters distributed
in a ring. Each star cluster has a radius of $R_{A}=25\,$pc. The other 10
runs, Simulations ‘B…’ and ‘C…’, are run with $512\times 512\times 768$ cells,
giving a cell resolution of $\Delta x\simeq 19.5$ pc. These lower resolution
simulations explore the effects to the flow when the ring geometry, mass-
loading rate, or ambient circumgalactic gas density is varied. Within the
starburst volume, we inject energy and mass at rates
$\dot{E}_{\mathrm{T}}=\alpha\times 3.1\times
10^{41}\times(\dot{M}_{\mathrm{SFR}}/\mathrm{M_{\odot}\
yr^{-1}})\,[\mathrm{ergs\ s^{-1}}]$ and
$\dot{M}_{\mathrm{T}}=\beta\times\dot{M}_{\mathrm{SFR}}$, where $\alpha$,
$\beta$, and $\dot{M}_{\mathrm{SFR}}$ are the thermalization efficiency, mass-
loading rate, and star formation rate, respectively. The cut-off in the mass
and energy injection rates (see equations A1 and A2) at the ring edge yields a
characteristic solution similar to Chevalier & Clegg (1985) and Wang (1995),
such that a sonic point forms at the edge of the ring volume.

The Cholla simulations (Schneider & Robertson, 2015) use a piece-wise linear
interface reconstruction, an HHLC Riemann solver, and an unsplit Van Leer
integrator, and a cooling curve that is a piece-wise parabolic fit to a Cloudy
curve assuming collisional ionization equilibrium for a solar metallicity gas
(see Schneider & Robertson, 2018). We approximately account for background UV
heating through the implementation of a cooling floor at
$T_{\mathrm{floor}}=10^{4}\,$K. The cooling everywhere is assumed to be
optically-thin. Because we do not include a disk in these simulations, the
wind from the ring escapes unobstructed, even along the azimuthal plane. These
flows would be blocked by the gas disk (see Fig. 1) in a more realistic
simulation setup, which we save for future work. Simulations ‘A’,
‘$\texttt{A}_{\star}$’, ‘B…’ begin with a constant static background of
temperature $T_{\mathrm{bkgd.}}=10^{6}\,$K and density
$n_{\mathrm{bkgd.}}=10^{-5}\,\mathrm{cm^{-3}}$. Simulations ‘C…’ also start
with a static background $10^{6}\,$K, however the background density is varied
for each simulation (see Table 1). At the start of the simulation, energy and
mass is instantaneously deposited into the uniform (‘A’), or distributed
cluster (‘$\texttt{A}_{\star}$’), ring geometry.

Simulations ‘A’ and ‘$\texttt{A}_{\star}$’ produce radiative flows. As we
discuss in Section 3, the condition for bulk cooling is set primarily by the
mass-loading term $\beta$ and the physical width of the wind-driving region
$R_{A}$, as in the analogous spherical models (Thompson et al., 2016)). Models
with lower mass-loading and/or thicker ring widths (see Appendices C and D) do
not produce radiative flows, and a hot wind emerges in all directions.

Simulation Name | Ring Type | $\alpha$ | $\beta$ |  $R_{A}\,$[pc] |  $R_{B}\,$[kpc] | $\Delta x\,\mathrm{[pc]}$
---|---|---|---|---|---|---
A | uniform | 1.0 | 0.6 | 25 | 1.0 | 13.0
A⋆ | non-uniform | 1.0 | 0.6 | 25 | 1.0 | 13.0
B0 | uniform | 1.0 | 0.6 | 100 | 1.0 | 19.5
B1 | uniform | 1.0 | 0.6 | 100 | 0.5 | 19.5
B2 | uniform | 1.0 | 0.6 | 100 | 2.0 | 19.5
B3 | uniform | 1.0 | 0.6 | 50 | 1.0 | 19.5
B4 | uniform | 1.0 | 0.6 | 200 | 1.0 | 19.5
B5 | uniform | 1.0 | 1.1 | 100 | 1.0 | 19.5
B6 | uniform | 1.0 | 0.2 | 100 | 1.0 | 19.5
C0 | uniform | 1.0 | 0.1 | 100 | 1.0 | 19.5
C1 | uniform | 1.0 | 0.1 | 100 | 1.0 | 19.5
C2 | uniform | 1.0 | 0.1 | 100 | 1.0 | 19.5
Table 1: Summary of the simulation parameters. Simulations ‘A’ and
‘$\texttt{A}_{\star}$’ serve as the fiducial models of the study. Simulations
‘A’, ‘$\texttt{A}_{\star}$’, and ‘B…’ explore the dynamics and morphology of
the steady-state wind, when the initialized gas background at $t=0$ is already
blown out of the simulation domain (see Appendices B and C). Simulations ‘B…’
explore steady-state wind properties for values of the mass-loading parameter
$\beta$. Simulations ‘C0’, ‘C1’, and ‘C2’ explore the early interaction
between the flow and the initialized background gas, with varied background
densities of $n_{\mathrm{bkgd.}}=10^{-3},\ 10^{-5},\ \mathrm{and}\
10^{-7}\,\mathrm{cm^{-3}}$, respectively (see Appendix D).

##  3 Results

###  3.1 Uniform Ring

We first discuss the results for the uniform ring model from Simulation ‘A’
(see Table 1). In Fig. 2, we plot three different 2d slices of the gas number
density, temperature, velocity, and entropy at 20 Myr of Simulation ‘A’. Each
row displays a different 2d slice through the volume. The bottom row is a
slice through the origin in the $\hat{x}\hat{y}$-plane, which is a slice
through the ring plane. The middle row is a slice through the
$\hat{x}\hat{z}$-plane, showcasing the bi-conical hot outflow. The top row is
another 2d slice through the $\hat{x}\hat{z}$-plane, but at a height of
$z=1.875\,$kpc. As seen middle row panels of Fig. 5, the hot unimpeded flow
emerging away from the ring undergoes bulk cooling, while the inner wind
material shocks, producing a hot bi-conical flow with a flared-cylinder
geometry. These effects are apparent in 1d skewers of the unimpeded and
shocked flows, shown below. After an initial transient (discussed more and
shown in Fig. 6), the system quickly reaches a steady-state within the first
few Myr. Because the resolution is comparable to the ring width, there are
grid effects.

Figure 2: 2d slices of density, temperature, velocity, and entropy at 20 Myr
from Simulation ‘B0’. The top row shows the $\hat{x}\hat{y}$ plane slice at
$z=1.875\,$kpc, the middle row shows the $\hat{x}\hat{z}$ plane slice at
$y=0\,$kpc, and the bottom row shows the $\hat{x}\hat{y}$ plane at $z=0\,$kpc.
The starburst ring (bottom row) produces a multi-phase outflow (middle row)
where the hot phase is enveloped with a cooler phase (top row). The starburst
ring model quickly reaches a steady-state within the first few Myr after the
initial transient passes. The simulation size is $10\times 10\times
15\,\mathrm{kpc}$ ($\hat{x}\times\hat{y}\times\hat{z}$).

Fig. 3 shows number density, temperature, velocity, entropy, Mach number, and
pressure profiles for 1d $\hat{z}$ (red) and $\hat{x}$ (blue) skewers. The
$\hat{z}$ and $\hat{x}$ skewers pass through the origin at $x=y=0$ and
$y=z=0$, respectively. The $\hat{z}$ skewer goes along the minor axis, where
the flow is collimated. For this $\hat{x}$-skewer (blue line) the peaks, most
notably in the upper left and middle panels, correspond to the region where
energy and mass is injected into the ring. Looking at the blue line in the
upper middle panel and starting at one of the peaks, we see that the flow
begins to under-go cooling as it moves towards the origin and is then abruptly
shock heated back to $10^{7}\,$K. This compresses the flow, increasing the
density and pressure, while thermalizing its kinetic energy, pushing the flow
back into the subsonic $\cal{M}<1$ regime (see middle lower panel,
$\hat{x}$-skewer, blue line). The effects of radiative cooling are most
apparent in the entropy profile (bottom left panel). The hot collimated flow
along the minor axis ($\hat{z}$-skewer, red line) maintains a constant
entropy. Conversely, the unobstructed flow ($\hat{x}$-skewer, blue line)
undergoes strong cooling with a rapid decrease in entropy as a function of
distance. The distance at which the flow becomes marginally radiative is most
apparent in the entropy profile (bottom left panel of Fig. 3), where the
entropy gradually decreases at around $\pm 2.5\,$kpc, or about 1.75 kpc away
from the ring. Near the base of the outflow, the hot collimated flow along the
minor axis (red line) has shallower temperature, density, and velocity
gradients than a spherical flow with the same parameters, as expected from
Nguyen & Thompson (2021). The asymptotic velocities of both the collimated hot
flow and the unimpeded cool flow are:
$v_{\infty}=\sqrt{2}\,\dot{E}_{T}^{1/2}\dot{M}_{T}^{-1/2}\simeq
992\,(\alpha/\beta)^{1/2}\,\mathrm{km\ s^{-1}}\simeq
1280\,\mathrm{km\,s^{-1}}$. The temperature of the shocked gas reaches
equivalent temperatures of the uniform wind-driving ring (see the central
peaks, blue line in Fig. 3). We find the central temperature of the ring-hole
to be identical to that derived by Chevalier & Clegg (1985) for a spherical
wind-driving region (hereafter CC85). In CC85, $\dot{E}_{T}$ and
$\dot{M}_{T}$, fill a sphere, whereas here, it is the supersonic wind, which
has already left the wind-driving ring, feeding into the shocked volume that
supplies the energy and mass. In the center of the ring hole, the kinetic
energy is thermalized, such that (see Strickland & Heckman (2009))
$T_{0}=T_{\mathrm{ring}}=0.4\,\dot{E}_{T}/\dot{M}_{T}\,\mu
m_{p}/k_{B}=1.43\times 10^{7}(\alpha/\beta)\,\mathrm{K}=2.38\times 10^{7}\,K$,
where $\mu m_{p}$ and $k_{B}$ are the average particle mass and Boltzmann
constant, respectively.

The hot wind is non-spherical. In Fig. 4, we plot the calculated transverse
(with respect to the minor axis) areas for each height $z$ for Simulation ‘A’
(black x’s) and a functional fit to these values (green line).We model the
flow as a flared-cylinder (Breitschwerdt et al., 1991; Nguyen & Thompson,
2021) with $A(z)=A_{0}[1+(z/z_{b})^{\Delta}$]. The flared-cylinder function
maintains an approximately constant area $A_{0}$ up until height $z_{b}$ when
it transitions to a constant solid angle power-law expansion ($\Delta=2$ for
spherical flow). The best fit of the flared-cylinder has
$A_{0}=0.54\,\mathrm{kpc^{2}}$, $z_{b}=0.58\,$kpc, and $\Delta=2.14$. The term
$d\ln A/dz$ describes the areal divergence rate of the flow (e.g., $A=\Omega
r^{2}$ and $d\ln A/dr=2/r$ for constant $\Omega$). As the right plot in Fig. 4
shows, the flared-cylinder has a slower areal divergence rate than a spherical
flow for $z\lesssim 1.5\,$kpc, corresponding to the flow through the cylinder
of $A(z)$. However, it eventually matches closely with the spherical expansion
model, implying the flared-cylinder fit approaches a constant solid angle
expansion rate. Slower-than-spherical areal flows lead to shallower
temperature and density gradients (Nguyen & Thompson, 2021). The combination
of the kinetic energy thermalization and cylindrical flow geometry mean that
the emerging flow does not become radiative along the minor axis. This is in
contrast to the unobstructed flow, which becomes radiative (see entropy
profiles, Figs. 2 and 3). This occurs because the cooling rate is dependent on
the gas temperature, and at the sustained temperatures of $T\geq 10^{6}$, is
at least an order of magnitude less than that of the unimpeded flow when it is
at $10^{5}\,$K.

Figure 3: 1d $\hat{z}$-skewer (red lined-markers), 1d $\hat{x}$-skewer (blue
lined-markers), and the analytic solutions to the shocked gas temperature and
asymptotic velocity in orange. The $\hat{z}$ and $\hat{x}$ skewers pass
through the origin at $x=y=0$ and $y=z=0$, respectively. The $\hat{z}$ skewer
goes along the minor axis, in the direction of the hot collimated wind. Figure
4: Left: The transverse area of the hot outflow along the minor axis from
Simulation ‘A’ (black dots) and a fit using a flared-cylinder areal function
(green): $A(z)=A_{0}(1+(z/z_{b})^{\Delta})$. We find a best fit of
$z_{b}=0.58\,$kpc, $A_{0}=0.54\,\mathrm{kpc^{2}}$, and $\Delta=2.14$. Right
the areal divergence rate for the fit (green) and that of spherical expansion
(orange). The hot outflow undergoes slower-than-spherical divergence for
$z\leq 1.5\,$kpc. Figure 5: 2d slices of density, temperature, velocity, and
entropy at 40 Myr from Simulation ‘B0⋆’. The top row shows the
$\hat{x}\hat{y}$ plane slice at $z=1.875\,$kpc, the middle row shows the
$\hat{x}\hat{z}$ plane slice at $y=0\,$kpc, and the bottom row shows the
$\hat{x}\hat{y}$ plane at $z=0\,$kpc. The ring of star clusters (bottom row)
produces a multi-phase outflow (middle row) where the hot phase is enveloped
with a cooler phase (top row). The simulation size is $10\times 10\times
15\,\mathrm{kpc}$ ($\hat{x}\times\hat{y}\times\hat{z}$).

###  3.2 Non-uniform Ring

For comparison with the uniform ring in Simulation ‘A’, we contrast with
Simulation ‘$\texttt{A}_{\star}$’, where energy and mass is injected in 8
spherical star clusters distributed in a ring-like geometry (see Table 1).
Each cluster generates a wind that is identical to that predicted by the CC85
wind model, but with strong bulk radiative cooling outside of the individual
cluster volume (Wang, 1995; Silich et al., 2004; Thompson et al., 2016). To
make a direct comparison with the uniform ring model presented in Section 3.1,
we use the same total energy and mass-injection rates $\dot{E}_{T}$ and
$\dot{M}_{T}$, appropriate for a mass-loaded SN-driven flow and a constant
star formation rate. We thus neglect the time-dependent nature of cluster
evolution where massive star wind precede SNe, with different relative energy
and mass injection rates (Silich et al., 2004; Lochhaas et al., 2018). In Fig.
5, we plot three different 2d slices of the gas number density, temperature,
velocity magnitude, and entropy as in Fig. 2. As in Simulation ‘A’, after an
initial transient the system reaches a quasi steady-state. The flow
interaction between pairs of clusters is reminiscent of colliding winds from
massive star binary star systems (Girard & Willson, 1987; Kenny & Taylor,
2005).

The non-uniform ring develops multi-phase outflow with a starker temperature
contrast than the layered wind generated by the uniform ring (in contrast with
Fig. 2). This is due to the relative mass-loading within each star-cluster.
The distance at which a flow becomes radiative is critically dependent on this
volumetric mass-loading rate (Wang, 1995; Silich et al., 2004; Thompson et
al., 2016). With fixed $\dot{E}_{T}$ and $\dot{M}_{T}$, for a ring of $N$-star
clusters, the cooling radius, centered on an individual star cluster, is
$R_{\mathrm{cool}}=81.5\,\mathrm{pc}\,(\alpha^{2.13}/\beta^{2.92})(N/8)^{0.79}(\mu/0.6)^{2.13}(R_{A}/25\,\mathrm{pc})^{1.79}$.
For this simulation, we find $R_{\mathrm{cool}}\simeq 0.36\,$kpc. As shown in
Fig. 5, this appears to be the case. The inner ring flow briefly undergoes
rapid cooling to $10^{4}\,$K before it is shock heated back up to $10^{7}\,$K
at the ring’s center. Material from the cool outflow is continuously
incorporated and mixed with the hot phase at the “edges” of the cone.
Intermediate warm ($T\simeq 10^{5}$) temperature “spokes” are formed by the
cluster spacing. The hot flow near the top and bottom of the simulation begins
to approach a cool-warm temperature, as the cooling rate becomes more dominant
over the advection rate. Surprisingly, the cool phase has a higher velocity
than the hot phase. Because the resolution is comparable to the radius of an
individual cluster, there are grid effects. Our resolution tests suggest that
the qualitative features shown here are robust to changes in resolution.
Details of mixing between the hot and cooler phases may be affected by spatial
resolution.

In Fig. 6, we plot temperature and velocity slices for five different early
times, $t=0.5,\,1.0,\,1.5,\,2.5,\,\mathrm{and}\,4.5\,$Myr, for Simulation
‘$\texttt{A}_{\star}$’. The simulation domain at $t=0$ is filled with an
initial background of a warm, low density, gas ($n=10^{-5}\,\mathrm{cm^{-3}}$
and $T=10^{6}\,$K) which gets blown out of the simulation domain within
$10\,$Myr. During the initial blow out, the maximum speed of the material on
the computational grid is $\sim 4000\,\mathrm{km\,s^{-1}}$ at the leading edge
of the collimated shock (see $t=1.0\,\mathrm{Myr}$ in Fig. 6), which is
approximately three times the asymptotic velocity of the flow driven from the
individual clusters ($v\simeq 1300\,\mathrm{km\,s^{-1}}$) or two times the
non-collimated shock front ($\sim 2000\,\mathrm{km\,s^{-1}}$). These speeds
occur only during the initial transient and not in both the bubble structure
(right panel of Fig. 6) and the final steady-state structures (Figs. 2 and 5).
A high velocity transient is also seen in shock fronts at the wind collision
region of binary star systems (Pittard, 2009) as small filaments, rather than
the bubble-like ejections shown here. The high velocity bubble leaves the
simulation domain ($z=7.5\,$kpc) within $t\sim 5\,$Myr (see Section D).

Figure 6: Temperature and velocity slices for
$t=0.5,\,1.0,\,1.5,\,2.5,\,\mathrm{and}\,4.5\,$Myr, for Simulations
‘$\texttt{A}_{\star}$’. Before the wind reaches its final time-steady form, it
has a bubble-like morphology. The emergent bubble from the ring hole is fast,
with peak velocities of $4000\,\mathrm{km\,s^{-1}}$, compared to the
$2000\,\mathrm{km\,s^{-1}}$ shock front from the unobstructed, outer-ring,
flow. The asymptotic velocity of the steady-state configuration is
$v_{\infty}\simeq 1300\,\mathrm{km\,s^{-1}}$. We see that by 4.5 Myr, the
leading bubble becomes radiative and cools faster than the shock front. This
cooling fast bubble would be brighter than the surrounding shock in X-ray. The
initial background at $t=0$, consists of a galactic halo-like hot gaseous
background that gets blown out of the simulation domain within $\simeq$ 8 Myr.
Figure 7: The X-ray surface brightness projected along the $\hat{y}$ direction
which represents an ideal edge-on view of the ring system at 20 Myr for
Simulations ‘A’ and ‘A⋆’ (from Figures 2 and 5). The simulation size is
$10\times 10\times 15\,\mathrm{kpc}$ ($\hat{x}\times\hat{y}\times\hat{z}$).

In Fig. 7, we calculate the instantaneous X-ray surface brightness of the
steady-state winds of Simulations ‘A’ and ‘A⋆’ as
$\mathcal{S}_{X}^{\nu_{1},\nu_{2}}(x,z)=\int_{\nu_{1}}^{\nu_{2}}d\nu\,\int_{0}^{L_{y}}dy\,n(x,y,z)^{2}\Lambda(T(x,y,z),\nu)$,
where $L_{y}$ is the length of the simulation domain. This calculation assumes
no obscuration, collisional equilibrium ionization of the plasma, solar
metallicity, and optically thin-cooling. The energy band of the calculation is
carried out in $E_{\nu}=0.3-2\,$keV. In both models, the edge-on projected
surface brightness peaks at the wind driving rings, and is bright along the
collimated hot outflows.

##  4 Summary

In this work we study the outflows from a starburst ring. We explore
parameterized models with both uniform rings and star-clusters distributed as
a ring (Section 3). We consider a range of different ring geometries (Appendix
B), different mass-loading rates (Appendix C), and different background
densities (Appendix D). In analogy with previous calculations of spherical
outflows, unobstructed flows that are sufficiently mass-loaded (see Fig. 9 in
Appendix C), and/or have a sufficiently thin ring radius $R_{A}$ (see Fig. 8
in Appendix B), undergo rapid bulk cooling (Wang, 1995; Silich et al., 2004;
Thompson et al., 2016). Critically, the ring model proposed here creates a
layered, multi-temperature, outflow (see Figs. 2 and 5) in which the hot
($T\geq 10^{7}\,$K) phase is collimated and surrounded by the cool
($T=10^{4}\,$K) phase. This layered morphology is qualitatively similar to the
that observed observed in nearby starbursts such as M82 (Leroy et al., 2015),
NGC253 (Bolatto et al., 2013), and NGC 3079 (Hodges-Kluck et al., 2020). This
model also produces a fast ($v\geq 1000\,\mathrm{km\,s^{-1}}$) cool phase.

The Milky Way is observed to have $\gamma$-ray (Su et al., 2010) and X-ray
(Predehl et al., 2020) emitting bubbles extending several kpc above and below
the disk. While the origin of the bubbles is under debate, a popular model is
AGN feedback. Yang et al. (2022) showed that a single event of jet activity
from the central super-massive black hole can explain the bubbles. In Fig. 6,
we show that the steady-state ring-driven wind is preceded by a transient
structure that exhibits a bubble-like morphology. Here we emphasize that star-
forming rings may be able to reproduce the bubble morphology, however we do
not carry out a study with CMZ-like parameters and conditions for a detailed
comparison. The kinematics and structure of the bubbles depend on the density
of the ambient medium (see Appendix D). For background densities of
$n_{\mathrm{bkgrd}}=10^{-3},\ 10^{-5},\ \mathrm{and}\
10^{-7}\,\mathrm{cm^{-3}}$ (see Appendix D), we find that the bubbles can
reach velocities of $v>3,000-15,000\,\mathrm{km\,s^{-1}}$. These velocities
are within the regime of AGN-driven outflows (e.g., King & Pounds 2015).

The two models highlighted in Figs. 2 and 5 employ a rather large ring width
of $R_{A}=25\,$pc. The radius of a typical super-star cluster in M82 is much
smaller, at around $5\,$pc (Smith et al., 2006). A limitation of our work is
the spatial resolution ($\Delta x=13.0\,$pc). As a result, we cannot
accurately model individual realistic clusters with much smaller radii. This
limitation is apparent in Simulations ‘A’ and ‘$\texttt{A}_{\star}$’, where
grid effects are present (see Figs. 2 and 2). Silich et al. (2004) found that
radiative outflows should be generic for super star clusters with mass and
energy injection by stellar winds. We note that we do not model the time-
dependent evolution of clusters, where a wind phase precedes SN energy and
mass injection. To make the comparison between models ‘A’ and
‘$\texttt{A}_{\star}$’ most clear, we treat the properties of the energy and
mass injection as constant over the short time required to establish a steady-
state structure. A future study will include different evolutionary phases of
each star cluster (as in Schneider et al. (2020)).

Our simulations are simple parameterized models that neglect some pieces of
physics in an effort to sketch the main components of winds driven by star-
forming rings. Probably the most important for the qualitative character of
the outflow is that we neglect the gas disk that would surround any star-
forming ring. A future work will explore simulations with a disk included.
Because of the large velocities obtained in these simulations, we have also
neglected gravity. For models with lower thermalization efficiency and higher
mass-loading, the outflows will have lower velocities for which gravity may
become important. We have further used a simple approximation radiative
cooling, and have neglected non-equilibrium ionization and irradiation
effects, which is important to the ionization state of the outflow (Gray et
al., 2019; Sarkar et al., 2022).

DDN and TAT thank the OSU Galaxy/ISM Meeting for useful discussions. DDN
thanks Evan Schneider for guidance with running Cholla simulations. DDN thank
Ohio Supercomputing Center engineer, Zhi-Qiang Yuo, for technical assistance.
DDN and TAT are supported by National Science Foundation Grant #1516967 and
NASA ATP 80NSSC18K0526.

##  Appendix A The uniform starburst ring geometry and volumetric injection
rates

The physical system is specified by energy and mass-injection within a ring.
We choose to specify the starburst volume in a method similar to that in CC85.
In the spherical CC85 model, the problem can be viewed in one-dimension,
namely, the region of energy and mass-injection is satisfied by the bounds:
$r\leq R$, where $R$ is the radius of the wind-driving sphere (Chevalier &
Clegg, 1985). The ring requires additional variables. The controlling
parameters of the ring system are the distance from the origin to the center
of the ring, $R_{B}$, and the radius of the ring, $R_{A}$. For a given $R_{A}$
and $R_{B}$, we determine the minimum and maximum radii,
$R_{\mathrm{min}}(\Gamma)$ and $R_{\mathrm{max}}(\Gamma)$, which bisects outer
and inner ring surfaces, respectively. For specific line of sight angle,
$\Gamma$, the line of sight chord endpoints, which skewers the ring surface
twice, are derived to be
$R_{\mathrm{min}}=R_{\mathrm{B}}\cos(\Gamma)-(R_{\mathrm{A}}^{2}-R_{\mathrm{B}}^{2}/2+R_{\mathrm{B}}^{2}\cos(2\Gamma)/2)^{1/2}$
and
$R_{\mathrm{max}}=R_{\mathrm{B}}\cos(\Gamma)+(R_{\mathrm{A}}^{2}-R_{\mathrm{B}}^{2}/2+R_{\mathrm{B}}^{2}\cos(2\Gamma)/2)^{1/2}.$
The bounds, in spherical coordinates, corresponding to the ring volume are
derived to be $R_{\mathrm{min}}\leq r\leq R_{\mathrm{max}},\
0\leq\theta\leq\pi,\ 0\leq\varphi\leq 2\pi,$ where
$0\leq\Gamma\leq\Gamma_{\mathrm{crit}},\ \Gamma=\pi/2-\theta,\ \mathrm{and}\
\Gamma_{\mathrm{crit}}=\tan^{-1}(R_{A}/R_{B}).$ This is done for every sight
line $\Gamma(\theta)$ and then revolved over the azimuthal axis. We see that
for opening angle $\Gamma=0$, we recover the expectation that
$R_{\mathrm{min}}=R_{B}-R_{A}$ and $R_{\mathrm{max}}=R_{B}+R_{A}$. The uniform
volumetric energy and mass-injection rates are then:

$\displaystyle
q(r,\theta)=\begin{cases}\dot{M}_{\mathrm{hot}}/V,&R_{\mathrm{min}}\leq r\leq
R_{\mathrm{max}}\\\ 0,&R_{\mathrm{max}}<r\ \mathrm{or}\
r<R_{\mathrm{min}}\end{cases}$ (A1) $\displaystyle
Q(r,\theta)=\begin{cases}\dot{E}_{\mathrm{hot}}/V,&R_{\mathrm{min}}\leq r\leq
R_{\mathrm{max}}\\\ 0,&R_{\mathrm{max}}<r\ \mathrm{or}\
r<R_{\mathrm{min}}\end{cases}$ (A2)

where the normalizations for $\dot{E}_{T}$ and $\dot{M}_{T}$ are described in
Section 2, and $V$ is the volume of the starburst ring:

$V=2\pi^{2}R_{A}^{2}R_{B}.$ (A3)

##  Appendix B Varied Ring Geometries

###  B.1 Varied Uniform Ring Geometries

In Fig. 8, we show 2d temperature slices for Simulations ‘B0’, ‘B1’, ‘B2’,
‘B3’, and ‘B4’ (see Table 1). As we would expect, changing the ring hole
radius $R_{B}$ affects the width of the cylindrical hot outflow. Smaller
$R_{A}$ values lead to a smaller cooling radius, thus a more prominent
temperature contrast between the unobstructed wind and the hot thermalized
wind along the minor axis.

Figure 8: 2d temperature slices of Simulations ‘B0’, ‘B1’, ‘B2’, ‘B3’, ‘B4’.
We see that the different ring geometries produce different outflow
morphologies despite equivalent energy and mass-injection rates $\alpha$ and
$\beta$.

##  Appendix C Varied Mass-loading Rates

In Fig. 9, we show 2d density, temperature, velocity, and Mach number slices
of Simulations ‘A0’, ‘B0’, and ‘C0’. As detailed in Table 1, these simulations
have varied the SN mass-loading rates described by $\beta=1.2,\ 0.8,\
\mathrm{and}\ 0.2$, respectively. The ring geometry produces the layered wind
morphology with a temperature contrast dependent on the mass-loading. Higher
mass-loading leads to a stronger contrast. The non-radiative ring wind models
to quickly approach a steady-state. The highly mass-loaded wind model
approaches a quasi-steady state, with small thermal instability filaments
forming within the bi-cone.

Figure 9: 2d slices of density, temperature, velocity, and Mach number at 40
Myr for varied energy thermalization efficiencies within the starburst ring:
$\alpha=1.5,\ 2.0,\ \mathrm{and}\ 2.5$ (Simulations ‘C0’, ‘C1’, and ‘C2’, see
Table 1). All models have the same ring geometry of $R_{A}=0.1\,\mathrm{kpc}$
and $R_{C}=0.5\,\mathrm{kpc}$.

##  Appendix D Varied initialized ambient backgrounds

In Fig. 10 we show 2d velocity slices of Simulations ‘C0’, ‘C1’, and ‘C2’. As
detailed in Table 1, all simulations here have the same energy and mass-
injection rates and equivalent ring geometries. Here, we vary the initialized
background gas density to probe the initial interaction between the flow and
ambient background during the transient phase, before the steady-state wind is
established, to understand how fast velocities at the leading edge along the
minor axis depend on the density of the surrounding medium. Simulations ‘C0’,
‘C1’, and ‘C2’ are are initialized with constant background density of
$n_{\mathrm{bkgd.}}=10^{-3},\ 10^{-5},\ \mathrm{and}\
10^{-7}\,\mathrm{cm^{-3}}$. The temperature of the background is
$T_{\mathrm{bkgd.}}=10^{6}\,$K. For background densities of
$n_{\mathrm{bkgd.}}=10^{-3},\ 10^{-5},\ \mathrm{and}\
10^{-7}\,\mathrm{cm^{-3}}$, the respectively maximum velocity of the leading
bubble edge is approximately $5,000,\ 10,000,\ \mathrm{and}\
15,000\,\mathrm{km\,s^{-1}}$. As we would expect, expansion through a denser
medium leads to slower bubble velocities.

Figure 10: 2d velocity slices for Simulations ‘C0’, ‘C1’, and ‘C2’ (top,
middle, bottom row, respectively) for times $t=0.25,\ 0.5,\ 0.75,\ 1.0,\
1.25\,\mathrm{Myr}$. All simulations have equivalent wind-driving region
parameters of $\alpha=1$, $\beta=0.1$, $R_{A}=0.025\,$kpc, and $R_{B}=1\,$kpc.
Simulations ‘C0’, ‘C1’, and ‘C2’ are initialized with a background gas density
of $n_{\mathrm{bkgd.}}=10^{-3},\ 10^{-5},\ \mathrm{and}\
10^{-7}\,\mathrm{cm^{-3}}$, respectively.

## References

  * Armillotta et al. (2019) Armillotta, L., Krumholz, M. R., Di Teodoro, E. M., & McClure-Griffiths, N. M. 2019, MNRAS, 490, 4401, doi: 10.1093/mnras/stz2880
  * Arnaboldi et al. (1995) Arnaboldi, M., Capaccioli, M., Cappellaro, E., Held, E. V., & Koribalski, B. 1995, AJ, 110, 199, doi: 10.1086/117507
  * Athanassoula (1992) Athanassoula, E. 1992, MNRAS, 259, 328, doi: 10.1093/mnras/259.2.328
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Böker et al. (2008) Böker, T., Falcón-Barroso, J., Schinnerer, E., Knapen, J. H., & Ryder, S. 2008, AJ, 135, 479, doi: 10.1088/0004-6256/135/2/479
  * Bolatto et al. (2013) Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450, doi: 10.1038/nature12351
  * Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79
  * Buta & Crocker (1993) Buta, R., & Crocker, D. A. 1993, AJ, 105, 1344, doi: 10.1086/116514
  * Chevalier & Clegg (1985) Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44, doi: 10.1038/317044a0
  * Comerón et al. (2010) Comerón, S., Knapen, J. H., Beckman, J. E., et al. 2010, MNRAS, 402, 2462, doi: 10.1111/j.1365-2966.2009.16057.x
  * Cooper et al. (2009) Cooper, J. L., Bicknell, G. V., Sutherland , R. S., & Bland-Hawthorn, J. 2009, ApJ, 703, 330, doi: 10.1088/0004-637X/703/1/330
  * Elmegreen (1994) Elmegreen, B. G. 1994, ApJ, 425, L73, doi: 10.1086/187313
  * Girard & Willson (1987) Girard, T., & Willson, L. A. 1987, A&A, 183, 247
  * Gray et al. (2019) Gray, W. J., Oey, M. S., Silich, S., & Scannapieco, E. 2019, ApJ, 887, 161, doi: 10.3847/1538-4357/ab510d
  * Gronke & Oh (2018) Gronke, M., & Oh, S. P. 2018, MNRAS, 480, L111, doi: 10.1093/mnrasl/sly131
  * Henshaw et al. (2022) Henshaw, J. D., Barnes, A. T., Battersby, C., et al. 2022, arXiv e-prints, arXiv:2203.11223.  https://arxiv.org/abs/2203.11223
  * Hodges-Kluck et al. (2020) Hodges-Kluck, E. J., Yukita, M., Tanner, R., et al. 2020, ApJ, 903, 35, doi: 10.3847/1538-4357/abb884
  * Hughes et al. (1994) Hughes, D. H., Gear, W. K., & Robson, E. I. 1994, MNRAS, 270, 641, doi: 10.1093/mnras/270.3.641
  * Kenny & Taylor (2005) Kenny, H. T., & Taylor, A. R. 2005, ApJ, 619, 527, doi: 10.1086/426309
  * King & Pounds (2015) King, A., & Pounds, K. 2015, ARA&A, 53, 115, doi: 10.1146/annurev-astro-082214-122316
  * Leroy et al. (2015) Leroy, A. K., Walter, F., Martini, P., et al. 2015, ApJ, 814, 83, doi: 10.1088/0004-637X/814/2/83
  * Leroy et al. (2018) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2018, ApJ, 869, 126, doi: 10.3847/1538-4357/aaecd1
  * Leroy et al. (2021) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43, doi: 10.3847/1538-4365/ac17f3
  * Lochhaas et al. (2018) Lochhaas, C., Thompson, T. A., Quataert, E., & Weinberg, D. H. 2018, MNRAS, 481, 1873, doi: 10.1093/mnras/sty2421
  * Mazzuca et al. (2008) Mazzuca, L. M., Knapen, J. H., Veilleux, S., & Regan, M. W. 2008, ApJS, 174, 337, doi: 10.1086/522338
  * Moon et al. (2021) Moon, S., Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2021, ApJ, 914, 9, doi: 10.3847/1538-4357/abfa93
  * Nakai et al. (1987) Nakai, N., Hayashi, M., Handa, T., et al. 1987, PASJ, 39, 685
  * Nguyen & Thompson (2021) Nguyen, D. D., & Thompson, T. A. 2021, MNRAS, 508, 5310, doi: 10.1093/mnras/stab2910
  * Pittard (2009) Pittard, J. M. 2009, MNRAS, 396, 1743, doi: 10.1111/j.1365-2966.2009.14857.x
  * Predehl et al. (2020) Predehl, P., Sunyaev, R. A., Becker, W., et al. 2020, Nature, 588, 227, doi: 10.1038/s41586-020-2979-0
  * Sarkar et al. (2022) Sarkar, K. C., Sternberg, A., & Gnat, O. 2022, arXiv e-prints, arXiv:2203.15814.  https://arxiv.org/abs/2203.15814
  * Scannapieco & Brüggen (2015) Scannapieco, E., & Brüggen, M. 2015, ApJ, 805, 158, doi: 10.1088/0004-637X/805/2/158
  * Schneider et al. (2020) Schneider, E. E., Ostriker, E. C., Robertson, B. E., & Thompson, T. A. 2020, arXiv e-prints, arXiv:2002.10468.  https://arxiv.org/abs/2002.10468
  * Schneider & Robertson (2015) Schneider, E. E., & Robertson, B. E. 2015, ApJS, 217, 24, doi: 10.1088/0067-0049/217/2/24
  * Schneider & Robertson (2018) —. 2018, ApJ, 860, 135, doi: 10.3847/1538-4357/aac329
  * Schneider et al. (2018) Schneider, E. E., Robertson, B. E., & Thompson, T. A. 2018, ApJ, 862, 56, doi: 10.3847/1538-4357/aacce1
  * Shapiro et al. (2009) Shapiro, K. L., Genzel, R., Quataert, E., et al. 2009, ApJ, 701, 955, doi: 10.1088/0004-637X/701/2/955
  * Shopbell & Bland-Hawthorn (1998) Shopbell, P. L., & Bland-Hawthorn, J. 1998, ApJ, 493, 129, doi: 10.1086/305108
  * Silich et al. (2004) Silich, S., Tenorio-Tagle, G., & Rodríguez-González, A. 2004, ApJ, 610, 226, doi: 10.1086/421702
  * Smith et al. (2006) Smith, L. J., Westmoquette, M. S., Gallagher, J. S., et al. 2006, MNRAS, 370, 513, doi: 10.1111/j.1365-2966.2006.10507.x
  * Strickland & Heckman (2007) Strickland, D. K., & Heckman, T. M. 2007, ApJ, 658, 258, doi: 10.1086/511174
  * Strickland & Heckman (2009) —. 2009, ApJ, 697, 2030, doi: 10.1088/0004-637X/697/2/2030
  * Strickland et al. (2004) Strickland, D. K., Heckman, T. M., Colbert, E. J. M., Hoopes, C. G., & Weaver, K. A. 2004, ApJS, 151, 193, doi: 10.1086/382214
  * Su et al. (2010) Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, ApJ, 724, 1044, doi: 10.1088/0004-637X/724/2/1044
  * Tanner et al. (2016) Tanner, R., Cecil, G., & Heitsch, F. 2016, ApJ, 821, 7, doi: 10.3847/0004-637X/821/1/7
  * Thompson et al. (2016) Thompson, T. A., Quataert, E., Zhang, D., & Weinberg, D. H. 2016, MNRAS, 455, 1830, doi: 10.1093/mnras/stv2428
  * Tress et al. (2020) Tress, R. G., Sormani, M. C., Glover, S. C. O., et al. 2020, MNRAS, 499, 4455, doi: 10.1093/mnras/staa3120
  * Veilleux et al. (2005) Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769, doi: 10.1146/annurev.astro.43.072103.150610
  * Wang (1995) Wang, B. 1995, ApJ, 444, 590, doi: 10.1086/175633
  * Weiß et al. (2001) Weiß, A., Neininger, N., Hüttemeister, S., & Klein, U. 2001, A&A, 365, 571, doi: 10.1051/0004-6361:20000145
  * Westmoquette et al. (2009) Westmoquette, M. S., Gallagher, J. S., Smith, L. J., et al. 2009, ApJ, 706, 1571, doi: 10.1088/0004-637X/706/2/1571
  * Yang et al. (2022) Yang, H. Y. K., Ruszkowski, M., & Zweibel, E. G. 2022, Nature Astronomy, doi: 10.1038/s41550-022-01618-x
  * Zhang et al. (2018) Zhang, D., Davis, S. W., Jiang, Y.-F., & Stone, J. M. 2018, ApJ, 854, 110, doi: 10.3847/1538-4357/aaa8e4


# Constraining the Rotation Profile in a Low-Luminosity Subgiant with a
Surface Rotation Measurement

Tanner A. Wilson,1,2 Andrew R. Casey,1,2 Ilya Mandel,1,3 Warrick H. Ball,4
Earl P. Bellinger,5,6 and Guy Davies 4,6
1School of Physics & Astronomy, Monash University, Victoria, Australia
2Center of Excellence for Astrophysics in Three Dimensions (ASTRO-3D)
3The ARC Center of Excellence for Gravitational Wave Discovery – OzGrav,
Australia
4School of Physics and Astronomy, University of Birmingham, Edgbaston, UK
5Max Planck Institute for Astrophysics, Garching, Germany
6Stellar Astrophysics Centre (SAC), Department of Physics and Astronomy,
Aarhus University, Aarhus C, Denmark Corresponding Author E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Rotationally-induced mode splitting frequencies of low-luminosity subgiants
suggest that angular momentum transport mechanisms are 1-2 orders of magnitude
more efficient in these stars than predicted by theory. Constraints on the
rotation profile of low-luminosity subgiants could be used to identify the
dominant mechanism for angular momentum transport. We develop a forward model
for the rotation profile given observed rotational splittings, assuming a
step-like rotation profile. We identify a consistent degeneracy between the
position of the profile discontinuity and the surface rotation rate. We
perform mock experiments that show the discontinuity position can be better
constrained with a prior on the surface rotation rate, which is informed by
star spot modulations. We finally apply this approach to KIC 12508433, a well-
studied low-luminosity subgiant, as an example case. With the observed surface
rotation prior, we obtain a factor of two increase in precision of the
position of strong rotation gradient. We recover the literature values of the
core and surface rotation rates and find the highest support for a
discontinuity in the radiative zone. Auxiliary measurements of surface
rotation could substantially improve inferences on the rotation profile of
low-luminosity subgiants with already available data.

###### keywords:

asteroseismology — stars: rotation

††pubyear: 2023††pagerange: Constraining the Rotation Profile in a Low-
Luminosity Subgiant with a Surface Rotation Measurement–11

##  1 Introduction

All stars rotate. The rotation rate and distribution of angular momentum
throughout a star evolves with time. The effects of rotation on the structure
and evolution of a star are substantial (e.g. Heger, 1998; Maeder & Meynet,
2000), and accurate prescriptions of rotation in stellar models are important
to reproduce observations, particularly asteroseismic measurements of
oscillation modes.

Low-mass subgiants and low-luminosity red giants pulsate in mixed modes, which
are sensitive to structure in both the core and envelope. Mixed modes result
from the near-surface convection, which drives the oscillation modes to
amplitudes that are detectable in space-based photometry missions like CoRoT
(Baglin, 2003) and Kepler (Borucki et al., 2010). Rotation lifts the
degeneracy of oscillation modes of the same angular degree and azimuthal
order. The change to the frequency of rotationally split modes is related to
the rotation profile. For this reason, measuring rotational splittings allows
us to constrain the rotation profile.

Current measurements of rotational splittings place low precision constraints
on the core and surface rotation rates ($\approx$10% and 30% respectively
(e.g. Deheuvels et al., 2014; Fellay et al., 2021)) and have little capability
to constrain the shape of the rotation profile connecting the core and the
surface. These observations, however, present some notable results. Subgiants
demonstrate solid-body rotation early in their transition off the main-
sequence (MS), like their MS counterparts (Deheuvels et al., 2020; Noll et
al., 2021). However, later in their post-MS evolution, the cores of more
evolved subgiants and low-luminosity red giant branch (RGB) stars rotate much
faster than their envelope. The core-to-surface rotation ratio can grow to
$\approx$20 for stars leaving the subgiant phase (Deheuvels et al., 2014;
Gehan et al., 2018; Eggenberger et al., 2019).

Stellar models predict differential rotation between the core and the surface
to be two to three orders of magnitude greater than observations suggest
(Eggenberger et al., 2012; Cantiello et al., 2014). The surface rotation rates
of white dwarfs agree well with the core rotation rates of RGB stars (Gough,
2015; Hermes et al., 2017), suggesting an angular momentum transport mechanism
that is much more efficient throughout the first ascent of the RGB, and
nowhere else (Eggenberger et al., 2012; Marques et al., 2013; Ceillier et al.,
2013; Fuller et al., 2015; Spada et al., 2016; Ouazzani et al., 2018).

The angular momentum transport mechanism sets the rotation profile. The core-
to-surface rotation ratio and the position and strength of the gradient of
rotation rate can characterise the rotation profile. Fellay et al. (2021)
suggest that tighter constraints can be made on angular momentum transport
mechanisms through more precise measures of the core-to-surface rotation
ratios of post-MS stars (Deheuvels et al., 2014), and of the position and
strength of a rotation rate gradient (Di Mauro et al., 2018). For example, a
rotation profile with a constant rotation rate internal to the base of the
convective zone (BCZ), and a decreased rotation rate that is inversely
dependent on radius in the convective zone, could be indicative of angular
momentum transport from deep fossil magnetic fields (Gough & Thompson, 1990;
Kissin & Thompson, 2015; Takahashi & Langer, 2021). This results from
differential rotation being damped along poloidal field lines (Garaud, 2002;
Strugarek et al., 2011). On the other hand, a steep rotational gradient near
the H-burning shell of a subgiant would indicate turbulent angular momentum
transport. This could be in the form of internal gravity waves (Pinçon et al.,
2017), leading to localised shallow gradients in the profile (Charbonnel &
Talon, 2005) or through magneto-rotational instabilities which arise from
steep angular momentum gradients (Spada et al., 2016; Balbus & Hawley, 1994;
Arlt et al., 2003; Menou & Mer, 2006). The gradient of the rotation profile of
subgiants is not well constrained through current asteroseismic data
(Deheuvels et al., 2014).

In this work, we consider the constraints to the position of a steep
rotational gradient, where we show that $\ell=1,2$ rotational splittings may
be sufficient to make valuable inferences about the rotation profiles of low-
luminosity subgiants if a precise auxiliary measure of surface rotation rate
is available. We specifically investigate the impact of employing surface
rotation periods from photometric variability owing to stellar spots (e.g.
those measured in McQuillan et al., 2014; Garcia et al., 2014; Santos et al.,
2021). In adopting these values, we utilise a data set that overlaps with the
subset used to measure asteroseismic rotational splittings. We deem it
appropriate to employ both constraints simultaneously due to the distinct
methods of measuring these quantities.

In Section 2 we describe a forward model to infer rotation profile parameters
given observed rotational splittings, assuming a step rotation profile. We
perform tests using mock data generated by three hypothetical profiles to show
the differences in constraining the rotation profile with realistic
independent measures of surface rotation rate from stellar spot brightness
modulations. Finally, we perform inference on the observed rotational
splittings of KIC 12508433 with different priors and compare the constraints
on the rotation profile. The implications are discussed in Section 4, and
summarised in Section 5.

##  2 Method and results

###  2.1 Rotational splittings

Stellar oscillations can be decomposed into oscillating spherical harmonic
modes. Individual modes frequencies ($\nu_{n,\ell,m}$) are characterised by
their radial order ($n$), angular degree ($\ell$) and azimuthal order ($m$).
Low-luminosity subgiants have much longer rotation periods (of order 101 days
(Deheuvels et al., 2014)) in the fast rotating core than the average
oscillation period (on the order of hours (Aerts et al., 2010)). We can
therefore treat rotation as perturbative to the structure. The effect of
stellar rotation on oscillation mode frequencies can be approximated as
perturbations to the non-rotating mode frequencies from $m=0$ to $m=-\ell$ and
$m=\ell$. This is a widely employed approximation in the field of
asteroseismic inversions of rotation (e.g. Deheuvels et al., 2014; Deheuvels
et al., 2015; Fellay et al., 2021). For more detail on this approach see Unno
et al. (1989); Aerts et al. (2010). To first order, the rotationally split
oscillation frequencies are

$\displaystyle\nu_{n,\ell,m}$ $\displaystyle=$
$\displaystyle\nu_{n,\ell,0}+\delta\nu_{n,\ell,m}$ (1) $\displaystyle=$
$\displaystyle\nu_{n,\ell,0}+m\ \delta\nu_{n,\ell}$ (2)

where $\nu_{n,\ell,m}$ is the frequency of the $n,\ell,m$ mode,
$\nu_{n,\ell,0}$ is the frequency in the non-rotating case, and
$\delta\nu_{n,\ell,m}$ is the change in oscillation frequency due to rotation,
known as the rotational splitting. The difference between the $m=0$ and $m\neq
0$ mode frequency is the $m$’th multiple of the $\delta\nu_{n,\ell}$
rotational splitting. In practice we are usually only able to observe $\ell=1$
and $\ell=2$ modes with $m=1$ and $m=-1$ rotational splittings in low-
luminosity subgiants from photometric time series data (Benomar et al., 2013;
Deheuvels et al., 2014).

Rotational splittings are scaled averages of the rotation profile. The scaling
is different for each oscillation mode and is quantified using a so-called
rotational kernel. The rotational kernels are inherent to the thermodynamic
structure of a star (see Aerts et al., 2010, for a derivation of these
kernels). Assuming spherical symmetry, the $n,\ell$ rotational splitting is
given by

$\delta\nu_{n,\ell}(\Omega)=\beta_{n,\ell}\int^{R}_{0}K_{n,\ell}(r)\Omega(r)\,\text{d}r$
(3)

where $K_{n,\ell}$ is the rotational kernel of the $n,\ell$ mode (determined
from a stellar model), $\Omega(r)$ is the scaled average, with respect to the
polar axis, 1D rotation profile along the radial axis, $\beta_{n,\ell}$ is a
normalisation constant, and $R$ is the outermost radius of the star. The
rotational kernel, and thus scaled averaged frequency shift, changes with each
oscillation mode. Changes to the rotation profile therefore result in distinct
variances for each rotational splitting. Some rotation profiles are more
likely to result in measurable rotational splittings than others. As a result
we are able to use forward modelling to determine the set of likely rotation
profiles given some observed rotational splittings.

###  2.2 Forward model

A forward model requires a set of rotation kernels to predict rotational
splittings given some profile. In this work, we use rotational kernels of the
low-luminosity subgiant KIC 12508433, a well-studied asteroseismic target
(e.g., Deheuvels et al., 2014). Models of KIC 12508433 indicate that it is
early in its evolution off the MS, which is supported by a relatively low
core-to-surface rotation ratio. It is the earliest star known in its post-MS
evolution with evidence of differential rotation.

We were provided with a model of KIC 12508433 from Ball & Gizon (2017)
generated using the astero module of the Modules for Experiments in Stellar
Astrophysics (MESA) evolutionary code (r7624; Paxton et al., 2010; Paxton et
al., 2013, 2015, 2019). The model was found by simultaneously matching the
non-seismic properties of KIC 12508433 (Teff, $\log{g}$, and [Fe/H] in Table
1), the global seismic quantities - the frequency where peak power of the
Gaussian asteroseismic power envelope occurs, $\nu_{\text{max}}$, and the
frequency spacing between consecutive radial order modes with the same angular
degree, $\Delta\nu$ \- and the observed mode frequencies (see Tables 1 and 3
in Deheuvels et al., 2014) to those predicted by the model. Mode frequencies
were calculated using ADIPLS (Christensen-Dalsgaard, 2008), with combined
surface effect corrections to the frequencies (Ball & Gizon, 2014, 2017).

Most of the observed and best-fit model quantities agree to within 1-$\sigma$.
The exceptions are T${}_{\text{eff}}$ and $L$ where the agreement is within
2-$\sigma$. Our model T${}_{\text{eff}}$ is closer to the reported
T${}_{\text{eff}}$ from the infrared flux method (Casagrande et al., 2010)
(5302 $\pm$ 124 K from Deheuvels et al. (2014)), but fully consistent with the
spectroscopic effective temperature. $L$ is not included in the $\chi^{2}$
fitting of the model ($M$, $R$ and $L$ are outputs of the model) and the value
from Deheuvels et al. (2014) was estimated from scaling relations. This level
of deviation varies with surface modelling assumptions (Ball & Gizon, 2017),
is consistent with previous works (e.g. Deheuvels et al., 2014; Li et al.,
2020), and while rotational inversions of subgiants can be dependent on model
uncertainties (Schunker et al., 2016) this discrepancy is unlikely to affect
the results of this work.

We used this model of KIC 12508433 to calculate the rotational kernels for all
observable rotational splittings. Radial positions with both high kernel
magnitude and inter-kernel variance (represented by the standard deviation of
kernel magnitudes, shown in red in Figure 1) are most sensitive to the
rotation profile. For KIC 12508433 these regions are near the core and at
$r/R$ > $0.8$.

Figure 1: Rotational kernels for the best-fitting model of KIC 12508433. The red curve shows three times the local standard deviation of the set of observed kernels (kernel deviation). Regions with large kernel deviation are expected to be sensitive to the rotation profile following forward modelling. In this model these regions can be seen in the H-burning shell core ($r/R<0.05$) and nearing the surface ($r/R$ > $0.8$). The inset shows the same kernels with two modes highlighted: the $\ell=1,\ n=9$ mode (blue) is sensitive to core rotation and the $\ell=2,\ n=6$ (orange) is sensitive to surface rotation. | KIC 12508433 | Model
---|---|---
$M\,(M_{\odot}$) | $1.20~{}\pm~{}0.16$ | $1.293$
$R\,(R_{\odot})$ | $2.20~{}\pm~{}0.10$ | $2.277$
$L\,(L_{\odot}$) | $3.25~{}\pm~{}0.45$ | $4.065$
$\log{g}(\ \mathrm{cm}\mathrm{s}^{-2})$ | $3.83~{}\pm~{}0.04$ | $3.834$
$T_{\mathrm{eff,spec}}\,(\mathrm{K})$ | $5248~{}\pm~{}130$ | $5434$
$[\mathrm{Fe/H}]$ (dex) | $0.25~{}\pm~{}0.23$ | $0.06$
Table 1: Measured properties of KIC 12508433 (Deheuvels et al., 2014), and
those of the best-fitting stellar model from which the rotation kernels are
generated.

Deheuvels et al. (2014) concluded that it is difficult to distinguish between
smooth and discontinuous rotation profiles of low-luminosity subgiants using
current data. Consequently, we assume a simple form of a rotation profile
$\Omega(r)$ using three parameters of the following form:

$\Omega(r)=\left\\{\begin{array}[]{ll}\Omega_{c}&\quad r/R\leq p\\\
\Omega_{s}&\quad r/R>p\end{array}\right.$ (4)

where $\Omega_{c}$ and $\Omega_{s}$ describe the core and surface rotation
rates respectively, and $p$ is the position of the step (in units of $r/R$).
We discuss the implications of the assumed step-function form of the rotation
profile further in Section 4. Initially we will assume weak uninformed,
uniform priors on the core and surface rotation rate and a uniform prior on
$p$:

$\displaystyle p$ $\displaystyle\sim$
$\displaystyle\mathcal{U}\left(0,1\right)$ (5) $\displaystyle\Omega_{s}/2\pi$
$\displaystyle\sim$ $\displaystyle\mathcal{U}\left(0,600\right)\ \mathrm{nHz}$
(6) $\displaystyle\Omega_{c}/2\pi$ $\displaystyle\sim$
$\displaystyle\mathcal{U}\left(0,1000\right)\ \mathrm{nHz}\quad$ (7)

where $\mathcal{U}\left(x,y\right)$ denotes a uniform prior between x and y.
We calculate the expected rotational splitting frequencies $\delta\nu_{n,l}$
for all observable $\\{n,l\\}$ modes given a model $\Omega(r)$ and the
rotational kernels ${K}_{n,\ell}(r)$ using Eq. 3, and assume the observed
splitting frequencies are normally distributed with a log-likelihood

$\ln\mathcal{L}(\boldsymbol{\delta\nu}|\Omega_{c},\Omega_{r},p,\mathbf{{K}(r)},\boldsymbol{\sigma_{\delta\nu}})\propto-\frac{1}{2}\sum_{n,l}\left(\frac{\delta\nu_{n,l}-\delta\nu_{\mathrm{obs}\
n,l}}{\sigma_{\delta\nu_{n,l}}}\right)^{2}$ (8)

where $\delta\nu_{\mathrm{obs}\ n,l}$ is the observed rotational splitting
frequency for mode $n$ and $l$, and its associated uncertainty is
$\sigma_{\delta\nu_{n,l}}$. The vector-valued symbols
$\boldsymbol{\delta\nu},\mathbf{{K}(r)},$ and
$\boldsymbol{\sigma_{\delta\nu}}$ indicate that the log-likelihood depends on
the corresponding values for all of the observed rotational splittings. We
constructed this model using PyMC3 (Salvatier et al., 2016) and used the ‘No
U-Turn Sampler’ (Hoffman & Gelman, 2011) to draw samples from the posterior.

##  3 Results

###  3.1 Mock data experiments with three hypothetical rotation profiles

We begin by generating mock data with our forward model to test the impact of
independent measures of surface rotation rate. We chose three rotation
profiles with extreme differences in the position of a strong rotational
gradient, which represent hypothetical angular momentum transport mechanisms
that could result in step-like rotation profiles. The three step positions
are: in the H-burning shell (purple); in the radiative zone (blue); and at the
base of the convective zone (red) as shown in Figure 2. The mock rotation
profiles are motivated by, but not representative of, various angular momentum
transport processes. The BCZ step rotation profile (red) is a signature of
angular momentum transport by fossil magnetic fields, which results in solid
body rotation in the radiative region and inverse rotation rate on radius in
the convective region (Kissin & Thompson, 2015; Takahashi & Langer, 2021). The
H-burning step rotation profile (purple) is indicative of turbulent angular
momentum transport through internal gravity waves (Pinçon et al., 2017) or
magneto-rotational instabilities (Spada et al., 2016; Balbus & Hawley, 1994;
Arlt et al., 2003; Menou & Mer, 2006) which result in a strong gradient in
rotation rate close to the core. The radiative zone step profile (blue)
corresponds to delocalized angular momentum transport from the core into the
radiative zone and is not indicative of a specific angular momentum transport
process.

The core and surface rotation rates will be realistically different for each
profile/angular momentum transport process. To account for this for each
profile, we fixed the step position and evaluated the log-likelihood (given
the observed splittings of KIC 12508433) at each combination of
$(\Omega_{c},\Omega_{s})$ and set the rotation rates to those with the maximum
log-likelihood. The resulting rotation profiles are shown in Figure 2:
radiative-zone step (blue: $p=0.2$), a BCZ step (red: $p=0.5$), and the
H-burning shell step (purple: $p=0.05$). The rotational splitting frequencies
of these profiles were then calculated using (Eq. 3). We adopted uncertainties
on those expected values given the precision of mode frequencies measured in
KIC 12508433 (Table 3 of Deheuvels et al. (2014)). This provides us with three
mock data sets to consider the rotation profile’s effect on the observaations.

Figure 2: Three rotation profiles used in the mock data experiments. These
profiles represent the extreme consequences of different angular momentum
transport mechanisms in low-luminosity post-MS stars. Cross-hatching
represents the convective surface region and diagonal-hatching represents the
H-burning shell. See Section 3.1 for descriptions of each profile. Figure 3:
Normalised posterior density following sampling for each mock rotational
splitting experiment. From left to right the normalised posterior density of
each of the parameters of the rotation profile model are shown: surface
rotation rate ($\Omega_{s}/2\pi$), core rotation rate ($\Omega_{c}/2\pi$), and
position of rotational gradient ($p$). Thin coloured dashed lines are samples
with no prior on the surface rotation rate, thin dotted black lines correspond
to the the introduced prior on surface rotation rate and thick coloured lines
correspond to samples when the informed surface rotation prior was introduced.
Vertical black lines correspond to the input values for each of the rotation
profile parameters used to generate the mock rotational splittings. From top
to bottom each row displays the result of sampling a different set of
rotational splittings with the same colours as the rotation profiles used to
generate the mock data in Figure 2: the H-burning shell step (purple),
radiative zone step (blue) and base of convective zone step (red).

We treated the mock data generated by each step profile as if it were real
data and performed inference using the model and sampler described in Section
2.2. We performed inference twice: first with a flat prior on surface rotation
and then with a Gaussian prior on surface rotation with mean equal to the
injected surface rotation rate, and a standard deviation 10% of the mean
value. Here we have chosen 10% as representative of the average uncertainty on
state-of-the-art measurements of main-sequence and subgiant stellar rotation
from photometric variation (Santos et al., 2021). We discuss this choice
further in Section 4. We drew 20,000 posterior samples in each case. The 1D
marginalised posterior samples with a flat and informed prior are compared in
Figure 3. The full posteriors are shown in the Appendix (Figures 6 to 11).

When sampling with a uniform prior on surface rotation (Figures 6 to 8), we
identify multi-modality and a strong correlation between $\Omega_{s}$ and $p$:
smaller $\Omega_{s}$ values tend to coincide with higher $p$ values. The
inferred step position $p$ is very uncertain in all three cases, showing
nearly uniform probability throughout the domain.

The impact of an auxiliary surface rotation rate measure on the posterior
distributions is evident when comparing the normalised posterior density when
using an uninformed prior (thin dashed) to the informed prior (thick solid) in
Figure 3. When the informed prior is introduced, the degeneracy between
surface rotation rate and $p$ is broken and the surface rotation rate and $p$
are better recovered. In the H Shell step (purple) and radiative step (blue)
experiments, the prior has collapsed support for $p>0.4$. The introduction of
the prior does not, however, allow us to differentiate between rotation
profiles deeper in the star. We find that for profiles where $p\lesssim 0.2$,
multi-modality remains and the introduction of the prior on $\Omega_{s}$
increased support for $p$ closer to the core (and closer to the true value).
Indeed the introduction of the informed prior for the H shell and radiative
step experiments has constrained the posterior on $p$, but the $1\sigma$ range
on these values overlaps significantly. We could not differentiate between
these profiles using this method and state-the-art data. This is not the case
for the BCZ mock data experiment. Introducing the informed surface rotation
prior allows us to place more significant constraints on $p$. The posterior on
$p$ was flat with an uninformed prior and the BCZ step. In contrast, with a
surface rotation rate prior, the posterior is now a single peak centred at the
injected value with a percentage standard deviation of the median value of
about $18\%$. This illustrates the constraining power on $p$ of the
independent measures of $\Omega_{s}$ in specific circumstances.

###  3.2 KIC 12508433

Our mock data experiments confirm that a realistic measure of surface rotation
rate, treated as a prior during inference can better constrain the internal
rotation profile. We now apply this method to observed rotational splittings
of KIC 12508433 (Deheuvels et al., 2014). First with a uniform prior on
surface rotation between $0$ and $600$ nHz, and then with an independent
surface rotation rate (measured by Garcia et al., 2014, from photometry)
treated as a Gaussian prior on the surface rotation rate with mean 172 nHz and
standard deviation 21 nHz. The posteriors are shown in Figures 4 and 5. The
best-fitting values and credible intervals of the rotation profile parameters
for each prior are shown in Table 2 and are compared to the results of
optimally localised average (OLA) inversions performed in Deheuvels et al.
(2014).

Figure 4: Posterior distributions of the core rotation rate $\Omega_{c}$, the
surface rotation rate $\Omega_{s}$, and the discontinuity position $p$ (Eq. 4)
given the observed $l=\\{1,2\\}$ rotational splittings of KIC 12508433 and
assuming a rotation profile with a step function. Figure 5: Posteriors on step
profile parameters given the rotational splittings of KIC 12508433 and
informed prior on surface rotation rate from (Garcia et al., 2014).

|  $\Omega_{c}/2\pi$ (nHz) |  $\Omega_{s}/2\pi$ (nHz) |  $p$ ($r/R$)
---|---|---|---
Deheuvels et al. (2014) (OLA) | $532\pm 79$ | $213\pm 26$ | …
Sampling - Flat $\Omega_{s}$ Prior |  $524$${}^{+36}_{-34}$ |  $141$${}^{+51}_{-79}$ |  $0.45$${}^{+0.22}_{-0.26}$
Sampling - Informed $\Omega_{s}$ Prior |  $525$${}^{+39}_{-35}$ |  $175$${}^{+19}_{-21}$ |  $0.29$${}^{+0.13}_{-0.12}$

Table 2: best-fitting rotation profile parameters given observed $\ell=1$ and
$2$ rotational splittings of KIC 12508433 from optimally localised averages
(OLA) (Deheuvels et al., 2014) and forward modelling with flat and informed
($\Omega_{s}/2\pi=172\pm 21$ nHz) priors. OLA does not place constraints on
the rotation profile outside the g- and p-mode cavities in the core and at the
star’s surface, respectively. $p$ is more precisely measured with the
introduction of the informed prior when forward modelling.

With a uniformed prior, the median of the posterior on the core and surface
rotation rates and corresponding 16th and 84th percentile credible intervals
are $\Omega_{c}/2\pi=524^{+36}_{-34}\,\mathrm{nHz}$ and
$\Omega_{s}/2\pi={141}^{+51}_{-79}\,\mathrm{nHz}$, with the maximum a
posteriori or MAP point at $\Omega_{c}/2\pi\approx 520$ nHz and
$\Omega_{s}/2\pi\approx 170$ nHz. With the informed surface rotation rate
prior the 16th and 84th percentile credible intervals are
$\Omega_{c}/2\pi=525^{+39}_{-35}\,\mathrm{nHz}$ and
$\Omega_{s}/2\pi={175}^{+19}_{-21}\,\mathrm{nHz}$. The inferences agree with
each other, as well as the inversions performed in Deheuvels et al. (2014).
Those authors use OLA to estimate the mean rotation rate in the g-mode cavity
to be $\langle\Omega_{g}\rangle/2\pi=532\pm 79\,\mathrm{nHz}$ and the mean
rotation rate in the p-mode cavity to be $\langle\Omega_{p}\rangle/2\pi=213\pm
26\,\mathrm{nHz}$. Taking their $\Omega_{g}$ and $\Omega_{p}$ as analogues for
$\Omega_{c}$ and $\Omega_{s}$, respectively, both sets of quantities agree
within about one joint standard deviation, without considering that the mode
of the $\Omega_{s}$ posteriors are closer to $\langle\Omega_{p}\rangle$ than
the median of $\Omega_{s}$.

Like in the mock data experiments, the independent surface rotation measure
removed support for discontinuities close to the BCZ (comparing Figures 4 and
5), but two peaks remain: $p\approx 0.25$ (closest to the median posterior
value and a peak in the posterior when considering a flat prior on
$\Omega_{s}$) and a much smaller peak closer to the core. Rotation profiles
with a discontinuity located $p>0.4$ for KIC 12508433 is strongly disfavoured,
and a discontinuity in the radiative region $~{}0.2\ r/R$ is favoured.
Comparing the confidence intervals when applying a flat and informed prior, we
obtain a factor of approximately two increase in measurement accuracy on $p$.

##  4 Discussion

We find that a surface rotation rate measured from star spot modulations can
substantially improve inferences on the rotation profiles of low-luminosity
subgiants. This takes advantage of a degeneracy between $\Omega_{s}$ and $p$,
found by performing inference using mock experiments with different
discontinuity positions.

Including an independent measure of surface rotation as a prior in the
analysis of KIC 12508433 yields increased support for step-like profiles with
rotational gradients in the radiative zone and decreased support for step-like
rotation profiles with $p\gtrsim 0.4$.

In our mock data tests, we proposed that a step-like rotation profile with a
rotational gradient at $0.5$ r/R (BCZ profile) is motivated by, but not
representative of, angular momentum transport by fossil magnetic fields.
However, the simplified step rotation profile does not accurately recreate the
inverse dependency of rotation rate with radius beyond the BCZ indicative of
this angular momentum transport process. Further, the signature rotation
profile becomes most apparent in the red giant branch (Takahashi & Langer,
2021; Fellay et al., 2021). Therefore, we are not proposing a lack of support
for eliminating fossil magnetic fields as a possible explanation for the
subgiant angular momentum transport problem.

The auxiliary measure of surface rotation we adopt as a surface rotation rate
prior requires some discussion. The ’data’ we use are rotational splitting
frequencies measured from a Fourier transform of a photometric time series.
The measured surface rotation rate from Garcia et al. (2014) also uses some
subset of the same photometric time series data. The two sets of measurements
are not strictly statistically independent, such that there is a risk that we
are ‘stealing information’. The two approaches do use subsets of the same data
set, but the methods are sufficiently distinct that we consider it reasonable
to use the surface rotation rate here.

Santos et al. (2021) suggest that the average percentage uncertainty on
photometric surface rotation measurements is $\sim$10% for main-sequence and
subgiant stars. The precision varies with temperature, and peaks around 6000 K
(the Kraft break). Those authors also concluded that only about 20% of
subgiants with long and precise photometric observations from Kepler have
measurable rotation periods. Subgiants may express smaller, shorter-lived
active regions, which, when combined with long rotation periods, would present
a smaller observable fraction and less precise measures of the surface
rotation rate when compared to main-sequence counterparts. On the other hand,
observing rotational splittings in subgiants requires a long observation time,
which may lend itself to a higher fraction with measurable surface rotation
rates.

It’s difficult to precisely measure the surface rotation rate without using
photometry. Rotational velocity ($v\sin i$) only requires a single high-
resolution spectrum, but requires an estimate of the inclination angle and the
star’s radius, which is quickly changing in the post-MS and is strongly model
dependent.

We assumed a step-like (discontinuous) function for the rotation profile,
consistent with other works (e.g., Fellay et al., 2021). Realistically this
may not reflect the true shape of the rotation profile of any low-luminosity
subgiant. Asteroseismic forward modelling of the stellar rotation profile is
an ill-posed problem. The observed rotational splittings alone do not uniquely
determine the shape of a star’s rotation profile. An assumption of the
rotation profile’s shape is necessary for the rotational splittings to
constrain the range of possible rotation profile parameters. A much more
flexible model than the step-like rotation profile could be employed. However,
the inferences must be interpreted in the context of the realistic evolution
of angular momentum transport. Ultimately, our work exemplifies that
introducing an informed surface rotation prior further reduces the potential
rotation profiles already constrained by asteroseismic data. Combining state-
of-the-art asteroseismic and surface rotation data is a practical method to
obtain more substantial constraints to the evolution of rotation from already
available data.

In this work, we take KIC 12508433 as an example of a low-luminosity giant
with asteroseismic observations. It is not yet clear whether the degeneracy
between surface rotation and step position generalises across all subgiants of
different masses, metallicities, and ages. Very few subgiants are known to
have mode measurements suitable for this kind of inference. Further work is
needed.

##  5 Conclusions

We investigate the impact of independent surface rotation rate priors on
inference of the rotation profile of sugiant stars. We perform forward
modelling of the rotation profile given observations of rotational splittings
assuming a step function rotation profile. Under these assumptions, we
identify a degeneracy between the surface rotation rate and the position of
the strong gradient. Mock experiments show that the introduction of a prior on
surface rotation breaks the degeneracy, allowing us to place stronger
constraints on the position of the gradient.

We introduce an independent measure of surface rotation from star spot
photometric modulations as a prior when performing inference on the rotation
profile of KIC 12508433. We find that introducing the independent measure of
surface rotation increases the measurement precision on the position of the
rotational gradient by a factor of two and that the data strongly disfavours
discontinuities outside $r/R=0.4$ in KIC 12508433. Including auxiliary surface
rotation measurements when performing asteroseismic forward modelling of the
rotation profile is a simple way of extending what we can learn from each star
with existing data.

## Acknowledgements

We thank the anonymous referee for their thorough review of this work and for
their helpful suggestions. A. R. C. is supported in part by the Australian
Research Council through a Discovery Early Career Researcher Award
(DE190100656). Parts of this research were supported by the Australian
Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions
(ASTRO 3D), through project number CE170100013. I. M. acknowledges support
from the Australian Research Council Centre of Excellence for Gravitational
Wave Discovery (OzGrav), through project number CE17010004. I. M. is a
recipient of the Australian Research Council Future Fellowship FT190100574. W.
H. B. thanks the UK Science and Technology Facilities Council (STFC) for
support under grant ST/R0023297/1. This work has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (CartographY GA. 804752). Funding for the
Stellar Astrophysics Centre is provided by The Danish National Research
Foundation (Grant agreement no.: DNRF106).

## Data Availability

The data and models underlying this article are available upon request to the
corresponding author.

## References

  * Aerts et al. (2010) Aerts C., Christensen-Dalsgaard J., Kurtz D. W., 2010, Asteroseismology, doi:10.1007/978-1-4020-5803-5. , https://ui.adsabs.harvard.edu/abs/2010aste.book.....A
  * Arlt et al. (2003) Arlt R., Hollerbach R., Rüdiger G., 2003, Astronomy & Astrophysics, 401, 1087
  * Baglin (2003) Baglin A., 2003, Advances in Space Research, 31, 345
  * Balbus & Hawley (1994) Balbus S., Hawley J., 1994, Monthly Notices of The Royal Astronomical Society - MON NOTIC ROY ASTRON SOC, 266
  * Ball & Gizon (2014) Ball W. H., Gizon L., 2014, Astronomy & Astrophysics, 568, A123
  * Ball & Gizon (2017) Ball W. H., Gizon L., 2017, Astronomy & Astrophysics, 600, A128
  * Benomar et al. (2013) Benomar O., et al., 2013, The Astrophysical Journal, 767, 158
  * Borucki et al. (2010) Borucki W. J., et al., 2010, Science, 327, 977
  * Cantiello et al. (2014) Cantiello M., Mankovich C., Bildsten L., Christensen-Dalsgaard J., Paxton B., 2014, The Astrophysical Journal, 788, 93
  * Casagrande et al. (2010) Casagrande L., Ramirez I., Melendez J., Bessell M., Asplund M., 2010, Astronomy and Astrophysics, 512, A54
  * Ceillier et al. (2013) Ceillier T., Eggenberger P., García R. A., Mathis S., 2013, Astronomy and Astrophysics, 555, A54
  * Charbonnel & Talon (2005) Charbonnel C., Talon S., 2005, Science, 309, 2189
  * Christensen-Dalsgaard (2008) Christensen-Dalsgaard J., 2008, Astrophysics and Space Science, 316, 113
  * Deheuvels et al. (2014) Deheuvels S., et al., 2014, Astronomy and Astrophysics, 564, A27
  * Deheuvels et al. (2015) Deheuvels S., Ballot J., Beck P. G., Mosser B., Østensen R., García R. A., Goupil M. J., 2015, Astronomy & Astrophysics, 580, A96
  * Deheuvels et al. (2020) Deheuvels S., Ballot J., Eggenberger P., Spada F., Noll A., Hartogh J. W. d., 2020, Astronomy & Astrophysics, 641, A117
  * Di Mauro et al. (2018) Di Mauro M. P., Ventura R., Corsaro E., Moura B. L. D., 2018, The Astrophysical Journal, 862, 9
  * Eggenberger et al. (2012) Eggenberger P., Montalbán J., Miglio A., 2012, Astronomy & Astrophysics, 544, L4
  * Eggenberger et al. (2019) Eggenberger P., et al., 2019, Astronomy & Astrophysics, 621, A66
  * Fellay et al. (2021) Fellay L., Buldgen G., Eggenberger P., Khan S., Salmon S. J. A. J., Miglio A., Montalbán J., 2021, Astronomy & Astrophysics, 654, A133
  * Fuller et al. (2015) Fuller J., Cantiello M., Stello D., Garcia R. A., Bildsten L., 2015, Science, 350, 423
  * Garaud (2002) Garaud P., 2002, Monthly Notices of the Royal Astronomical Society, 335, 707
  * Garcia et al. (2014) Garcia R. A., et al., 2014, Astronomy & Astrophysics, 572, A34
  * Gehan et al. (2018) Gehan C., Mosser B., Michel E., Samadi R., Kallinger T., 2018, Astronomy & Astrophysics, 616, A24
  * Gough (2015) Gough D. O., 2015, Space Science Reviews, 196, 15
  * Gough & Thompson (1990) Gough D. O., Thompson M. J., 1990, Monthly Notices of the Royal Astronomical Society, 242, 25
  * Heger (1998) Heger A., 1998, PhD thesis, Technical University of Munich
  * Hermes et al. (2017) Hermes J. J., Kawaler S. D., Bischoff-Kim A., Provencal J. L., Dunlap B. H., Clemens J. C., 2017, The Astrophysical Journal, 835, 277
  * Hoffman & Gelman (2011) Hoffman M. D., Gelman A., 2011, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, doi:10.48550/arXiv.1111.4246, http://arxiv.org/abs/1111.4246
  * Kissin & Thompson (2015) Kissin Y., Thompson C., 2015, The Astrophysical Journal, 808, 35
  * Li et al. (2020) Li Y., Bedding T. R., Li T., Bi S., Stello D., Zhou Y., White T. R., 2020, Monthly Notices of the Royal Astronomical Society, 495, 2363
  * Maeder & Meynet (2000) Maeder A., Meynet G., 2000, Annual Review of Astronomy and Astrophysics, 38, 143
  * Marques et al. (2013) Marques J. P., et al., 2013, Astronomy & Astrophysics, 549, A74
  * McQuillan et al. (2014) McQuillan A., Mazeh T., Aigrain S., 2014, The Astrophysical Journal Supplement Series, 211, 24
  * Menou & Mer (2006) Menou K., Mer J. L., 2006, The Astrophysical Journal, 650, 1208
  * Noll et al. (2021) Noll A., Deheuvels S., Ballot J., 2021, Astronomy & Astrophysics, 647
  * Ouazzani et al. (2018) Ouazzani R.-M., Marques J. P., Goupil M.-J., Christophe S., Antoci V., Salmon S. J. A. J., 2018, arXiv e-prints, 1801, arXiv:1801.09228
  * Paxton et al. (2010) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2010, The Astrophysical Journal Supplement Series, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, The Astrophysical Journal Supplement Series, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, The Astrophysical Journal Supplement Series, 220, 15
  * Paxton et al. (2019) Paxton B., et al., 2019, The Astrophysical Journal Supplement Series, 243, 10
  * Pinçon et al. (2017) Pinçon C., Belkacem K., Goupil M. J., Marques J. P., 2017, Astronomy & Astrophysics, 605, A31
  * Salvatier et al. (2016) Salvatier J., Wiecki T. V., Fonnesbeck C., 2016, PeerJ Computer Science, 2, e55
  * Santos et al. (2021) Santos A. R. G., Breton S. N., Mathur S., García R. A., 2021, The Astrophysical Journal Supplement Series, 255, 17
  * Schunker et al. (2016) Schunker H., Schou J., Ball W. H., Nielsen M. B., Gizon L., 2016, Astronomy & Astrophysics, 586, A79
  * Spada et al. (2016) Spada F., Gellert M., Arlt R., Deheuvels S., 2016, Astronomy & Astrophysics, 589, A23
  * Strugarek et al. (2011) Strugarek A., Brun S., Zahn J., 2011, Astronomy & Astrophysics - ASTRON ASTROPHYS, 532
  * Takahashi & Langer (2021) Takahashi K., Langer N., 2021, Astronomy & Astrophysics, 646, A19
  * Unno et al. (1989) Unno W., Osaki Y., Ando H., Saio H., Shibahashi H., 1989, Nonradial oscillations of stars. https://ui.adsabs.harvard.edu/abs/1989nos..book.....U

##  Appendix A Sampling Results

Here we provide the posteriors following sampling of each set of rotational
splittings.

Figure 6:  Posterior distributions using mock data generated with a step
function aligned with the H burning shell ($r/R=0.05$, purple profile in
Figure 2). True values are indicated in purple. There is considerable multi-
modality and degeneracy present. Figure 7: Posterior distributions using mock
data generated with a step function in the radiative region ($r/R=0.2$, blue
profile in Figure 2) and realistic uncertainties. True values in blue. Figure
8: Posterior distributions using mock data generated with a step function at
the BCZ ($r/R=0.5$; red profile in Figure 2), and realistic uncertainties.
True values in red. Figure 9: Posterior distributions using mock data
generated with a step profile at the g-mode cavity ($r/R=0.05$; purple profile
in Figure 2), with realistic uncertainties, and a 10% prior on surface
rotation $\Omega_{s}$. There is still degeneracy between $p$ and the rotation
parameters (e.g., Figure 6), but the prior has collapsed all other modes.
Figure 10: Posterior distributions using mock data generated with a step
profile in the radiative region ($r/R=0.20$; blue profile in Figure 2), with
realistic uncertainties, and a 10% prior on surface rotation $\Omega_{s}$
(compare with Figure 7). Figure 11: Posterior distributions using mock data
generated with a step profile at the BCZ ($r/R=0.50$; red profile in Figure
2), with realistic uncertainties, and a 10% prior on surface rotation
$\Omega_{s}$ (compare with Figure 8).