\section{Introduction}
\label{sec.intro}
Theoretical studies of spiral instabilities in galaxies
\citep[{\it e.g.}][]{To81, BLLT, SK91} have almost always assumed that bulge
and halo components can be treated as inert matter that contributes to
the central attraction, but does not otherwise affect the dynamics of
the disc component.  This assumption is often made in simulations also
\citep[{\it e.g.}][]{GKC12, DVH13, SC14, Baba15}.  Exceptionally, \citet{Mark76}
argued that spiral modes would be amplified by the loss of angular
momentum to the halo, which he predicted would take place mostly at
corotation of the spiral mode.  While many simulations of discs
embedded in live halos have been performed in subsequent years, there
does not appear to have been a quantitative test of Mark's prediction.

It has long been established that strong bars experience dynamical
friction from responsive halos \citep[{\it e.g.}][]{Sell80, Wein85, DS00}.
It has also been found that bar instabilites are more vigorous in
responsive halos \citep[{\it e.g.}][]{Athan02, SaNa13, BeSe16} than in rigid
ones, and that the more rapid growth begins at tiny amplitude
\citep{Sell16}.  This last result suggests that a rigid halo might
also be an inadequate approximation for studies of spiral
instabilities.

This paper therefore reports a suite of idealized simulations of discs
in halos to determine the effect on spiral modes of subsitituting a
responsive halo for a rigid one.

\section{Technique}
\label{sec.methods}
When using simulations to make quantitative measurements of the growth
rates of instabilities, it is important to set up the initial model
with considerable care in order that it be close to a settled
equilibrium.  Even quite mild adjustments to the initial model will
mask, and perhaps seed, early growth of instabilities, and the
instability will be that of the relaxed model and not that of the mass
distribution originally intended.  Equilibrium models of a disc
embedded in a halo present a particular challenge, since few
equilibrium distribution functions (DFs) are known for the halo
component when a massive disc is present.  In this paper we present
two separate equilibrium models and provoke a spiral instability in
the disc component of each.  We then compare the evolution of the
instability in both the fully-self consistent models having reponsive
halos and cases in which the halos are frozen.

\subsection{Numerical method}
The particles in our simulations move in a 3D volume that is spanned
by two separate grids; a cylindical polar mesh and a much larger
spherical grid. The gravitational field is calculated at grid points
and interpolated to the position of each particle.  The disc particles
are initially assigned to the polar grid, while the spherical grid is
used for the bulge and halo particles.  Naturally, all particles are
attracted by all others at every step.  A full description of our
numerical procedures is given in the on-line manual \citep{Se14} and
the code itself is available for download.

Tables~\ref{tab.PKpars} and \ref{tab.DBHpars} give the values of the
numerical parameters for most of the simulations presented in this
paper.  Since the gravitational field is a convolution of the mass
density with a Green function that is most efficiently computed by
Fourier transforms, it is easy to restrict the sectoral harmonics that
contribute to the field when using a polar grid.  In most simulations
we report here, non-axisymmetric forces arising from the particles are
confined to just the $m=3$ sectoral harmonic for the PK model
(\S\ref{sec.PKmodel}) and $m =2$ for the DBH model
(\S\ref{sec.DBHmodel}).  Aspherical components of forces computed via
the spherical grid are similarly confined to the spherical harmonic $l
= 3$ or 2.

\subsection{Other details}
We make measurements of the development of non-axisymmetric
disturbances from the azimuthal Fourier coefficients of the mass
distribution on the polar grid, generally combining them to report a
weighted average of the relative amplitude over a limited radial
range.  We also compute logarithmic spiral transforms of the disc
particle distribution, which is defined for particles of unequal mass
as
\begin{equation}
A(m,\tan\gamma,t) = {\sum_{j=1}^N \mu_j\exp[im(\phi_j + \tan\gamma \ln R_j)]
\over \sum_{j=1}^N \mu_j},
\label{eq.logspi}
\end{equation}
where $(R_j,\phi_j)$ are the cylindrical polar coordinates of the
$j$-th particle of mass $\mu_j$ at time $t$, and $\gamma$, the
complement to the spiral pitch angle, is the angle between the radius
vector and the tangent to an $m$-arm logarithmic spiral, with positive
values for trailing spirals.

In order to measure mode frequencies, including growth rates, we fit
modes to these data using the procedure described by \citet{SA86}.
The perturbed surface density of a mode is the real part of
\begin{equation}
\delta\Sigma(R,\phi,t) = {\cal A}_m(R)e^{i(m\phi - \omega t)},
\label{eq.mode}
\end{equation}
where the frequency $\omega = m\Omega_p + i\beta$ is complex with
$\beta$ being the growth rate.  The complex function ${\cal A}_m(R)$,
which is independent of time, describes the radial variation of
amplitude and phase of the mode.

\citet{SA86} showed that the mode frequency, especialy the growth
rate, could be estimated more precisely from quiet start simulations
in which noise was artificially suppressed by placing particles around
rings and filtering out the higher sectoral harmonics, allowing a
longer period of exponential growth.  As this trick is more difficult
to accomplish in 3D multi-component models, we here beat down shot
noise simply by employing large numbers of particles.

\begin{table}
\caption{Numerical parameters for the PK model}
\label{tab.PKpars}
\begin{tabular}{@{}ll}
Polar grid size & 85 $\times$ 128 $\times$ 125 \\
Grid scaling & $a= 10$ grid units \\
Vertical spacing & $\delta z = 0.02a$ \\
Active sectoral harmonics & $m = 0$, 3 \\
Softening length & $a/10$ \\
Spherical grid & 201 shells \\
Active spherical harmonics & $l = 0$, 3 \\
Disc mass fraction & $f_d = 1/3$ \\
Number of disc particles & $10^8$ \\
Number of halo particles & $10^8$ \\
Basic time-step & $(a^3/GM)^{1/2}/40$ \\
Time step zones & 3 \\
\end{tabular}
\end{table}

\section{Plummer-Kuzmin model}
\label{sec.PKmodel}
A composite disc and halo model that is very nearly in equilibrium can
be constructed using a Kuzmin disc embedded within a Plummer sphere.
The potential of an axisymmetric Kuzmin disc is
\begin{equation}
\Phi_{\rm K}(R,z) = -G M \left[R^2 + (a + |z|)^2\right]^{-1/2},
\label{eq.phitot}
\end{equation}
\citep[][eq.~2.68a]{BT08}, where $M$ is the mass of the disc and $a$
is a length scale.  That of a Plummer sphere, also of mass $M$ and
core radius $a$, is $\Phi_{\rm P}(r) = -GM[r^2 + a^2]^{-1/2}$, which
is identical in the $z=0$ plane and differs only slightly for $|z| >
0$.  Thus, a composite model of two superposed concentric components
having equal length scales $a$, with disc mass $f_dM$ and sphere mass
$(1-f_d)M$, will have the same potential in the mid-plane.  We
describe this as the PK model, and it is natural to adopt units such
that $G = M = a = 1$.  Though hardly a realistic galaxy model, it has
the following advantages for our study.

\citet{Kaln76} gave a family of equilibrium DFs having a parameter
$m_K$ that determines the degree of random motion in the razor-thin
Kuzmin disc model.  Also \citet{Dejo87} derived a family of DFs having
a parameter $q$ that determines the shape of the velocity ellipsoid in
the Plummer sphere.  Neither of these DFs would yield perfect
equilibrium in a composite model, however. The Plummer sphere would be
affected by the slight difference of the disc contribution to the
potential away from the mid-plane, which we find is small enough to
ignore.  More significantly, the disc needs to be thickened to avoid
local disc instabilities.  A disc lacking random motion is Jeans
unstable \citep{To64}, while a razor-thin disc having a finite radial
velocity dispersion suffers from buckling instabilities \citep{To66,
  Arak85}.  Since thickening the disc weakens the central attraction
in the mid-plane, as does gravity softening, radial balance requires a
compensating additional central force that is the difference between
that expected from the total mid-plane potential and the numerically
determined attraction from the particles in the simulation.  The
radial variation of the extra central force needed in the mid-plane is
computed at the start, and is applied as a spherically symmetric term
that is held fixed throughout the simulation.  It can be thought of as
an additional rigid halo, although the implied mass is just a few
percent of that of the disc.

The principal advantage of pursuing these oversimplified galaxy models
is that we can employ anisotropic DFs of the Plummer sphere in our
experiments.  We wish to determine the influence of a responsive halo
on the spiral instabilities within the disc, and this model enables us
to vary the shape of halo velocity ellipsoid from radial bais, through
isotropy to azimuthal bias, as was also exploited by \citet{Sell16} in
his study of bar instabilities.  Furthermore, we are able to set up
good equilibria for both the disc and halo by drawing particles from
the known DFs, using the method described by \citet{DS00}.  We use the
1D Jeans equation to set up the vertical balance of the disc
particles, which yields an acceptable equilibrium when the radial
velocity dispersion is modest.  The resulting model is very close to
overall equilibrium, and the mass profiles of both disc and halo
remain indistinguishable from the initial set up during the entire
linear growth phase of the instability.

Since the bar instability in a disc is made even more vigorous by a
live halo, spiral instabilities would be overwhelmed by the bar mode
unless the disc mass fraction were very low, and then the spirals
would be weak and multi-armed.  We therefore restrict non-axisymmetric
disturbance forces in these simulations to the $m = 3$ sectoral
harmonic in order to inhibit the bar mode, while allowing a single
spiral mode.

We provoke a spiral instability at a predictable radius by creating a
groove in the angular momentum profile of the disc \citep{SK91}.  In
this case, we multiply the DF given by \citet{Kaln76}, which is a
function of the two classical integrals $E$ and $L_z$, by the
Lorentzian factor
\begin{equation}
\left[ 1 - {D  w_j^2 \over ( L_z - L_z^* )^2 + w_j^2 } \right],
\label{eq.groove}
\end{equation}
and select equal mass disc particles from this modified DF.  In
expression (\ref{eq.groove}), $D$ is the relative depth of the
groove, $w_j$ is the width parameter, and $L_z^*$ is the angular
momentum of the groove centre.  Note that this is a groove in the
distribution of guiding centres and epicyclic blurring implies a
shallower and broader dent in the surface density profile.

For our simulations, we choose $f_d = 1/3$ and the disc DF parameter
$m_K = 30$, which yields $Q \simeq 2$ in the disc centre and
decreasing gently to $Q\simeq1.2$ in the outer disc.  The groove
parameters are: $L_z^* = 1.4$, $w_j= 0.04$ in our adopted units, and
$D=0.4$.  The halo DF parameters in the three models with live
halos are $q=0$ for an isotropic velocity distribution, $q=2$ for
maximum radial bias, and $q=-15$ for a strong azimuthal bias.  Note
that in all cases the central parts of the halo are closely isotropic,
while anisotropies increase with radius beyond the core.

\begin{figure}
\includegraphics[width=.95\hsize,angle=0]{PK-grow.ps}
\begin{center}
\includegraphics[width=.92\hsize,angle=0]{PK-spct.ps}
\includegraphics[width=.8\hsize,angle=0]{PK-mode.ps}
\end{center}
\caption{Top panel: The weighted average of the $m=3$ relative
  overdensity in the indicated radial range over which the spiral mode
  is strongest.  The four separate simulations differ only in the halo
  DF as indicated by the line colours, and the black line is for a rigid
  halo.  Middle panel: Power spectrum of the $m=3$ component from the
  isotropic halo model over the time interval $61 \leq t \leq 180$
  indicating a single dominant mode.  The solid curve marks
  $m\Omega_c$ and the dashed curves $m\Omega_c \pm \kappa$.  Bottom
  panel: The mode fitted to the data from the isotropic halo model
  over the time interval $80 \leq t \leq 180$.  The corotation radius
  is marked by the full-drawn circle and the Lindblad resonances by
  the dotted circles.}
\label{fig.PK-res}
\end{figure}

\subsection{Results}
\label{sec.results}
We have conducted four simulations of the PK model, each employing
100M particles in the disc component.  They differ in the shape of
velocity ellipsoid in the live halo component, which also has 100M
particles and in one case the halo was replaced by a rigid
unresponsive mass of the same density.  The groove provoked an $m=3$
spiral instability in all four cases.  The power spectrum in the
middle panel of Figure \ref{fig.PK-res} indicates that the isotropic
halo model supported a single dominant mode with corotation at
$R\simeq 2.6$; power spectra from the other three simulations were
very similar.  The upper panel shows the growth of the $m=3$
disturbance amplitude in the disc component in each PK model.  The
instability took slightly different lengths of time to emerge from the
noise, but over the interval $100 \leq t \leq 180$ the slope of each
line in this log-linear plot, which is the growth rate of the mode, is
remarkably similar.  This finding is already in strong contrast with
the behaviour reported for bar instabilities by \citet{Sell16}, where
growth rates in models with different halo velocity ellipsoids were
between 2 and 5 times greater than the bar mode in a disc with rigid
halo.

\begin{table}
\caption{Mode frequencies fitted to the PK models}
\label{tab.PKmodes}
\begin{tabular}{@{}lcc}
Halo  & $m\Omega_p$ & $\beta$ \\
\hline
Rigid &  $0.663\pm0.002$ &  $0.039\pm0.004$ \\
Isotropic & $0.659\pm0.001$ & $0.040\pm0.002$ \\
Azimuthal bias &  $0.653\pm0.002$ & $0.045\pm0.004$ \\
Radial bias & $0.660\pm0.001$ & $0.041\pm0.002$ \\
\end{tabular}
\end{table}

We have estimated the mode frequencies in all four models by the
method described in \citet{SA86}, fitting eq.~(\ref{eq.mode}) to data
from the simulations over selected time ranges.  We fit both the
disturbance density and logarithmic spiral transforms, selecting
differing time ranges, radial ranges, and numbers of values of
$\tan\gamma$.  We use these different fits to obtain estimates of the
uncertainties in the fitted frequencies.  The best fit values are
presented in Table~\ref{tab.PKmodes}, with the quoted uncertainties
bracketing the full range of credible fitted values in each case.  The
mode shape drawn in bottom panel of Figure \ref{fig.PK-res} is the
real part of best fit function ${\cal A}_m(R)$ (eq.~\ref{eq.mode}),
and the contours are of only the positive relative overdensity; the
principal resonances are marked by circles.

Only the model with the strong azimuthal bias has a higher growth rate
than the others, though barely so.  This finding is supported by the
blue line in the upper panel of Figure \ref{fig.PK-res}, which may
have a visibly steeper slope over the time interval $100 < t <180$
than the other three cases.  But it seems that the mode grows at
essentially the same rate in the rigid halo as in the isotropic and
radially biased halos, and scarcely any faster in the azimuthally
biased halo.

A supporting response from the halo would most likely come from those
orbits that are near circular, directly rotating, and oriented close
to the disc plane, since they would effectively augment the disc, as
\citet{Sell16} argued for bar modes.  As the number of such orbits
should be greatest when the halo velocity ellipsoid has a strong
azimuthal bias, a significant supporting response seems physically
reasonable, but even in this case the growth rate is enhanced by no
more than $\sim10$\%.  However, the trefoil spiral in these models may
not induce as strong a halo supporting response as would a
bi-symmetric spiral, as has been established for bar modes.

\begin{figure}
\includegraphics[width=.9\hsize,angle=0]{DBH-rc.ps}
\caption{The initial rotation curve of the disc-bulge-halo model
  (solid curve) measured from the particles.  The disc contribution is
  given by the dotted line, while the dot-dash lines indicate those of
  the bulge and compressed halo.}
\label{fig.DBH-rc}
\end{figure}

\section{Disc-bulge-halo model}
\label{sec.DBHmodel}
In order to study the behaviour in a somewhat more realistic model,
and to extend it to bisymmetric spirals, we adopt an exponential disc
model embedded within a halo and a central bulge.  The bulge is
required both to create the steeper inner rise in the rotation curve
typical of a large spiral galaxy, and also to prevent bar
instabilities \citep{To81} while allowing $m=2$ spirals.

\begin{itemize}
\item The exponential disc has the surface density
\begin{equation}
\Sigma(R) = {M_d \over 2\pi R_d^2} \exp(-R/R_d),
\end{equation}
which we taper to zero over the radial range $4.5 < R/R_d < 5$, and
$M_d$ is the mass of the notional infinite disc.  The vertical density
profile of the disc is Gaussian with a scale $0.1R_d$.

\item We use a dense Plummer sphere for the bulge, which has a mass of
  $0.1M_d$ and core radius $a=0.1R_d$, and adopt the isotropic DF
  given by \citet{Dejo87}.  The bulge is dense enough that it
  dominates the central attraction in the inner part of our galaxy
  model, and the analytic DF is close to equilibrium despite the
  presence of the disc and halo.

\item We employ a \citet{Hern90} model for the halo that has the
  density profile
\begin{equation}
\rho(r) = {M_h b \over 2\pi r ( b + r )^3},
\end{equation}
where $M_h$ is the total mass integrated to infinity, and $b$ is a
length scale.  We choose $M_h = 5M_d$ and $b=4R_d$.  Naturally, the
isotropic DF that Hernquist derived for this isolated mass
distribution would not be in equilibrium when the disc and bulge are
added, so we use the adiabatic compression procedure described by
\citet{SM05} that starts from the isotropic DF given by Hernquist and
uses the invariance of both the radial and azimuthal actions to
compute a revised density profile and DF as extra mass is inserted.
The revised DF has a slight radial bias.  We also apply an outer
cutoff to the selected particles that excludes any with enough energy
ever to reach $r>10b$, which causes the halo density to taper smoothly
to zero at that radius.
\end{itemize}

\begin{table}
\caption{Numerical parameters for DBH model}
\label{tab.DBHpars}
\begin{tabular}{@{}ll}
Polar grid size & 85 $\times$ 128 $\times$ 125 \\
Grid scaling & $R_d= 10$ grid units \\
Vertical spacing & $\delta z = 0.02R_d$ \\
Active sectoral harmonics & $m = 0$, 2 \\
Softening length & $R_d/10$ \\
Spherical grid & 501 shells \\
Active spherical harmonics & $l = 0$, 2 \\
Number of disc particles & $10^8$ \\
Number of halo particles & $10^8$ \\
Number of bulge particles & $10^7$ \\
Basic time-step & $(R_d^3/GM)^{1/2}/320$ \\
Time step zones & 6 \\
\end{tabular}
\end{table}

The resulting rotation curve of this model is illustrated in
Figure~\ref{fig.DBH-rc}, where we have adopted units such that $G =
M_d = R_d = 1$.  We also adopt a constant value of $Q=1.5$ at all
radii to determine the radial velocity dispersion of the disc
particles.  Although the disc is quite massive, the high value of the
epicyclic frequency, $\kappa$, near the centre in particular implies
disc random velocities are modest everywhere, and the Jeans equations
in the epicyclic approximation \citep{BT08} yield an excellent
equilibrium.

\begin{figure}
\includegraphics[width=.95\hsize,angle=0]{DBH-grow.ps}
\begin{center}
\includegraphics[width=.92\hsize,angle=0]{DBH-spct.ps}
\includegraphics[width=.8\hsize,angle=0]{DBH-mode.ps}
\end{center}
\caption{Top panel: The weighted average of the $m=2$ relative
  overdensity in the indicated radial range over which the spiral mode
  is strongest.  The six simulations of this DBH model have rigid
  (dotted lines) and live halos (solid lines) and $L_z^* = 2$ (red),
  $L_z^* = 1.5$ (green), and $L_z^* = 2.5$ (blue). The blue lines have
  been shifted horizontally to show the similarity of slopes without
  overlapping with the others.  Middle panel: Power spectrum of the
  $m=3$ component from the live halo simulation with $L_z^* = 2$ over
  the time interval $30 \leq t \leq 80$ indicating a single dominant
  mode.  Bottom panel: The mode fitted to the data from the live halo
  model over the time interval $30 \leq t \leq 90$.  The lines in the
  middle panel and the circles in the bottom panel have the same
  meanings as in Fig.~\ref{fig.PK-res}.}
\label{fig.DBH-res}
\end{figure}

Again we wish to seed a spiral instability at a predictable radius,
but require a different strategy in this case because we do not select
disc particles from a DF.  Instead we create a groove in the disc
surface density profile by multiplying the mass of every disc particle
by the same expression as above (eq.~\ref{eq.groove}).  As before,
this sharp groove in angular momentum is blurred by epicycle motions.

We run three cases placing the groove centre at $L_z^* = 1.5$, 2, and
2.5.  In each case we choose, $w_j= 0.04$ in our adopted units, and
$D=0.6$, and run both live and rigid halos and bulges.  We employ 100M
particles in each of the disc and halo and 10M in the bulge.  The
other numerical parameters are given in Table~\ref{tab.DBHpars}.

\begin{table}
\caption{Mode frequencies fitted to the DBH models}
\label{tab.DBHmodes}
\begin{tabular}{@{}lccc}
Halo & $L_z^*$ & $m\Omega_p$ & $\beta$ \\
\hline
Rigid & 1.5 &  $1.109\pm0.001$ &  $0.093\pm0.001$ \\
Live & 1.5 & $1.108\pm0.004$ & $0.095\pm0.001$ \\
Rigid & 2 &  $0.855\pm0.001$ &  $0.065\pm0.002$ \\
Live & 2 & $0.841\pm0.001$ & $0.069\pm0.003$ \\
Rigid & 2.5 &  $0.674\pm0.001$ &  $0.049\pm0.001$ \\
Live & 2.5 & $0.667\pm0.003$ & $0.051\pm0.001$ \\
\end{tabular}
\end{table}

\subsection{Results}
The middle panel of Figure~\ref{fig.DBH-res} gives the power spectrum
in the model with the $L_z^*=2$ and live halo while the upper panel
indicates the growth of the $m=2$ disturbance in all six simulations
with this galaxy model; those having a responsive halo and bulge and
marked with solid lines and those for which both spherical components
are rigid have dotted lines.  As for the PK models, differences
between the live and rigid halos for all three diferent groove radii
are tiny -- the bi-symmetric spiral mode appears to grow at roughly
the same rate in the live as in the rigid halo in each case.
Naturally, the mode pattern speed and growth rate decrease as the
groove radius is increased because the dynamical clock runs more
slowly.  The fitted frequencies in the six runs are listed in
Table~\ref{tab.DBHmodes}, from which it may be seen that the fitted
growth rates in the live halo runs are fractionally higher than in the
rigid halo cases, but the differences are all within the estimated
uncertainties.

The mode shape for the $L_z^*=2$ case is illustrated in the bottom
panel of Figure \ref{fig.DBH-res}.  The pronounced kink across
corotation reflects the mechanism of the groove mode in a cool disc
\citep{SK91}.  The mode shape is very similar for the rigid halo case
and, apart from a different spatial scale, for the other two groove
radii.

\begin{figure}
\includegraphics[width=.95\hsize,angle=0]{DBH-torque.ps}
\caption{The time evolution of the square-root of the torque between
  the disc and halo for the $L_z^*=2 case$.  Notice the quantitative
  scaling similarity between the solid red line in the top panel of
  Fig.~\ref{fig.DBH-res} and that in this figure, indicating that the
  torque varies as the square of the mode amplitude.}
\label{fig.DBH-torque}
\end{figure}

The live halo causes other differences, apart from a marginal increase
in growth rate.  For example in the $L_z^*=2$ case, the saturation
amplitude of the mode is some 20\% higher in the live halo than in the
rigid.  Also a small fraction, $\sim 0.2$\%, of the disc angular
momentum is taken up by the halo, and a tiny fraction, $\sim 0.004$\%,
by the bulge.  Though always small, the torque on the halo
(Figure~\ref{fig.DBH-torque}) increases approximately as the square of
the mode amplitude, {\it i.e.}\ at twice the growth rate, until it saturates,
which is the scaling expected in dynamical friction as a perturber is
decelerated by its wake \citep{TW84, Sell06, BT08}.

It is noteworthy that we could not find any evidence for a bar
instability in this model, in agreement with the linear theory
prediction for a model having a dense centre \citep{To81}.  However,
inner Lindblad resonance damping of disturbances is a prediction of
small amplitude perturbation theory and the strong bisymmetric spiral
extends over a broad radial range, yet it still did not lead to
non-linear trapping into a bar in the inner disc, at least for the
duration of the simulations.

\section{Discussion}
It may seem strange that bar instabilities in responsive halos grow
several times faster than in the equivalent rigid halo \citep{Athan02,
  Sell16}, while a live halo has a tiny effect on even
bisymmetric spiral modes.  Note that this difference appears in the
linear growth phase of the instabilities and therefore the fact that
bar modes may ultimately rise to larger amplitude than do spirals
cannot account for this particular discrepancy.

Another significant difference between the two cases is the second
order dependence of the torque from the spiral instability on the
halo, which contrasts with the first order torque reported by
\citet[][his fig.~2]{Sell16} for the strong halo enhancement of bar
instabilities.

The self-consistency requirement of a linear instability is that the
disturbance density results from orbit deflections caused by the
disturbance potential.  Thus, when a live halo has a large effect on
the growth rate of the bar mode, its supporting response must be an
integral part of the instability of the disc-plus-halo system.
Furthermore, the linear growth of the disc-halo torque indicates that
the instability is driven, in part at least, by angular momentum
loss from the disc.

In contrast, we report here that a live halo has little effect on
spiral modes, indicating that the halo takes no significant part in
this type of instability in the disc-halo system.  Dynamical friction
on the disc disturbance as it moves through the non-rotating halo
is inevitable, but is weak and second order.

So why do we find a negligible halo response to spiral modes, but had
previously reported an enthusiastic supporting response to bar modes?
We suggest the following reason: \citet{LB79} demonstrated that orbits
of arbitrary eccentricity that are subject to the gravitational
influence of a weak bar generally will be repelled from alignment with
the bar {\em except} in the inner parts of galaxy models having a
quasi-uniform core, {\it i.e.}\ where the rotation curve rises roughly
linearly.  The bar unstable models simulated by both \citet{Athan02}
and by \citet{Sell16} had slowly rising rotation curves, and therefore
the prograde halo orbits near the disc plane in the inner parts of
those models would be attracted into co-alignment with the growing bar
perturbation in the disc, creating a strong supporting response from
the halo that enhances the growth rate.  Where the rotation curve is
flat or declining, the response of higher angular momentum halo
particles does not reinforce the perturbation, and disturbances in the
disc are little affected by a responsive halo.  Thus an aligning halo
response for the bisymmetric spiral in the DBH model would not be
expected because the central bulge eliminates an extensive region
where the potential is quasi-harmonic.  We have reworked Lynden-Bell's
analysis for the trefoil spiral in the PK model, finding that the key
gradient $\partial \Omega_i /\partial L_z|_{J_f}$ indicates
anti-alignment in the region of corotation for the groove mode, and
therefore no supporting response from halo orbits of any eccentricity.
Here $\Omega_i = \Omega_\phi - \Omega_R/3$ and the fast action $J_f =
L_z/3 + J_R$; see \citet{LB79} for a fuller explanation.  This
argument would seem to account both for previously published results
for the bar mode and for the contrasting finding for spiral modes in
this paper.

\section{Conclusions}
\label{sec.concl}
We have constructed two equilibrium models of moderately heavy discs
embedded in live halos.  We created grooves in both discs in order to
provoke strong spiral instabilities and have compared the growth rates
of the spiral modes in simulations that employed both rigid and live
halos.  We cannot exclude that the spiral growth rate was increased by
the halo response in both models; there was a hint of a higher growth
rate for the $m=2$ spiral mode for all three grooves in the live halo
case in the DBH model and a slightly more significant enhanced growth
rate for the $m=3$ mode in the PK model when the halo velocity
ellipsoid was strongly azimuthally biased, but the growth rate
differences were small.  Since halos in real galaxies are not expected
to have azimuthally biased velocity ellipsoids, this mild boost to
spiral growth rates, if real, is unlikely to be of relevance to the
development of spirals in galaxies.

We did find evidence for additional differences between the live and
rigid halos in the DBH models.  The halo was weakly torqued by
dynamical friction from the spiral mode, but the loss of angular
momentum from the disc was only 0.2\% by the end of the simulation
and much less during the linear growth phase.  This exchange may
have allowed the final spiral amplitude to be about 20\% greater
in the live halo case.

The results reported here justify the usual assumption made in
theoretical work and in simulations that a frozen halo is an adequate
approximation when considering the dynamics of spiral formation.

\section*{Acknowledgements}
The author thanks Scott Tremaine for some very helpful comments on a
draft of this paper, Ray Carlberg for useful correspondence, and an
anonymous referee for a careful reading of the paper.  JAS
acknowledges the continuing hospitality of Steward Observatory.

\section*{Data availability}
The data from the simulations reported here can be made available
on request.  The simulation code can be downloaded from
{\tt http://www.physics.rutgers.edu/galaxy}


\section{Introduction}

Improvements in deep, wide photometric imaging in recent years has allowed a rapid increase in the number of observed, faint dwarf galaxies (see \citealt{Simon2019} for a review), establishing and characterising a population of ultra-faint dwarf (UFD) galaxies ($\rm L_{\odot} \lesssim 10^{5} \, \rm L_{\star}$). The low stellar masses of UFDs and their shallow potential wells make them particularly suited to constrain how star formation is regulated by feedback mechanisms within the interstellar medium (ISM; see \citealt{Somerville2015, Naab2017} for reviews). Achieving this, however, requires robust theoretical predictions, both resolving the ISM processes within such small galaxies while sampling their possible cosmological histories and environments. These imperatives place high demands on computing power, and has only recently been met by cosmological hydrodynamical simulations of galaxy formation (e.g. \citealt{Jeon2017, Maccio2017, Revaz2018, Munshi2019, Wheeler2019, Agertz2020EDGE, Applebaum2021}).

Following this effort, discrepancies between individual predictions and the observed population of UFDs have become apparent. In particular, the faint-end of the mass--metallicity relation has proven a challenging and promising regime, due to its acute sensitivity to the modelling of star formation and feedback processes within a UFD (\citealt{Agertz2020EDGE}). In isolation, a dwarf's metallicity is set by (i) its overall metal production following stellar evolution and (ii) its ability to retain metals while regulating star formation in a shallow potential well. A detailed account of such ISM processes remains complex however, due to uncertainties in associating metal yields to massive stars at low metallicities (e.g. \citealt{Buck2021, Muley2021}) and to modelling the coupling between supernova and photoionization feedback to a dwarf's ISM (\citealt{Smith2019, Agertz2020EDGE, Smith2020PhotoRT}). These uncertainties propagate into inconsistent predictions for the chemical properties of field UFDs (see discussions in \citealt{Agertz2020EDGE,Applebaum2021}), with current models typically underenriching the faintest galaxies and struggling to explain the observed `plateau' of constant $\rm [Fe/H]$ at the lowest stellar masses (\citealt{Kirby2013, Simon2019}).

Despite their sensitivity to internal ISM processes, these discrepancies might not be uniquely related to star formation and feedback within the dwarfs. The known population of UFDs is primarily observed within the virial radius of the Milky Way, thus also requiring a quantification of environmental processing (e.g. \citealt{Buck2019, Applebaum2021}). Tides from a massive host will preferentially strip the metal-poor outskirts, providing a mechanism to lower a dwarf's stellar mass at nearly constant $\rm [Fe/H]$ and potentially populating the faint-end plateau. Alternatively, outflows from a nearby host might pollute and pre-enrich a dwarf's ISM at early times, leading to higher metallicities at a given stellar mass than expected in isolation. Although such environmental effects likely play a role in shaping the properties of Milky-Way satellites, their detailed contribution and importance remains obscured by their coupling with uncertainties in modelling astrophysical processes before infall (\citealt{Applebaum2021}).

Furthermore, observations of UFDs are operating at the limits of current instrumental capabilities -- deeper, more complete observations of the metal-poor outskirts of UFDs have revised and lowered determinations of dwarf metallicities (e.g. \citealt{Chiti2021}), providing a distinct way to ease tension with simulated objects. The next generation of instruments provides us with a unique opportunity to tell apart the respective roles of these mechanisms. Forthcoming wide and deep photometric surveys (e.g. the Vera Rubin Observatory, the Euclid and Nancy Roman Space Telescopes) are likely to start discovering isolated UFDs in the field (\citealt{Mutlu-Pakdil2021, Borlaff2022}). Establishing the metallicities and kinematics of such isolated objects, either through high-resolution spectroscopy of single faint stars (with e.g. the Thirty-Meter and Giant-Magellan Telescopes) or efficient scanning thanks to multiplexed spectrographs (e.g. \citealt{McConnachie2016}), will enable a clean separation on the respective roles of environmental and internal effects.

To further pursue model exploration and prepare for these prospects, we investigate in this work the sensitivity of dwarf metallicities to physical assumptions in computing the available feedback and metal budget, to which the initial mass function (IMF), i.e. the mass spectrum of a stellar population, is a key input. The IMF is generally observed to be universal within the solar neighbourhood (\citealt{Bastian2010}), although explaining this lack of variations with local interstellar conditions remains challenging (see \citealt{Kroupa2013} for a review). Extragalactic observations have in fact repeatedly suggested the existence of systematic IMF variations in massive elliptical galaxies (e.g. \citealt{vanDokkum2010, Cappellari2012, Ferreras2013, Martin-Navarro2015}, see \citealt{Smith2020IMFETGs} for a review) and in dwarf galaxies (\citealt{Dabringhausen2009, Dabringhausen2012, Geha2013, Gennaro2018}). They provide an empirical picture in which metallicity is a key driver of IMF variations (\citealt{Geha2013, Martin-Navarro2015, Gennaro2018, Liang2021}), suggesting an IMF becoming top heavy in metal-poor, faint dwarfs.

None the less, there remains significant debate in the interpretation and significance of this correlation. Within faint dwarf galaxies and their old stellar populations, massive stars are long gone, thus only allowing a measurement of the low-to-intermediate-mass IMF slope. Inferring the high-mass slope can then be affected by finite observational depth (\citealt{El-Badry2017}, although see \citealt{Filion2020} for recent deeper efforts), and requires assuming a fixed universal, functional shape that might also vary (e.g. \citealt{Yan2020}). Furthermore, such observations are typically obtained from integrated observations across several stellar populations and dwarf galaxies. Theoretical models of the interstellar medium suggest metallicity to be a critical parameter for IMF variations due to its role in regulating the efficiency of cloud cooling and subsequent fragmentation (\citealt{Marks2012, Hopkins2013, Chabrier2014, Sharda2022BottomHeavyIMF}), but how to link, if at all, variations of the IMF over galactic scales to those on giant molecular cloud and stellar cluster scales, remains a challenge (e.g. \citealt{PflammAltenburg2009IMF, Kroupa2013, Dib2018, Jerabkova2018, Guszejnov2019}).

The acute sensitivity of UFDs to galaxy formation physics, thus, provides a chance to discriminate between models and clarify this picture -- at the low metallicities typical of UFDs ($\rm [Fe/H]\approx-2$), a metallicity-varying IMF would strongly affect the metal and feedback budget from massive stars, motivating a quantification of its impact on observables and in particular the mass--metallicity relation. In this work, we thus extend ongoing efforts assessing how IMF variations affect galactic properties (e.g. \citealt{Guszejnov2017, Barber2018, Gutcke2019}) to the previously unexplored faintest galaxies. To this end, we perform cosmological, zoom simulations of a field low-mass galaxy as part of the EDGE project (Engineering Dwarfs at Galaxy formation’s Edge; \citealt{Rey2019UFDScatter, Rey2020, Rey2022EDGEHI}; \citealt{Agertz2020EDGE, Orkney2021}). Each simulation includes hydrodynamics, a comprehensive galaxy formation model (\citealt{Agertz2013,Agertz2020EDGE}; \citealt{Agertz2015}) and sufficient resolution ($m_{\mathrm{DM}} = 960 \, \rm M_{\odot}$, 3 pc) to resolve individual supernova explosions within the dwarf's ISM. We complement this setup with a new implementation of a metallicity-dependent IMF according to \citet{Geha2013} (see also \citealt{Kroupa2002} and \citealt{Marks2012}), which we describe in Section~\ref{METHOD}. We then show how a metallicity-dependent IMF boosts both feedback and metal production in the faintest dwarfs, leading to a reduction in their final stellar masses while maintaining significant iron enrichment (Section~\ref{RESULTS}). We conclude in Section~\ref{CONCLUSIONS}.

\section{Methods} \label{METHOD}

As we aim to quantify and isolate the impact of a varying IMF, we focus on a single UFD with a final dynamical halo mass of $\approx 8 \, \times \, 10^{8} \ \rm M_{\odot}$. The host dark matter halo evolved with our baseline galaxy formation model (hereafter `Fiducial') was first presented in \citet{Agertz2020EDGE}. We briefly summarize the main features of the `Fiducial' model in Section~\ref{sec:Agertz2020}, before describing how we complement it with an implementation of IMF variations in Section~\ref{sec:imfmethods}.

\subsection{Fiducial numerical setup} \label{sec:Agertz2020}

We target an isolated dark matter halo (no neighbours more massive than it within five virial radii) embedded within a cosmic underdensity to minimize the chances of environmental interactions (see \citealt{Agertz2020EDGE} for further details and \citealt{Orkney2021} for a visual). We construct zoomed, cosmological initial conditions for this dark matter halo using the \textsc{genetIC} software (\citealt{Stopyra2021}) with dark matter particles of mass $m_{\mathrm{DM}} = 945 \, \rm M_{\odot}$ and baryonic resolution $\Omega_{\mathrm{b}} / \Omega_{\mathrm{m}} \, m_{\mathrm{DM}} = 161 \, \rm M_{\odot} $. We follow the evolution of stars, dark matter and gas using the adaptive-mesh-refinement code \textsc{ramses} (\citealt{Teyssier2002}), achieving a maximum spatial resolution of 3 pc across the dwarf's ISM. We account for the cooling of gas out of equilibrium (\citealt{Rosdahl2013}), and the formation of stars using a Schmidt-like law with a star formation efficiency of $10\%$ (\citealt{Agertz2020EDGE}). Star particles are formed through stochastic sampling on a cell-by-cell basis at every simulation fine step (\citealt{Rasera2006}) with a birth mass of $300\, \rm M_{\odot}$, sufficient to ensure a complete sampling of the canonical IMF in each stellar population (e.g. \citealt{Smith2021IMF}). We note that complete sampling for each individual population is not guaranteed once the IMF is allowed to vary (Section~\ref{sec:imfmethods}), but we show in Appendix~\ref{app:imfsampling} that the IMF averaged over the multiple stellar populations of a UFD's star formation history is well-sampled, thus ensuring the robustness of the total, integrated feedback and metal production budgets. As we wish to compare results with our `Fiducial'  model, we keep this parameter fixed in this study, and plan to explore in future work the importance of IMF sampling (see also discussion in Section~\ref{CONCLUSIONS}).

We track the injection of momentum, energy, iron, and oxygen from stellar winds, Type II and Type Ia supernovae in the dwarf's ISM (SNe II, SNe Ia, see \citealt{Agertz2013, Agertz2015, Agertz2020EDGE} for details). Our resolution allows us to directly inject thermal energy associated with supernovae and self-consistently follow the build-up of momentum by solving the hydrodynamics equations, thus, strongly reducing uncertainties associated with modelling supernova feedback (\citealt{Kim2015, Martizzi2015, Smith2019}). In addition, some of our simulations model photoionization feedback with radiative transfer, following the evolution of radiation discretized into six photon groups to form a low-resolution spectrum, and computing its advection across cell boundaries using \textsc{ramses-rt} (\citealt{Rosdahl2013,Rosdahl2015RAMSESRT}, see \citealt{Agertz2020EDGE} for further information). We record the evolution of two metals, iron and oxygen, and we initialize oxygen at $10^{-3} Z_{\odot}$ to account for unresolved enrichment from Population-III stars. We model reionization as a homogeneous time-dependent, ultraviolet (UV) background as in \citet{Faucher-Giguere2009}.

\subsection{Implementing a metallicity-dependent initial mass function} \label{sec:imfmethods}

We complement the `Fiducial'  model with an IMF which varies depending on metallicity. Each stellar particle is treated as a single stellar population (SSP) with an associated two-part Kroupa IMF (\citealt{Kroupa2001}):
 \begin{equation} \label{IMF}
    \Phi(M)
    =
    A
    \begin{cases}
        \begin{array}{lll}
        C_1 M^{-1.3}   & \mathrm{for} & 0.1 \leq M < 0.5 \ \rm M_{\odot} \, ,\\
        C_2 M^{-\alpha} & \mathrm{for} & 0.5 \leq M < 100 \ \rm M_{\odot} \, , \\
        \end{array}
    \end{cases}
\end{equation}
where $\Phi$ is the number of stars of mass $M$, $A$ is a normalization constant and $C_1, C_2$ ensure continuity between the two mass regimes. For our `Fiducial'  model with a fixed IMF, the high-mass slope is set to $\alpha = 2.35$.

In this work, we assume that the IMF varies with metallicity using an observationally motivated relation (\citealt{Geha2013}, see also \citealt{Kroupa2002}), given as
\begin{equation}\label{Variation}
    \alpha =
        \begin{cases}
        0.5 \times \mathrm{[Fe/H]}_i + 2.35 & \text{if $\mathrm{[Fe/H]}_{i} > -3$} \, , \\
        0.85 & \text{if $\mathrm{[Fe/H]}_{i} \leq -3$} \, , \\
        \end{cases}
\end{equation}
where ${\rm [Fe/H]}_i$ is the metallicity of star particle $i$. Equation~\ref{Variation} follows an empirical determination of IMF variations in UFDs (\citealt{Geha2013}, figure 5), which we extrapolate from the last observational point at ${\rm [Fe/H]}=-2.5$ to ${\rm [Fe/H]}=-3.0$\footnote{We impose a floor to the extrapolation at the lowest metallicities, to prevent extreme stellar populations in the first generation of stars initialized with vanishing ${\rm [Fe/H]}$. We verified that such events represent at most $\sim 3$ per cent of the total number of stellar populations, following rapid enrichment of the ISM to ${\rm [Fe/H]}\geq-3$ by the first few SNe II.} and renormalize to match our fiducial IMF (\citealt{Kroupa2001}) at solar metallicity. We stress that IMF variations driven through ISM metallicity remain theoretically motivated (e.g. \citealt{Sharda2022BottomHeavyIMF}) although debated, suggesting future implementations of additional variation channels through molecular cloud density (e.g. \citealt{Marks2012}) or ISM turbulence (e.g. \citealt{Hopkins2013}) are motivated.

A direct consequence of Equation~\ref{Variation} is that the number of massive stars in each stellar population increases as metallicity decreases ($\alpha$ decreases and the IMF becomes more top heavy). Quantitatively, the mass fraction of massive stars ($M > 8 \, \rm M_{\odot}$) at solar metallicity is 19 per cent, compared to 94 per cent at ${\rm [Fe/H]}=-3$. Modifying the number of massive stars affects the energy, metal and radiation budget available in each stellar population, which we quantify now.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/fig1.pdf}

    \caption{Modifications to the feedback and metal budget for a stellar population with an IMF varying with metallicity. As the high-mass slope decreases (top axis) with decreasing metallicities (bottom axis), the IMF becomes more top-heavy, boosting the number of available SNe II (red) and bolometric luminosity (black dashed). In parallel, oxygen (yellow) and iron (brown) production is enhanced due to the increased number of massive stars. We cap the metallicity evolution for the most metal-poor populations ($\rm [Fe/H] \leq -3$, grey shading).}
    \label{fig1}
\end{figure}

Figure~\ref{fig1} shows the fractional change in the number of SNe II in each stellar population, $N_{\mathrm{SNeII}}$ (red), compared to a fiducial, fixed IMF (grey). In the range most relevant to UFD galaxies ($ -3 \leq{\rm [Fe/H]}\leq-1$, which encompasses the complete metallicity distribution within our simulated galaxies), the number of SNe II is  boosted by a roughly constant factor of $2.5$ compared to the fiducial case. This directly translates into an increased energy budget for feedback through $E_{\rm SNII} = N_{\mathrm{SNeII}} \, \Tilde{E}_{\mathrm{SNeII}}$ as, following our fiducial model, we assume all individual supernovae events to release the same energy $\Tilde{E}_{\mathrm{SNeII}} = 10^{51} \rm \, erg$. We leave an exploration of progenitor-dependent explosion budgets (e.g. \citealt{Sukhbold2016}) as future work.

We note a peak, followed by a slight lessening, in the increase of $N_{\mathrm{SNeII}}$ at the lowest metallicities. This feature is due to our assumption that SNe II are sampled from progenitors with masses ranging from 8 to $40 \, \rm M_{\odot}$ within their parent stellar particle -- a specific IMF slope hence maximizes the number of massive stars within this mass range (here $\alpha\approx1.45$ at $\rm [Fe/H] \approx -1.8$). A more top-heavy IMF then overproduces supermassive stars that are not accounted as SNe II. Given the remaining uncertainties in the mass upper limit for progenitors of SNe II (see e.g. \citealt{Janka2012, Janka2016} and references therein), we leave a detailed quantification of its importance in determining the feedback budget of a dwarf to a future study.

We further update the metal budget by computing progenitor-mass dependent yields for each SNe II (\citealt{Kim2014}, equations 5 and 6) and show the total mass of iron (brown) and oxygen (yellow) released by SNe II in Figure~\ref{fig1} as the IMF is varying. Following the increase in $N_{\mathrm{SNeII}}$, the overall metal production is boosted compared to the fiducial, baseline scenario. Progenitor-mass-dependent yields, however, introduce different weightings for each individual metal as SNe II originating from massive progenitors contribute exponentially more than low-mass events. Oxygen production is thus boosted by up to a factor $\approx 7.5$, while iron production can be increased up to a factor $\approx 4$.

Finally, we account for modifications to the photon budget introduced by a metallicity-dependent IMF, as radiative feedback is a key process to regulating star formation in UFDs (\citealt{Agertz2020EDGE, Smith2020PhotoRT}). To capture the early, photoionizing phase during which the total luminosity of each SSP is UV-dominated, we rescale uniformly the flux in each fiducial photon group for $\sim 40$ Myr -- the time after which all massive stars tracked in the model have turned supernova. This is necessary to ensure the correct bolometric, and hence UV, luminosity for each SSP. We compute a library of IMF-weighted, bolometric luminosities for zero-age $300\, \rm M_{\odot}$ SSPs using \textsc{starburst99} (SB99; \citealt{Leitherer1999}), scanning through the IMFs slopes expected for ISM metallicities up to solar in our model.\footnote{Although \textsc{starburst99} is only calibrated down to $0.05 \, Z_{\odot}$, the metallicity dependency of stellar evolution in this regime is subdominant compared to changes to the number of massive stars as we vary the IMF.} For each IMF mapping onto a metallicity, we then compute the boost in luminosity compared to the canonical slope, which we linearly interpolate over the metallicity range to obtain the uniform rescaling of photon groups (black dashed in Figure~\ref{fig1}). As expected from an increased number of massive stars with more top-heavy IMFs, the integrated bolometric luminosity of each SSP increases with decreasing metallicity, by up to a factor of $\sim6.5$ compared to the fiducial case.

We stress that modifying the bolometric luminosity only provides a first approximation to the available photon budget with top-heavy IMFs -- in principle, each photon group should be rescaled independently to account for changes in the spectral energy distribution. Our model is nonetheless a lower limit for the UV flux, as top-heavier IMFs would shift the peak of the SSP spectrum towards shorter wavelengths, further increasing the energy of photoionization feedback. We leave an implementation of such schemes as future work.

We identify the dark matter halo of galaxies using the \textsc{hop} halo finder (\citealt{Eisenstein1998}) and extract stellar properties using the \textsc{pynbody} (\citealt{Pontzen2013}) library. We derive stellar spectra for each IMF associated with a stellar particle using SB99, and perform mock observations of the galaxy by computing the radiative transfer through astronomical dust using a method similar to the \textsc{sunrise} code (\citealt{Jonsson2006}) and a dust attenuation curve from \citet{Li2001}. We pick a random line of sight to compute surface brightness radial profiles and determine the half-light radius, $r_{\rm1/2}$, of each simulated galaxy. Since observational abundances are determined from bright stars in the inner galaxy (e.g. \citealt{Simon2007, Kirby2010}; S. Kim et al., in preperation), we select all stars within $2 \, r_{\rm1/2}$ to determine the galaxy's mean iron metallicity, $\rm [Fe/H]$ (\citealt{Agertz2020EDGE}, eq. 5), the 1D equivalent of velocity dispersion, $\sigma_{\star}$ and total absolute V-band magnitude, $M_{\rm V}$ ({\citealt{Agertz2020EDGE}, sec. 3.4).

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/fig2.pdf}

    \caption{Stellar mass growth of each UFD as we resimulate their same dark matter host halo with varying physical models. Allowing the IMF to vary with metallicity (solid red) enhances the number of SNe II (Figure~\ref{fig1}), causing a reduction in the stellar mass of the dwarf compared to the fiducial scenario (solid turquoise). This effect combines with the increased effectiveness of photoionization feedback under a varying IMF (dashed red), further reducing the final stellar mass of our ultra-faint by a dex compared to the fiducial scenario accounting for radiative feedback (dashed turquoise).}
    \label{fig2}
\end{figure}

\section{Results} \label{RESULTS}
Having accounted for variations in energy, metal, and photon budget with metallicity, we are now equipped to assess how IMF variations affect the observable properties of an UFD galaxy.

\subsection{Stellar masses and sizes with a varying IMF} \label{RESULTS:masses}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/fig3.pdf}
    \caption{Absolute, total V-band magnitudes and half-light radii of our simulated galaxies compared to a compilation of observed, Local Group dwarf galaxies (grey points). As the budget of feedback increases through a metallicity-dependent IMF (turquoise to red), the total magnitude of our simulated ultra-faint increases, mirroring its stellar mass (Figure~\ref{fig2}). All galaxies, however, have similar sizes ($r_{\rm1/2} \sim 150 \, \mathrm{pc}$), yielding both overall fainter and more diffuse dwarfs (diagonal grey lines show constant central surface brightnesses). Our most extensive model including a varying IMF and increased radiative feedback (red triangle) lacks observed counterparts -- its properties place it below the detection capabilities of current instruments but future photometric surveys (e.g. the Vera Rubin Observatory, the Euclid and Nancy Roman Space Telescopes) will enable us to verify its physical nature.}
    \label{fig3}
\end{figure}

In Figure~\ref{fig2}, we compare the stellar mass assembly of our simulated UFDs, evolved with and without a metallicity-dependent IMF. We select all stars at $z=0$ encompassed within a sphere defined by the virial radius and centered around the galaxy, extract their birth times and accumulate their final masses over time.\footnote{The final masses of each star particle are, thus, computed after their lifetime's mass-loss from supernovae and winds, which can lead to total galactic stellar masses lower than an individual birth mass early-on.} This provides an archaeological record of the stellar mass across all progenitors, as would be observed from a dwarf in the Local Group today.

We first note that star formation proceeds up to $z\sim4$ in all simulated versions of our UFD and is quenched permanently after this time. This quenching is imposed by cosmic reionization which, by $z\sim 8-6$, has heated the intergalactic medium up to $T\sim 10^{4}~\rm K$. This increased pressure support prevents gas accretion into the potential well of small dwarf galaxies (final dynamical masses $\lesssim 10^9\, \rm M_{\odot}$, e.g. \citealt{Efstathiou1992, Hoeft2006, Okamoto2008, Noh2014}). Self-shielded gas within the galaxy can allow star formation to proceed up to $z\sim4$ (\citealt{Susa2004, Onorbe2015, Agertz2020EDGE}), but the gas reservoir of the galaxy is never replenished, eventually leading to the definitive cessation of star formation.

We now compare the UFD evolved with the `Fiducial'  model detailed in \citet{Agertz2020EDGE} to our new implementation allowing for a metallicity-dependent IMF (`Fiducial + IMF(Z)'). A varying IMF leads to systematically lower stellar masses at all times, and a reduction in the final stellar mass by a factor of $\sim2.5$ compared to the `Fiducial'  case. This reduction is naturally explained by the increased number of massive stars with a metallicity-varying IMF -- all stellar particles within our ultra-faint have $\rm [Fe/H] \leq -1.1$ at $z=0$, and hence see their number of SNe II at least doubled (Figure~\ref{fig1}).

Turning to photoionisation feedback and its combined effect with a varying IMF, we first introduce the `Fiducial + RT'  version of the UFD from \citet{Agertz2020EDGE}, recovering that accounting for additional radiative feedback from massive stars lowers the final stellar mass of the galaxy by a factor of $\approx 5.4$ compared to the `Fiducial'  model. Further accounting for the increased photoionizing budget from a varying IMF (recall Figure~\ref{fig1}), here denoted `Fiducial + IMF(Z) + RT', results in a much stronger suppression of up to 2 dex compared to the `Fiducial'  model. By increasing the available supernova and radiative feedback budget, a metallicity-varying IMF strongly reduces star formation in UFDs before cosmic reionization permanently shuts down their gas supply. It therefore predicts the existence of low mass, field ultra-faints, here with $\rm M_{\star} \sim 10^{3} \, \rm M_{\odot}$ and dynamical masses $\approx 8 \, \times \, 10^{8} \, \rm M_{\odot}$ at $z=0$.

Figure~\ref{fig3} then compares the sizes and total absolute V-band magnitudes of our simulated objects to a compilation of known, satellite dwarf galaxies in the Local Group (grey points; \citealt{McConnachie2012, Kirby2014, Homma2019, Torrealba2019, Simon2019}). As expected from the evolution of the final stellar masses, the total V-band magnitude increases as we account for a varying IMF and radiative feedback. All models, however, produce an UFD with a similar half-light radius of $r_{\rm1/2} \approx 150-200 \, \mathrm{pc}$, and we further verified similar 1D stellar velocity dispersions ($\sigma_{\star} \approx 5-6 \, \rm km/s$). We therefore reaffirm that re-simulating the same dark matter history with variations to the feedback budget yield similar dynamical properties for UFDs (\citealt{Agertz2020EDGE}), as these are primarily set by the mass and growth history of the host dark matter halo (e.g. \citealt{Rey2019UFDScatter}).

A direct consequence of reducing stellar mass at fixed size is to create both overall fainter, and lower surface brightness dwarfs (constant central surface brightness lines are shown in grey in Figure~\ref{fig3}). All simulations are compatible with observed systems, except for the model featuring both IMF variations and radiative transfer (red triangle), which creates a more diffuse ultra-faint than currently observed.

This mismatch is best explained by the current observational difficulties of detecting objects at such low surface brightnesses. Given its total magnitude and size, this simulated ultra-faint would currently be undetectable within the virial radius of the Milky Way by a modern imaging survey (e.g. DES; \citealt{Drlica-Wagner2020}, figure 5), let alone in the field where our simulated objects reside. Another source of uncertainty in this comparison lies in that observed faint dwarfs all reside within the Milky Way. Their sizes and stellar contents could thus potentially be processed by their host's tidal environment during infall. However, very few, if any, UFDs appear to be tidally limited today based on their orbits (\citealt{Simon2018}), their stellar kinematics lack tidal processing signatures (\citealt{Zoutendijk2021TidesAndCores}), while the identification of member stars out to large radii in at least one case further argues against significant tidal stripping (\citealt{Chiti2021}). It, thus, remains unclear whether, and by how much, tidal processing affects the observables of currently known UFDs, suggesting our comparison is likely to be valid to leading order. Furthermore, dwarf galaxy sizes could also be underestimated at the faintest end, due to finite observational depth in the outskirts of galaxies (e.g \citealt{Chiti2021}). We plan in future work to quantify the importance of this effect when comparing simulated and observed objects (S. Kim et al., in preperation).

Such faint and diffuse systems beyond observational reach have been predicted to arise through different physical mechanisms, such as tidal stripping in a Milky-Way environment in which observed systems reside (e.g. \citealt{Applebaum2021}), from an altogether different modelling of feedback and galaxy formation physics (e.g. \citealt{Wheeler2019}), or from the diversity of possible assembly histories in field UFDs (\citealt{Rey2019UFDScatter}) -- a metallicity-dependent IMF therefore provides a distinct and new avenue to create faint UFDs beyond the reach of current surveys.

Our results thus suggest IMF-variations as a plausible framework to interpret photometric observations of UFDs, although it remains to be quantified whether systematic variations across every dwarf's ISM would be consistent with the overall observed population. Such quantification requires us to both (i) sample multiple cosmological histories at a given halo mass to quantify how IMF variations affect the scatter in stellar masses and sizes (see also \citealt{Rey2019UFDScatter}), and (ii) sample alternate halo masses to pinpoint the response of observables as the IMF tends towards its canonical shape at higher masses and metallicities. Our study provides a strong motivation for future studies to pursue this assessment, in particular to establish which levels of IMF variations are compatible with the more complete census of faint dwarfs' sizes and absolute magnitudes that will be obtained by upcoming photometric surveys such as the Vera Rubin Observatory and Nancy Roman telescope. We leave this to future work and focus now on the imprint of IMF variations on spectroscopic properties of faint dwarfs.

\subsection{Iron enrichment and the mass--metallicity relation} \label{RESULTS:metallicity}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/fig4.pdf}

    \caption{Mean iron metallicities and absolute V-band magnitudes of our simulated galaxies. Compared to the fiducial model (turquoise circle), radiative feedback (turquoise triangle) regulates star formation less explosively, lowering the final stellar mass while retaining similar metal amounts. Allowing for IMF variations (red circle) yields comparable outcomes, as the increased number of supernovae is compensated by an increased metal production (Figure~\ref{fig1}). Combining these two effects (red triangle) amplifies this trend, leading to a faint but significantly enriched field ultra-faint within the observed metallicity `plateau' from Local Group galaxies (both confirmed and candidates reported as grey points). Enhancing the strength of supernova feedback (yellow) also produces fainter systems but drastically reduces their mean metallicities, discrepant with observations due to overly ejective regulation of star formation.}
    \label{fig4}
\end{figure}

The mass--metallicity relation is a sensitive probe of alterations to the feedback, outflow and metal budgets in a dwarf galaxy (\citealt{Agertz2020EDGE}), all of which are modified when introducing a metallicity-dependent IMF (Figure~\ref{fig1}, see also e.g. \citealt{Koeppen2007}). We quantify the impact of a varying IMF on the metallicity content of our simulated dwarfs by showing the mass weighted, average, iron metallicity as a function of their absolute, total V-band magnitude in Figure~\ref{fig4}.

We start by comparing the impact of IMF variations without accounting for photoionization feedback from massive stars (circles). As discussed in Section~\ref{RESULTS:masses}, allowing for a metallicity-dependent IMF increases the number of SNe II, reducing the galaxy's final stellar mass and making it fainter overall. However, the overall production of iron by a stellar population is also significantly boosted when allowing for a varying IMF, resulting in an increased \rm [Fe/H] \ from $-2.1$ to $-1.8$. A metallicity-dependent IMF therefore provides a way to decrease an ultra-faint's stellar mass, while maintaining a similar metal content within the galaxy by increasing metal production.

In this first scenario, star formation is regulated by repeated blowouts following the explosions of supernovae in the dwarf's ISM. Photoionization from massive stars, however, provides an alternative regulation mechanism, keeping the dwarf's ISM warmer, decreasing the clustering of supernovae and in turn the strength of galactic outflows (\citealt{Agertz2020EDGE, Smith2020PhotoRT}). This allows for greater metal retention within the dwarf galaxy, thus also resulting in a lowering of the galaxy's stellar mass at roughly fixed \rm [Fe/H]\ (turquoise triangle).

These two effects combine further when accounting for both photoionization feedback and a metallicity-dependent IMF (red triangle). The overproduction of massive stars increases the effectiveness of radiative feedback by boosting the early UV luminosity of each stellar population (Figure~\ref{fig1}). This results in an ever-fainter dwarf galaxy, but with a metallicity only decreased from --2.1 to --2.5 compared with the fiducial case. Thus, an increasingly top-heavy IMF in metal-poor ISMs lowers the stellar masses of UFDs by up to 2 dex, without significantly lowering their metal content by increasing both metal production per population and retention within the ISM.

The significance of this statement is emphasized by comparing our simulated galaxies with an observed sample of dwarf galaxies around our Milky Way (grey points; \citealt{Kirby2013, Simon2019} showing both candidates and confirmed UFDs). An apparent feature in the observed mass--metallicity relation is the plateau in metallicity at $\rm [Fe/H] \approx -2.5$ for faint dwarf galaxies with $M_{\rm V} \gtrsim -7.5$. This plateau has proven to be a significant predictive challenge for current numerical simulations of UFDs (\citealt{Agertz2020EDGE}), with most studies underenriching dwarf galaxies in iron at the faintest end both when modelling faint dwarfs in the field and within a Milky-Way environment (\citealt{Maccio2017, Wheeler2019, Applebaum2021}). The robustness of the faint-end plateau remains uncertain due to the small number of confirmed observations (\citealt{Simon2019}), but its ability to discriminate between galaxy formation models makes it a particularly promising prospect in the near future.

We further show how such enrichment patterns cannot be reproduced by remaining uncertainties in supernova feedback. Although our resolution allows us to accurately follow the explosion of individual supernovae, pinpointing their clustering and subsequent coupling to the dwarf's ISM remains challenging (\citealt{Smith2019, Agertz2020EDGE, Smith2020PhotoRT}). We show re-simulations of the `Fiducial'  galaxy artificially increasing the strength of each supernovae event 10-fold, 100-fold and turning off all feedback channels (yellow symbols, see also \citealt{Agertz2020EDGE}). As the strength of supernova feedback increases, blowouts of gas become more efficient at regulating star formation, reducing the final stellar mass of the dwarf. This more ejective regulation mode, however, steeply lowers the retained metal content, resulting in underenriched dwarfs at a given stellar mass.

Our study therefore reaffirms the faint end of the mass--metallicity relation as an acute test of galaxy formation and ISM physics. We demonstrate that observationally motivated IMF variations (\citealt{Geha2013}) provide a self-consistent mechanism to enrich isolated, field UFD galaxies to the observed metallicity `plateau', strongly motivating further assessments of the role played by IMF variations in shaping the observables of the dwarf galaxy population.

In particular, quantifying whether IMF variations as considered in this work would stay consistent with the observed mass--metallicity relation at higher masses should be a key focus for future studies. As stellar masses and ISM metallicities increase, the IMF will smoothly transition towards its canonical value at solar metallicity to remain consistent with Milky-Way determinations. This potentially provides a natural explanation for the transition from a plateau to a slope in the mass--metallicity relation, but might also over-enrich dwarfs at higher masses compared to current data. Pinpointing which levels of IMF variations are allowed across the classical and faint dwarf population will thus be key to constrain our understanding of star formation. Another interesting signature to differentiate a varying from a canonical IMF could be to target the significant overproduction of $\alpha$ elements due to IMF variations (recall Figure~\ref{fig1}, gold). As an example, the average $\rm [O/Fe]$ abundance is increased from 0.75 to 0.92 when complementing the `Fiducial + RT'  model with IMF variations. Similar $\alpha$--enrichment patterns have been reported in local UFDs of similar masses (e.g. \citealt{Vargas2013, Frebel2014, Simon2019}), but abundance data is limited to a handful of bright stars within each observed dwarf and the galaxy-to-galaxy and star-to-star scatter makes quantitative comparisons currently challenging.

Finally, implementing a varying IMF at these low metallicities produces rather extreme stellar populations, some only consisting of a few massive stars. We show in Appendix~\ref{app:imfsampling} that the integrated metal and feedback budget over a dwarf's history is robust to the mass and phase--space sampling of the IMF, but the impact of sampling noise on the instantaneous coupling between feedback and the local ISM remains to be quantified with strongly top-heavy IMFs (see also \citealt{Kroupa2013, Su2018, Grudic2019, Applebaum2020, Smith2021IMF} for discussions with a canonical IMF). In future work, we plan on coupling our varying-IMF implementation with the ability to resolve energy and metal injection from individual massive stars in a galaxy's ISM (e.g. \citealt{Emerick2019, Andersson2020, Gutcke2021Model}). This will allow us to quantify the robustness of our predictions to IMF-sampling uncertainties, while allowing for a more robust comparison with spectroscopic observations focusing on limited numbers of individual bright stars.

\section{Discussion and conclusions} \label{CONCLUSIONS}
We have studied how a varying IMF with a metallicity-dependent slope affects the observable properties, such as mass and metallicity, of a faint dwarf galaxy. We performed zoom-in, cosmological simulations capable of resolving individual supernovae explosions ($m_{\mathrm{DM}} = 960 \, \rm M_{\odot}$, $3 \, \mathrm{pc}$), re-simulating a $\sim 8 \, \times \, 10^{8} \, \rm M_{\odot}$ host dark matter halo with varying physical models. We complemented the EDGE galaxy formation model (\citealt{Agertz2020EDGE}) with a new implementation of a metallicity-dependent IMF as suggested by observational reports (\citealt{Geha2013}), self-consistently modelling the resulting alterations to the energy, momentum, radiation and metal budget within the dwarf as the IMF varies.

Under this framework, the IMF becomes top heavy in stellar populations born with lower gas metallicities. In the metallicity range relevant to UFDs ($ -3 \leq{\rm [Fe/H]}\leq-1$), the subsequent increase in the fraction of massive stars per population boosts the overall number of supernovae and UV luminosity (Figure~\ref{fig1}), enhancing the effectiveness of supernova and photoionization feedback at regulating star formation.

This increase results in a hundred-fold suppression of the final stellar mass of the galaxy at fixed halo mass today compared to a canonical IMF (Figure~\ref{fig2}). Furthermore, because we re-simulate the same halo and growth history every time, the ultra-faint's final size does not significantly vary with each re-simulation. We therefore demonstrate that such IMF variations provide a new mechanism to create very low-surface-brightness, field dwarfs (Figure~\ref{fig3}), which will become observable with the next-generation of wide and deep photometric surveys such as the Vera Rubin Observatory.

In parallel to boosting the available feedback budget, a top-heavy IMF at low ISM metallicities enhances the production of metals in each stellar population (Figure~\ref{fig1}). The net result is a reduction of a dwarf's stellar mass, while maintaining nearly constant iron content (Figure~\ref{fig4}). Our implementation of IMF variations therefore provides a mechanism to produce faint, yet enriched, UFDs. This allows us to populate the faint-end ($M_{\rm V} \gtrsim - 6$) of the observed mass--metallicity relation, which has remained a significant challenge for modern, numerical simulations of UFDs (e.g. \citealt{Wheeler2019,Agertz2020EDGE, Applebaum2021}).

Our study, thus, establishes a metallicity-varying IMF as a plausible and compelling framework to interpret observations of UFDs, strongly motivating further studies. A key line of future work will be to extend our findings to multiple objects and formation histories, in order to quantify whether this change in physical model remains compatible with the overall observed population of faint and higher mass, classical dwarfs. Furthermore, our results reaffirm the faint end of the mass--metallicity relation as an acute probe of galaxy formation models. This sensitivity encourages further model testing in the ultra-faint regime, for example, extending our approach to a physically motivated ansatz for how the IMF varies (e.g. with local ISM density, pressure or turbulence; \citealt{Bonnell2006, Marks2012, Hopkins2013}), alternate IMF shapes (e.g. \citealt{Yan2020}) or different empirically motivated scalings (e.g. with stellar velocity dispersion; \citealt{Cappellari2012}). Undertaking this quantification will be key to derive constraints from current and coming photometric and spectroscopic surveys, particularly in light of the prospect to detect isolated, field UFDs in the near-future. This will offer a promising avenue to distinguish ISM-based alterations, like a varying IMF or additional feedback channels (e.g. \citealt{Agertz2020EDGE, Smith2020PhotoRT}), from environmental processing by a Milky-Way host (e.g. \citealt{Applebaum2021}).

\section*{Acknowledgements}
We thank Tobias Buck and the anonymous referees for their comments that helped clarify this manuscript. MP, MR, EA, ASO, OA and FR acknowledge support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grants 2014-5791 and 2019-04659) and the Royal Physiographic Society of Lund. MR is further supported by the Beecroft Fellowship funded by Adrian Beecroft. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 818085 GMGalaxies. AP is further supported by the Royal Society. This work was performed in part using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This work was also supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s890. Simulations were performed in part using computational resources at LUNARC, the centre for scientific and technical computing at Lund University, on the Swedish National Infrastructure for Computing (SNIC) allocation 2018/3-649, as well as allocation LU 2018/2-28 thanks to financial support from the Royal Physiographic Society of Lund. The authors acknowledge the use of the UCL Grace High Performance Computing Facility, the UCL Hypatia facility provided by the UCL Cosmoparticle Initiative, the Surrey Eureka supercomputer facility, and their associated support services.

\section*{Data Availability}
The data underlying this article will be shared on reasonable request to the corresponding author.

\bibliographystyle{mnras}