galaxies: evolution — galaxies: star formation — galaxies: stellar content

#  HINOTORI I: The Nature of Rejuvenation Galaxies

Takumi S. Tanaka 11affiliation: Department of Astronomy, Graduate School of
Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan 22affiliation: Kavli Institute for the Physics and Mathematics of the
Universe, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
33affiliation: Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan ∗ Kazuhiro Shimasaku
11affiliationmark:  44affiliation: Research Center for the Early Universe,
Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan Sandro Tacchella 55affiliation: Kavli Institute for
Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
66affiliation: Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge CB3 0HE, UK Makoto Ando 11affiliationmark:  Kei Ito
11affiliationmark:  Hassen M. Yesuf 22affiliationmark:  77affiliation: Kavli
Institute for Astronomy and Astrophysics, Peking University, Beijing 100871,
People’s Republic of China Suin Matsui 11affiliationmark:
<EMAIL_ADDRESS>
###### Abstract

We present the Hinotori (star formation History INvestigatiOn TO find
RejuvenatIon) project to reveal the nature of rejuvenation galaxies (RGs),
which are galaxies that restarted their star formation after being quiescent.
As the first step of Hinotori, we construct the largest RG sample with 1071
sources. We select these RGs from 8857 MaNGA (Mapping Nearby Galaxies at APO)
survey galaxies by reconstructing their star formation histories with
Prospector spectral energy distribution fitting code. Both optical spectral
data and UV to IR photometric data are used for the fitting. Using mock data,
we confirm that our method can detect weak rejuvenation events that form only
about 0.1% of the total stellar mass with high completeness. The RGs account
for $\sim 10\%$ of the whole sample, and rejuvenation events contribute on
average only about 0.1% of the total stellar mass in those galaxies but 17% of
the cosmic-star formation rate density today. Our RGs have a similar mass
distribution to quiescent galaxies (QGs). However, the morphology of the RGs
is more disk-like than QGs, suggesting that rejuvenation may occur selectively
in disk-like QGs. Our results also suggest the possibility of multiple-time
rejuvenation events in a single galaxy. Further spatially resolved analyses of
integral field unit data and radio observations and comparisons to simulations
are needed to identify the mechanism and the role of rejuvenation in galaxy
evolution.

##  1 Introduction

Galaxies segregate into two distinct groups in the star-formation rate (SFR) -
stellar mass ($M^{*}$) diagram: “star-forming galaxies” (SFGs) with active
star formation and “quiescent galaxies” (QGs) with low-to-no star formation
(e.g., [Renzini et al. (2015), Feldmann (2017)]). Because the fraction of QGs
increases with cosmic time (e.g., [Ilbert et al. (2013), Muzzin et al. (2013),
Davidzon et al. (2017), Weaver et al. (2022)]), a widely-accepted simple
galaxy evolution scenario is that SFGs stop star formation (quenching) and
become QGs (e.g., [Faber et al. (2007), Peng et al. (2010)]). Galaxies in
between SFGs and QGs are called “green valley” galaxies (GVs). GVs are widely
thought to be galaxies in transition from SFGs to QGs (e.g., [Salim (2014),
Angthopo et al. (2020)]).

Recent progress in the spectral energy distribution (SED) fitting technique
([Walcher et al. (2011)] for review) has enabled us to search for and analyze
galaxies with unique star-formation histories (SFHs), such as “rejuvenating
galaxies”, which have been overlooked in previous studies. Rejuvenating
galaxies have resumed their star formation after being quiescent (e.g.,
[Trayford et al. (2016), Cleland and McGee (2021), Zhang et al. (2022)]). We
may need to update the simple evolution scenario from SFGs to QGs to explain
rejuvenating galaxies.

Some studies (e.g., [Chauke et al. (2019), Mancini et al. (2019), Tacchella et
al. (2022)]) have also focused on “rejuvenated galaxies”. The difference
between “rejuvenating galaxies” and “rejuvenated galaxies” is that the latter
have terminated their resumed star-formation and have been quiescent again.
Since both populations have experienced rejuvenation events, this paper deals
with both populations and collectively refers to them as “rejuvenation
galaxies” (RGs). We do not use the abbreviation “RG” when distinguishing
between “rejuvenated galaxies” and “rejuvenating galaxies”.

RGs have been studied since the 2000s. Early studies, which are limited to the
low-$z$ universe, have selected early-type galaxies (ETGs) with recent star
formation as RGs by using UV detection ([Kaviraj et al. (2007), Donas et al.
(2007), Schawinski et al. (2007)]), spectral and photometric features ([Treu
et al. (2005), Cleland and McGee (2021)]), or a stellar population synthesis
([Thomas et al. (2010)]). The fraction of RGs (or ETGs with recent star
formation) in previous samples ranges from $10\%$ to $30\%$ depending on the
selection method and the properties of the parent samples. Recent developments
in SED fitting methods and spectroscopic surveys in the intermediate-$z$
universe enabled selecting $0.6\lesssim z\lesssim 2$ RGs based on their SFHs
or stellar ages reconstructed by SED fitting ([Belli et al. (2017), Gobat et
al. (2017), Carnall et al. (2019), Chauke et al. (2019), Mancini et al.
(2019), Akhshik et al. (2021), Tacchella et al. (2022), Paspaliaris et al.
(2023)]).

The mechanism and cause of rejuvenation and the role of rejuvenation in galaxy
evolution are major open questions. On the cause of rejuvenation, some studies
have suggested that rejuvenation is related to mergers. For example, Kaviraj
et al. (2009) have used a simulation to show that a minor merger can explain
the recent star formation of ETGs. Other studies have also suggested that gas
accretion onto elliptical galaxies after merging may cause rejuvenation and
enables rejuvenated ellipticals to acquire a disk and turn into red spirals or
S0 galaxies (e.g., Mapelli et al. (2015); Diaz et al. (2018); Hao et al.
(2019); Rathore et al. (2022)). Paspaliaris et al. (2023) have argued that HI
gas, as found in some ETGs by Thom et al. (2012), could cause rejuvenation.
Note, however, that Park et al. (2022) have suggested that red disk-like QGs
cannot be described with the above merger-rejuvenation scenario. Besides,
Martin-Navarro et al. (2022) have reported that many type-I active galactic
nuclei (AGN) had experienced rejuvenation, suggesting the relationship between
rejuvenation events and AGN activities.

Understanding the rejuvenation mechanism will deepen our knowledge of the star
formation activity in a galaxy. In particular, we can consider RGs as former
QGs that could not maintain a quiescent state. Comparing RGs with QGs enables
us to discuss the quenching process, especially the conditions needed to
maintain a quiescent state.

The role of rejuvenation in galaxy evolution also remains to be evaluated.
Chauke et al. (2019) and Tacchella et al. (2022) have reported that, at most,
only about 10% of the total mass is formed in a rejuvenation event. However,
we can detect only recent rejuvenation events because the time resolution of
the reconstructed star formation history of a galaxy decreases with look-back
time. Thus, we may underestimate the contributions of rejuvenation events to
the total star-formation activities. The relationship between RGs and GVs is
also relevant to the role of rejuvenation. GVs are widely thought to be in the
quenching phase. However, this simple view needs to be modified if a
significant fraction of GVs are RGs. Furthermore, Mancini et al. (2019)
reported that the bending of the star-forming main sequence (SFMS) in the
high-mass end is caused by RGs.

Despite the importance of RGs, previous studies have failed to discuss the
nature of RGs statistically. This is mainly because the previous studies are
based on either small samples of RGs or biased parent samples (e.g., RGs are
selected from only QGs). For example, Tacchella et al. (2022) have found that
RGs reside in massive dark halos based on nine RGs selected from 161 QGs.
Cleland and McGee (2021) have detected no significant environmental dependence
of rejuvenation fraction in a sample of $\simeq 350$ RGs. However, these
results need to be confirmed with a larger RG sample and with a more inclusive
parent sample not limited to QGs. In addition, different definitions of RGs
among the previous studies have also made it difficult to compare their
results with each other and to compare them with simulation results (e.g.,
Kaviraj et al. (2009); Trayford et al. (2016); Pandya et al. (2017); Nelson et
al. (2018); Behroozi et al. (2019); Alarcon et al. (2022)).

We launch the Hinotori111star-formation History INvestigatiOn TO find
RejuvenatIon. Hinotori means phoenix in Japanese. project to solve those
problems. Hinotori aims to understand the nature of RGs using well-constrained
SFHs. Selecting RGs with reliable SFHs is the most direct selection method. As
the first step of Hinotori, we select RGs from a large ($N\sim 10^{4}$) parent
sample of $z\simeq 0$ galaxies covering all morphological types. To
reconstruct the SFH for individual galaxies, we run the Bayesian SED fitting
library Prospector (Leja et al., 2017; Johnson et al., 2021) on their
spectroscopic and UV to MIR photometric data. In this paper, we present the
most extensive catalog to date of SFH-selected RGs and make the first secure
statistical analysis of RGs.

Section 2 describes the data and how to make the input data for Prospector.
Section 3 describes the SED fitting method and the selection method of RGs. In
Sections 4 and 5, we show the results from a mock sample and real data,
respectively. We discuss the nature of RGs in Section 6 and present a
conclusion and prospects for future RG studies in Section 7. Throughout this
work, we assume the cosmological parameters from the WMAP-9 (Hinshaw et al.,
2013).

##  2 Data

###  2.1 MaNGA survey

We construct a parent sample from the sample of the MaNGA survey (Bundy et
al., 2015), an integral field unit survey in the SDSS-IV (Blanton et al.,
2017). This survey was carried out with the SDSS 2.5-meter telescope (Gunn et
al., 2006) equipped with a fiber feed system (Drory et al., 2015) and the
Baryon Oscillation Spectroscopic Survey (BOSS) spectrograph (Smee et al.,
2013) covering 3600 – 10300 Å at $R\sim 2200$. The target galaxies of the
MaNGA survey were selected from the NASA-Sloan atlas (NSA v1_0_1). See Yan et
al. (2016) for the survey design. As detailed in Wake et al. (2017), the MaNGA
sample consists of the following three samples selected with redshifts,
magnitudes, and colors to achieve a uniform number density distribution of
$M_{i}$ and additional targets selected for specific sciences.

  * •

Primary Sample: Galaxies whose $1.5R_{e}$ area is covered by the IFU.

  * •

Secondary Sample: Galaxies whose $2.5R_{e}$ area is covered by the IFU.

  * •

Color-Enhanced Sample: Selected to fill in the region on the ${\rm NUV}-i$
versus $M_{i}$ color-magnitude diagram poorly sampled in the Primary and
Secondary Samples. Galaxies in this region include low-luminosity red
galaxies, high-luminosity blue galaxies, and GVs. The IFU covers their
$1.5R_{e}$, the same as the Primary Sample. The Primary Sample and the Color-
Enhanced Sample are combined as the Primary+ Sample.

  * •

Ancillary Targets222https://www.sdss4.org/dr17/manga/manga-target-
selection/ancillary-targets/: Additionally selected 977 targets, such as
luminous AGN, galaxies in voids, pairs and mergers, dwarf galaxies, blight
cluster galaxies, and so on.

These galaxies are distributed over $5\times 10^{8}\leq M^{*}/M_{\odot}\leq
3\times 10^{11}$ at $0.01\leq z\leq 0.15$. In this study, we use the IFU data
of 11,273 datacubes in the SDSS DR17 (Abdurro’uf et al., 2022) reduced with
the version 3.1.1 MaNGA Data Reduction Pipeline (DRP: Law et al. (2016,
2021)), which performs flux calibration, sky subtraction, correction of the
Galactic dust attenuation, and coadding individual exposures to make the final
3D datacubes of each galaxy.

Assuming that the RG fraction, the number ratio of RGs to total galaxies, is
$f_{\rm RG}\sim 10\%$ as suggested by Nelson et al. (2018); Behroozi et al.
(2019); Chauke et al. (2019); Tacchella et al. (2022), we expect that the
entire MaNGA sample includes $N\sim 10^{3}$ RGs, enabling a reliable
statistical analysis. Besides, because all selected RGs have IFU spectroscopic
data, we can also analyze the rejuvenation mechanism with spatially resolved
characteristics.

###  2.2 Removing AGN, mergers, and pair galaxies

In this study, we do not use galaxies hosting active galactic nuclei (AGN)
because Prospector has no AGN template other than dust torus emission in IR
wavelengths (Leja et al., 2018). Comerford et al. (2020) published an AGN
catalog selected from the 8th MaNGA Product Launch (MPL-8) sample. However,
since this sample ($N=6261$ galaxies) is smaller than the DR17 sample we use,
we extend the selection method used by Comerford et al. (2020) to the DR17
MaNGA sample to exclude AGN. In total, we remove 406 AGN from the whole
sample.

Merger and pair galaxies are also excluded using the MANGA TARGET3
information, a flag for ancillary targets, because it is difficult to create a
1D spectrum of a target galaxy by subtracting the contribution from its
companion(s) in the 1D-flattening process described in section 2.4. We remove
119 galaxies that have either TARGET3_PAIR_SIM, TARGET3_PAIR_RECENTER,
TARGET3_PAIR_2IFU, or TARGET3_PAIR_ENLARGE flag.

###  2.3 Crossmatching with photometric catalogs

For SED fitting, we use both photometric data from UV to MIR and IFU-flattened
1D spectroscopic data (see Section 2.4) to obtain more accurate SFHs than the
existing SED fitting results, such as Pipe3D (Sánchez et al., 2016a, b, 2018)
and FIREFLY (Goddard et al., 2017; Parikh et al., 2018). For UV (FUV and NUV)
photometry from Galaxy Evolution Explorer (GALEX) and optical ($ugriz$)
photometry from the SDSS DR8, we use Sérsic fluxes summarized in the NSA and
correct them for Galactic extinction using Schlegel et al. (1998). We also use
IR photometry (W1, W2, W3, and W4) from the unWISE catalog (Schlafly et al.,
2019), which is deeper than the ALLWISE catalog due to using coadd images and
a forced photometry method (Lang et al., 2016). Furthermore, because the
forced photometry is based on galaxy positions and profiles from the SDSS, the
unWISE photometry is more consistent with the NSA photometry described above.
By crossmatching MaNGA galaxies with the unWISE catalog with a maximum
separation of 10 arcsec, we select 8857 sources with valid photometry (i.e.,
without missing values) as the main parent sample. The $S/N$ distribution for
each band and the redshift distribution of the sample are shown in figure 1.
We find that crossmatching with the unWISE catalog does not affect the
redshift distribution significantly.

Figure 1:  Summary of the photometry of the main parent sample. Histograms
show the $\log S/N$ distribution of our sample in each photometric band. The
first row represents FUV, NUV, $u$, $g$, $r$, $i$, and $z$, respectively, from
left to right. The second row represents W1, W2, W3, and W4, respectively,
from left to right. The median $S/N$ value is shown at the top of each panel.
The bottom right panel shows the redshift distribution. The blue line shows
the distribution of MaNGA galaxies matched with the unWISE catalog, i.e., our
parent sample, while the red line shows the distribution of MaNGA galaxies not
matched with the unWISE catalog. The p-value from the Kolmogorov–Smirnov (KS)
test is shown in the panel.

###  2.4 1D flattening of datacubes

To transform 3D IFU data (with two spatial and one spectral dimensions) into a
form that can be used for SED fitting, for each galaxy, we sum the spectra of
all spatial positions to create a 1D spectrum in the following manner.

First, based on MANGA_DRP3PIXMASK information, we exclude invalid data with
FORESTAR or DONOTUSE mask. FORESTAR means that the spaxel is contaminated by a
foreground star while DONOTUSE means that the data of the spaxel are
unreliable for some reason. We do not use a spaxel if it has even one
wavelength bin masked with either flag. The mask_1 panel in figure 2 shows the
positions removed by this operation in dark blue, and the mask_2 panel shows
the wavelengths removed by this operation (the value of 0 indicates that the
wavelength has been removed).

Then, we remove spaxels that either have an invalid value at least one
wavelength bin or are located near the edge of the field of view, i.e., at a
distance from the center being greater than 0.8 times the maximum distance of
all spaxels. The locations removed by this operation are shown in dark blue in
figure 2, mask_3.

In addition, we mask three wavelength ranges [5500 Å, 5600 Å], [5850 Å, 5950
Å], and [6250 Å, 6450 Å] because luminous airglow emission lines contaminate.
These wavelength ranges are shown in dark gray in the spectral and $S/N$
wavelength dependence diagrams in figure 2.

Finally, we coadd the spaxels not masked by the above operations to produce
1D-spectrum data. We calculate the errors in 1D spectra by assuming a normal
distribution of errors for each spaxel. An example coadded 1D-spectrum and its
$S/N$ are shown in the second and third lines in figure 2.

As mentioned in section 2.1, the size of the MaNGA field of view varies
depending on the sample type. The field of view covers $1.5R_{e}$ for the
Primary Sample and $2.5R_{e}$ for the Secondary Sample, where $R_{e}$ is the
effective radius. However, we do not perform aperture correction according to
sample type. This is because Prospector has the function of photospectral
correction as explained in section 3.1.2. Note that if there is non-negligible
radial dependence in spectra, e.g., strong star formation only outside
$1.5R_{e}$, there may be a bias between the Primary Sample and the Secondary
Sample.

Figure 2:  An example of our 1D flattening procedure. The first row shows the
0th moment of an example galaxy and the distributions of masks used in the
1D-flattening procedure. The second and third rows show a coadded 1D spectrum
and its $S/N$ distribution, respectively, with the red (blue) line indicating
the result with (without) masking.

###  2.5 Value-added catalogs

The MaNGA sample has various value-added
catalogs333https://www.sdss4.org/dr17/data_access/value-added-catalogs/ (VAC).
Another advantage of using MaNGA data set is those VACs, and we can discuss
the characteristics of RGs in more detail and more easily with such catalogs.

In this study, we use two VACs. First, we use “MaNGA PyMorph DR17 photometric
catalog” (Domínguez Sánchez et al., 2022) in section 5.4 to discuss the
morphology of the selected RGs. Second, we use “MaNGA Pipe3D value-added
catalog: Spatially resolved and integrated properties of galaxies for DR17”
(Pipe3D, Sánchez et al. (2022)) in appendix to validate the results of our SED
fitting. We describe the detail of each catalog and the used values in the
individual sections.

##  3 Methodology

###  3.1 SED fitting

SED fitting is a technique of fitting observed data (photometry and/or
spectroscopy) with model spectra (for a review, see Walcher et al. (2011)).
Through SED fitting, we can estimate basic parameters of galaxies such as
stellar mass ($M^{*}$), dust attenuation ($\tau$), and star formation history
(SFH).

####  3.1.1 Prospector

Prospector (Leja et al., 2017; Johnson et al., 2021) is a Baysian-based SED
fitting code with which it is possible to fit with a flexible non-parametric
SFH (e.g., Cid Fernandes et al. (2005); Ocvirk et al. (2006); Tojeiro et al.
(2007); Dye (2008)) by dividing look-back time into several age bins and
assuming a constant SFR in each bin.

With high-quality input data, non-parametric SFH fitting can reproduce the
true SFH more accurately than parametric SFH methods because of fewer
assumptions on the SFH (e.g., Leja et al. (2019)). Lower et al. (2020) have
suggested that a non-parametric SFH can also estimate $M^{*}$ more accurately.
However, we should note that non-parametric SFH estimation still has prior
dependence as Leja et al. (2019) have demonstrated. Prospector can fit
spectroscopic and photometric data together to constrain parameters more
tightly than fitting only photometric data (Johnson et al., 2021).

In this study, we input a flattened 1D spectrum and photometric data of each
galaxy to Prospector and estimate a non-parametric SFH to select RGs. Below we
briefly describe our parameter setting (for a detailed description of
Prospector, see Leja et al. (2017); Johnson et al. (2021)).

####  3.1.2 Prospector settings

Stellar population: Prospector uses the Flexible Stellar Population Synthesis
(FSPS, Conroy et al. (2009); Conroy and Gunn (2010)) to generate model SEDs of
galaxies. We assume a Chabrier (2003) initial mass function. We set stellar
metallicity, total stellar mass, and the velocity dispersion of stellar
components as free parameters and assume that their priors are uniform in
$\log\left(Z^{*}/Z_{\odot}\right)$ between $-2$ and 0.19, in
$\log\left(M^{*}/M_{\odot}\right)$ between $10^{8}$ and $10^{13}$, and uniform
between $10~{}{\rm km/s}$ and $300~{}{\rm km/s}$.

Nebular emission: We use nebular emission templates in Prospector that are
based on CLOUDY (Ferland et al., 1998, 2013) and described in Byler et al.
(2017). We assume that gas-phase metallicity is equal to stellar metallicity.
We set the ionization parameter $\log U$ as a free parameter with a uniform
prior between $-4$ and $-1$.

Dust emission: We also use dust emission templates with fixed parameters. We
do not use AGN dust emission templates because we have excluded AGN-hosting
galaxies in the sample selection (section 2.2).

Dust attenuation: We assume the two-component dust attenuation model by
Charlot and Fall (2000). This model has two types of attenuation. One is the
attenuation only for birth clouds, $\tau_{\rm dust,1}\left(\lambda\right)$,
calculated as:

$\tau_{\rm dust,1}\left(\lambda\right)=\hat{\tau}_{\rm
dust,1}\left(\frac{\lambda}{5500~{}{\rm\AA}}\right)^{-1}.$ (1)

The other is the attenuation for the entire galaxy, $\tau_{\rm
dust,2}\left(\lambda\right)$. We use the equation of Noll et al. (2009):

$\tau_{\rm dust,2}\left(\lambda\right)=\hat{\tau}_{\rm
dust,2}\left(k\left(\lambda\right)+D\left(\lambda\right)\right)\left(\frac{\lambda}{5500~{}{\rm\AA}}\right)^{n},$
(2)

where $k\left(\lambda\right)$ is Calzetti et al. (2000) attenuation curve and
$D\left(\lambda\right)$ is the Lorentzian-like Drude profile parameterizing
the UV bump. We use the $D\left(\lambda\right)$ definition by Kriek & Conroy
(2013). This attenuation model has three parameters; $\hat{\tau}_{\rm
dust,2}$, the ratio of $\hat{\tau}_{\rm dust,1}$ to $\hat{\tau}_{\rm dust,2}$,
and the index $n$. We set all three parameters as free parameters and assume
their priors to be uniform between $0$ and $4$ for $\hat{\tau}_{\rm dust,2}$,
uniform between $1.0$ and $0.4$ for $n$, and a normal distribution with
$\mu=1.0$ and $\sigma=0.3$ truncated at $0.0$ and $2.0$ for $\hat{\tau}_{\rm
dust,1}/\hat{\tau}_{\rm dust,2}$.

SFH: We fix the number of age bins for non-parametric SFH fitting to eight and
calculate the SFR for each bin. Leja et al. (2019) have examined the
dependence of the number of age bins on the reconstructed results and found
that there is no significant dependence if the number of age bins is greater
than 4. We also fix the recent five age bins to [$0$ yr, $10^{7.5}$ yr],
[$10^{7.5}$ yr, $10^{8}$ yr], [$10^{8}$ yr, $10^{8.33}$ yr], [$10^{8.33}$ yr,
$10^{8.67}$ yr], and [$10^{8.67}$ yr, $10^{9}$ yr] in lookback time $t_{l}$,
and set the remaining three age bins by splitting equally $\log t_{l}$ in the
logarithmic interval from $10^{9}$ yr to $0.95t_{H}$.

The typical priors of non-parametric SFHs used in Prospector are the
continuity prior (Leja et al., 2019; Johnson et al., 2021) and the Dirichlet
prior (Leja et al., 2017, 2019). The continuity prior is for the ratio of SFRs
in adjacent age bins, and the Student’s-t distribution is widely used. For the
Dirichlet prior, we assume a Dirichlet distribution (Leja et al., 2017, 2018)
for the fractional sSFR in each age bin. Leja et al. (2019) have examined the
prior dependence of SFHs by applying Prospector on mock data, finding that
both priors can reconstruct the shape of the assumed SFH in mock data.
However, they have also found that the continuity and Dirichlet priors tend to
reconstruct the SFH “dispersively” and “concentratedly”, respectively. In this
study, we use both priors and compare the results. We have eight age bins, and
total stellar mass is also a free parameter. Thus, the number of free
parameters for the SFH is seven for both priors: seven SFR ratios between the
adjacent bins for the continuity prior and seven Dirichlet parameters.

Spectrophotometric calibration: MaNGA data are already flux-calibrated by the
DRP. However, for any galaxy, the area of IFU spectroscopy is different from
photometric apertures. We apply the spectrophotometric calibration model in
Prospector to do aperture correction. This model fits a polynomial calibration
function to scale the input observed spectrum to model spectra so that its
pseudo photometry is consistent with the observed multiband photometry
(Johnson et al. (2021) for detail). We confirm that this calibration model
works well with mock data fitting results.

Emission-line marginalization: Because Prospector assumes all nebular emission
lines to originate from star formation, our SED fitting will fail to
reconstruct the spectrum around an observed emission line if the line is
contaminated by a LIER, an AGN, or shock heating. To avoid this problem, we
marginalize the amplitude of each emission line in each fitting step. We fit
for the velocity dispersion of gas, $\sigma_{\rm gas}$, assuming a uniform
prior between $10~{}{\rm km/s}$ and $300~{}{\rm km/s}$. We also confirm that
introducing line marginalization can improve the model spectrum around a
complex single-line feature consisting of narrow emission and broad absorption
lines.

Likelihood calculation and the outlier model: Finally, we calculate the
likelihood $\mathcal{L}$ as the product of the likelihood for the photometric
data $\mathcal{L}_{\rm phot}$ and the spectrum $\mathcal{L}_{\rm spec}$, as

$\ln\mathcal{L}=\ln\mathcal{L}_{\rm phot}+\left(1-f_{\rm
out}\right)\ln\mathcal{L}_{\rm spec}+f_{\rm out}\ln\mathcal{L}_{\rm out},$ (3)

where $\mathcal{L}_{\rm out}$ is the penalty from the outlier model described
below (equation 5). Here, $\mathcal{L}$ is calculated on the assumption of
Gaussian uncertainties, as

$\ln\mathcal{L}\left(f,m,\sigma\right)=-\sum_{i}^{N}\frac{\left(f_{i}-m_{i}\right)^{2}}{2\sigma_{i}^{2}},$
(4)

where $f$, $\sigma$, and $m$ are the observed flux density, its uncertainty,
and the model flux density, respectively.

We apply the outlier model (Hogg et al., 2010; Johnson et al., 2021) to ignore
outlier data not fitted well by the model spectrum (e.g., cosmic rays, airglow
lines, and complex line features). On the assumption that a certain percentage
$f_{\rm out}$ of data points are outliers, this model adds a penalty term in
the likelihood calculation. The penalty term is calculated by multiplying the
uncertainty by $s_{\rm out}$ for all data points, as

$\ln\mathcal{L}_{\rm
out}\left(f,m,\sigma\right)=-\sum_{i}^{N}\frac{\left(f_{i}-m_{i}\right)^{2}}{2\left(s_{\rm
out}\sigma_{i}\right)^{2}}.$ (5)

With this penalty term, the model becomes less sensitive to the spectral data.
In this study, we applied $s_{\rm out}=50$.

Table 1 summarizes the free parameters and their priors. This parameter
setting is determined by referring to the run time, the distribution of each
parameter obtained from several test runs, the range of each parameter in the
Pipe3D and the MaNGA targeting catalog, and settings in some previous studies
with Prospector (Tacchella et al., 2022).

As an example, figure 3 shows the observed and modeled (best-fit) spectra and
photometry of a galaxy at $z=0.044$ (plateifu:10216-6102). We calculate
$\chi^{2}$ values for spectroscopic data ($\chi^{2}_{\rm spec}$) and
photometric data ($\chi^{2}_{\rm photo}$) separately. We find that the median
of $\chi^{2}/N_{\rm data}$ ($N_{\rm data}$ is the number of data points) for
the whole sample is $\chi^{2}_{\rm spec}/N_{\rm
data,spec}=1.15^{+1.82}_{-0.62}$ and $\chi^{2}_{\rm photo}/N_{\rm
data,photo}=2.44^{+3.85}_{-1.51}$. While the median $\chi^{2}_{\rm
spec}/N_{\rm data,spec}$ value is reasonable ($\simeq 1$), we find that
objects whose spectra around the short-wavelength end are of poor quality have
large $\chi^{2}_{\rm spec}/N_{\rm data,spec}$. Indeed, the $\chi^{2}_{\rm
spec}/N_{\rm data,spec}$ of those objects decreases close to unity if the data
around the short-wavelength end are excluded. On the other hand, the large
median $\chi^{2}_{\rm photo}/N_{\rm data,photo}$ value is primarily because
the WISE data (especially W4) of most objects are not reproduced well. This
problem may be resolved by adding new free parameters to the dust model.

Figure 3:  The top (bottom) panel compares photometric (spectroscopic) data
with the best-fit model for a galaxy at $z=0.044$ (plateifu:10216-6102), with
the lower plot in each panel showing $\chi=\left(f_{\nu,{\rm obs}}-f_{\nu,{\rm
model}}\right)/\sigma$.

The joint posterior distribution plot (corner plot) of six key parameters for
the same example galaxy as figure 3 is shown in figure 4. This galaxy is one
of the rejuvenation galaxies selected with the method described in section
3.2. We find that basic parameters such as $M^{*}$, $Z^{*}$ and
$\hat{\tau}_{\rm dust,2}$ are strongly constrained. The estimated SFH
distribution shows that this galaxy experienced rejuvenation at
$t_{l}=10^{8}~{}{\rm yr}$ with ${\rm SFR}_{\rm 32~{}Myr}\sim 10M_{\odot}/{\rm
yr^{-1}}$, and the resumed star formation ended in the most recent time bin.
As the strong negative correlation between $\hat{\tau}_{\rm dust,2}$ and
$\log\left(Z^{*}/Z_{\odot}\right)$ suggests, dust-metallicity degeneracy can
also be seen in our results. However, the $1\sigma$ range of $\hat{\tau}_{\rm
dust,2}$ and $\log\left(Z^{*}/Z_{\odot}\right)$ is much smaller than the
median value of each parameter, suggesting that our method can constrain these
parameters reasonably well. Besides, $\hat{\tau}_{\rm dust,2}$ and
$\log\left(Z^{*}/Z_{\odot}\right)$ are not strongly correlated with the SFRs
in the recent three bins that are very sensitive to the selection of RGs (see
section 3.2 for our selection method). These results support our RG selection
method.

Figure 4:  Joint posterior distribution plot (corner plot) of $M^{*}$,
$Z^{*}$, $\hat{\tau}_{\rm dust,2}$, and SFRs in the recent three age bins
([$0$ yr, $10^{7.5}$ yr], [$10^{7.5}$ yr, $10^{8}$ yr], and [$10^{8}$ yr,
$10^{8.33}$ yr]) for the example galaxy in figure 3 (plateifu:10216-6102,
$z=0.044$). The diagonal panels show the posterior distribution of each
parameter with black and gray dashed lines indicating the median and the
$1\sigma$ confidence level, and the panels below each diagonal panel show
correlations with the remaining parameters. The upper right panel shows the
estimated SFH of this galaxy; the solid black line indicates the median SFR,
with a gray shaded region corresponding to the $1\sigma$ confidence level of
the fit.

We also check whether Prospector can select true RGs under this setting by
running it on mock data (section 4) and comparing the results for the real
data with existing catalogs (appendix).

Table 1: Free parameters and their priors for Prospector fitting parameter |  Description |  prior
---|---|---
$\log\left(M^{*}/M_{\odot}\right)$ |  total stellar mass |  Uniform: min=8, max=13
$\log\left(Z^{*}/Z_{\odot}\right)$ |  stellar metallicity |  Uniform: min=-2, max=0.19
$\log U$ |  ionization parameter |  Uniform: min=-4, max=-1
$n$ |  power-law index of the dust attenuation curve of the diﬀuse dust |  Uniform: min=-1, max=0.4
$\hat{\tau}_{\rm dust,2}$ |  optical depth of the diffuse dust attenuation |  Uniform: min=0, max=4
$\hat{\tau}_{\rm dust,1}$ |  optical depth of the birth-cloud dust attenuation |  Clipped normal in $\hat{\tau}_{\rm dust,1}/\hat{\tau}_{\rm dust,2}$: min=0, max=2, $\mu=1$, $\sigma=0.3$
$\sigma_{*}$ |  velocity dispersion of the stellar component |  Uniform: min=$10~{}{\rm km/s}$, max=$300~{}{\rm km/s}$
$\sigma_{\rm gas}$ |  velocity dispersion of gas |  Uniform: min=$10~{}{\rm km/s}$, max=$300~{}{\rm km/s}$
$f_{\rm out}$ |  outlier fraction of spectral data points |  Uniform: min=$10^{-5}$, max=$0.5$
SFR ratios |  $\Delta\log{\rm SFR}$ between adjacent age bins (only when using the continuity prior) |  Student’s-t distribution with $\sigma=0.3$ and $\nu=2$
Dirichlet parameters |  weight of SFRs in each age bin (only when using the Dirichlet prior) |  Dirichlet distribution with $\alpha_{D}=1$

###  3.2 Selection method

As introduced in section 1, RGs are galaxies that have restarted star
formation after being quenched. We select RGs using reconstructed SFHs.

For each galaxy, first, we calculate $\mathcal{D}\left(z\right)$, the fraction
of the mass that will be formed by the given redshift $z$ if the star
formation rate is constant with ${\rm SFR}\left(z\right)$, as

$\displaystyle\mathcal{D}\left(z\right)$ $\displaystyle={\rm
sSFR}\left(z\right)\times t_{H}\left(z\right)$ (6) $\displaystyle=\frac{{\rm
SFR}\left(z\right)}{M^{*}\left(z\right)}\times t_{H}\left(z\right),$

where $t_{H}(z)$ is the age of the universe at $z$. Then, we classify the
galaxy in each age bin into an SFG, GV, or QG using
$\mathcal{D}\left(z\right)$ to track the evolution of type with cosmic time.
We determine the threshold $\mathcal{D}$ for classification by fitting the
$\log\mathcal{D}$ distribution with a three-component mixture Gaussian
profile. Figure 5 shows the fitting result of the $\log\mathcal{D}$
distribution in the latest age-bin, i.e., the bin from $t=10^{7.5}\ {\rm yr}$
to the observed time for the continuum prior. We have also confirmed that the
Dirichlet prior gives almost the same results. We fix the SFG/GV and GV/QG
boundaries at $\log\mathcal{D}=-1.2$ and $-1.8$, respectively, based on the
$1\sigma$ values of the peaks corresponding to SFGs and QGs.

Because the location of the SFMS depends on $z$ (e.g., Noeske et al. (2007);
Elbaz et al. (2007); Whitaker et al. (2012); Speagle et al. (2014); Renzini et
al. (2015); Tomczak et al. (2016); Pearson et al. (2018); Popesso et al.
(2023)), the threshold $\mathcal{D}$ value can vary with redshift. For
example, Tacchella et al. (2022) have used $1/3$ and $1/20$ for the
$\mathcal{D}$ thresholds ($-0.48$ and $-1.3$ for $\log\mathcal{D}$ thresholds)
with a $z\sim 0.8$ sample. However, because our RG selection only focuses on
the most recent $10^{8}$ years much shorter than the timescale of the change
of the SFMS location, we can ignore the redshift dependence of the
$\mathcal{D}$ threshold value in this study.

Figure 5:  The colored histogram is the distribution of $\mathcal{D}$ values
in the most recent age-bin obtained with the continuity prior. The solid black
line indicates the result of mixture Gaussian fitting to the histogram, with
three gray dashed lines being the three Gaussian components. The best-fit
values of each Gaussian are written in gray. Red dotted lines indicate the
threshold $\mathcal{D}$ values and the histogram is colored according to the
classification (orange, green, and cyan for QGs, GVs, and SFGs).

Finally, we define an RG as a galaxy that returned to a GV or SFG from a QG
within the recent $1\times 10^{8}~{}{\rm yr}$. This means $t_{\rm rej}$, the
lookback time when a rejuvenation started, is either $10^{8}~{}{\rm yr}$ or
$10^{7.5}~{}{\rm yr}$ in this study. In other words, we define an RG as
follows,

  * •

$t_{\rm rej}=10^{7.5}~{}{\rm yr}$
$\left[0,10^{7.5}~{}{\rm yr}\right]$: GV or SFG
$\left[10^{7.5}~{}{\rm yr},10^{8.0}~{}{\rm yr}\right]$: QG

  * •

$t_{\rm rej}=10^{8.0}~{}{\rm yr}$
$\left[0,10^{7.5}~{}{\rm yr}\right]$: QG, GV or SFG
$\left[10^{7.5}~{}{\rm yr},10^{8.0}~{}{\rm yr}\right]$: GV or SFG
$\left[10^{8.0}~{}{\rm yr},10^{8.33}~{}{\rm yr}\right]$: QG

It is conceivable that galaxies restart star formation before $10^{8}~{}{\rm
yr}$. However, the age-bin setting in our Prospector run is equally spaced in
logarithmic space; thus, returning to a GV or SFG in an age-bin older than
$10^{8}{\rm yr}$ is a different-timescale event from the above definition. For
this reason, we limit $t_{\rm rej}$ to $10^{7.5}$ and $10^{8.0}$ and focus on
the rejuvenation in the recent $\sim 100~{}{\rm Myr}$.

We apply the above selection method to chain data obtained from Prospector
runs and calculate $p_{\rm RG}$, the probability that the galaxy is an RG, as

$p_{\rm RG}=\frac{n_{\rm chain,RG}}{n_{\rm chain}},$ (7)

where $n_{\rm chain}$ is the length of all chains and $n_{\rm chain,RG}$ is
the number of chains classified as an RG. As mentioned above, our RG
definition has two patterns in terms of $t_{\rm rej}$, and we calculate
$p_{\rm RG}$ as the total RG probability; as the sum of the probability of
each $t_{\rm rej}$ pattern to select galaxies more likely to have rejuvenated
regardless of their $t_{\rm rej}$. We define the $t_{\rm rej}$ of a given
galaxy as the timing with a higher $p_{\rm RG}$.

Each galaxy has two $p_{\rm RG}$ values: one with the continuum prior, $p_{\rm
RG,con}$, and the other with the Dirichlet prior, $p_{\rm RG,Dir}$. We select
galaxies with $p_{\rm RG,con}>0.8$ and $p_{\rm RG,Dir}>0.5$ as RGs. We
determine these criteria with the SED fitting results of our parent sample and
the mock data as described in sections 4.3.

Figure 6 shows examples of the specific star-formation history (sSFH) for the
selected RGs. Note that our method depends on the classification into
QG/GV/SFG and hence that a galaxy with a clearly increasing sSFR may have a
relatively low $p_{\rm RG}$ because its sSFR before rejuvenation is just below
the threshold between QGs and GVs (e.g., 8078-6104 in figure 6 column 5 row
3).

Figure 6:  Examples of the sSFHs of selected high-$p_{\rm RG}$ RGs with
$p_{\rm RG}=1$ (rows 1, 2, and 3) and $p_{\rm RG}\neq 1$ (rows 4, 5, and 6,
where $p_{\rm RG}$ is given in each panel). Lines 1 and 2 are for objects with
$\log\left(t_{\rm rej}/{\rm Myr}\right)=7.5$ while lines 3 and 4 for
$\log\left(t_{\rm rej}/{\rm Myr}\right)=8$, where $t_{\rm rej}$ is indicated
by a purple dashed line. In each panel, the solid black line indicates the
median sSFR, with a gray shaded region corresponding to the $1\sigma$
confidence level of the fit. The gray dashed lines in each panel show the
SFG/GV/QG thresholds.

##  4 Verification of our RG selection Method

We run Prospector with mock data to verify that our method works well. Since
SED fitting with a non-parametric SFH has been pointed out to have prior
dependence (Leja et al., 2019), we also discuss prior dependence.

###  4.1 Making mock data

We make mock data at $z=0.1$ with Prospector. Since we also use Prospector in
the SED fitting (see section 3.1 for the fitting process), we cannot evaluate
the validity of the model spectra in this validation. The purpose of using a
mock sample is only to discuss whether our fitting process can reconstruct the
correct SFH and select RGs even with noise.

We construct two mock samples, a single-peak-SFH sample and a double-peak-SFH
sample, based on the delayed-tau SFH, one of the popular parametric SFHs.

Double-peak-SFH mock sample: To examine if our selection method can select RGs
and distinguish RGs from non-rejuvenation galaxies (nRGs) with RG-like SFHs,
we make mock data with double-peak SFHs that are a superposition of two
delayed-tau models as follows,

$\displaystyle{\rm SFH}\left(t\right)=C_{1}$
$\displaystyle\left(t_{1}-t_{l}\right)\exp\left(-\frac{t_{1}-t_{l}}{\tau_{1}}\right)+$
$\displaystyle\begin{cases}C_{2}\left(t_{2}-t_{l}\right)\exp\left(-\frac{t_{2}-t_{l}}{\tau_{2}}\right)&\left(t_{l}>t_{1}\right)\\\
0&\left(t_{l}>t_{2}\right),\end{cases}$ (8)

where $\tau_{1}$ ($\tau_{2}$) is the time scale of the first (second) delayed
tau model, $t_{1}$ ($t_{2}$ [$<t_{1}$]) is the lookback time when the first
(second) delayed tau star formation starts, and $C_{1}$ ($C_{2}$) is a
coefficient of the first (second) delayed tau model. We calculate $C_{1}$ and
$C_{2}$ from $f$, the ratio of the formed mass in the second delayed tau SFH
to the total $M^{*}$. We fix $t_{1}$ as $0.95t_{H}$, and $M^{*}$ as
$10^{11}M_{\odot}$. In equation (8), the first term shows the typical galaxy
evolution scenario of quenching after active star formation. The second term
corresponds to a rejuvenation or pseudo-rejuvenation event. Note that the
first term has a negative value in $t_{l}>t_{1}$; thus, we calculate the SFH
only for $t_{l}>t_{1}$.

To examine various RGs and pseudo RGs, we generate 144 SFHs as a combination
of the following parameters,

  * •

$f$: $10^{-1}$, $10^{-2}$, $10^{-3}$, $10^{-4}$

  * •

$\tau_{1}$: 1 Gyr, 1.5 Gyr, 2 Gyr

  * •

$\left(t_{H}-t_{2}\right)/\tau_{2}$: 1, 10, 50

  * •

$t_{H}-t_{2}$: $10^{7}$ yr, $10^{7.5}$ yr, $10^{8}$ yr, $10^{8.5}$ yr.

The specific star-formation histories (sSFHs) of the mock data are shown in
figure 7. Note that although these mock galaxies have two star-formation
peaks, not all mock galaxies are RGs. For example, if we apply our definition
introduced in section 3.2, some galaxies were SFGs or GVs just before the
start of the second star formation, and these are not classified as RGs.
Besides, we should also note that some galaxies with $f\sim 0.1$ have
extremely high SFRs such as $\log\left({\rm SFR}/M_{\odot}\ {\rm
yr^{-1}}\right)\gtrsim 3$.

Figure 7:  The sSFHs of the double-peak-SFH mock sample. Rows 1, 2, 3, and 4
are for $f=0.0001$, $0.001$, $0.01$, and $0.1$, respectively, while lines 1,
2, and 3 for $\left(t_{H}-t_{2}\right)/\tau_{2}=1$, $10$, and $50$,
respectively. Different colors of lines correspond to different combinations
of $\tau_{1}$ and $t_{H}-t_{2}$. Gray dashed lines show the boundaries of
SFG/GV/QG classification.

We calculate mock SEDs with the above modeled SFHs. We add dust and nebular
emissions to the stellar continuum. We fix $\log\left(Z/Z_{\odot}\right)$ at
$-0.5$ and $\hat{\tau}_{2}$ at $0.3$, which are similar to the estimated
values for our galaxy sample as shown in appendix. We assume a Chabrier (2003)
initial mass function as in the case of real data. Then, we use the mock SEDs
to calculate mock photometric and spectral data by adding noise following the
$S/N$ values shown in the table 2. These $S/N$ are mainly based on the median
$S/N$ for the main sample, but we fix the mock $S/N$ to $20$ for $S/N>20$ and
$5$ for $S/N>5$ WISE photometry to account for possible systematic errors,
such as aperture correction errors. Thus, our mock data are more conservative
than the real data in terms of $S/N$.

Table 2: Median $S/N$ for the main sample and the assumed $S/N$ in making mock observation data. | GALEX | SDSS | WISE | MaNGA
---|---|---|---|---
| FUV | NUV | $u$ | $g$ | $r$ | $i$ | $z$ | W1 | W2 | W3 | W4 | blue | red
Median | 5.1 | 17.6 | 42.9 | 290.1 | 357.6 | 341.3 | 167.6 | 485.7 | 155.7 | 17.8 | 4.5 | 23.0 | 46.6
Mock | 5.1 | 17.6 | 20 | 20 | 20 | 20 | 20 | 5 | 5 | 5 | 4.5 | 20 | 20

Single-peak-SFH mock sample: To discuss the possibility that galaxies without
rejuvenation or a second star-formation event are erroneously selected, we
also generate mock data with single peak SFHs using the single delayed-tau
model as,

$\displaystyle{\rm SFH}\left(t\right)=C_{1}$
$\displaystyle\left(t_{1}-t_{l}\right)\exp\left(-\frac{t_{1}-t_{l}}{\tau_{1}}\right).$
(9)

This time, we also change $t_{1}$: the lookback time of starting star
formation. We generate 25 SFHs as a combination of the following parameters,

  * •

$\tau_{1}$: $2.00\times 10^{8}$ yr, $3.56\times 10^{8}$ yr, $6.32\times
10^{8}$ yr, $1.12\times 10^{9}$ yr, $2.00\times 10^{9}$ yr

  * •

$t_{H}-t_{1}$: $1.25\times 10^{9}$ yr, $2.19\times 10^{9}$ yr, $3.84\times
10^{9}$ yr, $6.74\times 10^{9}$ yr, $1.18\times 10^{10}$ yr,

where $\tau_{1}$ is equally sampled in log space between $2\times 10^{8}$ yr
and $2\times 10^{9}$ yr, and $t_{1}-t_{l}$ is equally sampled in log space
between $0.1t_{H}$ and $0.95t_{H}$.

The sSFHs of the single-peak-SFH mock sample are shown in figure 8. They have
a wide range of current sSFRs over the SFR, GV, and QG regimes and have never
experienced a second active star formation, including rejuvenation.

Figure 8:  The sSFHs of the single-peak-SFH mock sample with
$\tau_{1}=2.00\times 10^{8}$ yr, $3.56\times 10^{8}$ yr, $6.32\times 10^{8}$
yr, $1.12\times 10^{9}$ yr, and $2.00\times 10^{9}$ yr from left to right.
Different colors of lines indicate different $t_{H}-t_{1}$. Gray dashed lines
show the boundaries of SFG/GV/QG classification.

In the calculation of mock SEDs, we make the same assumptions for the
metallicity, dust extinction, and IMF as for the double-peak-SFH mock sample.

###  4.2 Basic parameters

Using the settings described in section 3.1, we run Prospector on the mock
galaxies and estimate their SFH and other parameters.

We compare the estimated parameters with the assumed ones in the mock data in
figure 9. Figures 9 (a), (b), and (c) shows the distribution of $M^{*}$, ${\rm
Z}$, and $\hat{\tau}_{\rm 2,dust}$, respectively. We find that Prospector
reproduces the values of these parameters well. We also find that there are no
large differences between the results for the two priors.

Next, we test how well the assumed SFHs are reconstructed by calculating mean
SFRs in three bins: $t_{l}<100~{}{\rm Myr}$, $100~{}{\rm
Myr}<t_{l}<1000~{}{\rm Myr}$, and $t_{l}>1000~{}{\rm Myr}$ (figure 9 (d), (e),
and (f), respectively). Figure 9 (d) compares ${\rm SFR_{100~{}Myr}}$, the SFR
in $t_{l}<100~{}{\rm Myr}$, showing a strong correlation between the mock and
estimated values. For galaxies with low assumed SFRs ($\log{\rm
SFR_{100~{}Myr}}\lesssim-0.5$), the continuity prior reproduces the mock
values better than the Dirichlet prior. However, both priors overestimate the
SFRs by about 0.3-0.4 dex.

Figure 9 (e) compares ${\rm SFR_{100-1000~{}Myr}}$, the SFR in $100~{}{\rm
Myr}<t_{l}<1000~{}{\rm Myr}$, finding a similarly strong correlation. For
galaxies with medium assumed SFRs ($\log{\rm SFR_{100-1000~{}Myr}}\sim 0$),
the continuity prior often underestimates ${\rm SFR_{100-1000~{}Myr}}$. On the
other hand, for galaxies with low assumed SFRs ($\log{\rm
SFR_{100-1000~{}Myr}}\sim-2$), the Dirichlet prior often overestimates ${\rm
SFR_{100-1000~{}Myr}}$.

Figure 9 (f) compares ${\rm SFR_{before~{}1000~{}Myr}}$, the SFR in
$t_{l}>1000~{}{\rm Myr}$, and finds a weak correlation in both priors. This
study does not focus on past star formations. However, because we define RGs
with sSFHs, underestimation of the old stellar population may cause an
overestimation of the current sSFR and affect our selection. We discuss this
problem in section 4.3.

The distribution of estimated metallicities (figure 9 (b)) shows a small peak
with lower metallicities ($\log\left({\rm Z^{*}}/{\rm
Z_{\odot}}\right)\lesssim-0.8$) in addition to the primary peak around the
assumed value ($\log\left({\rm Z^{*}}/{\rm Z_{\odot}}\right)\sim-0.5$). The
galaxies around the small peak have underestimated masses
($\log\left(M^{*}/M_{\odot}\right)\lesssim 10.7$). The double-peak-SFH mock
sample have 16 galaxies with $\log\left({\rm Z^{*}}/{\rm
Z_{\odot}}\right)<-0.8$, and 15 of them ($\sim 94\%$) have
$t_{H}-t_{2}=10^{7}~{}{\rm yr}$. Besides, ten and six (about 63 and 38%) have
$f=0.1$ and $0.01$, respectively. The ${\rm SFR_{before~{}1000~{}Myr}}$ of
these galaxies are greatly underestimated. It is likely that Prospector cannot
reconstruct old star formation due to an extremely large contribution from
recent star formation, i.e., the young stellar population. Thus, Prospector
may underestimate the stellar mass and metallicity of galaxies with very
recent and active resumed star formation. However, since many $f=0.1$ SFHs
have extremely high SFRs over $10^{3}M_{\odot}~{}{\rm yr}^{-1}$ for low-$z$
galaxies, this underestimation is not a serious problem in applying to real
data. However, this underestimation may need to be considered to search for
RGs in the high-z universe because they have low stellar masses and high SFRs.

Figure 9:  Panels (a), (b), and (c) show the histograms of estimated total
stellar masses, metallicities, and optical depths of the diffuse dust
attenuation, respectively, for the continuity prior in purple and the
Dirichlet prior in green. The assumed values (the true answers) are indicated
by gray dashed lines. Panels (d), (e), and (f) compare estimated SFRs with
mock SFRs within 100 Myr, between 100 Myr and 1000 Myr, and before 1000 Myr,
respectively, for the continuity prior in purple and the Dirichlet prior in
green. Spearman’s correlation coefficient $\rho$ for each prior is shown in
the upper left corner of each panel. The solid black lines show the $y=x$
relation. We perform fitting with $y=x+\Delta$ only for panel (d). Generally,
Our method reconstructs each parameter well, except for SFR in old age (${\rm
SFR_{before~{}1000~{}Myr}}$).

To evaluate how well Prospector can detect rejuvenation events, we calculate
$\mu$, the ratio of the formed mass in the recent 100 Myr to the total stellar
mass. Note that $\mu$ is different from $f$ because we calculate $\mu$ by
dividing the formed mass in the recent 100 Myr by the total $M^{*}$, while
calculating $f$ by dividing the formed mass in the second delayed-tau model by
the total $M^{*}$. Although estimating $f$ from reconstructed SFHs is
challenging, $\mu$ can be easily calculated. Figure 10 compares estimated
$\mu$ with assumed ones. We find that Prospector tends to overestimate $\mu$
about 0.5 dex. This tendency is mainly due to the overestimation of ${\rm
SFR_{100~{}Myr}}$ (about 0.4 dex higher).

Figure 10:  Comparison of $\mu_{\rm mock}$ calculated from mock SFHs with
$\mu_{\rm est}$ obtained from reconstructed SFHs with the continuity prior
(panel [a]) and the Dirichlet prior ([b]). Diamonds are for single-peak SFH
galaxies, and stars for double-peak-SFH galaxies. Yellow and purple colors
mean nRGs and RGs, respectively. Spearman’s correlation coefficient $\rho$ is
given in the upper left corner of each panel. The solid black lines show the
$y=x$ relation. We perform fitting with $y=x+\Delta$, and the value of
$\Delta$ and its 68% error is given in each panel. Our method can select only
$\mu\sim 10^{-3}$ double-peak-SFH objects.

###  4.3 RG selection

The correlation between $p_{\rm RG,con}$ and $p_{\rm RG,Dir}$ for mock and
real galaxies is shown in figure 11 (a). We classify mock galaxies into RGs
and nRGs in the same manner as for real galaxies using reconstructed SFHs.
These mock RGs and nRGs are shown in purple and orange, respectively, in
figure 11 (a). We find that almost all the mock RGs have $p_{\rm
RG,con}>p_{\rm RG,Dir}$.

For real galaxies, a high Spearman’s correlation coefficient of $\rho\simeq
0.79$ is obtained. However, the value drops to $\rho\simeq 0.38$ if limited to
sources with intermediate $p_{\rm RG}$ ($0.1<p_{\rm RG,con}<0.9$), suggesting
that the high coefficient value for all galaxies is caused by very RG-like
($0.9<p_{\rm RG}$) and very non-RG-like ($p_{\rm RG}<0.1$) galaxies. We note
that both priors are relatively consistent with each other in that they can
correctly select both very RG-like and very non-RG-like objects.

Despite the high correlation coefficient, real galaxies (and mock galaxies)
are not evenly distributed around the equality line; objects with medium
$p_{\rm RG,con}$ values tend to be lower than the line. This distribution may
be explained as follows. As seen in section 4.2, the continuity prior tends to
give lower ${\rm SFR}_{\rm 100-1000~{}Myr}$ and higher ${\rm SFR}_{\rm
before~{}100~{}Myr}$ than the Dirichlet prior, especially for low SFR objects.
This means that the continuity prior tends to reconstruct more V-shaped, or
more RG-like, SFHs than the Dirichlet prior. Thus, $p_{\rm RG,con}$ is likely
to be higher than $p_{\rm RG,Dir}$ especially for low to medium $p_{\rm
RG,con}$ objects, which are often located around the boundary of QGs and GVs.

Considering these trends in the mock and main parent samples, we set up three
regions (A), (B), and (C) in the $p_{\rm RG,Dir}$ \- $p_{\rm RG,con}$ plane as
candidate regions for RG selection and discuss which region selects RGs best.
Regions (A), (B), and (C) are defined as,

  * •

(A) :$p_{\rm RG,Dir}>0.5$ and $p_{\rm RG,con}>0.8$,

  * •

(B) :$p_{\rm RG,Dir}<0.4$ and $p_{\rm RG,con}>0.6$,

  * •

(C) :$p_{\rm RG,Dir}<0.2$ and $p_{\rm RG,con}<0.2$.

The SFHs and sSFHs of mock RGs and mock nRGs in each region are shown in
figure 11 (b). The parameters of the single-peak-SFH mock sample in each
region are also shown in table 3.

Figure 11:  (a) Comparison of $p_{\rm RG,con}$ and $p_{\rm RG,Dir}$. Diamonds,
stars, and dots indicate single-peak-SFH mock galaxies, double-peak-SFH mock
galaxies, and real data. Yellow and purple colors mean nRGs and RGs,
respectively. The dashed black lines show the $y=x$ relation. Regions (A),
(B), and (C) are shaded in purple, gray, and yellow, respectively. Region (A)
can select $f=0.001$ RGs with high completeness. (b) The mock SFHs and sSFHs
of mock galaxies in each region. RGs and nRGs are separately plotted. In the
sSFH plots, the gray dashed lines indicate the SFG/GV/QG boundaries.

In region (A), we have 12 mock galaxies, and six of them are classified as RGs
with mock SFHs. All six RGs have $f=0.001$ and $t_{H}-t_{2}=10^{7}~{}{\rm
yr},10^{7.5}~{}{\rm yr}$. Similarly, all the mock nRGs in region (A) also have
$f=0.001$ and $t_{H}-t_{2}=10^{7}~{}{\rm yr}$. From these results, galaxies in
region (A) are likely to have experienced a very recent second SF with $f\sim
0.001$.

Then, in region (B), we have 18 mock galaxies, and five of them are classified
as RGs with mock SFHs. As far as we can judge from the mock data, the
contamination rate in region (B) is 72%, which is higher than in region (A)
(50%). Region (B) has more mock galaxies with $f=0.01$ and
$t_{H}-t_{2}=10^{7.5}~{}{\rm yr}$ than region (A).

Finally, in region (C), we have 106 mock galaxies, and 21 of them are
classified as RGs with mock SFHs. The contamination rate is 80%, the highest
in the three regions. Most single-peaked-SFH mock galaxies ($22/25$) are
distributed in region (C), meaning that our method can correctly diagnose
single-peaked-SFH galaxies as nRGs. We find that mock RGs with
$f=0.0001,0.01,$ and $0.1$ are concentrated in region (C) while there are no
mock RGs with $f=0.001$ there. Besides, most of the mock RGs with
$t_{H}-t_{2}\gtrsim 10^{8}~{}{\rm yr}$ are in region (C). Although our method
focuses on rejuvenation within $10^{8}~{}{\rm yr}$, rejuvenation events
starting from $t_{l}=10^{8}~{}{\rm yr}$ may not be detected in our method.

Based on these results, we choose the region for RG selection. If we use
region (A) to select RGs, the contamination rate of RG selection in mock data
is 50% ($6/12$). However, the contaminants, i.e., nRGs in region (A), also
have a recent second star formation with $f=0.001$. Besides, region (A) does
not include single-peaked-SFH mock galaxies. Therefore, region (A) can select
RGs and RG-like secondary-star-formation galaxies with $f\sim 0.001$. Since
the parameter distribution of mock galaxies should not be the same as the
distribution of real galaxies, we cannot calculate the exact contamination
rate. We can only affirm the possibility of selecting such contaminants.

Next, the completeness of RG selection using region (A) is found to be 17%
($6/36$). However, similar to the above discussion, we cannot know the exact
completeness with our mock data. For example, as shown in figure 11(b), the
mock SFHs in regions (B) or (C) have extremely high SFRs like ${\rm
SFR}\gtrsim 10^{2}M_{\odot}~{}{\rm yr^{-1}}$. In the low-$z$ Universe,
galaxies with such extremely high SFRs are extremely rare, and the settings of
mock data may not be realistic. If we focus only on galaxies with $f=0.001$,
the contamination is $67\%$ ($6/9$). The results with mock data suggest that
using region (A), we can detect weak rejuvenation events with high
completeness, although maybe missing galaxies with strong rejuvenation. From
these results, we decide to use region (A).

As for contamination, region (A) may select galaxies with a second star
formation whose progenitors are not QGs but GVs or SFGs. As for completeness,
region (A) can miss galaxies that have experienced strong rejuvenation with
$f\gtrsim 0.01$ or weak rejuvenation with $f\sim 0.0001$. However, region (A)
can select RGs with $f\sim 0.001$ with high completeness.

We use the results with the continuity prior in the remaining part of this
paper. This is mainly because the continuity prior has a stronger correlation
in ${\rm SFR_{100~{}Myr}}$ (figure 9 (d)) and in $\mu$ (figure 10). Besides,
$p_{\rm RG,con}$ tends to be higher for $f=0.001$ mock RGs than $p_{\rm
RG,Dir}$, suggesting that the continuity prior is better in reproducing
$f=0.001$ rejuvenation events. Note that our selection method uses both
$p_{\rm RG,con}$ and $p_{\rm RG,Dir}$ as introduced in section 3.2 and
discussed above.

The low completeness of our selection method comes from difficulties in
estimating non-parametric SFHs. For large-$f$ galaxies ($f\gtrsim 0.01$), the
light from the young stellar population overwhelms that from the old stellar
population, i.e., the outshining problem (c.f., Maraston et al. (2010); Sorba
& Sawicki (2018); Giménez-Arteaga et al. (2023); Narayanan et al. (2023)). In
contrast, for small-$f$ galaxies ($f\sim 0.0001$), the light from the young
stellar population is buried by that from the old stellar population. It will
be important to reduce this uncertainty as much as possible by improving the
method, e.g., improving the SED fitting technique and searching for indices to
find RGs easily. The latter will become especially important in the future as
data are increasing rapidly. Because the completeness of selecting $f=0.001$
mock RGs is high, we focus on $f\sim 0.001$ RGs in this study. We will discuss
the effects of this incompleteness in section 6.

Again, note that because the mock data are generated with Prospector, the
comparison in this section cannot demonstrate the validity of model spectra.
To further demonstrate the validity of our SED fitting results with
Prospector, we also checked that the basic results, such as $M^{*}$ and SFR,
estimated with Prospector are strongly correlated with an existing SED fitting
result, the Pipe3D DR17 catalog (Sánchez et al., 2022). The detail of this
comparison is introduced in appendix.

Table 3: The number of mock galaxies with each mock parameter in each region. |  | $f$ | $\tau_{1}$ | $\left(t_{H}-t_{2}\right)/\tau_{2}$ | $t_{H}-t_{2}$
---|---|---|---|---|---
|  | 0.0001 | 0.001 | 0.01 | 0.1 | 1 Gyr | 1.5 Gyr | 2 Gyr | 1 | 10 | 50 |  $10^{7}$ yr |  $10^{7.5}$ yr |  $10^{8}$ yr |  $10^{8.5}$ yr
All | All | 36 | 36 | 36 | 36 | 48 | 48 | 48 | 48 | 48 | 48 | 36 | 36 | 36 | 36
A | 0 | 12 | 0 | 0 | 6 | 3 | 3 | 4 | 4 | 4 | 9 | 3 | 0 | 0
(0%) | (33%) | (0%) | (0%) | (13%) | (6.3%) | (6.3%) | (8.3%) | (8.3%) | (8.3%) | (25%) | (8.3%) | (0%) | (0%)
B | 0 | 2 | 10 | 6 | 5 | 8 | 5 | 2 | 9 | 7 | 3 | 14 | 1 | 0
(0%) | (5.6%) | (28%) | (17%) | (10%) | (17%) | (10%) | (4.2%) | (19%) | (15%) | (8.3%) | (39%) | (2.8%) | (0%)
C | 36 | 17 | 23 | 30 | 33 | 33 | 39 | 38 | 32 | 35 | 23 | 15 | 32 | 35
(100%) | (47%) | (64%) | (83%) | (69%) | (69%) | (81%) | (79%) | (67%) | (73%) | (64%) | (42%) | (89%) | (97%)
RG | All | 9 | 9 | 9 | 9 | 36 | 0 | 0 | 12 | 12 | 12 | 12 | 12 | 12 | 0
A | 0 | 6 | 0 | 0 | 6 | 0 | 0 | 2 | 2 | 2 | 3 | 3 | 0 | 0
(0%) | (67%) | (0%) | (0%) | (17%) | - | - | (17%) | (17%) | (17%) | (25%) | (25%) | (0%) | -
B | 0 | 1 | 2 | 2 | 5 | 0 | 0 | 0 | 3 | 2 | 0 | 4 | 1 | 0
(0%) | (13%) | (22%) | (22%) | (14%) | - | - | (0%) | (25%) | (17%) | (0%) | (33%) | (8.3%) | -
C | 9 | 0 | 5 | 7 | 21 | 0 | 0 | 8 | 6 | 7 | 9 | 4 | 8 | 0
(100%) | (0%) | (56%) | (78%) | (58%) | - | - | (67%) | (50%) | (58%) | (75%) | (33%) | (67%) | -

##  5 Results

###  5.1 RG fraction

Our method described in section 3 selects 1071 RGs out of the $N=8857$ parent
sample; thus, the RG fraction calculated by simply dividing the number of RGs
by the number of all galaxies is about 12%. Considering the volume weight $w$
provided in the MaNGA targeting catalog (Wake et al., 2017) for the Primary+
sample and the secondary sample (ESRWEIGHT), we calculate the intrinsic RG
fraction, $f_{\rm RG}$, as,

$f_{\rm RG}=\frac{\sum_{i\in{\rm RG}}w_{i}}{\sum_{i\in{\rm All}}w_{i}},$ (10)

and obtain $f_{\rm RG}\simeq 8.1\%$. Furthermore, when restricted to galaxies
with $f\sim 0.001$ and considering the completeness discussed in section 4.3,
the true RG fraction is $f_{\rm RG}\sim 12\%$. Note, however, that due to the
distribution differences between mock galaxies and real galaxies, it is
difficult to estimate the true RG fraction using our results and completeness
from the mock data, especially when considering RGs with all $f$ values.

As shown in section 2.2, we have excluded 520 objects as AGN or merging
galaxies from the parent sample. If these galaxies have higher $f_{\rm RG}$
than our parent sample, the $f_{\rm RG}$, including these galaxies, will also
be higher than the above value.

###  5.2 Images and spectra

Figure 12 compares images of the selected RGs constructed from MaNGA datacubes
with HSC-SSP PDR3 (Aihara et al. (2022)) $gri$ composite images. HSC images
are deeper and cover a larger area than MaNGA datacubes. We find that although
RGs were QGs before rejuvenation, many RGs have a disk structure. The
morphologies of RGs will be analyzed in detail in section 5.4.

Figure 12:  Cutout images of 30 RGs included in the HSC-SSP PDR3, with their
plateifu given at the top. The left side of each panel is a $gri$ image
constructed from the MaNGA datacube and the right side is a $1^{\prime}\times
1^{\prime}$ HSC $gri$ image. Each image is scaled with the arcsinh method
(Lupton et al., 2004). The yellow line in each panel indicates ten arcsecs.
For each object, $p_{\rm RG}$, $t_{\rm rej}$, and the current state
(QG/GV/SFG) are shown below the plateifu, from left to right, with the
letters’ color corresponding to the current state: orange for QG, yellow-green
for GV, and cyan for SFG.

The upper panels of figure 13 show stacked normalized spectra of current QGs,
GVs, and SFGs, separated by whether they are RGs or not; i.e., comparison of
rejuvenated QGs (RQGs) with non-rejuvenation QGs (nRQGs), rejuvenating GVs
(RGVs) with non-rejuvenation GVs (nRGVs), and rejuvenating SFGs (RSFGs) with
non-rejuvenation SFGs (nRSFGs). Then, to compare the spectra of RGs with those
of non-rejuvenation galaxies (nRGs), we calculate a normalized residual
spectrum, $\Delta F_{\lambda}\left(\lambda\right)$, as,

$\Delta F_{\lambda}\left(\lambda\right)=\frac{F_{\lambda,{\rm
RG}}\left(\lambda\right)-F_{\lambda,{\rm
nRG}}\left(\lambda\right)}{F_{\lambda,{\rm nRG}}\left(\lambda\right)},$ (11)

as shown in the lower panels of figure 13. Additionally, we examine the
distribution of galaxies on the well-studied ${\rm Dn4000}$ versus the ${\rm
H\delta}$ equivalent widths (EW) plane in the lowest panels. We calculate
${\rm Dn4000}$ as follows,

${\rm Dn4000}=\frac{\int_{4000~{}{\rm\AA}}^{4100~{}{\rm\AA}}f_{\nu}~{}{\rm
d}\lambda}{\int_{3850~{}{\rm\AA}}^{3950~{}{\rm\AA}}f_{\nu}~{}{\rm d}\lambda}.$
(12)

Firstly, it is found that the spectral difference between RGs and non-
rejuvenation galaxies is less than the $1\sigma$ range, consistent with the
previous results from observations by Chauke et al. (2019) and simulated
spectra by Zhang et al. (2022). As also mentioned in Zhang et al. (2022), this
result suggests that it is difficult to select RGs without an SED fitting
method.

Focusing on QGs, panel (a), $\Delta F_{\lambda}$ is close to 0 for wavelengths
longer than $4000~{}{\rm\AA}$, but increases to around 0.1 for shorter
wavelengths. This indicates that RQGs have a slightly weaker $4000~{}{\rm\AA}$
break than nRQGs. Similarly, in the ${\rm Dn4000}$ versus ${\rm H\delta}$ EW
plane, RQGs have similar EWs to nRQGs but smaller ${\rm Dn4000}$. This is
probably due to the presence of a young stellar population, with a smaller
${\rm Dn4000}$ value, formed in the rejuvenation event. Also seen is a strong
excess in the ${\rm H\alpha}$ emission line, which also indicates recent star
formation.

In the case of GVs, panel (b), $\Delta F_{\lambda}$ is found to decrease with
decreasing wavelength. This is because RGVs retain a larger population of old
stars formed before quenching than nRGVs. RGVs have a slightly stronger
$4000~{}{\rm\AA}$ break than nRGVs, but not as pronounced as in QGs. This
tendency is also seen in the ${\rm Dn4000}$ versus ${\rm H\delta}$ EW plane,
i.e., there is no significant difference in the distributions of RGVs and
nRGVs. Looking at $\Delta F_{\lambda}$ for each line, H recombination lines
(Balmer lines) have larger $\Delta F_{\lambda}$ than the continuum. This is
likely because RGVs have stronger Balmer emission lines due to active recent
star formation. In contrast, neutral metal lines like Nai $\lambda\lambda
5885,5889$ (D1, D2 lines) and Mgi$\lambda 5175$ have lower $\Delta
F_{\lambda}$ than the continuum, probably because of an aging stellar
population formed before quenching.

In the case of SFGs, panel (c), we can see a stronger decreasing trend of
$\Delta F_{\lambda}$ than in GVs. This is because RSFGs have a larger amount
of old stellar population formed before quenching than nRSFGs. Similar to GVs,
RSFGs have stronger Balmer emission lines and absorption lines of neutral
metals than nRSFGs, suggesting the presence of recent star formation and an
older stellar population before quenching.

Figure 13:  (Row 1) Stacked normalized spectra
$F_{\lambda}\left(\lambda\right)$ of QGs (panel [a]), GVs ([b]), and SFGs
([c]). The purple and black lines show the spectra of those selected and not
selected as RGs, respectively, with shaded areas corresponding to the
$1\sigma$ confidence intervals. The wavelengths of several strong Balmer lines
and neutral metal lines are indicated with cyan and gray dashed lines,
respectively. The Balmer lines plotted are from H$\alpha$ to H$\zeta$, from
left to right. The metal lines plotted are Nai $\lambda\lambda 5889,5885$ (D2
and D1 line), Mgi $\lambda 5175$, Cai $\lambda 4307$ (G line), Caii$\lambda
3968$ (H line), and Caii$\lambda 3933$ (K line). (Row 2) Differences in the
stacked normalized spectra $\Delta
F_{\lambda}\left(\lambda\right)=\left\\{F_{\lambda,{\rm
RG}}\left(\lambda\right)-F_{\lambda,{\rm
not~{}RG}}\left(\lambda\right)\right\\}/F_{\lambda,{\rm
not~{}RG}}\left(\lambda\right)$ of QGs, GVs, and SFGs (the order is the same
with row 1), where the shaded areas correspond to the $1\sigma$ confidence
intervals. Similar to row 1, the wavelengths of Balmer lines and neutral metal
lines are indicated. (Row 3) Magnified views around $4000~{}{\rm\AA}$. (Row 4)
Relationship between H$\delta$ EW and Dn4000. Gray and purple dots indicate
nRGs and RGs, respectively.

###  5.3 $M^{*}$ and SFR distributions

Figure 14 (a) shows that the distribution of RGs in the SFR-$M^{*}$ plane is
biased toward high $M^{*}$ and SFR values, with their median SFR being higher
than those of nRSFGs and nRGVs. We also find that the median $M^{*}$ of RGs is
comparable to that of nRGVs but lower than that of nRQGs. Note that although
RGs have a similar distribution to nRGVs in the 1D distributions of SFR and
$M^{*}$, the distribution of RGs in the SFR-$M^{*}$ plane is closer to that of
SFGs than of nRGVs, as indicated by the black contours in figure 14 (a).

Figure 14 (a) also gives an insight into the path in the SFR-$M^{*}$ plane
during rejuvenation events. Since we focus on recent rejuvenation events
within $10^{8}$ yr (by our definition), even if the sSFR increases to about
$10^{-9}~{}{\rm yr^{-1}}$ due to rejuvenation, the ratio of the formed mass
during the rejuvenation to the total stellar mass is only $10^{-1}$. This
means that galaxies move upward almost vertically in the SFR-$M^{*}$ plane by
rejuvenation. This result is consistent with Chauke et al. (2019).

Figure 14:  (a) The distribution of RGs (purple diamonds) in the SFR-$M^{*}$
plane, with calculated trajectories during the rejuvenation overplotted as
solid purple lines (only for randomly selected 100 RGs with $t_{\rm
rej}=10^{7.5}~{}{\rm yr}$). The dots show all galaxies in our sample, and
black solid contours indicate the 1, 2, and 3 $\sigma$ levels of their
distribution. The two black dashed lines correspond to the boundaries of
SFG/GV/QG classification at the median redshift of all galaxies. The purple,
orange, green, and cyan histograms in the subplots are for RGs, not
rejuvenation QGs, not rejuvenation GVs, and not rejuvenation SFGs,
respectively. In each subplot, the median of each distribution is indicated by
a short bar of the same color as the histogram. (b) The $f_{\rm RG}$
distribution in the SFR-$M^{*}$ plane. The two black dashed lines correspond
to the boundaries of SFG/GV/QG classification at the median redshift of the
parent sample. In both panels, we use SFRs in the recent 32 Myr (the most
recent age bin).

Figure 14 (b) shows the distribution of RG fractions in the SFR versus $M^{*}$
plane. The RG fraction is found to be higher for GVs and massive
($M^{*}>10^{11}M_{\odot}$) SFGs. The RG fraction calculated with equation (10)
for only galaxies with $M^{*}>10^{11}M_{\odot}$ is about 15%, indicating that
RGs are more abundant on the high mass side. These results are consistent with
previous studies, which compared rejuvenated galaxies with QGs and found that
RGs have a smaller mass than QGs (e.g., Chauke et al. (2019); Tacchella et al.
(2022)).

###  5.4 Morphologies

In order to compare morphologies between RGs and other galaxies, we show in
figure 15 the distributions of the $r$-band Sérsic index $n_{r}$, the $r$-band
effective radius $R_{e,r}$, the bulge-to-total flux ratio $B/T$, and the ratio
of rotational velocity $V$ to velocity dispersion $\sigma$, $V/\sigma$. For
$n_{r}$, $R_{e,r}$, and $B/T$, we use the PyMorph catalog (Domínguez Sánchez
et al., 2022), which summarizes photometric morphological parameters obtained
from $g$, $r$, and $i$ images of MaNGA DR17 galaxies. This PyMorph catalog has
two types of results, one fit by a Sérsic profile and the other by a
Sérsic+exponential profile, with FLAG_FIT parameter indicating which result is
reliable: FLAG_FIT=0 indicates that both fits are reliable; FLAG_FIT=1
indicates that only the Sérsic fit is successful while the Sérsic+exponential
profile fit may not be reliable; FLAG_FIT=2 indicates that only the
Sérsic+exponential fit is successful while the Sérsic profile fit may not be
reliable; FLAG_FIT=3 indicates both fits may be unreliable. Using this
parameter, we construct a different sample depending on the morphological
parameter to be examined. We use only galaxies with FLAG_FIT=1 when discussing
$n_{r}$, only galaxies with FLAG_FIT=0 or 2 when discussing $B/T$, and only
galaxies with FLAG_FIT=0, 1, or 2 when discussing $r_{e,r}$. For the kinematic
parameters, $V$ and $\sigma$, we use the DAPall catalog (Westfall et al.
(2019); Belfiore et al. (2019)).

First, we focus on the photometric morphology parameters. From figure 15 (a)
and (b), we find that RGs have a very different morphology from nRQGs (the
$p$-value of Kolmogorov-Smirnov test is $p\ll 0.05$); RGs show more a disk-
like morphology (lower $n_{r}$ and $B/T$) than nRQGs and a more elliptical-
like morphology (higher $n_{r}$) than nRSFGs. From figure 15 (c), we find that
RGs tend to have a larger radius than nRSFGs and a smaller radius than nRQGs.
Similar to the $n_{r}$ and $B/T$ distributions, RGs have almost the same
$R_{e,r}$ distribution as nRGVs. Based on the above results, we conclude that
the RGs have a photometric parameter distribution between nRQGs and nRSFGs,
similar to nRGVs.

We can also see similar trends in the kinematic morphology parameters. Figure
15 (d) shows that RGs have an intermediate distribution of $V/\sigma$ between
nRSFGs and nRQGs. Particularly, RGs have a more rotation-supported morphology
than nRQGs, even though the mass distribution of nRQGs is not significantly
different from that of the RGs. RGs also show a similar $V/\sigma$
distribution to nRGVs, but nRGVs have a slightly more rotational-supported
than RGs.

Figure 15:  The relationship between stellar mass and four indices of
morphology: Sérsic index in the $r$ band $n_{r}$ (panel [a]), bulge to total
flux ratio $B/T$ ([b]), effective radius in the $r$ band $r_{e,r}$ ([c]), the
ratio of the rotation velocity $V$ to the velocity dispersion $\sigma$ ([d]).
Plotted in panels (a), (b), and (c) are only galaxies with FLAG_FIT=1,
FLAG_FIT=0 or 2, and FLAG_FIT=0, 1, or 2. In each panel, purple, orange,
green, and cyan colors indicate RGs, nRQGs, nRGVs, and nRSFGs, respectively,
where dots are individual objects, contours are the 0.5, 1, and 2$\sigma$
levels of the distribution, and histograms are the distributions of individual
parameters, with medians indicated by short bars. Shown in the top right of
each $y$-axis histogram are the $p$-values of the Kolmogorov–Smirnov test
between the RG sample and the three non-rejuvenation galaxy samples (nQGs,
nGVs, nSFGs).

###  5.5 Environment

To examine the environmental dependence of RGs, we cross-match our parent
sample with the group catalog of Tempel et al. (2012) that applied a modified
friends-of-friends method to the SDSS DR8, obtaining 7585 galaxies with
environmental information. We consider those with the richness $n_{\rm
rich}=1$ to be in the field or void environment, while those with $n_{\rm
rich}>4$ to be in the cluster environment.

Figure 16 compares among galaxy types the distribution of normalized
environmental densities within $r=1h^{-1}~{}{\rm Mpc}$ that are calculated by
Tempel et al. (2012) from galaxy luminosity density fields in the same manner
as Liivamägi et al. (2012). We find that the density distribution of RGs is
intermediate between those of QGs and SFGs, being similar to GVs. This trend
is similar to the morphology trend found in section 5.4 and consistent with
previous studies such as Schawinski et al. (2007); Thomas et al. (2010);
Chauke et al. (2019) in that RGs tend to be in lower-density environments than
QGs.

Figure 16:  Same as Figure 15, but for stellar mass versus the scaled
environmental density within $1h^{-1}~{}{\rm Mpc}$ from each galaxy.
Individual galaxies are not plotted.

Next, we draw a phase space diagram (Bertschinger (1985)) for galaxies in the
cluster environment ($n_{\rm rich}>4$), with the clustocentric radius $r_{\rm
sep}$ normalized by virial radius $r_{\rm vir}$ being the horizontal axis and
the clustocentric velocity in the line-of-sight direction $v_{\rm LOS}$
normalized the velocity dispersion in the cluster $\sigma_{\rm cl}$ being the
vertical axis, using the group catalog by Tempel et al. (2012). Following Rhee
et al. (2017), we classify galaxies into five regions, A, B, C, D, and E,
according to their positions on the phase-space diagram as shown in figure 17.
Rhee et al. (2017) have reported that region A is dominated by interlopers and
galaxies that have fallen into the cluster recently. In contrast, region E is
dominated by galaxies that have been in the cluster for a long time.

Table 4 shows $f_{\rm RG}$ for the five regions and the field. We find that
$f_{\rm RG}$ is higher in region A and the field than in regions B-D and E. We
also calculate $f_{\rm RG/nRQG}$, the ratio of the number of RGs to the number
of nRQGs. As RGs were QGs before rejuvenation, we consider $f_{\rm RG/nRQG}$
as an indicator of the fraction of QGs that caused rejuvenation, i.e., how
likely rejuvenation occurs in QGs. $f_{\rm RG/nRQG}$ is found to be highest in
the field and decreases toward the cluster center in the order of regions A,
B-D, and E. This suggests that rejuvenation is more likely to occur in
galaxies in the field and galaxies that have just entered a cluster than in
galaxies that have been in a cluster for a long time.

Table 4: The number of each type of galaxies for each environment. The error in the fraction only includes the Poisson error.  |  | Number | Fraction
---|---|---|---
|  | All | RG | nRSFG | nRGV | nRQG |  $f_{\rm RG}$: RG/All |  $f_{\rm RG/nRQG}$; RG/nRQG
All Environment | 6743 | 731 | 2831 | 376 | 2805 | $10.8\pm 0.4\%$ | $26.1\pm 1.1\%$
Field ($N_{\rm rich}=1$) | 2063 | 221 | 1179 | 122 | 541 | $10.7\pm 0.8\%$ | $40.9\pm 3.3\%$
Cluster ($N_{\rm rich}>4$) | A | 1056 | 129 | 362 | 65 | 500 | $12.2\pm 1.1\%$ | $25.8\pm 2.5\%$
B-D | 1294 | 102 | 326 | 69 | 797 | $7.9\pm 0.8\%$ | $12.8\pm 1.3\%$
E | 272 | 20 | 38 | 13 | 201 | $7.4\pm 1.7\%$ | $10.0\pm 2.3\%$
Figure 17:  The phase-space diagram of our galaxies in the cluster
environment. Following Rhee et al. (2017), we divide the cluster environment
into five regions (A to E) whose boundaries are shown by solid gray lines.
Dots indicate all galaxies, color-coded according to the region they reside
in. Purple stars are RGs, and purple contours are the 1 and $2\sigma$ levels
of their distribution.

##  6 Discussion

###  6.1 RGs are important players in galaxy evolution

One of the key questions in the study of rejuvenation is the importance of RGs
in galaxy evolution, especially how much of the total stellar mass was formed
by rejuvenation and whether or not rejuvenation has a significant contribution
to cosmic star formation. Chauke et al. (2019) have reported that about 10% of
the total stellar mass of rejuvenated galaxies is formed in the rejuvenation
event and that RGs contribute only about 0.3% of the cosmic SFR density
(CSFRD) over $z\sim 0.7-1.5$. Tacchella et al. (2022) have reported that, at
most, only about 10% of the total mass is formed in rejuvenation events.

To address this question, first, we calculate $C_{j}$, the contribution of
each current galaxy population $j$ ($j=$nRSFG, nRGV, nRQG, and RG) to the
CSFRD in the last 100 Myr, the timescale of rejuvenation in our definition.
Here, $C_{j}$ is defined as the ratio of the total $M^{*}$ formed in $j$-type
galaxies in the last 100 Myr to the total $M^{*}$ formed in all galaxies in
the last 100 Myr,

$C_{j}=\frac{\sum_{i\in j}w_{i}M^{*}_{i,{\rm 100~{}Myr}}}{\sum_{i\in{\rm
All}}w_{i}M^{*}_{i,{\rm 100~{}Myr}}},$

where $M^{*}_{i,{\rm 100~{}Myr}}$ is the formed stellar mass in a galaxy $i$
in the last 100 Myr and $w_{i}$ is the volume weight of this galaxy given in
the MaNGA catalog. The mass dependence of $C_{j}$ is shown in the upper figure
18. The sum of $C_{\rm RG}$ for all $M^{*}$ bins is about $17\%$, meaning that
RGs contribute about $17\%$ to the CSFRD in the recent 100 Myr, about two
orders of magnitude larger than the value obtained by Chauke et al. (2019) for
rejuvenated galaxies at $z\sim 0.8$. If the detail of rejuvenation events is
similar at $z\sim 1$ and $z\sim 0$, this result implies that rejuvenation
events have higher importance at lower $z$. However, this difference may be
explained by the difference in the definition of rejuvenation and the
difference in the parent sample (they search for RGs only from QGs). It may
also be attributed to the $z$ dependence of the CSFRD because the CSFRD at
$z\sim 1$ is about one order of magnitude higher than at $z\sim 0$ (e.g.,
Madau & Dickinson (2014)). Note also that we cannot rule out the possibility
that $C_{j}$ is overestimated or underestimated in this study due to the
contamination or completeness discussed in section 4.3.

Next, we calculate $\mu$, the ratio of the formed stellar mass in the last 100
Myr, $M^{*}_{\rm 100~{}Myr}$, to the total stellar mass $M^{*}$ for each RG
(same as the discussion in section 4.2). The lower panel of figure 18 shows
$\mu$ as a function of $M^{*}$. We find that rejuvenation events increase
$M^{*}$ only by a few percent at most, with $\sim 0.1\%$ on average. This
result roughly agrees with Akhshik et al. (2021) reporting that $0.5\%$ of
$M^{*}$ is formed in 100-Myr-timescale rejuvenation events at $z=1.883$. This
agreement may suggest that the detail of rejuvenation events does not change
with redshift.

Note that the values of $\mu$ obtained here are consistent with the result for
mock data presented in section 4.2 that our method can select galaxies with
$\mu\sim 10^{-3}$ with high completeness. This consistency suggests that our
RG sample is not significantly contaminated. On the other hand, we cannot rule
out the possibility that $\mu\sim 10^{-3}$ is obtained just because our method
fails to select RGs with other $\mu$ values. To do a more solid statistical
discussion, it will be important to evaluate contamination and completeness
using a larger mock sample covering wider parameter ranges or maybe
cosmological simulation data and to develop a more robust selection method.

The very small increase in $M^{*}$ by rejuvenation obtained here is primarily
because the rejuvenation timescale, 100-Myr, is too short to increase $M^{*}$
significantly. However, if galaxies experienced rejuvenation events many times
in the past, rejuvenation may be an important channel of galaxy mass growth.
Although we cannot directly detect such past rejuvenation events for our
sample because of insufficient time resolution of the reconstructed SFHs, we
discuss the possibility of multiple rejuvenation events in section 6.2.

Figure 18:  Upper panel. Distribution of $C_{j}$ as a function of $M^{*}$ for
four types of galaxies coded by color. The horizontal dotted line indicates
$C_{\rm RG}$, the contribution of RGs averaged over all masses. Lower panel.
$\mu$ (purple solid line) and its $1\sigma$ confidence level (purple shaded
area) as a function of $M^{*}$.

###  6.2 Multiple rejuvenation scenario

Whether rejuvenation can occur more than once in a galaxy’s lifetime is still
unclear. Nelson et al. (2018) have found that $\sim 10\%$ and $\sim 1\%$ of
massive passive galaxies in the TNG100 simulation have experienced one and
more than one rejuvenation event, respectively. In this subsection, we roughly
estimate $N_{\rm rej}$, the number of times galaxies have experienced
rejuvenation, using $f_{\rm RG}$ obtained in section 5.1.

We simply assume that $f_{\rm RG}$ is $\sim 8\%$ (the value obtained in
section 5.1) over $z\sim 0-1$. This assumption comes from the fact that
previous studies of $z\lesssim 1$ galaxies have obtained similar $f_{\rm RG}$
values (Tacchella et al., 2022; Chauke et al., 2019). We also estimate the
timescale of rejuvenation $\tau_{\rm rej}$, the time from restarting star
formation to stopping the resumed star formation, to be $\tau_{\rm
rej}\left(10^{8}~{}{\rm yr}-10^{7.5}~{}{\rm yr}\right)/0.21\sim 320~{}{\rm
Myr}$ because about 21% of the RGs with $t_{\rm rej}=10^{8}~{}{\rm yr}$
returned to QGs at $10^{7.5}~{}{\rm yr}$. Note that this calculation
implicitly assumes a constant probability that the resumed star formation will
end within a certain time.

With these assumption and estimate, the averaged $N_{\rm rej}$ for all
galaxies between $z\sim 0-1$ can be calculated as $N_{\rm
rej}=\left(t|_{z=1}-t|_{z=0}\right)f_{\rm RG}/\tau_{\rm rej}\sim 2.1$,
suggesting that galaxies have experienced rejuvenation twice at $z\lesssim 1$
on average. Because $f_{\rm RG}$ increases with stellar mass as shown in
figure 14, massive galaxies can have larger $N_{\rm rej}$. For example, the
$N_{\rm rej}$ of $M^{*}>10^{11}M_{\odot}$ galaxies, with $f_{\rm RG}\sim
18\%$, can be as high as $N_{\rm rej}\sim 4.4$, meaning that such massive
galaxies are likely to have experienced multiple rejuvenation events.

As shown in figure 14, the path of a rejuvenation event in the SFR-$M^{*}$
plane is almost a vertical ascent. Considering that the time scale of
rejuvenation is $\tau_{\rm rej}\sim 0.3~{}{\rm Gyr}$ and that rejuvenation may
occur many times, essentially all rejuvenation events have to end up with re-
quenching. Therefore, we infer that massive galaxies are frequently moving up
and down on the massive side of the SFR-$M^{*}$ plane.

To verify whether rejuvenation can occur many times, we use the reconstructed
SFHs of current QGs to calculate the time from the last quenching $t_{\rm
lq}$, the time since the galaxy last entered the QG regime from the SFG or GV
regime. In this verification, we consider two patterns of quenching, “long
quenching” and “mini quenching”.

In the “long quenching”, it is assumed that galaxies never (at least for a
longer time than $\tau_{\rm rej}$) rejuvenate once quenched. This corresponds
to the popular galaxy evolution scenario in that SFGs become QGs. For
simplicity, we assume that the probability of “long quenching” for SFGs does
not change with cosmic time. This assumption is reasonable because the
cosmological evolution of quenching possibility can be ignored for our
focusing timescale; e.g., if an observed object has $z\sim 0.03$, even the
second-to-last age bin only reaches about $z\sim 0.2$, and the quenching rate
at $z\sim 0.2$ is about only two times higher than that at $z\sim 0.03$
according to Peng et al. (2010).

In contrast, “mini quenching” is temporary quenching that is terminated by
rejuvenation. Recently, Dome et al. (2023); Looser et al. (2023b) have
reported that high-$z$ (low-mass) galaxies can experience mini-quenching. Dome
et al. (2023) have found the duration of a mini-quenching event to be $20-40$
Myr at $z\sim 7$, being comparable to the dynamical time of galaxies at this
redshift. Because the dynamical time scales as $t_{\rm
dyn}\propto\rho^{-1/2}\propto\left(1+z\right)^{-3/2}$, we expect the duration
of a mini-quenching event of local galaxies to be $\sim$400-900 Myr, which is
indeed similar to the rejuvenation timescale $\tau_{\rm rej}$ of our RGs. If
galaxies can have multiple rejuvenation events, quenching other than the last
mini quenching (i.e., the long quenching and the previous mini quenchings)
will not be reflected in $t_{\rm lq}$ because $t_{\rm lq}$ only focuses on the
last quenching event. Therefore, the distribution of $t_{\rm lq}$ is expected
to have an excess on a shorter time scale than $\tau_{\rm rej}$.

In the histogram of $t_{\rm lq}$ and the number normalized by the duration of
those time bins shown in figure 19, we indeed find an excess with $t_{\rm
lq}\lesssim\tau_{\rm rej}$. This excess can be explained by considering that
the information on mini-quenching that occurred at $t_{l}\gtrsim\tau_{\rm
rej}$ is overwritten by rejuvenation and the next mini-quenching and hence
that only mini-quenching at $t_{l}\lesssim\tau_{\rm rej}$ contributes to the
histogram. However, this interpretation is inconsistent with the fact that
rejuvenation occurs at $z>0$ as well, e.g., $z\sim 0.8$ (Chauke et al., 2019;
Tacchella et al., 2022) and $z=1.88$ (Akhshik et al., 2021). We conclude that
galaxies, especially massive ones, have likely experienced multiple
rejuvenation events.

Figure 19:  The distribution of $t_{\rm lq}$. The upper panel plots the number
of RGs falling into each $t_{\rm lq}$ bin while the lower panel plots a scaled
number, i.e., the number divided by the bin width. Vertical black dashed lines
indicate $t_{\rm lq}=\tau_{\rm rej}$.

###  6.3 Why are they disk-like?

As shown in figures 15 (a) and (c), RGs have a more disk-like or rotational-
supported morphology distribution than nRQGs, even though RGs were QGs before
rejuvenation by definition. In this section, we discuss the reason for this
morphological difference between RGs and the nRQGs. Here, we propose the
following two hypotheses to explain the disk-like structure of RGs.

  1. 1.

Selective rejuvenation (SR)
As high as one-third $(=1259/3806)$ of the QGs in our sample have disk-like
morphology, i.e., $V/\sigma>1$. Furthermore, it is well known that some spiral
galaxies are anemic spirals without star formation (e.g., van den Bergh
(1976); Shimakawa et al. (2022)). Fudamoto et al. (2022) have found the
reddest spiral galaxies at $z\sim 1-3$ from JWST/NIRCam images and reported
that one of them is well-described as a passive galaxy. If only disk-like QGs,
or passive spiral galaxies, can rejuvenate, and if their morphology does not
change during rejuvenation, then RGs will have a disk-like morphology. This is
“selective” rejuvenation, where only disk-like QGs can rejuvenate.

  2. 2.

Disk-forming rejuvenation (DR)
The disk-like morphology distribution can be explained if elliptical-like QGs
form a disk during rejuvenation. This scenario does not prohibit the presence
of RGs that evolved from disk-like QGs.

We argue that the DR scenario is unlikely because of two serious problems.
First, it is extremely difficult to change morphology in a short period of
$\tau_{\rm rej}\sim 320~{}{\rm Myr}$. Although Diaz et al. (2018) claim that
elliptical galaxies can form a disk by merging with gas-rich satellite
galaxies, it takes a Gyr time scale, much longer than $\tau_{\rm rej}$.
Besides, as shown in figure 14, galaxies increase the stellar mass only by
$\sim 0.1\%$ in a single rejuvenation event. It is challenging to make a
noticeable disk with such a small amount of mass increase. Second, there are
significant differences in morphologies between RGs and nRQGs (figure 15).
These differences, especially in kinetics morphology $V/\sigma$, may not be
explained by simple star formation in the outer part of galaxies due to
rejuvenation. For these reasons, we conclude that the DR scenario is unlikely.

We then test the SR scenario using the $M^{*}$ and $n_{{\rm
dens},1h^{-1}~{}{\rm Mpc}}$ plane. As shown in figure 14, $M^{*}$ does not
change significantly before and after rejuvenation. Because $n_{{\rm
dens},1h^{-1}~{}{\rm Mpc}}$ should be almost unchanged before and after
rejuvenation, the spatial distribution of RGs and their progenitors in this
plane should be the same. By comparing the distribution of RGs, disk-like
nRQGs, and elliptical nRQGs, we can test the SR scenario.

In the SR scenario, only disk-like QGs are the progenitors of RGs. Because not
all disk-like QGs may be able to rejuvenate, the distribution of disk-like QGs
does not necessarily have to be the same as of RGs. However, the distribution
of disk-like QGs must cover that of RGs. In contrast, elliptical QGs can take
a significantly different distribution from RGs.

Figure 20 shows the $M^{*}$ and $n_{{\rm dens},1h^{-1}~{}{\rm Mpc}}$
distributions of RGs, disk-like nRQGs, and elliptical nRQGs. Here, we define
galaxies with the same $V/\sigma$ distribution with RGs as disk-like nRQGs and
construct a sample of disk-like nRQGs by selecting an nRQG with the nearest
$V/\sigma$ for each RG. On the other hand, we define galaxies with
$V/\sigma<0.5$ as elliptical nRQGs. We find that all RGs are within the area
where disk-like nRQGs are found, while a significant fraction of RGs are
outside the area where elliptical nRQGs are found. This result is consistent
with the SR scenario. Note that this result does not rule out the existence of
RGs that evolved from elliptical-like QGs. An interesting feature in this
figure is that the distribution of disk-like nRQGs has a secondary peak around
$\log\left(M^{*}/M_{\odot}\right)\sim 10$ and $\log\left(n_{{\rm
dens},1h^{-1}~{}{\rm Mpc}}\right)\sim 2.5$ where no RGs exist. This may imply
that low-mass sources in high-density regions cannot rejuvenate, consistent
with the discussion in section 5.5. In a subsequent paper, we will further
discuss the mechanism of rejuvenation, including the conditions in which
rejuvenation occurs.

Figure 20:  Same as Figure 16, but for RGs (purple), disk-like nRQGs (orange),
and elliptical nRQGs (red).

Finally, we propose further tests for the two scenarios focusing on the
mechanisms expected in these scenarios. One possible mechanism with the SR
scenario is that the remaining gas in disk-like QGs is used for rejuvenation.
In this case, disk-like QGs must have a higher molecular gas fraction $f_{\rm
gas}$ than normal QGs, and the star formation in RGs must be driven by an
increase in star formation efficiency (SFE). Another possible mechanism in the
SR scenario is that gas-poor disk-like QGs obtain gas to restart star
formation. In this case, we expect lower or similar $f_{\rm gas}$ in disk-like
QGs. In any case, the SR scenario can be tested using a combination of $f_{\rm
gas}$ and SFE data. Additionally, the two scenarios can be verified by
analyzing the spatially resolved SFH and kinetics of gases in RGs; e.g.,
evidence of gas inflows suggests the DR scenario. Thus, a specially-resolved
discussion of $f_{\rm gas}$, SFE, and kinetics using radio observation like
ALMA and IFU data from MaNGA will help us constrain the mechanism and verify
the SR and DR scenarios. These hypotheses can also be tested by using
simulation data to examine the progenitors of the RGs.

###  6.4 The contribution of AGN

In this study, we exclude 520 objects in section 2.2 as AGN or merging
galaxies. This number amounts to only $6\%$ of the size of the parent sample
and also smaller than the number of selected RGs. Therefore, AGN activities as
detected in MaNGA data are unlikely to be the primary driver of rejuvenation.
Note, however, that the parent sample may have some undetected AGNs and
merging sources. It is also possible that past AGN activities or mergers that
are not detected from current observations could trigger rejuvenation because
rejuvenation started $10^{7.5}~{}{\rm yr}$ or $10^{8}~{}{\rm yr}$ ago.

##  7 Conclusion

In this paper, we have presented the Hinotori project that aims at
understanding the nature of RGs using the SFH. As the first step of this
project, we have performed a flexible SED fitting analysis using Prospector on
8857 MaNGA galaxies with photometric and spectroscopic data and constructed a
sample of 1071 RGs, the largest ever SFH-selected RG sample. Our main results
are summarized as follows:

  * •

Tests with mock data show that our selection method can select $f\sim 0.001$
RGs with high completeness ($\sim 67\%$).

  * •

The RG fraction is $f_{\rm RG}=8.1\%$ after volume-weight correction. The
selected RGs contribute $\sim 20\%$ of the CSFRD in the recent 100 Myr.

  * •

From a rough estimate of $N_{\rm rej}$, we find that single galaxies can
rejuvenate multiple times. Multiple rejuvenations may significantly contribute
to the mass assembly of galaxies, while a single rejuvenation increases the
stellar mass only $\sim 0.1\%$.

  * •

RGs have more disk-like morphology than nRQGs, suggesting the “selective
rejuvenation” scenario where disk-like QGs rejuvenate more likely than
elliptical QGs.

The next goal of Hinotori is to specify the mechanism of rejuvenation. We plan
to analyze the IFU data of the selected RGs in the future. Using simulation
data is also an effective way to constrain the rejuvenation mechanism.
Furthermore, radio observations of molecular gas with, e.g., ALMA will allow
us to discuss the star formation activity in RGs in more detail. We also plan
to examine the relationship between AGN activity and rejuvenation (c.f.,
Martin-Navarro et al. (2022)).

Recent studies have found high-$z$ low-$M^{*}$ QGs ($M^{*}\sim 5\times
10^{8}M_{\odot}$ at $z\sim 7.3$, Looser et al. (2023a)) and post-starburst
galaxies ($M^{*}\sim 4\times 10^{7}M_{\odot}$ at $z\sim 5.2$, Strait et al.
(2023)) in JWST data. These findings may imply that high-$z$ low-mass QGs
evolve into more massive galaxies at lower $z$ through repeated rejuvenation
events (c.f., Looser et al. (2023a); Strait et al. (2023); Dome et al.
(2023)). To evaluate the importance and the role of rejuvenation in galaxy
evolution, we have to explore RGs in a wide redshift range. Besides, the
$z$-dependence of $f_{\rm RG}$ will place stronger constraints on the
possibility of multiple rejuvenations. High-$z$ RG exploration needs a deep
infrared spectroscopic survey with such as JWST/NIRSpec, TAO/SWIMS, Roman/WFI,
and Euclid/NISP. The development of a low-cost RG selection method (e.g.,
Zhang et al. (2022)) is another effective way to select RG candidates for
follow-up spectroscopic observation and advance RG research. As our results
with mock data suggest, the outshining problem is serious when we reconstruct
the SFHs of galaxies with strong rejuvenation. Spatial resolved (or pixel-by-
pixel) SED fitting may be able to solve this problem as Sorba & Sawicki (2018)
and Giménez-Arteaga et al. (2023) suggested.

We thank Keiya Hirashima and Kazunari Iwasaki for giving us useful advice for
using a high-performance computing system.

This work is partly supported by the 2022 Undergraduate Research Abroad in
Science Program of School of Science, the University of Tokyo. We thank
Roberto Maiolino and other members of the Galaxy Formation and Evolution Group
at Kavli Institute for Cosmology, Cambridge (KICC) for fruitful discussion in
group meetings and with a cup of coffee.

This research is partially supported by Initiative on Promotion of
Supercomputing for Young or Women Researchers, Supercomputing Division,
Information Technology Center, The University of Tokyo. Numerical computations
were in part carried out on Small Parallel Computers at Center for
Computational Astrophysics, the National Astronomical Observatory of Japan,
iDark cluster at Kavli Institute for the Physics and Mathematics of the
Universe, The University of Tokyo, and the Fujitsu PRIMERGY CX400M1/CX2550M5
(Oakbridge-CX) at the Information Technology Center, The University of Tokyo.

TT is supported by Forefront Physics and Mathematics Program to Drive
Transformation (FoPM), a World-leading Innovative Graduate Study (WINGS)
Program, the University of Tokyo. H.Yesuf was supported by JSPS KAKENHI Grant
Number JP22K14072 and the Research Fund for International Young Scientists of
NSFC (11950410492).

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS acknowledges support and resources from the
Center for High-Performance Computing at the University of Utah. The SDSS web
site is www.sdss4.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration including the Brazilian Participation
Group, the Carnegie Institution for Science, Carnegie Mellon University,
Center for Astrophysics — Harvard & Smithsonian (CfA), the Chilean
Participation Group, the French Participation Group, Instituto de Astrofísica
de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and
Mathematics of the Universe (IPMU) / University of Tokyo, the Korean
Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut
für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA
Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-
Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatório Nacional / MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical
communities of Japan and Taiwan, and Princeton University. The HSC
instrumentation and software were developed by the National Astronomical
Observatory of Japan (NAOJ), the Kavli Institute for the Physics and
Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High
Energy Accelerator Research Organization (KEK), the Academia Sinica Institute
for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University.
Funding was contributed by the FIRST program from the Japanese Cabinet Office,
the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the
Japan Society for the Promotion of Science (JSPS), Japan Science and
Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK,
ASIAA, and Princeton University.

This paper is based on data collected at the Subaru Telescope and retrieved
from the HSC data archive system, which is operated by the Subaru Telescope
and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out
with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ. We
are honored and grateful for the opportunity of observing the Universe from
Maunakea, which has the cultural, historical and natural significance in
Hawaii.

## Comparison of Prospector’s outputs with existing SED fitting results

Several studies have estimated galaxy parameters for the MaNGA sample with SED
fitting, although their results cannot be used to select RGs because they are
based on parametric SFHs. In this subsection, to further demonstrate the
validity of our SED fitting results with Prospector, we compare the estimated
basic parameters with existing SED fitting results in the Pipe3D DR17 catalog
(Sánchez et al., 2022). This catalog summarizes the results of a stellar
population and ionized gas analysis of the MaNGA data set using the pyPipe3D
pipeline (Lacerda et al., 2022). Note that the true values of individual
parameters, like assumed mock values in section 4.1, are unavailable and hence
that we cannot know which SED fitting results are closest to the truth.
However, we expect our results to be generally consistent with the values in
the Pipe3D catalog.

Figure 21:  Comparison of Prospector-estimated SFRs and sSFRs with the Pipe3D
catalog values. In lines 1 and 2 (lnes 3 and 4), we compare estimated SFRs
(sSFRs) in the most recent bin ($t<32~{}{\rm Myr}$) with Pipe3D SFRs (sSFRs)
in recent $32~{}{\rm Myr}$ calculated from simple stellar population analysis
and ${\rm H\alpha}$ fluxes, respectively. Rows 1, 2, and 3 are for current
QGs, GVs, and SFGs, respectively, classified with sSFRs in the most recent age
bin estimated by Prospector. Spearman’s correlation coefficient $\rho$ is
shown in the upper left corner of each panel. The solid black lines show the
$y=x$ relation. The dashed lines in $\rho>0.8$ panels show the results of
$y=x+\Delta$ fit.  Figure 22:  Same as figure 21, but in lines 1, 2, and 3, we
compare estimated $M^{*}$, $A_{\rm V}$, and $Z$ with Pipe3D values.

The Pipe3D catalog has two types of SFR; log_SFR_ssp obtained from simple
stellar population analysis and log_SFR_Ha obtained from the Ha flux. We use
both SFRs in the Pipe3D catalog in the comparison. Lines 1 and 2 of figure 21
show that the SFRs estimated with Prospector strongly correlate with those
from the Pipe3D catalog for GVs and SFGs. Indeed, Spearman’s correlation
coefficient calculated for all galaxies is $\rho\simeq 0.75$ for SSP-based
SFRs and $0.88$ for H$\alpha$-based ones. We note, however, that the SFRs by
Prospector are systematically lower by about 0.6 dex. In contrast, the
correlation is very weak for QGs, especially those with ${\rm SFR}\lesssim
10^{-1}M_{\odot}~{}{\rm yr^{-1}}$.

Line 1 of figure 22 shows that the stellar masses by Prospector strongly
correlate with those in the Pipe3D catalog (log_Mass), with a Spearman’s
correlation coefficient calculated for all types of galaxies of $\rho\simeq
0.96$, although Prospector tends to give slightly higher masses for QGs and
GVs. We also find that SFGs with higher sSFR in Prospector tend to have lower
$M^{*}$ than in the Pipe3D.

In our RG selection, we classify a galaxy into an SFG, GV, or QG at each of
the eight age bins using the sSFR and $t_{H}$ of the bin. Because the redshift
range of galaxies in the parent sample is very narrow ($0.01\leq z\leq 0.15$),
the differences in $t_{H}$ at each age bin among the galaxies are less than
$\sim 15\%$. Therefore, comparing Prospector’s sSFRs in the most recent age
bin with the Pipe3D’s values provides a sensitive test for the validity of our
RG selection.

Lines 3 and 4 of figure 21 compare sSFRs for each galaxy type. For the whole
sample, we obtain a high Spearman’s correlation coefficient of $\rho\simeq
0.84$ (vs. SSP-based sSFR) and $0.91$ (H$\alpha$-based sSFR). However, as seen
from the figure, the correlation in a given galaxy type is not so strong,
being weaker than that for SFR and $M^{*}$. This is mainly because the
classification of SFGs, GVs, and QGs is based on sSFR estimated by Prospector;
thus, the distribution in each panel has a cutoff at an sSFR.

Line 2 of figure 22 shows that the $A_{\rm V}$ estimated with Prospector are
weakly correlated with those from the Pipe3D catalog (Av_ssp_Re) with a
Spearman’s correlation coefficient (for all galaxies) of $\rho\simeq 0.48$.
The correlation is weakest for QGs. Similarly, we find a very weak correlation
for $Z$ (ZH_LW_Re_fit in Pipe3D) in line 3 of figure 22.

In summary, the values of SFR, $M^{*}$, and sSFR by Prospector correlate
strongly with those from the Pipe3D catalog despite the presence of systematic
offsets. Therefore, we conclude that our results are reliable enough to select
RGs. The criteria for RG selection are based solely on Prospector’ outputs.
This means that the systematic offsets from Pipe3D’s values will not
significantly affect our classification.

## References

  * Abdurro’uf et al. (2022) Abdurro’uf et al. 2022, ApJS, 259, 35
  * Aihara et al. (2022) Aihara, H., et al. 2022, PASJ, 74, 247
  * Akhshik et al. (2021) Akhshik, M., et al. 2021, ApJ, 907, L8
  * Alarcon et al. (2022) Alarcon, A., Hearin, A. P., Becker, M. R., & Chaves-Montero, J. 2022, MNRAS, 518, 562
  * Angthopo et al. (2020) Angthopo, J., Ferreras, I.,& Silk, J. 2020, MNRAS, 495, 2720
  * Behroozi et al. (2019) Behroozi, P., Wechsler R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143
  * Belfiore et al. (2019) Belfiore, F., et al. 2019, AJ, 158, 160
  * Belli et al. (2017) Belli, S., et al. 2017, ApJ, 841, L6
  * Bertschinger (1985) Bertschinger, E. 1985, ApJS, 58, 39
  * Blanton et al. (2017) Blanton, M. R., et al. 2017, AJ, 154, 28
  * Bundy et al. (2015) Bundy, K., et al. 2015, ApJ, 798, 7
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., Kinney, A. L., Koornneef, J., & Storchi-Bergmann, T. 2000, ApJ, 533, 682
  * Carnall et al. (2019) Carnall, A. C., et al. 2019, MNRAS, 490, 417
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Charlot and Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718
  * Chauke et al. (2019) Chauke, P., et al. 2019, ApJ, 877, 48
  * Cid Fernandes et al. (2005) Cid Fernandes, R., et al. 2005, MNRAS, 358, 363
  * Cleland and McGee (2021) Cleland, C., & McGee, D. L. 2021, MNRAS, 500, 590
  * Comerford et al. (2020) Comerford, J. M., et al. 2020, AJ, 901, 159
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, AJ, 699, 486
  * Conroy and Gunn (2010) Conroy, C., & Gunn, J. E. 2010, AJ, 712, 833
  * Davidzon et al. (2017) Davidzon, I., et al. 2017, A&A, 605, A70
  * Diaz et al. (2018) Diaz, J., Bekki, K., Forbes, D. A., Couch, W. J., Drinkwater, M. J., & Deeley, S. 2018, MNRAS, 477, 2030
  * Dome et al. (2023) Dome, T., Tacchella, S., Fialkov, A., Dekel, A., Ginzburg, O., Lapiner, S., & Looser, T. J. 2023, arXiv e-prints, arXiv:2305.07066
  * Domínguez Sánchez et al. (2022) Domínguez Sánchez, H., Margalef, B., Bernardi, M., & Huertas-Company, M. 2022, MNRAS, 509, 4024
  * Donas et al. (2007) Donas, J., et al. 2007, ApJS, 173, 597
  * Drory et al. (2015) Drory, N., et al. 2015, AJ, 149, 77
  * Dye (2008) Dye, S. 2008, MNRAS, 389, 1293
  * Elbaz et al. (2007) Elbaz, D., et al. 2007, A&A, 468, 33
  * Faber et al. (2007) Faber, S. M., et al. 2007, ApJ, 665, 265
  * Feldmann (2017) Feldmann, R. 2017, MNRASLett., 470, L59
  * Ferland et al. (1998) Ferland, G. J., Korista, K. T., Verner, D. A., Ferguson, J. W., Kingdon, J. B., & Verner, E. M. 1998, PASP, 110, 761
  * Ferland et al. (2013) Ferland, G. J., et al. 2013, Revista Mexicana de Astronomía y Astrofísica, 49, 137
  * Fudamoto et al. (2022) Fudamoto, Y., Inoue, A. K., & Sugahara, Y. 2022, ApJ, 938, L24
  * Giménez-Arteaga et al. (2023) Giménez-Arteaga, C., et al. 2023, ApJ, 948, 126
  * Gobat et al. (2017) Gobat, R., et al. 2017, A&A, 599, A95
  * Goddard et al. (2017) Goddard, D., et al. 2017, MNRAS, 466, 4731
  * Gunn et al. (2006) Gunn, J. E., et al. 2006, AJ, 131, 2332
  * Hao et al. (2019) Hao, C-N., Shi, Y., Chen, Y., Xia, X., Gu, Q., Guo, R., Yu, X., & Li, S. 2019, ApJ, 883, L36
  * Hinshaw et al. (2013) Hinshaw, G., et al. 2013, ApJS, 208, 19
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv e-prints, arXiv:1008.4686
  * Ilbert et al. (2013) Ilbvert O., et al. 2013, A&A, 556, 19
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22
  * Kaviraj et al. (2007) Kaviraj, S., et al. 2007, ApJS, 173, 619
  * Kaviraj et al. (2009) Kaviraj, S., Peirani, S., Khochfar, S., Silk, J., & Kay, S. 2009, MNRAS, 394, 1713
  * Kriek & Conroy (2013) Kriek, M., & Conroy, C. 2013, ApJ, 775, L16
  * Lacerda et al. (2022) Lacerda, E. A. D., Sánchez, S. F., Mejía-Narváez, A., Camps-Fariña, A., Espinosa-Ponce, C., Barrera-Ballesteros, J. K., Ibarra-Medel, H., & Lugo-Aranda, A. Z. 2022, NewA, 97, 101895
  * Lang et al. (2016) Lang, D., Hogg, D. W., & Schlegel, D. J. 2016, AJ, 151, 36
  * Law et al. (2016) Law, D. R., et al. 2016, AJ, 152, 83
  * Law et al. (2021) Law, D. R., et al. 2021, AJ, 161, 52
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, 837, 170
  * Leja et al. (2018) Leja, J., Johnson, B. D., Conroy, C., & van Dokkum, P. 2018, ApJ, 854, 62
  * Leja et al. (2019) Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3
  * Liivamägi et al. (2012) Liivamägi, L. J., Tempel, E., and Saar, E. 2012, A&A, 539, 80
  * Looser et al. (2023a) Looser, J. T., et al. 2023a, arXiv e-prints, arXiv:2302.14155
  * Looser et al. (2023b) Looser, J. T., et al. 2023b, arXiv e-prints, arXiv:2306.02470
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., Johnson, B. D., Conroy, C., & Davé, R. 2020, ApJ, 904, 33
  * Lupton et al. (2004) Lupton, R., Blanton, M. R., Fekete, G., Hogg, D. W., O’Mullane, W., Szalay, A., & Wherry, N. 2004, The Publications of the Astronomical Society of the Pacific, 116, 133
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * Mancini et al. (2019) Mancini, C., et al. 2019, MNRAS, 489, 1265
  * Mapelli et al. (2015) Mapelli, M., Rampazzo, R., & Marino, A. 2015, A&A, 575, 16
  * Maraston et al. (2010) Maraston, C., Pforr, J., Renzini, A., Daddi, E., Dickinson, M., Cimatti, A., & Tonini, C. 2010, MNRAS, 407, 830
  * Martin-Navarro et al. (2022) Martín-Navarro, I., Shankar, F., & Mezcua, M. 2022, MNRAS Lett., 513, 10
  * Muzzin et al. (2013) Muzzin, A., et al. 2013, ApJ, 777, 18
  * Narayanan et al. (2023) Narayanan, D., et al., 2023, arXiv e-prints, arXiv:2306.10118
  * Nelson et al. (2018) Nelson, D., et al. 2018, MNRAS, 475, 624
  * Noeske et al. (2007) Noeske, K. G., et al. 2007, ApJ, 660, 43
  * Noll et al. (2009) Noll, S., Burgarella, D., Giovannoli, E., Buat, V., Marcillac, D., & Muñoz-Mateos, J. C. 2009, A&A, 507, 1793
  * Ocvirk et al. (2006) Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut,E. 2006, MNRAS, 365, 46
  * Pandya et al. (2017) Pandya, V., et al. 2017, MNRAS, 472, 2054
  * Parikh et al. (2018) Parikh, T., et al. 2018, MNRAS, 477, 3954
  * Park et al. (2022) Park, M., et al. 2022, MNRAS, 515,213
  * Paspaliaris et al. (2023) Paspaliaris, E. D., Xilouris, E. M., Nersesian, A., Bianchi, S., Georgantopoulos, I., Masoura, V. A., Magdis, G. E., & Plionis, M. 2023, A&A, 669, A11
  * Pearson et al. (2018) Pearson, W. J., et al. 2018, A&A, 615, 146
  * Peng et al. (2010) Peng, Y-j, et al. 2010, ApJ, 721, 193
  * Popesso et al. (2023) Popesso, P., et al. 2023, MNRAS, 519, 1526
  * Rathore et al. (2022) Rathore, H., Kumar, K., Mishra, P. K., Wadadekar, Y., & Bait, O. 2022, MNRAS, 513, 389
  * Renzini et al. (2015) Renzini, A.,& Peng, Y. 2015, ApJ, 801, L29
  * Rhee et al. (2017) Rhee, J., Smith, R., Choi, H., Yi, S. K., Jaffé, Y., Candlish, G., & Sánchez-Jánssen, R. 2017, ApJ, 843, 128
  * Sánchez et al. (2016a) Sánchez, S. F., et al. 2016a, RMxAA, 52, 21
  * Sánchez et al. (2016b) Sánchez, S. F., et al. 2016b, RMxAA, 52, 171
  * Sánchez et al. (2018) Sánchez, S. F., et al. 2018, RevMex, 54, 217
  * Sánchez et al. (2022) Sánchez, S. F. et al. 2022, ApJS, 262, 36
  * Salim (2014) Salim, S. 2014, Serbian Astronomical Journal, 189, 1
  * Schawinski et al. (2007) Schawinski, K., et al. 2007, ApJS, 173, 512
  * Schlafly et al. (2019) Schlafly, E. F., Meisner, A. M.,& Green, G. M. 2019, ApJS, 240, 30
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Shimakawa et al. (2022) Shimakawa, R., Tanaka, M., Bottrell, C., Wu, P.-F., Chang, Y.-Y., Toba, Y., & Ali, S. 2022, PASJ, 74, 612
  * Smee et al. (2013) Smee, S. A., et al. 2013, AJ, 146, 32
  * Sorba & Sawicki (2018) Sorba, R., & Sawicki, M. 2018, MNRAS, 476, 1532
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15
  * Strait et al. (2023) Strait, V., et al., 2023, ApJ, 949, L23
  * Tacchella et al. (2022) Tacchella, S., et al. 2022, ApJ, 926, 134
  * Tempel et al. (2012) Tempel, E., Tago, E., & Liivamägi, L. J. 2012, A&A, 540, A106
  * Thom et al. (2012) Thom, C., et al. 2012, ApJ, 758, L41
  * Thomas et al. (2010) Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775
  * Tojeiro et al. (2007) Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, MNRAS, 381, 1252
  * Tomczak et al. (2016) Tomczak, A. R., et al. 2016, ApJ, 817, 118
  * Trayford et al. (2016) Trayford, J. W., Theuns, T., Bower, R. G., Crain, R. A., Lagos, C. d. P., Schaller, M., & Schaye, J. 2016, MNRAS, 460, 3925
  * Treu et al. (2005) Treu, T., et al. 2005 ApJ, 633, 174
  * van den Bergh (1976) van den Bergh, S. 1976, ApJ, 206, 883
  * Wake et al. (2017) Wake, D. A., et al. 2017, AJ, 154, 86
  * Walcher et al. (2011) Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Astrophysics and Space Science, 331, 1
  * Weaver et al. (2022) Weaver, J. R., et al. 2022, arXiv e-prints, arXiv:2212.02512
  * Westfall et al. (2019) Westfall, K. B., et al. 2019, AJ, 158, 231
  * Whitaker et al. (2012) Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, 29
  * Yan et al. (2016) Yan, R., et al. 2016, AJ, 152, 197
  * Zhang et al. (2022) Zhang, J., Li, Y., Leja, J., Whitaker, K. E., Nersesian, A., Bezanson, R., & van der Wel, A. 2022, ApJ, 952, 6