JWST Reveals Widespread AGN-Driven Neutral Gas Outflows in Massive z ~ 2 Galaxies

We use deep JWST/NIRSpec R~1000 slit spectra of 113 galaxies at 1.7<z<3.5, selected from the mass-complete Blue Jay survey, to investigate the prevalence and typical properties of neutral gas outflows at cosmic noon. We detect excess Na I D absorption (beyond the stellar contribution) in 46% of massive galaxies ($\log$ M$_*$/M$_\odot>$ 10), with similar incidence rates in star-forming and quenching systems. Half of the absorption profiles are blueshifted by at least 100 km/s, providing unambiguous evidence for neutral gas outflows. Galaxies with strong Na I D absorption are distinguished by enhanced emission line ratios consistent with AGN ionization. We conservatively measure mass outflow rates of 3 - 100 $M_\odot$ yr$^{-1}$; comparable to or exceeding ionized gas outflow rates measured for galaxies at similar stellar mass and redshift. The outflows from the quenching systems (log(sSFR)[yr$^{-1}$] $\lesssim$ -10) have mass loading factors of 4 - 360, and the energy and momentum outflow rates exceed the expected injection rates from supernova explosions, suggesting that these galaxies could possibly be caught in a rapid blowout phase powered by the AGN. Our findings suggest that AGN-driven ejection of cold gas may be a dominant mechanism for fast quenching of star formation at z~2.


INTRODUCTION
Determining the physical mechanism(s) responsible for quenching star-formation in massive galaxies is key to our understanding of galaxy evolution.Cosmological simulations typically quench massive galaxies via feedback from active galactic nuclei (AGN) which both expels cold gas from galaxies and heats halo gas, preventing it from cooling and being re-accreted to replenish the reservoir of fuel for star-formation (e.g.Di Matteo et al. 2005;Springel & Hernquist 2005;Bower et al. 2006;Croton et al. 2006;Hopkins et al. 2006;Somerville et al. 2008;Erb 2015;Beckmann et al. 2017).However, definitive observational evidence for a link between AGN feedback ★ Contact e-mail: rdavies@swin.edu.au and star-formation quenching has yet to be established (see Harrison 2017 and references therein).
Over the past decade, large galaxy surveys have significantly improved our understanding of outflows during the peak epoch of starformation and black-hole growth at  ∼ 1 -3, when feedback is expected to be most active.It is now well established that outflows are ubiquitous in massive star-forming galaxies at this epoch (e.g.Shapley et al. 2003;Weiner et al. 2009;Rubin et al. 2014;Harrison et al. 2016;Förster Schreiber et al. 2019).The strongest outflows are powerful enough to rapidly suppress star-formation in their host galaxies (e.g.Cano-Díaz et al. 2012;Cresci et al. 2015;Carniani et al. 2016;Kakkad et al. 2016;Davies et al. 2020), but these are generally associated with the most luminous AGN which are present in a small fraction of massive galaxies at any given time.It re-mains unclear whether outflows driven by more typical AGN are capable of quenching star-formation in their host galaxies.Measurements based on optical emission lines (tracing ionized gas) suggest that most outflows remove gas less rapidly than it is consumed by star-formation (e.g.Harrison et al. 2016;Leung et al. 2019;Förster Schreiber et al. 2019), whilst UV absorption line measurements (tracing neutral gas) suggest that the mass outflow rates are comparable to the star-formation rates of the host galaxies (e.g.Weiner et al. 2009;Kornei et al. 2012).
Observations based on a single gas phase provide a very incomplete picture of outflows which contain gas at a range of temperatures and densities including hot (10 6−7 K) X-ray emitting gas, warm (10 4−5 K) ionized gas, cool (100 K) neutral gas and cold (10 K) molecular gas (e.g.Strickland et al. 2004;Feruglio et al. 2010;Leroy et al. 2015;Krieger et al. 2019).Despite significant observational advances, the vast majority of outflows have only been observed in one gas phase, and constraining the total mass of gas ejected by outflows remains very challenging.State-of-the-art simulations of star-formation driven outflows predict that the majority of the outflowing mass is carried in the neutral and molecular phases (e.g Kim et al. 2020) (although it remains unclear whether cool clouds survive to large galactocentric radii or are shredded by the hot wind; e.g.Schneider et al. 2020;Fielding & Bryan 2022).Outflowing molecular gas is notoriously difficult to detect (see Veilleux et al. 2020, and references therein) but is often found to carry much more mass than the ionized phase (e.g.Fluetsch et al. 2019;Vayner et al. 2017;Brusa et al. 2018;Herrera-Camus et al. 2019).Cool neutral gas in outflows is commonly probed using low ionization rest-frame far-UV absorption lines.However, it is difficult to detect the far-UV continuum of massive, dusty AGN host galaxies, and UV absorption line measurements at high redshift are generally restricted to the strongest transitions which are often saturated, providing only lower limits on the outflowing mass (see Veilleux et al. 2020, and references therein).
An alternative tracer of neutral outflows is the resonant Na i d  5891,5897Å doublet.With a first ionization potential of 5.1 eV, Na i exists primarily in neutral regions where it is shielded by significant columns of gas and dust (e.g.Savage & Sembach 1996;Baron et al. 2020).Due to its location in the rest-frame optical spectrum, observations of Na i d already exist for many thousands of nearby galaxies.A few percent of local massive star-forming galaxies show blueshifted Na i d absorption indicative of neutral gas outflows (e.g.Nedelchev et al. 2019;Avery et al. 2022).In galaxies with both neutral and ionized outflows, the neutral outflow rates are 10 -100 times larger (e.g.Roberts-Borsani 2020; Avery et al. 2022;Baron et al. 2022), confirming that ionized gas likely represents a small fraction of the total mass budgets of typical nearby outflows.The Na i d absorption originates on spatial scales ≲ 10 kpc (e.g.Martin 2006;Rupke & Veilleux 2015;Rupke et al. 2017;Baron et al. 2020;Roberts-Borsani 2020;Avery et al. 2022;Rubin et al. 2022), indicating that it traces recently launched outflows rather than gas in the circumgalactic medium.Although Na i d is easily accessible in the local Universe, it has been significantly harder to detect at  ∼ 1 -3 where the line shifts into the observed near-infrared.
The unprecedented infrared sensitivity of JWST enables the detection of Na i d in distant galaxies, providing a new probe of neutral outflows in the early Universe.Initial observations have already revealed Na i d absorption tracing neutral outflows in three AGN host galaxies at  ∼ 2 -3, of which one is a quasar (Cresci et al. 2023;Veilleux et al. 2023) and two are post-starburst galaxies (Belli et al. 2023;D'Eugenio et al. 2023).These outflows are also detected in ionized gas emission lines, enabling direct comparisons of the mass outflow rates in different gas phases.Focusing on the post-starburst galaxies, both Belli et al. (2023) andD'Eugenio et al. (2023) find that the neutral mass outflow rates are about two orders of magnitude larger than the ionized outflow rates and exceed the current star formation rates (SFRs) of the host galaxies.The rapid ejection of cold gas by powerful AGN-driven outflows may have led to the recent, fast quenching of star-formation in these galaxies.
The detection of strong neutral gas outflows in two post-starburst AGN host galaxies provides tantalizing evidence that ejective AGN feedback may be an important mechanism for quenching massive galaxies at cosmic noon.However, it is unclear whether these objects are representative of the overall galaxy population.In this paper, we characterize the incidence and typical properties of neutral outflows across the galaxy population using 113 galaxies at 1.7 <  < 3.5 from the mass-selected Blue Jay survey.We discuss the sample and observations in Section 2, present the census of Na i d absorption in Section 3 and examine the neutral outflow properties in Section 4. We discuss the connection between neutral outflows, AGN activity and star-formation quenching in Section 5 and present our conclusions in Section 6.

Blue Jay
This work is based on observations from the JWST Cycle 1 program Blue Jay (GO 1810; PI Belli).The NIRSpec micro-shutter assembly (MSA; Ferruit et al. 2022;Rawle et al. 2022) was used to obtain R ≃ 1000 spectra of 151 galaxies spread over two masks in the COSMOS field.Four of these galaxies are filler targets at  ∼ 6, and the remaining 147 galaxies form a mass-selected sample (9 < log( * / ⊙ ) < 11.5) at cosmic noon (1.7 <  < 3.5).All galaxies were observed using the three medium resolution gratings (G140M, G235M and G395M) with exposure times of 13h, 3.2h and 1.6h respectively.A slitlet made of at least 2 MSA shutters was placed on each target and we employed a 2-point A-B nodding pattern along the slit.The data were reduced using a modified version of the JWST Science Calibration Pipeline v1.10.1, and version 1093 of the Calibration Reference Data System.Master background subtraction was performed using a spectrum measured from dedicated background slits and galaxy 1D spectra were optimally extracted (Horne 1986).The spectrum extraction failed for 6 galaxies which are excluded from our sample.Full details of the Blue Jay sample selection, observations and data reduction will be provided in the survey paper (Belli et al., in prep).
The individual grating spectra were combined to produce wide spectra covering rest-frame wavelengths of at least 3000Å -1.2m for all galaxies.Figure 1 shows the spectrum of COSMOS-10245 at  = 1.81 over rest-frame wavelengths of 3800 -6700Å.The wavelength coverage is not continuous due to gaps between the NIRSpec detectors, and we excluded 7 galaxies for which Na i d falls within a detector gap.Finally, we excluded 21 galaxies for which no spectroscopic redshift could be determined due to an absence of identifiable emission or absorption line features.Our final sample consists of 113 galaxies and includes COSMOS-11142, the post-starburst galaxy analysed by Belli et al. (2023).

Stellar Population Fitting
The goal of this paper is to search for neutral gas outflows traced by interstellar Na i d absorption.However, Na i d absorption can also originate in stellar atmospheres and is particularly prominent in latetype stars (e.g.O'Connell 1976;Peterson 1976;Carter et al. 1986;Alloin & Bica 1989;Worthey 1998).Therefore, it is imperative to accurately remove the stellar absorption contribution prior to our analysis.
We model the stellar continuum using Prospector, a Bayesian stellar population inference code designed to simultaneously fit photometry and spectroscopy spanning UV to mid-IR wavelengths (Johnson et al. 2021).We adopt the synthetic stellar population library fsps (Conroy et al. 2009;Conroy & Gunn 2010), the mist isochrones (Choi et al. 2016), and the Chabrier initial mass function.The stellar metallicity is free to vary and individual elemental abundances are assumed to be solar scaled.We note that varying the stellar [Na/Fe] ratio within a reasonable range does not significantly impact our results (see Section 3.1).We adopt a non-parametric starformation history with 14 bins spaced logarithmically in time except for the lowest age bin which is placed at 30 Myr.Prospector accounts for dust absorption and re-emission which are assumed to be in energy balance.The dust absorption model consists of a primary component which applies to all stars and follows the Kriek & Conroy (2013) attenuation curve, as well as a multiplicative term representing extra attenuation towards young stars (with ages < 10 Myr).The Prospector model also includes a multiplicative 'jitter' term that scales the measurement errors to better represent the statistical fluctuations in the data, as well as a polynomial distortion term that corrects for shape mismatches between the spectra and the stellar templates resulting from imperfect flux calibration and/or slit losses.
We use Prospector to fit the JWST spectra along with publicly available HST/ACS+WFC3 (Skelton et al. 2014;Momcheva et al. 2016) and Spitzer/IRAC (Laigle et al. 2016) photometry.The NIR-Spec observations cover many age-sensitive spectral features including the 4000Å break and the Balmer absorption series (see Figure 1), providing strong constraints on the stellar population properties.During the fitting, we mask prominent emission lines as well as the Na i d and Ca ii h + k absorption lines which can have significant contributions from interstellar gas.Full details of the Prospector fitting will be provided in Park et al. (in prep).
The orange curve in Figure 1 shows the best-fit stellar continuum model for COSMOS-10245.The model provides a very good fit to the well-detected Balmer absorption series and enables us to accurately quantify the stellar contribution to the observed Na i d absorption.In this case, the observed Na i d absorption is significantly stronger than expected from the stellar continuum alone.Unless otherwise noted, all subsequent references to Na i d absorption refer to absorption in excess of the stellar contribution.
Prospector outputs probability distribution functions which are used to calculate the best-fit values (median) and corresponding uncertainties (16th-84th percentile range) for all model parameters.These include the distortion polynomial (which we use to produce flux-calibrated spectra), the best-fit jitter term (used to scale the measurement errors), the stellar mass (M * ) and the non-parametric star-formation history.The SFRs reported in this paper refer to the SFR in the youngest age bin (averaged over the last 30 Myr), but similar results are obtained using SFRs averaged over 100 Myr or SFRs computed from the H emission line luminosity.

Emission and absorption line fitting
We fit the spectrum of each Blue Jay galaxy over the wavelength region between 3800 -6700Å, including contributions from the stellar continuum ( * ,Prospector ), ionized gas emission lines ( gas ) and excess Na i d absorption ( NaD,excess ): The emission and absorption line components must be fit simultaneously because the He i  5876Å emission line falls in close proximity to Na i d (see Figure 1).The emission and absorption line models are convolved with the wavelength-dependent NIRSpec line spread function prior to fitting1 .The emission line model includes the two strongest Balmer lines (H and H), the two strongest forbidden line doublets within the fitted wavelength region ([O iii]  4959,5007Å and [N ii]  6548,6583Å), and He i  5876Å.Outflowing gas can produce redshifted Na i d emission due to resonant scattering off gas in the receding side of the outflow (e.g.Prochaska et al. 2011), but this has only been observed in a handful of objects, (e.g.Rupke & Veilleux 2015;Perna et al. 2019;Baron et al. 2020Baron et al. , 2022;;Sun et al. 2023), and we do not find clear evidence for Na i d emission in any of our spectra.
We initially fit each emission line with a single Gaussian profile, constraining all lines to have the same velocity offset and dispersion.A single Gaussian component is sufficient to explain the vast majority of observed line profiles given the relatively low spectral resolution of the observations.Two galaxies show complex Balmer and forbidden line profiles that are not well represented by a single kinematic component, and for these objects we add an extra Gaussian component to all 7 emission lines.Two other galaxies show prominent AGN broad line region emission which is modelled as a broad Gaussian component in the H, H and He i lines.
Interstellar Na i d absorption is parametrized using the standard partial covering model (Rupke et al. 2005b): Here,   is the covering fraction of the absorbing gas against the background continuum source, and   () and   () are the optical depth profiles of the blue (Na i d  5891Å) and red (Na i d  5897Å) doublet lines, respectively.We assume that the optical depth has a Gaussian velocity distribution: The optical depth at the centre of the blue line ( 0, ) is fixed to be twice the optical depth at the centre of the red line ( 0, ), reflecting the known doublet ratio. 2 We fit each Na i d line using a single but note that the true resolution is notably higher than this for compact sources (e.g. de Graaff et al. 2023).As a consequence, the measured velocity dispersions represent lower limits on their true values. 2 The equivalent width ratio is not fixed and varies between 2 in the optically thin regime and 1 in the optically thick regime.This is because the curve of growth representing the relationship between optical depth and equivalent width is non-linear.
Gaussian velocity distribution, which is sufficient to describe the observed absorption profiles in all cases.The spectral resolution is comparable to the velocity separation between the doublet lines, making it very difficult to resolve multi-component velocity structure.
The Na i d absorption kinematics are allowed to vary independently of the emission line kinematics.
The absorption optical depth and covering fraction can become degenerate when the Na i d doublet is blended (e.g.Rupke et al. 2005b; see discussion in Appendix A).To obtain accurate constraints on the parameter uncertainties and degeneracies, we perform the fitting using emcee (Foreman-Mackey et al. 2013), an Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble sampler.The walkers are initialised in small regions around the best-fit values obtained from preliminary least squares fitting.
Similarly to the Prospector parameters, the best-fit emission and absorption line parameters represent the medians of the emcee posterior distributions and the error bars reflect the 16th -84th percentile ranges.The best-fit emission line and Na i d absorption line profiles for COSMOS-10245 are shown by the magenta and blue curves in Figure 1, respectively.

Incidence
We visually inspect the fits for all Blue Jay galaxies and identify sources with significant interstellar Na i d absorption.For the first time, we report evidence for widespread Na i d absorption in  ∼ 2 galaxies, as shown in Figure 2. The absorption profiles are grouped into four categories based on the velocity shift of the Na i d absortion feature (see Section 3.2).Within each category, each pair of panels represents an individual galaxy and shows the observed spectrum over the region covering ± 85Å around the Na i d doublet (top, black), the best-fit continuum-only (orange) and continuum + line (magenta) models, and the residual spectrum after removing the continuum and line emission (bottom, grey).From our initial sample of 113 galaxies, 30 galaxies (27%) have Na i d absorption much stronger than expected from the stellar populations alone.The MCMC posterior distributions confirm that the excess absorption is detected at ≥ 3 significance in all cases.
The main panel of Figure 3 shows how the Na i d detections (colored squares) are distributed as a function of stellar mass ( * ) and specific SFR (sSFR).Galaxies without Na i d absorption are shown with filled grey circles.It is clear that interstellar Na i d absorption is detected almost exclusively in massive galaxies.27/59 (46%) of the log( * / ⊙ ) > 10 galaxies with spectroscopic redshift measurements show Na i d absorption, whereas the detection fraction in lower mass galaxies is 6%.Interestingly, the detections are spread almost uniformly over four orders of magnitude in sSFR, from highly star-forming galaxies to quenching systems.The prevalence of interstellar Na i d absorption in the mass-selected Blue Jay sample indicates that large neutral gas reservoirs are prevalent in massive  ∼ 2 galaxies.
The incidence of interstellar Na i d absorption is not strongly dependent on the assumed sodium abundance.We generate stellar absorption profiles for a range of sodium abundances using the alf code (Conroy & van Dokkum 2012;Conroy et al. 2018), and find that extremely sodium-enhanced stellar populations ([Na/Fe] = +0.6)can produce excess absorption with a rest-frame equivalent width of up to 1.2Å.This is weaker than all our observed excess absorption profiles which have measured equivalent widths (after removing the stellar contribution) of 1.5 -11.4 Å (median 4.4Å; see Table 1).

Redshift (inflow/ interaction) Type 1 AGN (broad HeI)
Figure 2. Spectral cutouts covering 5805 -5975Å for the 30 galaxies with Na i d absorption significantly exceeding the stellar contribution.Each pair of panels represents an individual galaxy and shows the observed spectrum (top, black), the best-fit continuum-only (orange) and continuum + line (magenta) models, and the residual spectrum after removing the continuum and line emission (bottom, grey).Numbers in the lower panels match the sample IDs in Table 1.The -axes are scaled independently for each panel.The dotted vertical lines show the systemic wavelengths of the blue and red Na i d lines in the galaxy rest frames, and the blue and red symbols in the lower (residual) panels show the centroids of the best-fit line profiles (centered at the 50th percentile velocity with horizontal error bars indicating the 16th -84th percentile confidence intervals).The absorption profiles are grouped into four categories: blueshifted (top), systemic (middle), Type 1 AGN with broad He i emission (bottom left) and redshifted (bottom right).Some absorption profiles appear shifted in velocity but are classified as systemic because the 68% confidence interval encompasses zero velocity.

Stacking to search for weaker NaD absorption
The detection of interstellar Na i d absorption requires a robust measurement of the stellar continuum, so our Na i d sample may be biased towards the brightest continuum sources.We search for interstellar absorption in galaxies lacking individual Na i d detections by stacking them in two bins of stellar mass (above and below log( * / ⊙ ) = 10).The spectra are continuum-normalized prior to stacking and weighted by the average continuum signal-to-noise ratio within 150Å of the Na i d line.Unweighted stacks are noisier but lead to the same overall conclusions.
The bottom panels of Figure 3 show stacked spectra of galaxies with and without individual Na i d detections (black and grey, respectively).Even after stacking 51 galaxies, we do not find any evidence of excess Na i d absorption in low mass galaxies lacking individual Na i d detections.There are at least two possible reasons for this.Stacking does not reveal any additional absorption in low mass galaxies.There is some evidence for weak systemic absorption in high mass galaxies lacking individual Na i d detections.
Firstly, Na I has an ionization potential of 5.1 eV and therefore cannot exist in large quantities without a significant amount of dust shielding.In the local Universe there is a well-known relationship between Na i d absorption strength and dust attenuation (quantified by e.g. the  band line-of-sight attenuation   , the color excess E(B-V), or the Balmer decrement  (H)/  (H); e.g.Heckman et al. 2000;Veilleux et al. 2005;Chen et al. 2010;Veilleux et al. 2020;Avery et al. 2022, and this correlation is also seen in our sample (see left-hand panel of Figure 4).Most low mass galaxies likely do not have enough dust and metals to retain detectable amounts of neutral sodium.A second possibility is that we do not see strong excess Na i d absorption because it primarily traces AGN-driven outflows (see Sections 3.2.1 and 4), and these are rare in low mass galaxies (e.g.Genzel et al. 2014;Förster Schreiber et al. 2019;Leung et al. 2019).
There is some evidence for weak excess Na i d absorption in the stack of 32 massive galaxies lacking individual detections.This excess absorption falls at approximately zero velocity and could plausibly be explained by additional stellar absorption (see Section 3.2.1).Regardless, it is significantly weaker than the absorption seen in individually detected galaxies.The fact that we see strong interstellar Na i d absorption in 46% of massive galaxies whilst the remaining 54% show weak or no absorption suggests that the distribution of absorption strengths is not continuous, and may be more bimodal.

Link between strong Na i d absorption and galaxy properties
The presence of strong Na i d absorption in massive galaxies does not appear to be governed by dust properties: when restricting the sample to galaxies with log( * / ⊙ ) > 10, there is no significant difference between the   or E(B-V) distributions galaxies with and without Na i d absorption (see histograms in the upper left panel of Figure 4).We investigate whether the Na i d absorption strength may instead be related to galaxy inclination.Detections of strong excess Na i d absorption in the local Universe seem to be preferentially associated with outflows (e.g.Heckman et al. 2000;Rupke et al. 2005b;Cazzoli et al. 2016;Rupke et al. 2017).Blueshifted wind material is primarily observed in face-on galaxies (e. and the galaxy axis ratio measured from CANDELS HST imaging (van der Wel et al. 2012).However, this imaging covers rest-frame UV and blue-optical wavelengths where young stellar populations dominate.Ongoing JWST NIRCam and MIRI imaging of the COS-MOS field will enable more complete mapping of the galaxy stellar mass distributions and thus more accurate axis ratio measurements in the future.
Intriguingly, the most striking difference between the galaxies with and without Na i d absorption is their [N ii]/H ratios, which are discussed further in Section 4.1.

Origin of absorbing gas
The origin of the interstellar Na i d absorption can be determined by examining the velocity of the absorption relative to the galaxy systemic velocity (shown by the vertical dotted grey lines in Figure 2).In many cases, the excess Na i d absorption is clearly offset from the galaxy systemic redshift, with measured velocity shifts ranging from -680 to +400 km s −1 .We use the absorption velocity posterior probability distributions to classify the galaxies into three categories: blueshifted absorption (84th percentile velocity less than zero), redshifted absorption (16th percentile velocity greater than zero), and systemic absorption (consistent with zero).These classifications are shown in Figure 2 and determine the colors of the markers in Figure 3.The absorption velocities for the two broad line AGN are very sensitive to the fitting of the broad He i emission, so we cannot reliably classify them.Of the remaining 28 galaxies, 14 (50%) show blueshifted absorption (discussed in Section 3.2.1),11 (39%) are consistent with the systemic velocity (Section 3.2.2),and 3 (11%) show redshifted absorption (Section 3.2.3).

Outflowing gas
Half of the classifiable absorption profiles (14/28 or 50%) are blueshifted by at least 100 km s −1 (50th percentile velocity), which is an unambiguous sign of neutral gas outflows (e.g.Phillips 1993;Heckman et al. 2000;Rupke et al. 2002;Schwartz & Martin 2004;Martin 2005;Rupke et al. 2005b).The high fraction of outflows among the Na i d-detected galaxies is consistent with studies of local (ultra-)luminous infrared galaxies which find that excess Na i d absorption is preferentially blueshifted and associated with winds (e.g.Heckman et al. 2000;Rupke et al. 2005b;Cazzoli et al. 2016;Rupke et al. 2017).The overall incidence of Na i d outflows across the full mass range of our sample (log( * / ⊙ ) = 8.5 -11.5) is 12% (14/113).This is notably higher than the ∼ 1% incidence of neutral outflows among typical star-forming and AGN host galaxies with similar stellar masses at  ∼ 0 (e.g.Nedelchev et al. 2019;Avery et al. 2022), but comparable to the ∼ 20% incidence of outflows in local post-starburst galaxies (e.g.Sun et al. 2023).Compared to other measurements at  ∼ 2, the observed frequency of Na i d outflows is lower than the 20 -30% incidence of ionized outflows traced by optical emission lines among massive galaxies (e.g.Förster Schreiber et al. 2019), and significantly lower than the ∼ 100% incidence of neutral outflows traced by rest-frame UV absorption lines in UVbright galaxies (e.g.Steidel et al. 2010).We note that the Na i d outflow fraction may be underestimated by up to a factor of 2 at high stellar masses due to the low spectral resolution of the observations (R ∼ 1000).Galaxies with strong systemic absorption could have weaker outflow components that would not be detectable in our data (see Section 3.2.2).The low incidence of Na i d outflows in low mass galaxies is likely driven to a large degree by the requirement for dust shielding which prevents the detection of Na i d absorption in UV-bright galaxies (see also Avery et al. 2022).Rest-frame UV and Na i d absorption lines may probe outflows in almost entirely separate populations of galaxies.
The incidence of neutral outflows in our sample appears to be independent of star-formation activity.The outflow sources, indicated by the blue squares in Figure 3, are distributed over a wide range in sSFRs, extending all the way to the quenching galaxy regime.This is somewhat surprising given that in the local Universe, blueshifted Na i d absorption is preferentially found in highly star-forming systems (e.g.Rupke et al. 2005c;Chen et al. 2010;Concas et al. 2019;Roberts-Borsani & Saintonge 2019;Avery et al. 2022).Excess Na i d absorption has been found in many nearby massive quiescent galaxies (e.g.Carter et al. 1986;Alloin & Bica 1989;Worthey 1998;Thomas et al. 2003;Worthey et al. 2011;Concas et al. 2019;Roberts-Borsani & Saintonge 2019), but this absorption is typically consistent with the systemic velocity.Furthermore, the Na i d absorption strength in local early type galaxies is found to correlate with that of the Mg b triplet, which only arises in stellar atmospheres (e.g.Heckman et al. 2000;Rupke et al. 2002Rupke et al. , 2005b;;Alatalo et al. 2016).Therefore, the excess Na i d absorption is generally attributed to enhanced stellar absorption, perhaps due to elevated [Na/Fe] (e.g.O'Connell 1976;Peterson 1976;Parikh et al. 2018) and/or a bottom-heavy initial mass function (e.g.van Dokkum & Conroy 2010, 2012;Spiniello et al. 2012).
The excess Na i d absorption we detect in low sSFR galaxies at  ∼ 2 has distinctly different properties from the excess stellar absorption seen at  ∼ 0. Firstly, all sources classified as outflows have Na i d absorption blueshifted by at least 100 km s −1 , indicating that a non-systemic component is required to explain the observed absorption profile.Secondly, the observed Na i d absorption is too strong to explain with excess stellar absorption (see Section 3.1).Jafariyazani et al. (2020) similarly detected excess Na i d absorption in a lensed quiescent galaxy at  ∼ 2 and showed that it was too strong to be explained by enhanced Na abundance.Thirdly, the excess Na i d absorption is not associated with excess Mg b absorption.The righthand panel of Figure 4 shows stacked residual spectra of galaxies with neutral gas outflows, zoomed in on the region around the Mg b line.There is no evidence for significant excess absorption, suggesting that the stellar absorption contribution has been fully accounted for in the Prospector fitting.Furthermore, the outflow host galaxies in our sample fall above the Na i d -Mg b correlation observed in local early-type galaxies (Alatalo et al. 2016), indicating that the Na i d absorption is much stronger than expected based on the Mg b absorption.Finally, the Na i d absorption strength (quantified by the rest-frame equivalent width  (NaD)) is positively correlated with   (Figure 4, left), indicating that the excess Na i d absorption is interstellar in origin.The combination of these factors provides strong evidence that the blueshifted Na i d absorption observed in massive, low sSFR galaxies at  ∼ 2 traces neutral gas outflows.

Systemic absorption
A further 39% (11/28) of the classified Na i d absorption profiles have centroid velocities consistent with the galaxy systemic velocity.This means that most of the neutral gas follows the bulk motion of the galaxy i.e. it is likely to be located primarily in the interstellar medium (ISM).Massive, high redshift galaxies are known to harbour large cold gas reservoirs (see Tacconi et al. 2020, and references therein), and these could be responsible for producing the strong systemic absorption we observe.However, it is also possible that a non-negligible fraction of the absorption classified as systemic could arise from outflows.Firstly, we are only able to robustly identify outflows with velocity offsets exceeding ∼ 100 km s −1 due to the relatively low spectral resolution of the observations and the close relative proximity of the Na i d doublet lines.Secondly, our classification based on the average absorption velocity would not identify outflow components that are hidden underneath strong ISM absorption.Observations at higher spectral resolution (i.e. = 2700 with JWST/NIRSpec) would enable us to perform 2-component fitting and separate neutral gas in the ISM from outflowing material (see Section 4.3).

Infalling gas
The remaining 11% (3/28) of the classified Na i d absorption profiles are redshifted, with velocity offsets of 27, 255 and 404 km s −1 .We note that although the absorption towards COSMOS-19572 is formally classified as infalling, the absolute velocity offset of 27 km s −1 is small and comparable to the measurement error, and therefore this absorption could plausibly be systemic.The other two redshifted absorption sources have significantly larger velocity shifts, suggesting that there is neutral gas flowing towards these galaxies.The infalling material could originate in bulk flows within interacting systems or could be directly accreting onto the galaxies, providing cold gas that may sustain (or rejuvenate) star-formation.
Redshifted Na i d absorption has been observed in local galaxies.Roberts-Borsani & Saintonge (2019) find that redshifted absorption is prevalent in massive edge-on star-forming galaxies, consistent with a picture where the absorption traces gas accreting along the disk plane, potentially originating from galactic fountains.Roy et al. (2021) find evidence for neutral gas inflows in a substantial fraction of passive 'red geyser' galaxies, again likely originating from internal recycling and/or minor mergers (see also Cheung et al. 2016).
We investigate the likely origin of redshifted absorption in the Blue Jay galaxies by examining their star formation histories and morphologies (using HST/ACS+WFC3 imaging from the 3D-HST survey and JWST/NIRCam imaging from the PRIMER survey; GO 1837, PI Dunlop).COSMOS-19572 shows evidence for a nearby companion and could be an interacting system.This galaxy also shows strong N ii and S ii line emission which could plausibly trace shocks induced by tidal forces.High spatial resolution maps of emission line fluxes and kinematics would help to determine whether interactions are significantly impacting the dynamical state of this system.The remaining two galaxies (COSMOS-21452 and COSMOS-13174) do not show any evidence for multi-component structure, suggesting they may host neutral gas inflows.Interestingly, both galaxies are at the peak of their star-formation histories, suggesting that the high SFRs may be fuelled by ongoing cold gas accretion.

OUTFLOW PROPERTIES
We have reported the first evidence for widespread Na i d absorption in massive (log( * / ⊙ > 10) galaxies at cosmic noon, revealing that these galaxies have large neutral gas reservoirs.Approximately half of the detected absorption profiles are blueshifted, providing unambiguous evidence of neutral gas outflows.Other galaxies may have weaker outflows which are undetected at R ∼ 1000 due to the presence of strong ISM absorption.In this section, we investigate the properties and principal driving mechanisms of the detected neutral outflows.
The first clue regarding the driving mechanism comes from the outflow demographics: the Blue Jay neutral gas outflows are spread almost uniformly over more than four orders of magnitude in sSFR (see Figure 3).This is in tension with expectations for star-formation driven outflows, for which the incidence should increase with SFR.However, the incidence of AGN-driven ionized gas outflows at cosmic noon is observed to be independent of sSFR (e.g.Leung et al. 2019;Förster Schreiber et al. 2019), suggesting that the neutral gas outflows we detect may be AGN-driven.
We further explore the link between neutral outflows and AGN activity by examining galaxy emission line ratios (Section 4.1) and the outflow velocities, mass outflow rates and energetics (Sections 4.2,4.3 and 4.4,respectively).

Emission line ratios
Optical emission line ratios are valuable diagnostics of the principal power sources within galaxies (e.g.Baldry et al. 2008;Veilleux & Osterbrock 1987;Kewley et al. 2001;Kauffmann et al. 2003).Figure 5 shows that massive galaxies with strong interstellar Na i d absorption (black) have distinctly different emission line ratios from those without strong Na i d absorption (brown).The line ratios of individual galaxies are shown in small markers and the histograms show the line ratio distributions for the two populations.Large triangles show the ratios measured from median stacked profiles.We have verified that similar values are obtained from mean stacked profiles and by averaging the individual line ratio measurements.Galaxies lacking significant Na i d absorption typically lie in the star-forming and composite regions of the [N ii]/H vs. [O iii]/H diagnostic diagram (Baldwin et al. 1981;Kewley et al. 2001;Kauffmann et al. 2003) and fall close to the locus of  ∼ 2.3 star-forming galaxies from the MOSDEF survey (Shapley et al. 2015).In contrast, galaxies with detected Na i d absorption have significantly larger [N ii]/H ratios consistent with AGN host galaxies at similar redshifts (e.g.Coil et al. 2015).In the local Universe, elevated [N ii]/H ratios can alternatively trace shock-excitation in star-formation driven outflows (e.g.Sharp & Bland-Hawthorn 2010).However, star-formation driven outflows at  ∼ 2 typically do not show very elevated [N ii]/H ratios (e.g.Newman et al. 2012a;Davies et al. 2019;Freeman et al. 2019), perhaps because the detected line emission primarily originates from regions close to the galaxy disk where ionizing radiation from young stars dominates.Therefore, we hypothesize that strong Na i d absorption is preferentially associated with AGN activity.

Outflow Velocity
Observations of outflows in ionized, molecular and neutral gas have shown that AGN-driven outflows typically have more extreme velocities than star-formation-driven outflows, where velocities ≳ 1000 km s −1 are primarily associated with AGN activity (e.g.Sturm et al. 2011;Rupke & Veilleux 2013;Arribas et al. 2014;Cicone et al. 2014;Harrison et al. 2016;Förster Schreiber et al. 2019).We estimate the outflow velocities for the Blue Jay targets using the velocity offset  and dispersion  of the absorption profiles:  out = || + 2.The measured outflow velocities range from 200 -1100 km s −1 , with a median value of ∼ 500 km s −1 .The fastest of the Blue Jay outflows are more likely to be AGN-driven than star-formation-driven, but the velocity information is insufficient to determine the driving mechanisms of the more moderate velocity outflows.Comparison between the emission line ratios of massive galaxies (log(  * / ⊙ ) > 10) with significant Na i d absorption (black), massive galaxies without Na i d absorption (brown), and low mass galaxies (grey).Small circles represent line ratios measured for individual galaxies with all four diagnostic emission lines detected, and the histograms show the line ratio distributions for each class of galaxies.Large triangles show line ratios measured from median stacked profiles of all galaxies in each class.The solid black curve is the theoretical upper bound for pure star-formation (Kewley et al. 2001) and the dashed black curve is the empirical upper bound of the  ∼ 0 star-forming galaxy locus (Kauffmann et al. 2003).The pink curve shows the best-fit locus of  ∼ 2.3 star-forming galaxies from the MOSDEF survey (Shapley et al. 2015).

Calculations
We estimate the neutral gas outflow rates using the time-averaged shell model presented in Rupke et al. (2005b) and updated in Baron et al. (2022): where  Ω is the large-scale covering factor related to the opening angle of the wind, N(H i) is the hydrogen column density,  out is the outflow radius and  out is the outflow velocity.
The small-scale covering fraction   is obtained directly from the line fitting (see Equation 2).The measured values range from 0.1 -0.9 (median 0.3); similar to what has been found in the local Universe (e.g.Avery et al. 2022).We assume that the outflows cover 50% of the solid sphere (i.e. Ω = 0.5), consistent with the incidence of neutral outflows in local infrared galaxies (e.g.Rupke et al. 2005a).The geometry of neutral outflows at  ∼ 2 is very uncertain, but the fact that we detect outflows in ≳ 25% of massive galaxies suggests that the covering fraction cannot be much smaller than 0.25.We therefore consider the systematic uncertainty on  Ω to be a factor of 2 (i.e. Ω = 0.25 -1).
We calculate N(Na i) from the optical depth at the centre of the red Na i d line,  0, , using the relationship from Draine (2011): where  lu = 0.32 and  lu = 5897Å are the oscillator strength and restframe wavelength of the transition, respectively, and  is the Doppler parameter, equivalent to √ 2.By directly converting  0, to N(Na i), we are assuming that the observed absorption comes primarily from outflowing gas, with no significant contribution from gas in the ISM.The low spectral resolution of our observations means that we are unable to constrain multi-component fits allowing for contributions from both the ISM and outflows.Belli et al. (2023) found that ISM gas could account for up to 44% of the Ca ii k (and by extension Na i d) absorption in COSMOS-11142.It is unlikely that the ISM component contributes more than half of the observed absorption for sources classified as outflows, because the outflow component must dominate to produce the observed negative velocity shift.
We convert N(Na i) to N(H i) assuming Milky-Way-like Na abundance and dust depletion factors, and a 10% neutral fraction (see Rupke et al. 2005b).This neutral fraction is based on values measured towards Milky Way stars (Stokes 1978) and a cold extragalactic H i cloud (Stocke et al. 1991), and is likely to underestimate the ionization fraction in more extreme outflow environments.Baron et al. (2020) measured a 5% ionization fraction in a local AGN-driven outflow, which would increase the mass outflow rates by a factor of 2 compared to our calculations.
The radial extent of the outflowing neutral gas cannot be measured from our observations.We estimate the likely radial extent using size measurements of 1) neutral outflows in the local Universe and 2) ionized outflows at cosmic noon.Resolved studies of Na i d outflows in the local Universe suggest that they typically extend a few kiloparsecs, with measured sizes ranging from ∼ 1 -15 kpc (e.g.Martin 2006;Rupke & Veilleux 2015;Rupke et al. 2017;Baron et al. 2020;Roberts-Borsani et al. 2020;Avery et al. 2022).Similarly, Na i d absorption in quasar spectra is only observed within 15 kpc of galaxies (Rubin et al. 2022).The two Na i d outflows to have been spatially analyzed at cosmic noon have sizes ≤ 1 kpc (Cresci et al. 2023;Veilleux et al. 2023) and2.7 kpc (D'Eugenio et al. 2023).In comparison, ionized gas outflows at cosmic noon typically extend to at least the galaxy effective radius, on the order of a few kpc (e.g.Newman et al. 2012b;Davies et al. 2020;Belli et al. 2023).We conservatively adopt a 1 kpc extent and note that the mass outflow rates could be up to 10 times higher if the outflows are significantly larger than this.
When calculating the mass outflow rates we use the full Monte Carlo posterior probability distributions for , ,   and  0, .As mentioned in Section 2.3,   and  0, are degenerate because the Na i d doublet lines are blended in our observations.However, the mass outflow rate scales with the product of these two parameters (Equation 4), and the posterior probability distributions for the mass outflow rate are well constrained (see Appendix A for more details).

Results
The measured properties of the detected Na i d absorption profiles and the derived neutral gas masses and outflow rates are listed in Table 1.We measure outflow masses ranging from log( out / ⊙ ) = 7.0 -8.1 (median 7.6) and mass outflow rates spanning  out = 3 -100 M ⊙ yr −1 (median 17 M ⊙ yr −1 ).These are consistent with neutral outflow properties measured for star-forming and AGN host galaxies in the local Universe (e.g.Rupke et al. 2005a;Cazzoli et al. 2016;Baron et al. 2020;Roberts-Borsani 2020;Avery et al. 2022) as well as at  ∼ 2 (Perna et al. 2015;Cresci et al. 2023).The neutral gas outflow rates are also comparable to the ionized gas outflow rates measured for galaxies at similar stellar mass and redshift (e.g.Förster Schreiber et al. 2019).We emphasize that the estimates presented here are based on conservative assumptions for the outflow extent and Na ionization fraction, and the true outflow rates could plausibly be an order of magnitude larger.In the two  ∼ 2 -3 galaxies for which neutral and ionized mass outflow rates have been directly compared, the neutral outflow rates exceed the ionized outflow rates by approximately a factor of 100 (Belli et al. 2023;D'Eugenio et al. 2023).Our results emphasize that it is important to account for the neutral phase in order to paint a complete picture of ejective feedback.
Interestingly, we do not find any correlation between neutral gas outflow rate and galaxy SFR, as shown in the left-hand panel of Figure 6.As a consequence, the outflow mass loading factors (, defined as mass outflow rate divided by SFR) differ strongly between the star-forming and quenching populations (Figure 6, right).The outflows launched from the most actively star-forming galaxies (log(SFR)[M ⊙ yr −1 ] ≳ 0) have mass-loading factors of  ≲ 1, consistent with expectations for star-formation-driven outflows (e.g.Finlator & Davé 2008;Davé et al. 2011;Somerville & Davé 2015).In contrast, the mass loading factors for the lower SFR galaxies range from 4 -360.It is unlikely that energy injection by young stars could remove gas so much faster than the stars themselves are forming.However, many of the low SFR galaxies in our sample show strong Balmer absorption lines indicative of a recent rapid decline in SFR.We therefore investigate whether it is possible that the outflows were launched during a recent starburst phase, in which case the high mass loading factors could be an artefact of the time delay between the launch of the outflows and the SFR measurement (which is averaged over the last 30 Myr).The slowest outflow in our sample has a velocity of 210 km s −1 , and absorption line measurements have found that Na i d absorption is only observed within 15 kpc of galaxies (e.g.Rubin et al. 2022).This corresponds to a maximum reasonable outflow travel time of 70 Myr.We compute mass-loading factors using SFRs in different age bins from the Prospector fitting.Mass-loading factors of order unity are only found when using SFRs more than 100 Myr in the past.This is at least 50% longer than the maximum reasonable outflow travel time, suggesting that past starformation could not reasonably have powered these outflows and providing further evidence that they are driven by AGN activity.

Energy and momentum rates
Next, we investigate whether the current levels of star-formation and AGN activity in the Blue Jay galaxies are sufficient to explain the energetics of the neutral gas outflows.Figure 7 compares the kinetic energy and momentum rates of the outflows (  out and  out , respectively) with the rates of energy and momentum injection by supernovae and AGN.For the supernovae, we adopt the mechanical energy and momentum rate scalings from Veilleux et al. (2005) based on solar metallicity Starburst99 models (Leitherer et al. 1999):  SN = 7 × 10 41 × SFR[ ⊙ yr −1 ] erg/s,  SN = 5 × 10 33 × SFR[ ⊙ yr −1 ] dyne.The AGN bolometric luminosity is estimated from the [O iii] luminosity applying a bolometric correction factor of 600 (Netzer 2009).The [O iii] luminosity is corrected for extinction using the median   from the Prospector posterior probability distribution (including extra attenuation towards towards young stars), and the uncertainty on   is propa- emission in most of the outflow host galaxies is dominated by AGN activity (see Figure 5), and Netzer (2009) show that even for composite galaxies where more than half of the Balmer line emission is due to star-formation, the total [O iii] luminosity predicts the bolometric luminosity to within a factor of 2. Energy conserving AGN-driven outflows are expected to have kinetic energy rates equivalent to 5% of the AGN bolometric luminosity (see King & Pounds 2015, and references therein).The AGN momentum flux output is  AGN =  AGN /, but in energy-conserving outflows, the momentum outflow rate can be boosted by a factor of ∼ 5-20 due to entrainment of ISM gas in the wind (e.g.Faucher-Giguère & Quataert 2012).To account for this, we plot lines for outflow momentum rates equivalent to  AGN / and 20 ×  AGN /.
The top row of Figure 7 compares the outflow energetics with predictions for star-formation-driven outflows.In the most actively star-forming galaxies (log( SF )[erg s −1 ] ≳ 44), the energy injected by star-formation is likely sufficient to power the observed outflows.There is no correlation between  out and SFR, and as a consequence,  differs strongly between the star-forming and quenching populations.In high SFR galaxies (log(SFR)[ ⊙ yr −1 ] ≳ 0),  ≲ 1, consistent with expectations for self-regulating star formation feedback.However, the quenching galaxies have mass loading factors of 4 -360, suggesting that AGN activity is required to power the outflows in these systems.
However, in the lower SFR systems, the energy (momentum) injection rates are up to 60 (280) times larger than the predicted inputs.It is unlikely that this discrepancy can be explained by systematic uncertainties on the mass outflow rates.The shaded grey regions show the range of plausible values accounting for uncertainties on the ISM absorption contribution, ionization fraction and wind opening angle (see Section 4.3).We adopt a very conservative outflow size of  out = 1 kpc, and assuming a larger extent would only increase the outflow energy further above the energy and momentum injection by supernovae.This, together with the implausibly large mass-loading factors (see Section 4.3), provides strong evidence to suggest that the neutral gas outflows from the low SFR galaxies are unlikely to be powered by star-formation.
The bottom row of Figure 7 compares the outflow energetics with predictions for AGN-driven outflows.We see that in all cases, the AGN are powerful enough to drive the observed outflows.

DISCUSSION
Our investigation of the outflow driving mechanisms indicates that AGN activity likely plays a major role in powering the observed neutral gas outflows in massive  ∼ 2 galaxies.The incidence of neutral outflows is independent of (s)SFR (Figure 3).Galaxies with strong Na i d absorption show high [N ii]/H ratios consistent with AGN ionization, whereas galaxies without Na i d absorption show lower [N ii]/H ratios consistent with photoionization by young stars (Figure 5).Some outflows have velocities exceeding 500 km s −1 (Section 4.2).The case for AGN-driven outflows is particularly strong for low SFR galaxies where the mass loading factors range from 4 -360 (Figure 6) and the outflows are removing energy and momentum tens to hundreds of times faster than they can be injected by young stars (Section 7).
In summary, the neutral gas outflows in the Blue Jay sample are evenly distributed across star-forming and quenching galaxies, and AGN accretion appears to play a major role in driving these outflows.This is in contrast to the local Universe where the majority of neutral outflows are found in star-forming galaxies and are consistent with being star-formation-driven (e.g.Rupke et al. 2005c;Chen et al. 2010;Bae & Woo 2018;Concas et al. 2019;Nedelchev et al. 2019;Roberts-Borsani & Saintonge 2019;Avery et al. 2022).There is some evidence that outflows in local low SFR galaxies are preferentially associated with AGN activity (e.g.Concas et al. 2019;Roberts-Borsani & Saintonge 2019;Avery et al. 2022;Sun et al. 2023), consistent with our findings.The role of AGN in driving outflows may be enhanced in our sample because we are probing significantly brighter AGN: the median AGN luminosity of the Blue Jay outflow hosts (6 × 10 44 erg/s) is about two orders of magnitude higher than that of optically selected samples at  ∼ 0 (e.g.Roberts-Borsani & Saintonge 2019; Avery et al. 2022); consistent with the known redshift evolution in AGN luminosity (e.g.Rosario et al. 2012;Carraro et al. 2020).In the local Universe, neutral outflows from quasar host galaxies are faster and have higher mass outflow rates than outflows from star-forming galaxies (e.g.Veilleux et al. 2005;Rupke & Veilleux 2013;Cazzoli et al. 2016), suggesting that luminous AGN play a significant role in driving outflows at all redshifts.
Neutral outflows from low sSFR galaxies are much more prevalent at cosmic noon than in the local Universe.This may be because low sSFR galaxies at  ∼ 2 have had a lot less time to grow and quench than their  ∼ 0 counterparts, and as a result they have much younger stellar populations (e.g.Belli et al. 2019;Carnall et al. 2019;Tacchella et al. 2022).85% of the massive, low sSFR galaxies in the Blue Jay sample have light-weighted ages less than 1 Gyr (Park et al., in prep) and could therefore be post-starburst galaxies.The molecular gas reservoirs of post-starburst galaxies have been observed to decline with time since quenching (e.g.French et al. 2018;Bezanson et al. 2022;Baron et al. 2023), making it less likely to observe neutral gas outflows from the most evolved sources (although some older poststarburst galaxies do have detectable quantities of molecular gas; e.g.Rowlands et al. 2015).Both the incidence and velocity of Na i d outflows in local massive galaxies appear to decrease as galaxies age and quench (e.g.Sun et al. 2023).We do not find any evidence for a correlation between outflow velocity and either SFR or lightweighted age in our sample; however, it is likely that even if such a correlation was present, it would be masked due to the small sample size and relatively large measurement errors.3).The shaded region does not account for variations in outflow extent which could increase the outflow energy and momentum by up to an order of magnitude (moving the data points up).The outflows from the low SFR galaxies (log(L SF )[erg s −1 ] ≲ 44) are too powerful to be driven by star-formation, but could easily be powered by the observed AGN activity.
to be driven by powerful outflows, which have been observed in many post-starburst galaxies (both at  ∼ 0 and  ∼ 1; e.g.Tremonti et al. 2007;Davis et al. 2012;Alatalo 2015;Maltby et al. 2019;Baron et al. 2020Baron et al. , 2022)).The low sSFR outflow host galaxies in our  ∼ 2 sample may similarly trace a 'blowout' phase where strong AGN-driven outflows are ejecting large amounts of cold gas, leading to rapid quenching of star-formation.This picture is supported by detailed analyses of two post-starburst galaxies at  ∼ 2 -3 (Belli et al. 2023;D'Eugenio et al. 2023).Both galaxies experienced a burst of star-formation 0.3 -0.8 Gyr ago followed by a rapid decline in star-formation activity.These galaxies host powerful, AGN-driven neutral gas outflows that are ejecting cold gas 10 -100 times faster than it can be converted into stars.The rapid quenching of these systems may therefore be fully explained by ejective AGN feedback.
We have shown that similarly powerful neutral gas outflows are prevalent across the massive galaxy population at cosmic noon.The outflows from the quenching galaxies in our sample have massloading factors of 4 -360 (see Figure 6), consistent with the case studies above.Our work indicates that AGN-driven neutral gas outflows may represent a dominant avenue for fast quenching at  ∼ 2. Chemical evolution modelling of massive quiescent galaxies at  ∼ 1 has shown that mass-loading factors of order 10 are required to explain the stellar magnesium abundances, suggesting that these galaxies experienced powerful outflows prior to quenching (e.g.Leethochawalit et al. 2019;Zhuang et al. 2023).Many high redshift quiescent galaxies show emission line ratios consistent with AGN ionization (e.g.Belli et al. 2017;Newman et al. 2018;Belli et al. 2019;Park et al. 2022, Bugiani et al., in prep), providing additional evidence that AGN feedback plays a crucial role in quenching star-formation.
It is important to note that only a small fraction of the outflowing gas we detect may be able to escape the galaxy halos.The halo escape velocity is expected to be 2.5 -3 times the galaxy circular velocity  circ (e.g.Weiner et al. 2009;Swinbank et al. 2019).We are not able to measure circular velocities directly because we only have slit spectra, so we adopt  circ ∼ 300 km s −1 which is typical of AGN host galaxies at this redshift (Förster Schreiber et al. 2019).This corresponds to halo escape velocities of  esc,halo = 750 -900 km s −1 .Only 4/14 (29%) of the outflows in our sample exceed these velocities, suggesting that the majority of gas ejected in the outflows will remain in the halo and may eventually be re-accreted onto the galaxies.It is possible that the slowest outflows in our sample may not even escape from their host galaxies.However, it is difficult to robustly determine the fate of the outflowing gas for several reasons: 1) the errors on the outflow velocity measurements in Table 1 are quite large (the median error is ∼ 200 km s −1 ), 2) the measured velocity dispersions are likely under-estimated given the compact nature of our targets (see footnote 1), meaning that the outflow velocities are also underestimated, 3) the calculated outflow velocity depends on the assumed geometry and definition for  out , and 4) the escape velocity depends on the galactocentric radius which is unknown.Furthermore, Sun et al. (2023) found that outflows diminish with age, meaning that low SFR galaxies may have driven much faster outflows in the past.
Star-formation quenching likely involves a combination of gas removal and heating.Roy et al. (2021) found that a large fraction of radio-detected quiescent galaxies show infalling neutral gas probed by redshifted Na i d absorption, and the authors suggest that radio jets may be responsible for heating the accreted gas and preventing it from forming new stars.Slow outflows may similarly inhibit star-formation through turbulence or redistribution of gas within the galaxies (see also Luo et al. 2022;Sun et al. 2023).The powerful ejection of gas through outflows is crucial to explain the observed rapid quenching of galaxies in the early Universe (e.g.Belli et al. 2019;Park et al. 2022), whilst maintenance mode feedback is required to prevent rejuvenation and keep galaxies quiescent over long timescales.

SUMMARY AND CONCLUSIONS
We have used JWST/NIRSpec observations of 113 galaxies at 1.7 <  < 3.5 selected from the mass-complete Blue Jay survey to investigate the demographics and properties of neutral gas outflows, traced by Na i d absorption, at cosmic noon.Our observations have revealed for the first time that interstellar Na i d absorption is widespread in massive (log( * / ⊙ ) > 10) galaxies at  ∼ 2. Our main findings are as follows: • We detect interstellar Na i d absorption in 30/113 galaxies.The detections are almost exclusively associated with massive (log( * / ⊙ ) > 10) galaxies, for which the detection fraction is 46%.Lower mass galaxies likely have insufficient columns of gas and dust to shield Na i d against ionization.
• 50% of the Na i d absorption profiles are blueshifted by at least 100 km s −1 , providing unambiguous evidence for neutral gas outflows.These neutral outflows are observed across the entire massive galaxy population, with similar incidence rates in star-forming and quenching galaxies.
• 39% of the Na i d profiles are consistent with the galaxy systemic velocity.These primarily trace cool gas in the ISM, but may also have weaker underlying outflow components that are hidden at R ∼ 1000.
• 3 galaxies (11%) show redshifted absorption profiles indicative of infalling gas.Of these, one galaxy shows a complex morphology and emission-line ratios consistent with shock excitation, suggesting that the redshifted absorption may trace bulk flows of gas within an interacting system.The other two galaxies appear isolated and are at the peaks of their star-formation histories, suggesting that their star-formation may be fuelled by ongoing accretion of cool gas.
• Assuming a conservative outflow extent of 1 kpc, we compute neutral mass outflow rates of 3 -100  ⊙ yr −1 .These are comparable to or greater than ionized gas outflow rates previously reported for other galaxies with similar stellar masses and redshifts.Existing measurements of neutral gas outflow sizes range from 1 -15 kpc, so the true outflow extents and mass outflow rates from the Blue Jay galaxies could plausibly be up to an order of magnitude larger than we report.
• Multiple lines of evidence indicate that the outflows are likely to be AGN-driven.Galaxies with strong interstellar Na i d absorption have enhanced [N ii]/H ratios indicative of AGN activity.The outflow incidence does not depend on the level of star-formation activity.Star-formation cannot power the outflows from the low SFR galaxies, where the outflow mass loading factors range from 4 -360 and the energy and momentum outflow rates exceed the injection rates from supernovae by at least an order of magnitude.
• The presence of strong neutral outflows in quenching systems could indicate that they are undergoing a post-starburst 'blowout' phase powered by the AGN.The outflow velocities range from 200 -1100 km s −1 , albeit with large uncertainties (typical error 200 km s −1 ).It is difficult to robustly constrain the fate of the outflowing gas, but in most cases the estimated outflow velocities are lower than the expected halo escape velocities, suggesting that the bulk of the outflowing material will remain in the galaxy halos.Nevertheless, the strong AGN-driven ejection of cold gas provides a mechanism to explain the rapid quenching of star-formation in massive quiescent galaxies in the early Universe.Maintenance mode feedback (e.g. through radio jets) may also be required to prevent the re-accretion of cold gas and keep the galaxies quiescent.
Our results indicate that powerful, AGN-driven neutral gas outflows are prevalent across the massive galaxy population at  ∼ 2 and are likely to be a dominant channel for fast quenching at this epoch.

Figure 1 .
Figure 1.R ≃ 1000 NIRSpec MSA spectrum of COSMOS-10245 at  = 1.81 over rest-frame wavelengths of 3800 -6700Å (black).The orange curve shows the best-fit stellar continuum model from Prospector and the magenta and blue curves show the best-fit models for the strong emission lines and the Na i d doublet absorption, respectively, obtained from MCMC fitting.The inset in the bottom row shows a zoom-in on the region around the He I and Na i d lines.The observed Na i d absorption is significantly stronger than expected for the stellar component.

Figure 3 .
Figure3.Top: Distribution of galaxies in the  * -sSFR plane.Galaxies with detected Na i d absorption are shown as squares, where the colors reflect the velocity shift of the absorption as shown in Figure2.Filled grey circles indicate galaxies with no significant absorption, open circles indicate galaxies for which Na i d falls in a detector gap and 'x' symbols indicate galaxies with no spectroscopic redshift.We detect Na i d absorption in 46% of massive galaxies (log  * / ⊙ > 10), distributed almost uniformly over more than four orders of magnitude in sSFR.The panel on the right shows the typical  * and sSFR errors for galaxies in two sSFR bins divided at log(sSFR[Gyr −1 ]) = -10.Bottom: Stacked spectra of galaxies with and without individual Na i d detections (black and grey, respectively) in two stellar mass bins divided at log( * / ⊙ ) = 10.The two Na i d-detected galaxies with Type 1 AGN are omitted from the stacks.Stacking does not reveal any additional absorption in low mass galaxies.There is some evidence for weak systemic absorption in high mass galaxies lacking individual Na i d detections.

Figure 4 .
Figure 4. Left: Rest-frame equivalent width of excess Na i d absorption as a function of   for all Na i d-detected galaxies (black).The   values include both the primary attenuation component and the extra attenaution towards young stars.The   errors represent the 16th-84th percentile range from the Prospector posterior probability distributions.The Na i d absorption strength correlates with   , indicating that the absorption is interstellar in origin.Galaxies lacking Na i d detections are shown as brown circles with W(Na i d) artificially set to zero.The upper histograms show the   distributions of massive galaxies with and without detected Na i d absorption.The two distributions do not differ significantly.Right: Stacked residual spectra of outflow host galaxies, zoomed-in on the Mg b triplet.There is no significant residual Mg b absorption, indicating that the stellar absorption is fully accounted for in the Prospector fitting.
Figure5.Comparison between the emission line ratios of massive galaxies (log(  * / ⊙ ) > 10) with significant Na i d absorption (black), massive galaxies without Na i d absorption (brown), and low mass galaxies (grey).Small circles represent line ratios measured for individual galaxies with all four diagnostic emission lines detected, and the histograms show the line ratio distributions for each class of galaxies.Large triangles show line ratios measured from median stacked profiles of all galaxies in each class.The solid black curve is the theoretical upper bound for pure star-formation(Kewley et al. 2001) and the dashed black curve is the empirical upper bound of the  ∼ 0 star-forming galaxy locus(Kauffmann et al. 2003).The pink curve shows the best-fit locus of  ∼ 2.3 star-forming galaxies from the MOSDEF survey(Shapley et al. 2015).

Figure 6 .
Figure6.Mass outflow rate (  out , left) and mass loading factor (, right) as a function of SFR.There is no correlation between  out and SFR, and as a consequence,  differs strongly between the star-forming and quenching populations.In high SFR galaxies (log(SFR)[ ⊙ yr −1 ] ≳ 0),  ≲ 1, consistent with expectations for self-regulating star formation feedback.However, the quenching galaxies have mass loading factors of 4 -360, suggesting that AGN activity is required to power the outflows in these systems.

Figure 7 .
Figure 7. Measured energy (left) and momentum (right) rates of the outflows compared to the luminosity and momentum flux from young stars (top) and AGN (bottom).Dashed lines indicate the expected energy and momentum injection into the outflows.Shaded bands represent a factor of four uncertainty in either direction, accounting for potential ISM Na i d absorption and/or incorrect assumptions about the wind opening angle or ionization fraction (see discussion in Section 4.3).The shaded region does not account for variations in outflow extent which could increase the outflow energy and momentum by up to an order of magnitude (moving the data points up).The outflows from the low SFR galaxies (log(L SF )[erg s −1 ] ≲ 44) are too powerful to be driven by star-formation, but could easily be powered by the observed AGN activity.

Figure A1 .
Figure A1.Left: Outflow covering fraction   as a function of the measured maximum absorption depth.The dotted line indicates a 1:1 relation.The maximum absorption depth provides a lower limit on   .Right: Single and joint posterior probability distributions for the Na i d absorption parameters (velocity offset Δ, dispersion , covering fraction   , and optical depth  0, ) as well as the mass outflow rate  out for COSMOS-10565.  and  0, are degenerate but  out , which scales with the product of these two parameters, is well constrained.

Table 1 .
Properties of the detected Na i d absorption. 1) Sample ID, matching the labels in Figure2.
The rapid depletion of the molecular gas reservoirs is thought