<figure><img src="image_1.png" /></figure>

# The 2175 Å Dust Feature in Star-forming Galaxies at $1.3 \le z \le 1.8$: The Dependence on Stellar Mass and Specific Star Formation Rate

Daichi Kashino<sup>1,2</sup>①, Simon J. Lilly<sup>1</sup>①, John D. Silverman<sup>3,4</sup>①, Alvio Renzini<sup>5</sup>①, Emanuele Daddi<sup>6</sup>①, Sandro Bardelli<sup>7</sup>①, Olga Cucciati<sup>7</sup>①, Jeyhan S. Kartaltepe<sup>8</sup>①, Vincenzo Mainieri<sup>9</sup>①, Roser Pelló<sup>10</sup>①, Ying-jie Peng<sup>11,12</sup>, David B. Sanders<sup>13</sup>①, and Elena Zucca<sup>7</sup>①

<sup>1</sup>Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 27, CH-8093 Zürich, Switzerland; kashino.daichi@mbox.nagoya-u.ac.jp

<sup>2</sup>Institute for Advanced Research, Nagoya University, Furocho, Chikusa-ku, Nagoya, 464-8601, Japan

<sup>3</sup>Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo, Kashiwa, 277-8583 (Kavli IPMU, WPI), Japan

<sup>4</sup>Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan

<sup>5</sup>INAF-Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

<sup>6</sup>CEA, Irfu, DAP, AIM, Université Paris-Saclay, Université de Paris, CNRS, F-91191 Gif-sur-Yvette, France

<sup>7</sup>INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, I-40129, Bologna, Italy

<sup>8</sup>School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

<sup>9</sup>European Southern Observatory, Karl-Schwarzschild-Strasse 2, Garching bei München, Germany

<sup>10</sup>Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, F-13388 Marseille, France

<sup>11</sup>Kavli Institute for Astronomy and Astrophysics, Peking University, 5 Yiheyuan Road, Beijing 100871, People's Republic of China

<sup>12</sup>Department of Astronomy, School of Physics, Peking University, 5 Yiheyuan Road, Beijing 100871, People's Republic of China

<sup>13</sup>Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

Received 2020 September 29; revised 2021 January 19; accepted 2021 January 21; published 2021 March 19

## Abstract

We present direct spectroscopic measurements of the broad 2175 Å absorption feature in 505 star-forming main-sequence galaxies at $1.3 \le z \le 1.8$ using individual and stacked spectra from the zCOSMOS-deep survey. Significant 2175 Å excess absorption features of moderate strength are measured, especially in the composite spectra. The excess absorption is well described by a Drude profile. The bump amplitude expressed in units of $k(\lambda) = A(\lambda)/E(B-V)$, relative to the featureless Calzetti et al. law, has a range $B_k \approx 0.2-0.8$. The bump amplitude decreases with the specific star formation rate (sSFR), while it increases moderately with the stellar mass. However, a comparison with local “starburst” galaxies shows that the high-redshift main-sequence galaxies have stronger bump features, despite having a higher sSFR than the local sample. Plotting the bump strength against the $\Delta\text{logsSFR} \equiv \log(\text{SFR}/\text{SFR}_{\text{MS}})$ relative to the main sequence, however, brings the two samples into much better concordance. This may indicate that it is the recent star formation history of the galaxies that determines the bump strength through the destruction of small carbonaceous grains by supernovae and intense radiation fields coupled with the time delay of ~1 Gyr in the appearance of carbon-rich asymptotic giant branch stars.

Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Interstellar medium (847); Interstellar dust (836); Interstellar dust extinction (837)

## 1. Introduction

A broad excess in the extinction curve at a rest wavelength $\lambda_{\text{rest}} \approx 2175$ Å (often called the “ultraviolet (UV) bump”) is the strongest signature of dust in the interstellar medium (ISM). It has attracted attention since its discovery by Stecher (1965) as a unique probe of the nature of dust in galaxies (see Salim & Narayanan 2020 for the latest review). Direct measurements of the extinction curves toward individual stars are limited in the Milky Way (MW) and a few very nearby extragalactic objects. The strong UV bump is ubiquitously seen in the extinction curves of stars in the MW (Fitzpatrick & Massa 1986) and the Large Magellanic Cloud (LMC). The extinction curves in the LMC2 supershell region, which is close to 30 Doradus (a starbursting giant HII region in the LMC), shows a noticeably weaker bump than that in the average LMC extinction curve (Fitzpatrick 1986; Gordon et al. 2003). The 2175 Å feature was also found to be weaker in M31 (Bianchi et al. 1996) and almost absent in the Small Magellanic Cloud (SMC; Prevot et al. 1984). These observations indicate that the nature of interstellar dust grains depends significantly on the local environment.

On galaxy scales, the impact of dust on the integrated spectrum is described by an attenuation curve, which includes the effects due to a more complicated spatial arrangement of dust and sources

of light than the simple case of a single star behind a dust screen. Nearby starburst galaxies have been found to have no significant UV bump in their attenuation curves, like the SMC (Calzetti et al. 1994; Gordon et al. 1997; Calzetti et al. 2000). Using radiative transfer calculations, Witt & Gordon (2000) showed that the so-called Calzetti attenuation law could be reproduced by an SMC-like extinction curve combined with a clumpy shell configuration. Vijh et al. (2003) showed that dust attenuation in Lyman break galaxies at high redshifts ($z>2$) is also consistent with that predicted with dust that has SMC-like, rather than MW-like, characteristics. The Calzetti et al. (2000) attenuation curve with no bump feature has thus been commonly assumed for both local and distant star-forming galaxies.

More recent observations, however, have detected and even accurately measured the bump feature for low- and high-redshift star-forming galaxies based on both photometry and spectroscopy (Noll et al. 2007, 2009b; Buat et al. 2011, 2012; Scoville et al. 2015; Battisti et al. 2017; Salim et al. 2018; Battisti et al. 2020; Shivaei et al. 2020). Notable work by Noll et al. (2009b) used a statistical sample of ~200 star-forming galaxies at $1 < z < 2.5$ to show that at least 30% of the sources exhibit a significant attenuation bump feature in their spectra and to directly determine the UV bump profiles using stacks. Another approach has been to
fit a galaxy’s spectral energy distribution (SED) measured in multiple photometric bands with stellar population synthesis models while applying different attenuation laws in which UV bump features of different strengths are implemented (e.g., Buat et al. 2011; Salim et al. 2018). All of these studies, however, have indicated that the strengths of the UV bump feature in $z \gtrsim 1$ galaxies are typically much weaker than observed in the MW extinction curve, usually of intermediate strength between the LMC2 supershell and the SMC extinction curves, with the latter having no bump at all.

Measurement of the UV bump across a range of redshifts is important for understanding the formation of dust grains through cosmic time. The bump feature is known to be reproduced by resonant absorption by carbon $sp^2$ bonds that efficiently occurs on the surfaces of carbonaceous small grains (grain size $a \lesssim 10$ nm $\ll \lambda$; Gilra 1971; Draine & Lee 1984; Papoular & Papoular 2009). However, the carrier(s) of the bump feature is still under debate, and other potential species have also been proposed (e.g., Bradley et al. 2005). Regardless of the precise identification of the bump carrier, analysis of the likely evolution of the distribution of grain sizes in a galaxy has succeeded in reproducing the observed extinction curves (e.g., Asano et al. 2013, 2014; Hirashita 2015). This suggests that the bump strength is linked to the relative abundance of small dust grains. The main processes determining the grain size distribution include dust production by stellar sources (see Gall et al. 2011 for a review), grain growth through metal accretion (e.g., Dwek 1998; Inoue 2011), grain destruction in the hot ISM (Dwek & Scalo 1980; McKee 1989), coagulation (e.g., Ormel et al. 2009), shattering (e.g., Yan et al. 2004),<fnref n="14" /> and possible selective removal by galactic winds (e.g., Bekki et al. 2015). The efficiencies of many of these different processes will depend on the properties of the host galaxies. However, our knowledge of the link between the galaxy properties and the UV bump characteristics is currently limited, particularly at high redshifts.

In this work, we directly measure the UV bump profiles for a large sample of star-forming galaxies at $1.3 \le z \le 1.8$ and correlate the derived bump strength with the global properties of the galaxies, specifically the stellar mass and star formation rate (SFR). We utilize the rest-frame UV spectra that have been obtained by the VIsible Multi-Object Spectrograph (VIMOS) mounted on the Very Large Telescope (VLT) UT3 in the zCOSMOS-deep survey (Lilly et al. 2007; S. J. Lilly et al. 2021, in preparation).

The paper is organized as follows. Section 2 presents an overview of the observations and describes the sample selection. Section 3 describes our spectral analysis, including the SED fitting for the sample. The results are presented in Section 4. These are discussed in Section 5, and Section 6 provides a summary of the paper. This paper uses a standard flat cosmology ($h = 0.7$, $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$). Magnitudes are quoted on the AB system. A Chabrier (2003) initial mass function (IMF) is used throughout.

## 2. Data

We investigate the possible presence of the UV bump feature in dust attenuation curves using rest-frame UV spectra obtained through the zCOSMOS-deep project. This section gives a brief summary of the survey and describes the sample selection.

### 2.1. Observations

The zCOSMOS-deep redshift survey has observed $\sim 10^4$ galaxies in the central $\sim 1$ deg$^2$ of the COSMOS field (Scoville et al. 2007). Observations were carried out with the VIMOS (Le Fèvre et al. 2003) mounted on the 8 m VLT/UT3 telescope from 2005 to 2010. Observations used the low-resolution (LR) blue grism with 1"/0 slits, yielding a spectral resolution of $R \sim 180$ and a spectral range of $\approx 3600-6700$ Å. All masks are observed with the slits oriented north-south.

The selection of the targets was performed based on a then-current version of the COSMOS photometric catalog. All of the objects were color-selected to preferentially lie at high redshifts, mostly through a $BzK$ method (Daddi et al. 2004) with a nominal $K_{AB}$ magnitude cut at 23.5, supplemented by the ugr selection (Steidel et al. 2004). An additional blue magnitude selection $B_{AB} < 25.25$ was adopted for most objects. These selection criteria yield a set of star-forming galaxies that lie mostly in the redshift range $1.3 < z < 3.5$ (Lilly et al. 2007). In total, the survey field covers a $0.92 \times 0.91$ deg$^2$ region centered on R.A. $= 10^h00^m43^s$, decl. $= +02^\circ10'23"$. The central $0.60 \times 0.62$ deg$^2$ region is referred to as the full sampling area, where the sampling rate achieved is $\approx 57\%$, while it is 30% in the outer region.

Data reduction was carried out using the VIPGI software (Scodeggio et al. 2005). Redshift measurements were visually inspected using 2D and 1D reduced spectra by identifying multiple prominent spectral features. In general, we often rely on strong absorption lines, such as Si IV $\lambda$1393, 1402, C IV $\lambda$1548, 1550, Fe II $\lambda$1608, and Al II $\lambda$1670 at $z \sim 2-3$. The Ly$\alpha$ emission line and break are key features at $z > 2.1$. At lower redshifts ($z \sim 1.3-2.0$), some redder absorption lines are important for redshift identification, such as Fe II $\lambda$2344, 2374, 2382, 2586, 2600 and Mg II $\lambda$2796, 2803. The spectroscopic redshifts have been determined purely from the spectra, independently of a photometric redshift or other information about the objects. A typical velocity error is $\approx 420$ km s$^{-1}$ ($\sigma(z)/(1+z) \approx 0.0014$) in the redshifts including the uncertainties in the wavelength calibration. We excluded the flux data in the wavelength ranges that are highly impacted by sky lines (5536–5617, 5836–5976, 6189–6399 Å) and those impacted by zeroth-order contamination of individual spectra. We will describe further spectrophotometric correction for the slit loss in Section 3.4.

### 2.2. Sample Selection

We constructed the sample from the full catalog of the zCOSMOS-deep survey (S. J. Lilly et al. 2021, in preparation). The redshift range for the current analysis is limited to $1.3 \le z \le 1.8$ to ensure that the spectral window of the VIMOS LR blue grism covers the rest-frame wavelength interval around the 2175 Å feature. To spectrophotometrically calibrate the spectra and perform SED fitting, we limited the sample to having a clear photometric counterpart in the COSMOS2015 photometric catalog (Laigle et al. 2016) and excluded a handful of targets that do not have one. Galaxies detected in X-rays were also excluded using the column type in the COSMOS2015 catalog to remove possible active galactic nuclei from the sample.

For the majority of the sample, we use the spectroscopic redshift ($z_{spec}$) measurements from zCOSMOS-deep. We use all objects with a very secure zCOSMOS-deep redshift (class = 3 or

<footnote n="14">¹⁴ See, e.g., Hirashita (2015) for a brief summary of all of these processes on dust in different size regimes.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Flowchart of the sample selection. The consistency with $z_{phot}$ is examined based on a threshold of $|z_{phot} - z_{spec}|/(1 + z_{spec}) \le 0.1$ against the photometric redshift in the COSMOS2015 catalog.</figcaption></figure>

4) within the selection range. For those with confidence class = 2, we use only those that are consistent to within $|z_{phot} - z_{spec}|/((1+z_{spec}) \le 0.1$ of the photometric redshift in the COSMOS2015 catalog. The redshifts in these two categories are both estimated to be $\ge 99\%$ reliable (S. J. Lilly et al. 2021, in preparation).

For some sources without reliable zCOSMOS-deep redshifts, we use the $z_{spec}$ measurements from the FMOS-COSMOS survey (Silverman et al. 2015; Kashino et al. 2019), which is a near-infrared (IR) spectroscopic survey of star-forming galaxies and covers most of the zCOSMOS-deep survey field. These spectroscopic redshifts are based on the detection of the H$\alpha$ emission line and some additional rest-frame optical lines (e.g., [N II] $\lambda$6584, H$\beta$, and [O III] $\lambda$5007) in the H- and J-band medium-resolution ($R \sim 3000$) spectra. We accepted the FMOS $z_{spec}$ measurements if the quality flag (zFlag) is $\ge 3$ (which corresponds to detecting $\ge 1$ line at $\ge 5\sigma$) and $z_{spec}$ is consistent with $z_{phot}$. Figure 1 shows the flowchart of selecting the sample and determining which zCOSMOS-deep or FMOS-COSMOS measurement of $z_{spec}$ is adopted.

Finally, we excluded 17 sources for which the UV continuum is barely detected or that suffered from severe spectral contamination or presented strong broad emission lines.

The final sample consists of 505 galaxies, of which $z_{spec}$ comes from the FMOS-COSMOS catalog for 160 galaxies. Figure 2 shows the redshift distribution of our sample, where the subset of 345 galaxies whose $z_{spec}$ is based on the zCOSMOS-deep catalog shows a nearly homogeneous coverage of the redshift range of interest. The additional component of 160 sources based on the FMOS-COSMOS catalog is mostly limited to the range $1.4 \lesssim z \lesssim 1.7$. The median redshift of the entire sample is $\langle z \rangle_{med} = 1.556$. As shown in Section 3.3, the galaxies in our sample are typical main-sequence (MS) galaxies in the stellar mass range of $9.5 \lesssim \log(M_*/M_\odot) \lesssim 11$.

## 2.3. Local Sample

Our goal is to measure the excess absorption due to the UV bump with respect to the so-called Calzetti law that is often adopted as a baseline attenuation curve. We thus compiled

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spectroscopic redshift distribution of our zCOSMOS-deep sample of 505 star-forming galaxies at $1.3 \le z_{spec} \le 1.8$. The filled blue histogram indicates the subset ($N = 345$) for which the zCOSMOS-deep $z_{spec}$ measurements were adopted, and the orange hatched histogram indicates the additional sources ($N = 160$) whose $z_{spec}$ comes from the FMOS-COSMOS catalog.</figcaption></figure>

nearby ($\langle z \rangle = 0.012$) starburst galaxies from Calzetti et al. (1994), Calzetti (1997), and Calzetti et al. (2000). The Calzetti et al. series of papers utilized a total of 47 local “starburst” galaxies to infer the featureless attenuation curve.<fnref n="15" /> For 23 of these, Calzetti et al. (1994) measured the excess absorption in their spectra and showed that the bump feature is, on average, negligible compared to that seen in the MW and LMC extinction laws (see their Figure 13). We will utilize these measurements together with the available reddening $E(B-V)$ measurements (see Table 1 of Calzetti et al. 1994) for comparison with our own results.

For these local sources, we adopted the estimates of the stellar mass ($M_*$) and SFR inferred from Galaxy Evolution Explorer (GALEX) and Wide-field Infrared Survey Explorer (WISE) photometry by Leroy et al. (2019). The $M_*$ and SFR estimates are available for 28 sources. For the majority of these (25/28), for which both GALEX and WISE data are available, the SFR is estimated from a linear combination of the observed luminosities in the GALEX far-UV (FUV) and WISE4 (22 $\mu$m) bands, and the stellar mass is converted from WISE1 (3.4 $\mu$m) with the SFR/WISE1-dependent mass-to-luminosity ratio. For the remaining three with no GALEX data, the SFR is estimated only from WISE4, and the mass is converted from WISE1 with the WISE4-WISE1 color-dependent mass-to-luminosity ratio. The average MS and the scatter ($\approx 0.3$ dex) of all local galaxies in Leroy et al. (2019, see Equation (19)) are consistent with other literature results.

Throughout the paper, we utilize the 34 local sources that have either $M_*/SFR$ or bump strength measurements. Of these, both the $M_*/SFR$ and bump measurement are available for 17, only the $M_*/SFR$ estimates are available for 11, and only the bump measurement is available for six.

## 3. Methodology

In this section, we describe the method used to measure the 2175 Å bump feature in the spectra. We also describe how

<footnote n="15">¹⁵ The sample sizes in three Calzetti et al. (1994, 1997, and 2000) papers are, respectively, 39, 24, and 8. Excluding duplicates, the number of unique starburst galaxies is 47.</footnote>
other galaxy properties are derived. In the following subsec-
tions, we will first introduce the formulation of the attenuation
curves and the approach to determining the attenuation curve
from precisely spectrophotometrically calibrated spectra.

We will then move to analyzing the VIMOS spectra
themselves. Every single spectrum is carefully corrected for
slit losses in a way that will not introduce artificial effects that
may mimic the UV bump. After the correction, we then
characterize the shapes of the UV continua of the individual
galaxies and also employ stacking analysis to increase the
signal-to-noise ratios (S/Ns) and precisely measure the UV
bump profiles of aggregates of galaxies. Each step of the
analysis is described in detail below.

3.1. Dust Attenuation Prescription

As usual, we define the attenuation at a given wavelength in
magnitudes, A(λ), to be

$$
A(\lambda) = 2.5 \log \left( \frac{f_{\lambda}^{\text{int}}(\lambda)}{f_{\lambda}^{\text{obs}}(\lambda)} \right), \quad (1)
$$

where $f_{\lambda}^{\text{int}}$ and $f_{\lambda}^{\text{obs}}$ are the intrinsic flux density (before dust obscuration) and the observed flux density (after dust obscuration), respectively.

The wavelength dependence of dust attenuation is formu-
lated in different ways in the literature (see Salim et al. 2018 for
a brief summary). Since the effects of dust are multiplicative,
the observed A(λ) is normalized by some measure of the
overall amount of extinction present (such as E(B − V )) to
obtain the wavelength dependence of the extinction, indepen-
dent of the overall amount of extinction. Throughout the paper,
we will refer to k(λ) ≡ A(λ)/E(B − V ) as an extinction curve if
it is measured toward individual stars and an attenuation curve
if it is for the integrated spectra of galaxies.

In this work, we will assume that the dust attenuation curves $k(\lambda)$ of high-redshift star-forming galaxies have two components: one with a smooth dependence on wavelength, which is assumed to be common to all sources, and another “excess” component from the UV bump feature at $\lambda_{\text{rest}} \approx 2175$ Å, whose strength we wish to determine. For the smooth component, we adopt the Calzetti et al. (2000) attenuation curve, $k_{\text{Cal}}(\lambda)$, and thus denote the observed attenuation curve as a sum of $k_{\text{Cal}}(\lambda)$ and an excess ($k_{\text{bump}}(\lambda)$) due to the UV bump:

$$
k(\lambda) = k_{\text{Cal}}(\lambda) + k_{\text{bump}}(\lambda). \quad (2)
$$

For the bump component, we employ a Drude profile with a
parameterization of Fitzpatrick & Massa (1990),

$$
k_{\text{bump}}(\lambda) = c_3 D(x; x_0, \gamma), \quad (3)
$$

$$
D(x; x_0, \gamma) = \frac{x^2}{(x^2 - x_0^2)^2 + x^2 \gamma^2}, \quad (4)
$$

where $x$, $x_0$, and $\gamma$ are expressed in units of inverse wavelength; $x_0$ is the inverse of the central wavelength ($x_0 = \lambda_0^{-1}$), and $\gamma$ is the FWHM of the UV bump in $\lambda^{-1}$ space. The FWHM in wavelength is therefore given as $w_\lambda \approx \gamma/\lambda_0^2$. The parameter $c_3$ is related to the peak amplitude of the excess, $B_k$, as $B_k \equiv c_3/\gamma^2$, measured from the $k(\lambda)$ attenuation curve. The fraction of the total absorption at 2175 Å that is due to the bump feature, which is denoted as $f_{\text{bump}}$, is therefore given by $f_{\text{bump}} = B_k/(B_k + k_{\text{Cal}}(2175 \text{ Å}))$, where $k_{\text{Cal}}(2175 \text{ Å})$ is fixed to have a value of 8.48.

Observationally, we require the knowledge of E(B − V ) to
derive k(λ), and thus k<sub>bump</sub>(λ) and B<sub>k</sub>, from the “observed” A
(λ). We could instead directly measure the absolute level of the
excess absorption, A<sub>bump</sub>(λ) in A(λ). We therefore define B<sub>A</sub> as
the amplitude of the bump in the observed A(λ). By definition,
this means that B<sub>A</sub> = B<sub>k</sub> × E(B − V ).

It is important to note that the measurement of $B_k$ depends on the assumed shape of the baseline smooth attenuation curve over a wide range of wavelengths (i.e., $k_{\text{Cal}}(\lambda)$ in our case). For a given observed UV continuum, with a certain slope and bump signature, an application of a baseline curve that has a steeper rise toward the FUV, such as the SMC extinction curve, will result in a lower estimated $E(B-V)$ and thus require a higher $B_k$ to account for the same level of the absolute excess absorption seen in $A(\lambda)$. On the other hand, the estimate of $B_A$ does not depend much on the shape of the baseline curve.

Our primary goal is to measure the strength of the UV bump that manifests itself in the attenuation curves, i.e., $B_k$, as a function of galaxy properties. However, we will also show the corresponding $B_A$ measurements for reference, as this is the more fundamental observable quantity and largely independent of the assumptions about the shape of the smooth attenuation curve. We will later discuss the effects of possible changes in the overall UV slope of the attenuation curve in Section 5.4.

3.2. SED Fitting

In order to measure the UV bump in an individual spectrum or stack of spectra, we require a model for the underlying “intrinsic” spectrum. Our approach to determine this model is as follows. We first fit the observed photometric SED of the galaxies with stellar population models, including the effects of a featureless attenuation curve, i.e., $k_{\text{Cal}}(\lambda)$, without any bump feature. This SED fit excludes all of the photometric bands in the rest-frame region of the bump feature. If we then compare the observed spectrum with this fitted model SED through the region of the bump, we will be able to detect whether there is any excess absorption relative to the model spectrum, i.e., whether there is any excess absorption that can be attributed to a bump feature in the attenuation law.

One could, of course, have used the excluded photometric bands in the region of the bump for this purpose (Noll et al. 2009a; Buat et al. 2011, 2012), but by using the spectra, a high-resolution spectrum of the absorption bump can be obtained, in which, for instance, the effect of discrete absorption or emission lines in the galaxy, both of variable strength, can be easily isolated and the shape of the underlying bump feature precisely determined. This method of course requires a high degree of spectrophotometric consistency between the photometry and the spectra, and this will be constructed in Section 3.4 below.

In this section, we first describe our photometric SED fitting to
obtain the model spectra that are used to measure the bump
features. We employ the software LePhare (Arnouts et al. 1999;
Ilbert et al. 2006) with a template library containing synthetic
spectra generated using the population synthesis model of Bruzual
& Charlot (2003) assuming a Chabrier (2003) IMF. We considered
12 models, combining a constant star formation history (SFH) and
delayed SFHs ($t/\tau^2 e^{-t/\tau}$ with $\tau = 0.1-3$ Gyr) with two metalli-
cities ($Z=0.008$ and 0.004) applied. The contribution of the
emission lines (Ly$\alpha$, [O II] $\lambda3727$, H$\beta$, [O III] $\lambda\lambda$ 4959, 5007, and
H$\alpha$) in the different filters is included following the recipe
described in Ilbert et al. (2009). As explained in Section 3.1, we
considered a single featureless attenuation curve from Calzetti et al. (2000) so that none of the resulting model spectra can possibly contain a bump feature. We allowed the reddening value to vary among $E(B-V) = 0, 0.01, 0.02, 0.06, 0.08$, and $0.10-0.70$ in fine steps of 0.025. The effects of possible variations in the overall shape of the baseline attenuation curve will be discussed in Section 5.4.

We used photometry from the COSMOS2015 catalog (Laigle et al. 2016) measured within 30 broad-, intermediate-, and narrowband filters from GALEX near-UV to Spitzer/IRAC ch2 (4.5 µm), as listed in Table 3 of Laigle et al. (2016). For CFHT, Subaru, and UltraVISTA photometry, we used the fluxes measured in a $3''$ diameter aperture and applied the offsets provided in the catalog to convert them to the total fluxes.

As noted above, a key feature of the analysis is that we exclude all of the photometric bands whose rest-frame central wavelengths are within $1960 \AA \le \lambda_{cen}/(1+z) \le 2440 \AA$ to ensure that the SED fitting is not affected by the SED of the galaxy in the region of the $2175 \AA$ bump feature. The number of filters excluded is four or five, depending on the redshift. The resulting SED fit is therefore quite independent of whether the actual attenuation curve of the galaxy in question does or does not have a UV bump feature.

One complicating fact is the presence of several sharp spectral features in the galaxy spectra, arising from absorption and emission in the ISM, which are not included in the stellar templates. This was dealt with in the SED fitting by applying small corrections to the photometry of the filters whose rest-frame central wavelengths are in the range $1500-2900 \AA$. Applying this correction is important in order to obtain the best possible continuum model for comparison with the observed spectrum. These corrections were determined using a stacked spectrum of the entire sample (see Section 3.6) to estimate the effects of these spectral lines for each of the photometric bands at the redshift of a given galaxy. Most individual spectra indeed do not have high enough S/Ns to allow accurate determination of the effects of the ISM emission and absorption features. Therefore, we here ignore the possible galaxy-to-galaxy variations of the strengths of these features. Although these effects of the ISM features are small (<5%) in most photometric bands, some strong absorption lines on the red side (i.e., Fe II $\lambda2344, 2374, 2382, 2586, 2600$ and Mg II $\lambda2796, 2803$) could cause $\gtrsim 10\%$ reduction of the fluxes in the intermediate and/or narrow bands that sample these particular wavelengths.

For each galaxy, we then analyze the likelihood, $\propto \exp(-\chi^2/2)$, that is associated with every possible model in the grid of models. Specifically, we computed the marginalized probability distribution function (PDF) of the flux value, $f_{\lambda}^{\text{mod}}(\lambda)$, at each wavelength grid using the likelihoods and then derived the model spectrum as the median values of the marginalized PDFs. Figure 3 presents some examples of SED fits, demonstrating that the overall SEDs are well fitted by the models. In some cases, the UV bump feature is seen as a gap between the photometric fluxes and the model spectrum. To measure the bump in units of $k(\lambda)$, we also estimated the reddening $E(B-V)$ as the median value of the marginalized PDF. Throughout the paper, the bump measurements are made by using the PDF-weighted median spectra and $E(B-V)$. The results of the bump measurements are consistent within 1% with those obtained using the best-fit model spectra and the corresponding best-fit $E(B-V)$ values that provide the maximum likelihood. However, the use of the likelihoods enables us to account for the uncertainties

in the models, not only the errors in the observed spectra. We will describe the estimation of the errors in the $B_k/B_A$ measurements using the marginalized PDFs in Section 3.7.

### 3.3. Stellar Masses and SFRs

For estimating the stellar masses of the galaxies, we again rely on SED fitting with LePhare. Throughout the paper, the stellar mass represents the mass of living stars at the time of observation and is estimated as the median value from the marginalized PDF based on the likelihoods computed for all possible models.

In the last subsection, our purpose was to obtain the model spectra that fit to the UV continuum of the individual galaxies. However, it is well known that unconstrained SED fitting often leads to unrealistically young ages, and thus low stellar masses, for $z \gtrsim 1$ star-forming galaxies when the age is left as a free parameter.<fnref n="16" /> This is because the SED of such galaxies is dominated by the youngest stellar populations, which outshine the older ones (e.g., Maraston et al. 2010). It has been demonstrated that limiting the starting times of star formation in such SED fitting to $z \sim 3-5$ better recovers the stellar masses of high-$z$ star-forming galaxies (e.g., Pforr et al. 2012).

We therefore impose a lower limit on the age of a given galaxy to estimate the stellar mass. The fitting procedure and photometric data are the same as those used to derive the model spectra in the last subsection, but now we restrict the SFH to begin before a redshift of 3 (i.e., ZFORM_MIN = 3). The resulting stellar masses of the galaxies are, on average, $\approx 0.3$ dex larger than those derived with age as a free parameter. This level of the systematic offset is consistent with what has been found in previous work (Pforr et al. 2012).

The total SFR is estimated from the UV luminosity as follows. Dust attenuation is accounted for based on the UV slope $\beta_{UV}$ of the rest-frame UV continua, which is defined as $f_{\lambda} \propto \lambda^{\beta_{UV}}$. We measured the rest-frame FUV ($1600 \AA$) flux density and $\beta_{UV}$ by fitting a power-law function to the broad- and intermediate-band fluxes within $1200 \AA \le \lambda_{cen}/(1+z) \le 2700 \AA$ but excluding the bands around the $2175 \AA$ feature, i.e., $1960 \AA \le \lambda_{cen}/(1+z) \le 2440 \AA$. The slope $\beta_{UV}$ is converted to the attenuation at $1600 \AA$ following Calzetti et al. (2000):

$$A_{1600} = 4.85 + 2.31\beta_{UV}. \qquad (5)$$

The dust-corrected UV luminosity density, $L_{1600}$, at rest frame $1600 \AA$ is then converted to SFR using a relation from Daddi et al. (2004),

$$\text{SFR}_{\text{UV,corr}}(M_{\odot} \text{ yr}^{-1}) = \frac{L_{1600}(\text{erg s}^{-1} \text{ Hz}^{-1})}{1.7 \times 8.85 \times 10^{27}}, \qquad (6)$$

where a factor of $1/1.7$ is applied to convert from a Salpeter (1955) IMF to a Chabrier (2003) IMF.

Figure 4 summarizes the distribution of stellar mass, SFR, and specific SFR (sSFR; $= \text{SFR}/M_*$) for our sample. Our sample includes galaxies with $9.5 \lesssim \log(M_*/M_{\odot}) \lesssim 11.0$ and $0.72 \lesssim \log \text{SFR}_{\text{UV,corr}}[M_{\odot} \text{ yr}^{-1}] \lesssim 2.19$ (the 2.5th–97.5th percentiles) with median values of $\langle \log(M_*/M_{\odot}) \rangle_{\text{med}} = 10.18$ and $\langle \text{SFR}_{\text{UV,corr}} \rangle_{\text{med}} = 10^{1.43} M_{\odot} \text{ yr}^{-1}$. The sSFR ranges across $-9.2 \lesssim \log \text{sSFR}_{\text{UV,corr}}[\text{yr}^{-1}] \lesssim -8.4$ with a median $\langle \text{sSFR}_{\text{UV,corr}} \rangle_{\text{med}} = 10^{-8.78} \text{ yr}^{-1}$. Note that the subset that is

<footnote n="16">¹⁶ Note that here the age represents the time elapsed since the onset of the SFH, which is thus the age of the oldest stellar population present in the galaxy.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. Examples of SED fits with LePhare. The red lines indicate the PDF-weighted median model spectra with the 16th–84th percentiles (light red). Symbols indicate the photometric fluxes in the broad bands (black circles) and the medium and narrow bands (blue squares). The photometric bands between the vertical dotted green lines are excluded from SED fitting (open symbols) so that the possible excess absorption by the bump feature does not affect the SED fit. The vertical dashed lines mark the peak wavelength of the UV bump ($\lambda_{rest}$ = 2175 Å). The object ID and redshift are noted in the top left corner of each panel.</figcaption></figure>

based on the FMOS $z_{spec}$ measurements has compensated at some level for the lack of the massive dusty population. The FMOS subset is biased toward higher stellar masses (+0.2 dex) and SFRs (+0.07 dex) and lower sSFRs (–0.05 dex) relative to the remaining subset.

Figure 5 shows the distribution of the $M_*$ and UV-based dust-corrected SFR for the sample, together with the local starburst sample (Section 2.3). Our galaxies lie closely along a linear relation (cyan line in the figure)

$$ \log \text{SFR}_{\text{UV,corr}} = 1.25 + 0.94(\log M_*/M_{\odot} - 10) \quad (7) $$

with a standard deviation of 0.18 dex. For reference, the MS $M_*$–SFR relations from the literature are shown and are all converted to a Chabrier (2003) IMF if necessary. Our sample is in good agreement with the relation derived at $z = 1.5–2.0$ by Whitaker et al. (2014). The Whitaker et al. (2014) relation precisely matches in normalization at $M_* \sim 10^{10} M_{\odot}$ the relation at the median redshift ($z = 1.56$) based on the time-dependent parameterization derived by Speagle et al. (2014) while better recovering the low-mass end slope. The agreement of our sample with the literature MS relations demonstrates that our sample is representative of MS galaxies at these epochs. We will use the Whitaker et al. (2014) relation to compute the

offsets ($\Delta\log\text{SFR} \equiv \log\text{SFR}/\langle\text{sSFR}\rangle_{\text{MS}}$) from the MS for our sample in Section 5.1. For local sources, we constructed the MS relation combining the Speagle et al. (2014) relation at $z = 0.01$ ($M_* \ge 10^{9.66} M_{\odot}$) and the one from Renzini & Peng (2015; $<10^{9.66} M_{\odot}$) for better reproducing the low-mass slope.

It should be noted that our sample of galaxies at high redshift has higher sSFRs, by roughly 0.5 dex, than the local “starburst” sample. However, relative to the evolving MS, it is the local starburst sample that has the higher $\Delta\log$ sSFR, on average. One might surmise that the high-redshift galaxies have had a more steady recent SFH than the local starbursts, which have indeed likely experienced a recent short-lived “burst” in SFR. This distinction will become significant in our interpretation (Section 5.2).

A final point to note is that the dust correction to the UV luminosity can be quite large ($A_{1600} \gtrsim 2.5$ mag) in some galaxies in our sample. We have thus checked that the conclusions do not change even if we replace the dust-corrected UV-based SFRs (SFR$_{UV,corr}$) with the UV+IR-based SFRs (SFR$_{UV+IR}$; if available), which are the sum of the UV-based SFRs but not corrected for dust attenuation and the SFR based on the total IR luminosity. For the latter, we utilized the public “super-deblended” catalog presented in Jin et al. (2018). These two SFRs are in broad agreement, and no substantial
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Histograms of stellar mass, SFR, and sSFR for our sample, all showing the number count per 0.1 dex. The SFRs are based on dust-corrected UV luminosity. The blue histograms show the subset of 345 galaxies whose zspec comes from the zCOSMOS-deep catalog, while the orange areas correspond to the additional 160 sources for which the zspec measurement from the FMOS-COSMOS catalog was adopted. The vertical dashed lines indicate the median value of the entire sample, which is also noted in each panel with the 2.5th–97.5th percentiles.</figcaption></figure>

population of heavily obscured starbursting galaxies was found in our sample. We present SFRs using the IR photometric data and the corresponding results from the reanalysis in Appendix.

### 3.4. Flux Calibration of the VIMOS Spectra

In this study, we will measure the excess absorption based on the comparison between the observed zCOSMOS spectra (and stacks of spectra) and the model spectra that were obtained from the broad-, medium-, and narrowband SED fitting. It should be recalled that these SED fits neither included a UV bump feature in the assumed attenuation curve nor considered photometric filters that were in the vicinity of the wavelength of the UV bump feature. The comparison requires an accurate spectrophotometric calibration of the observed spectrum and, in particular, correction of the slit losses that are likely to depend on the wavelength. This correction is also based on the broad- and medium-band photometry and separately calculated for each of the individual spectra that are produced from the standard zCOSMOS-deep reduction process. These spectra do have a nominal flux calibration applied to all spectra based on standard star observations, but they do not include effects such

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Stellar mass vs. SFR based on the dust-corrected UV luminosity. The cyan solid line indicates the linear fit to the data points. Orange diamonds indicate local starburst galaxies (see Section 2.3). For reference, the MS relations from the literature are shown. The thin dotted lines indicate the relations at the median redshift of the sample ($z = 1.56$) and at $z = 0.01$ based on the time-dependent parameterization derived by Speagle et al. (2014). The purple dashed line indicates a relation derived within $z = 1.5–2.0$ by Whitaker et al. (2014), which will be used to calculate $\Delta\log$ sSFR for galaxies in our sample (see Section 5.1). The black dashed line indicates the MS relation for local sources, combining the formulations from Speagle et al. (2014; $M_* \ge 10^{9.66} M_\odot$) and Renzini & Peng (2015; $<10^{9.66} M_\odot$). Note that the local starburst galaxies have lower sSFRs than the high-redshift sample but greater $\Delta\log$ sSFRs, on average, relative to the evolving MS.</figcaption></figure>

as slit width, inaccuracies in slit centering, and the effects of atmospheric dispersion. It is clear that a key requirement is that this correction must not induce any spurious feature in the spectra that may mimic any excess absorption due to the UV bump.

In constructing this correction, we make use of ~5000 spectra available in zCOSMOS-deep, regardless of whether there was a successful measurement of the spectroscopic redshift and whether the object lies in the redshift range used in the present study. For each object, we use the available photometry in four broad bands (CFHT/MegaCam *u* and Subaru/Suprime-Cam *B*, *V*, and *r*) and eight intermediate bands (Subaru/Suprime-Cam IA427, IA464, IA484, IA505, IA527, IA574, IA624, and IA679), which lie within the spectral coverage of the VIMOS LR blue grism used in zCOSMOS-deep. As was done for SED fitting (Section 3.2), we adopted the photometric flux measurements from a 3" diameter aperture from the COSMOS2015 catalog (Laigle et al. 2016) and converted these to total fluxes with the appropriate offsets. We refer to the photometric flux in the *i*th filter as $F_i^{\text{phot}}$ and the corresponding central wavelength as $\lambda_i$. We then also compute the “pseudo” broad- and intermediate-band fluxes based on the VIMOS spectra in the same filter passbands using the filter transmission curves (see Figure 6). We refer to the fluxes measured on the input spectra using the *i*th filter curve as $F_i^{\text{spec}}$. The flux ratio $R_i^{s/p} \equiv F_i^{\text{spec}}/F_i^{\text{phot}}$ is therefore a measure of the slit loss at the different wavelengths, modulo the effects of observational noise in both the imaging and spectral measurements. Here $R_i^{s/p}$ is calculated for all filters, regardless of their rest-frame wavelength.
<figure><img src="image_7.png" /><figcaption>Figure 6. Altered flux ratios $Y_i$ defined as Equation (9). The circles indicate the median values of $Y_i$ from 4972 spectra at the central wavelength of each filter, whose transmission curve is shown in the bottom part of the panel with the corresponding color. The black solid and dashed lines indicate, respectively, the 16th–84th and 5th–95th percentiles. The third-order polynomial fit to the median values, shown by the red curve, is adopted as $\xi(\lambda)$ in Equation (8).</figcaption></figure>

Our goal here is to derive a smooth functional form of the spectrophotometric correction that adequately describes the flux ratios $R_i^{s/p}$ for each spectrum. We therefore assume that the correction for a single spectrum consists of (i) an overall constant normalization, (ii) a $\lambda$-dependent term whose shape is assumed to be common to all spectra but whose amplitude may vary from spectrum to spectrum, and (iii) an additional smooth $\lambda$-dependent term that is to be determined individually for each spectrum.

The additional $\lambda$-dependence should have the simplest form that describes any residual monotonic trend of the correction in each spectrum. After some tests, we chose the following functional form for the total correction for each spectrum, in which the third term is a simple power law in wavelength:

$$ \log R^{s/p}(\lambda) = p_1 + p_2 \xi(\lambda) + p_3 \log (\lambda/5000 \AA), \quad (8) $$

where $\xi(\lambda)$ is a function that is determined from (and then applied to) all of the spectra collectively, and the three parameters, $p_1$, $p_2$, and $p_3$, are determined for each spectrum individually by fitting the flux ratio $R_i^{s/p} = F_i^{\text{spec}}/F_i^{\text{phot}}$. The parameters $p_1$ and $p_2$ determine, respectively, the overall scaling and amplitude of $\xi(\lambda)$. We found that the inclusion of the third term, which is arguably rather arbitrary, significantly improved the representation of $R^{s/p}(\lambda)$.

We first need to derive the shape of the common $\lambda$-dependence, $\xi(\lambda)$. We define $Y_i$ by dividing each $R_i^{s/p}$ to remove the overall scaling and tilt as follows:

$$ Y_i = \log R_i^{s/p} - (a + b \log \lambda_i), \quad (9) $$

where the term $(a + b \log \lambda_i)$ is a linear fit to $\log \lambda_i - \log R_i^{s/p}$. We then compile the $Y_i$ of all of the individual spectra to derive the $\xi(\lambda)$ function.

Figure 6 shows the average $\langle Y_i \rangle$ as a function of wavelength. The symbols indicate the median values in the filters at their corresponding central wavelengths, and the black solid and dashed lines indicate the 16th–84th and 5th–95th percentiles of the range of $Y_i$. To compute these medians and percentiles, we limited the spectra to those with $S/N \ge 5$ in $F_i^{\text{phot}}$ while

imposing a quite relaxed limit ($S/N \ge 1$) on $F_i^{\text{spec}}$ to avoid biasing the spectra toward those with smaller flux losses. The average $\langle Y_i \rangle$ shows a moderately curved shape decreasing (i.e., representing increasing spectral losses) at both extremities of the spectral window. This median trend in $\langle Y_i \rangle$ is well described by a third-order polynomial function of $\log \lambda$ (red curve in the figure). This overall trend is then adopted as our $\xi(\lambda)$ function.

It is worth noting that the derived $\xi(\lambda)$ function may have relatively large uncertainties at the shortest wavelengths ($\lambda \lesssim 4000 \AA$) because, as shown in Figure 6, the flux loss becomes rapidly worse at $<4000 \AA$, and we have only the CFHT/MegaCam u band there. Even worse, the u filter has a blue wing that extends down to $3300 \AA$, largely beyond the wavelength coverage. This likely induces additional uncertainties in comparisons between $F_i^{\text{phot}}(u)$ and $F_i^{\text{spec}}(u)$. However, these uncertainties in the flux calibration have negligible effects on the bump measurements because the central wavelength of the bump (2175 Å) is observed at $>5000 \AA$. We will refer to this concern later when we present the stacked spectrum in Section 3.6.

Given $\xi(\lambda)$, we can then determine the $p_1$, $p_2$, and $p_3$ in Equation (8) for each spectrum and correct the individual spectra by dividing by the resulting $\log R^{s/p}(\lambda)$. Figure 7 shows that the pseudo broad- and intermediate-band fluxes recomputed in the flux-corrected spectra are now in excellent agreement with the photometric fluxes. The median values of the corrected $F_i^{\text{spec}}$-to-$F_i^{\text{phot}}$ ratios in each band are all within $\pm 0.012$ dex, presumably reflecting small residual errors in the photometric calibration of the photometry. The scatter at a given wavelength is 0.01–0.05 dex and larger, relatively, in the intermediate bands because of their generally lower $S/N$.

The lower panel of Figure 7 shows the same data but now shifted to the rest-frame wavelength. There is no systematic trend in the corrected flux ratios as a function of the rest-frame wavelength. The running medians with a window size of 200 Å are all within 0.01 dex over the entire wavelength range of interest, and, in particular, they are quite close to zero (i.e., corrected $F_i^{\text{spec}}/F_i^{\text{phot}} = 1$) around the rest-frame wavelength of the UV bump (2175 Å). Simply as a reference, we show in this lower panel the excess absorption due to the UV bump that will be observed in the stacked spectrum of the entire sample (Section 3.6). This comparison simply shows that the accuracy of the spectrophotometric flux (re)calibration is, by a good margin, sufficient to detect the UV bump signature or, equivalently, that any residual flux calibration uncertainties can have only a very small effect on the measurement of this feature.

In the remainder of the paper, the term “observed VIMOS spectra” will always refer to these accurately spectrophotometrically recalibrated spectra.

### 3.5. Characterization of the UV Continua

In this subsection, we characterize the UV continua of the galaxy spectra in the region of the 2175 Å feature. We follow the parameterization introduced by Noll & Pierini (2005) and used in their subsequent papers (Noll et al. 2007, 2009b). The first parameter $\gamma_{34}$ represents the difference $(\gamma_3 - \gamma_4)$ between the power-law slopes $\gamma_3$ and $\gamma_4$ measured in the wavelength ranges 1900–2175 and 2175–2500 Å, respectively. Thus, the parameter $\gamma_{34}$ is a measure of the strength of the 2175 Å feature; a larger negative $\gamma_{34}$ corresponds to a stronger feature. The second observational parameter is the overall slope $\beta_b$, an
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Upper panel: individual flux ratios, $F_i^{\text{spec}}/F_i^{\text{phot}}$, after correction for the sample of 505 galaxies. Symbols for each single spectrum are connected by a line. Different spectra are colored differently for visual purposes. The red solid, dashed, and dotted lines indicate, respectively, the median values, 16th–84th percentiles, and 5th–95th percentiles. Lower panel: same flux ratios but as a function of rest-frame wavelength. The red solid and dashed lines indicate running medians and 16th–84th percentiles with a window size of 200 Å. The horizontal and vertical lines indicate $F_i^{\text{spec}}/F_i^{\text{phot}} = 1$ and the center wavelength of the UV bump (2175 Å), respectively. The dotted-dashed curve indicates the typical level of flux absorption caused by the UV bump feature measured in the stack of the whole sample (see Section 3.6).</figcaption></figure>

indicator of reddening of the UV continuum, measured across the UV bump using the rest-frame wavelength interval 1750–2600 Å but excluding 1950–2400 Å. In measuring these from the actual observed spectra, we exclude any strong narrow spectral features.

These two parameters, $\gamma_{34}$ and $\beta_b$, will vary together as a function of reddening following a relation that is determined by the dust extinction/attenuation law. We simulated these quantities using a range of Bruzual & Charlot (2003) stellar synthesis templates with typical ranges of SFHs and ages and for two metallicities ($Z = 0.004$ and 0.008). We modified these artificial spectra by applying dust attenuation with four different dust extinction/attenuation laws: the Calzetti et al. (2000) law, which is used in our SED fitting (Section 3.2); the MW extinction law from Fitzpatrick & Massa (2007); the LMC2 supershell curve; and the SMC one (Gordon et al. 2003).

Figure 8 shows that the relationship between $\beta_b$ and $\gamma_{34}$ depends significantly on the dust attenuation law, with some

<figure><img src="image_9.png" /><figcaption>**Figure 8.** The $\beta_b$ vs. $\gamma_{34}$ diagrams for model spectra. Upper panel: colored symbols indicate $\beta_b$ and $\gamma_{34}$ computed from representative Bruzual & Charlot (2003) models (a delayed SFH with $\tau = 1$ Gyr, age of 0.50 Gyr, and two stellar metallicities) in steps of $\Delta E(B-V) = 0.05$. The symbol size increases with $E(B-V)$. Different colors and symbols correspond to different dust extinction/attenuation laws, and different line types are for different metallicities of the stellar synthesis models (see legend). Gray lines are for models within plausible ranges of SFH and age, which indicate the expected scatters of the $\beta_b-\gamma_{34}$ relations due to the galaxy stellar populations for each extinction law. Lower panel: $\beta_b$ vs. $\gamma_{34}$ for attenuation laws that consist of the baseline Calzetti et al. (2000) law and the bump excess of different peak amplitudes ($0 \le B_k \le 2$, in steps of 0.05). The Bruzual & Charlot (2003) model with a delayed SFH with $\tau = 1$ Gyr, age of 0.50 Gyr, and $Z = 0.008$ is used as a representative. Circles mark $E(B-V)$ in steps of 0.05.</figcaption></figure>

scatter caused by the variation in the input intrinsic UV continuum spectra, as shown by the thin gray lines. The $\beta_b-\gamma_{34}$ relations show variations of $\Delta\gamma_{34} \lesssim 0.5$ and $\Delta\beta_b \lesssim 0.2$ at fixed $E(B-V)$. For reference, we highlight the loci for the intrinsic UV continuum spectra of a delayed SFH with $\tau = 1$ Gyr at an age of 0.50 Gyr for each constant stellar metallicity.
<figure><img src="image_10.png" /></figure>

**Figure 9.** Composite VIMOS spectrum of the entire sample of 505 galaxies at $1.3 \le z \le 1.8$. Top panel: number of spectra that have been stacked at each wavelength grid. Middle panel: stacked spectrum (black line) in which some prominent absorption and emission features are identified as marked by color-coded labels: interstellar absorption features (orange), nebular emission lines (green), and fluorescence Fe II* lines (purple). The spectral region of the Rix et al. (2004) 1978 Å index is marked by triple-dotted-dashed lines. The narrow gray region along the spectrum indicates the noise level in the composite spectrum estimated by bootstrap resampling. The blue line indicates the composite of the PDF-weighted median spectra obtained from SED fitting in which dust extinction following the bump-free Calzetti et al. (2000) attenuation law is applied. Bottom panel: 2175 Å bump feature in the attenuation curve. The black line indicates the stack of the individual $k_{\text{bump}}(\lambda)$ computed by Equation (11) for the whole sample. The red line indicates the best-fit Drude profile fit of the bump component shown.

Figure 8 demonstrates that the presence of the 2175 Å feature could, in principle, be determined by measuring any $\gamma_{34} \lesssim -1$. The figure also shows that, if all galaxies in a given sample follow a single universal attenuation law, then measurements of ($\beta_b, \gamma_{34}$) for the sample enable us to constrain the strength of the 2175 Å feature.

Similarly, we simulated $\gamma_{34}$ and $\beta_b$ by applying synthetic attenuation laws that consist of the Calzetti et al. (2000) curve as a baseline, to which is added a bump component with different amplitudes (Equation (3); the Drude profile width parameter $\gamma$ is fixed to $0.64 \ \mu\text{m}^{-1}$). The lower panel of Figure 8 shows that the gradient of the $\beta_b-\gamma_{34}$ relation becomes steeper for higher $B_k$, as expected. The relations of $B_k \approx 1.2$ yield a $\beta_b-\gamma_{34}$ relation that is similar to the prediction for the empirical LMC2 supershell extinction law.

### 3.6. Stacking Analysis and Detection of the 2175 Å Bump Feature

To accurately determine the shape of the 2175 Å bump, we rely on a stacking analysis of subsamples of galaxies selected by various intrinsic properties, specifically stellar mass and SFR. We stack the observed VIMOS spectra, as well as the PDF-weighted median model spectra derived from SED fitting, as follows. We first transform all of the individual VIMOS and model spectra to the rest-frame wavelength based on their spectroscopic redshift. Both the VIMOS and model spectra are resampled to a common wavelength grid with a spacing of

$2.0 \ \text{Åpixel}^{-1}$. Each spectrum is then normalized by the average flux density within the rest-frame 1950–2200 Å region. The composite stacked spectra are then produced by taking the median value of the individual spectra at each wavelength grid while ignoring any spectral regions that are missing and/or contaminated by, for example, zeroth-order contamination. The associated errors are estimated using bootstrap resampling.

Figure 9 shows the composite stacked spectrum (middle panel) of the whole sample. Some prominent narrow spectral features are clearly identified in the stacked spectra, as marked by vertical lines. Note that, given the redshift range of our sample, the rest-frame wavelength range of 1550–2400 Å is covered by nearly all of the sources, while shorter and longer wavelengths are less fully represented within the sample. The top panel of the figure shows the number of spectra that have contributed to the stack at each wavelength. This begins to decrease at $\lambda > 2000 \ \text{Å}$ because of the removal of a significant number of spectra that are contaminated at a given wavelength grid by sky lines that fall on the red side of the observing window (see Section 2.1).

In Figure 9, the stacked spectrum is in good agreement with the stacked model spectrum at both the blue and red ends of the spectrum, either side of the wavelength interval of the bump feature (~1900–2500 Å). To be clear, the stacked VIMOS spectra are independent of SED fitting. On the other hand, the SED fitting does not include any information from the observed spectra except the minor corrections of the photometry for the
narrow spectral features (see Section 3.2). Therefore, the agreement is a consequence of the precise flux calibration constructed in Section 3.4. Conversely, however, there is a clear gap between the stacked and model spectra across $\lambda \sim 1900-2500$ Å. Recalling that the latter was constructed using a featureless Calzetti et al. (2000) attenuation law, we can conclude that this gap is the result of a broad absorption excess that is attributed to the usual UV bump feature.

The absolute excess absorption, $A_{\text{bump}}(\lambda)$, can be measured as

$$A_{\text{bump}}(\lambda) = 2.5 \log \left( \frac{f_{\lambda}^{\text{mod}}(\lambda)}{f_{\lambda}^{\text{obs}}(\lambda)} \right), \qquad (10)$$

where $f_{\lambda}^{\text{obs}}$ and $f_{\lambda}^{\text{mod}}$ are, respectively, the observed and model spectra constructed using the featureless attenuation curve. Correspondingly, the bump component in the actual observed attenuation curve is given as

$$k_{\text{bump}}(\lambda) = A_{\text{bump}}(\lambda)/E(B-V), \qquad (11)$$

where $E(B-V)$ is the reddening value derived from the SED fitting.

To measure the bump component, we therefore first calculate $A_{\text{bump}}(\lambda)$ or $k_{\text{bump}}(\lambda)$ for the individual galaxies and then coadd them, instead of computing it directly from the stacked spectra and, for $k_{\text{bump}}(\lambda)$, a representative value of $E(B-V)$ in stacking samples. By doing so, we mitigate the possible effects of the variation in $E(B-V)$ from one galaxy to another, though the conclusion does not depend on the detailed method of stacking. The associated errors in the stacks are estimated via bootstrap resampling. See Section 3.7 for error estimation including the uncertainties in the model spectra.

The bottom panel of Figure 9 shows a Drude profile fitted to the stacked $k_{\text{bump}}(\lambda)$. For fitting, we excluded the known strong spectral lines that are marked in the figure and the region of the Rix et al. (2004) 1978 Å index (1935–2020 Å) where significant absorption features arise from the blending of numerous Fe III transitions. We fixed the central wavelength of the Drude profile, $\lambda_0$, to 2175 Å. We will come back to this result in detail in Section 4.2.

It is worth commenting on the small disagreement between the composite observed and synthesized spectra (Figure 9, middle panel) and, correspondingly, the deviation in $k_{\text{bump}}(\lambda)$ from zero at $\lambda_{\text{rest}} < 1500$ Å (bottom panel). This deviation would plausibly be due to the relatively large uncertainties in the secondary flux calibration at the bluest wavelengths and thus hardly affects the bump measurements as mentioned in Section 3.4.

### 3.7. Uncertainties in the Bump Measurements

The errors in the measurements of the bump amplitudes, $B_A$ and $B_k$, are subject not only to the errors in the observed VIMOS spectra but also to the uncertainties in the models. For simplicity, we approximately estimate the errors on $B_A$ and $B_k$ by assuming that the squared errors can be expressed as the sum of two independent terms, which are attributed to either the errors in the observed spectra or the uncertainties in the models:

$$\sigma_{\text{tot}}^2 = \sigma_{\text{obs}}^2 + \sigma_{\text{mod}}^2. \qquad (12)$$

The former term is estimated through the Drude profile fitting based on the stack and the associated error spectrum but

ignoring the uncertainties in the models (i.e., $f_{\lambda}^{\text{mod}}$ and $E(B-V)$). The error spectrum is obtained via bootstrap resampling.

Turning to the latter term, it is the flux uncertainties in the model spectra around $\lambda_{\text{rest}} \approx 2175$ Å that affect the measurements of the bump amplitude. We therefore consider the mean flux value, $f_{2175}^{\text{mod}}$, obtained by averaging $f_{\lambda}^{\text{mod}}(\lambda)$ within $\lambda_{\text{rest}} = 2075-2275$ Å. For each galaxy, we then calculate the marginal PDF of $r_{2175} = f_{2175}^{\text{mod}}/\tilde{f}_{2175}^{\text{mod}}$, where $\tilde{f}_{2175}^{\text{mod}}$ is the value obtained from the PDF-weighted median spectra. Given the marginal PDFs of $r_{2175}$ for all objects, we can calculate the model-origin uncertainties on $B_A$ as

$$\sigma_{\text{mod}}^2(B_A) = \text{var}(\text{med}(2.5 \log r_{2175})). \qquad (13)$$

Here $\text{var}(x)$ denotes the expected variance of $x$, and $\text{med}(x)$ denotes the median of $x$ taken over all objects in a subsample of a stack. The value $r_{2175}$ of each object is a stochastic variable that follows the relevant PDF.

In the case of $B_k$, we need to include the PDF of $E(B-V)$ and can approximately express the value of $\sigma_{\text{mod}}^2(B_k)$ as

$$\sigma_{\text{mod}}^2(E_k) = \text{var}\left(\text{med}\left(\frac{\langle B_A \rangle + 2.5 \log r_{2175}}{E(B-V)}\right)\right). \qquad (14)$$

Here $r_{2175}$ and $E(B-V)$ are both stochastic variables, which are assumed to be independent of each other, and $\langle B_A \rangle$ is the best estimate obtained using the stack of the subsample in question. In practice, we estimated $\sigma_{\text{mod}}^2(B_k)$ and $\sigma_{\text{mod}}^2(B_A)$ by (1) generating random $r_{2175}$ and $E(B-V)$ values for each galaxy following the marginalized PDFs, (2) calculating the median values of $2.5 \log r_{2175}$ and $(\langle E_A \rangle + 2.5 \log r_{2175})/E(B-V)$ over the subsample, and (3) repeating this process $10^4$ times with different random values to estimate the variance of these quantities. We found that the contribution to the final errors on the $B_k$ measurements from the uncertainties in the models is roughly comparable to that from the errors in the observed VIMOS spectra, i.e., $\sigma_{\text{mod}} \sim \sigma_{\text{obs}}$.

## 4. Results

In this section, we first show the characteristics of the UV continuum of individual galaxies and then present results on the detailed shapes of the 2175 Å bump feature that are obtained when the spectra are stacked according to stellar mass and SFR.

### 4.1. The Shape of the UV Continua in Individual Galaxies

Figure 10 shows the measurements of $\beta_b$ and $\gamma_{34}$ (see Section 3.5) from the individual spectra as a function of stellar mass. It is clear that these parameters are both correlated with $M_*$.

The positive correlation between $M_*$ and $\beta_b$ is a reflection of the correlation between $M_*$ and the reddening as already reported in the literature (e.g., Pannella et al. 2015). The interpretation of $\gamma_{34}$ is more complicated because the excess attenuation due to the UV bump depends on both the bump strength implemented in the attenuation curve and the absolute level of attenuation.

Figure 10 shows that $\gamma_{34}$ decreases with $M_*$, but the correlation becomes unclear at $\log M_*/M_{\odot} \gtrsim 10.5$ with a possible increase of the scatter. The correlation also appears to vanish at lower masses ($\log M_*/M_{\odot} \lesssim 9.4$) as $\gamma_{34}$ gets close to the upper limit corresponding to the zero attenuation. Note that
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Upper panel: stellar mass vs. $\beta_b$. Lower panel: stellar mass vs. $\gamma_{34}$. Blue circles indicate 452 objects with $\sigma(\gamma_{34}) < 1.0$, and gray circles are for the remaining 53 objects. Red solid and dashed lines indicate the running medians and 16th–84th percentiles of $\gamma_{34}$.</figcaption></figure>

the errors on $\gamma_{34}$ and $\beta_b$ are subject just to the noise in the observed spectra and not to the uncertainties in the SED modeling because these values are purely measured on the observed individual spectra.

Figure 11 shows $\beta_b$ versus $\gamma_{34}$ together with the simulated relations for different extinction and attenuation laws (see Figure 8). Here we show as a reference a stellar population model with $Z=0.008$ and adopting a delayed SFH with $\tau=1$ Gyr at an age of 0.50 Gyr, which is the representative case shown by solid lines in Figure 8. The majority of the sample is located between the ones for the bump-free (i.e., the Calzetti et al. 2000 and SMC) dust laws and the LMC2 supershell extinction curve, with some outliers that are presumably due to measurement errors. The location occupied by the sample is very consistent with that found by Noll et al. (2009b, see their Figures 3 and 5).

There is a clear correlation between $\beta_b$ and $\gamma$. A linear fit to the relatively secure data points ($\sigma(\gamma_{34}) < 1.0$) yields

$$ \gamma_{34} = -2.19\beta_b - 4.16. \quad (15) $$

The best-fit relation closely passes the location of zero attenuation. The slope of the fit is between the relations for the SMC and LMC2 supershell extinction curves. Thus, we can

<figure><img src="image_12.png" /><figcaption>**Figure 11.** The $\beta_b$ vs. $\gamma_{34}$ for the sample of 505 galaxies. Blue circles indicate 452 objects with $\sigma(\gamma_{34}) < 1.0$, and gray circles are for the remaining 53 objects (same as in Figure 10). Red solid and dashed lines indicate the running medians and 16th–84th percentiles of $\gamma_{34}$ in bins of $\beta_b$. The purple dotted–dashed line indicates a linear fit (Equation (15)) to the secure data points (blue circles). For comparison, the gray lines indicate the predictions for a representative model ($Z=0.008$, delayed SFH with $\tau=1$ Gyr, and age of 0.50 Gyr) that are highlighted in Figure 8.</figcaption></figure>

state that, to first order, the dust attenuation of the sample can be represented by an attenuation curve with a moderate UV bump that can nearly reproduce the best-fit $\beta_b-\gamma_{34}$ relation. However, the substantial scatter remains (rms $\approx 1.15$ in $\gamma_{34}$) after accounting for the measurement errors, which clearly indicates the presence of intrinsic variations in the attenuation curves across the sample.

Given the threshold of $\gamma_{34} = -2$ adopted by Noll et al. (2009b), we find the fraction of the sample having a nominal $\gamma_{34} < -2$ to be 30% (27% if limited to those with $\sigma(\gamma_{34}) < 1.0$). This fraction is in very good agreement with that found by Noll et al. (2009b) at $1 < z < 2.5$.

## 4.2. The UV Bump Profiles in Stacked Spectra

In this subsection, we show the direct measurements of the UV bump profiles based on different stacked spectra. As already demonstrated in Figure 9, we found that the stacked $k_{\text{bump}}(\lambda)$ of the entire sample of 505 galaxies is well described by a Drude profile. The UV bump parameters of the fit that minimize the $\chi^2$ value are given in Table 1. The amplitude is found to be $B_k = 0.538 \pm 0.008$,<fnref n="17" /> which corresponds to a fractional bump absorption of $f_{\text{bump}} = B_k/(B_k + k_{\text{Cal}}(2175\AA)) = 0.060$. The width $\gamma = 0.473 \pm 0.007 \ \mu\text{m}^{-1}$ (or $224 \pm 3 \AA$) is in agreement with that measured for $1 < z < 2.5$ galaxies with a similar approach by Noll et al. (2009b). The fiducial results have been obtained by fixing the central wavelength $\lambda_0$ to 2175 Å. If it is treated as a free parameter, we found $\lambda_0 = 2167 \AA$, while the amplitude $B_k$ and width $\gamma$ change very little.

We show in Figure 12 the stacked spectra and $k_{\text{bump}}(\lambda)$ in bins of the nominal value of $\gamma_{34}$. As expected, the stacked

<footnote n="17">¹⁷ The errors on the amplitude parameters ($B_k$ and $B_\lambda$) denote the 1$\sigma$ uncertainties, including the contribution of the uncertainties in the SED models.</footnote>
<table><caption>**Table 1**
UV Bump Parameters<sup>a</sup></caption><thead><tr><th>Sample</th><th>N</th><th>⟨γ<sub>34</sub>⟩<sub>med</sub><sup>b</sup></th><th>⟨M<sub>*</sub>⟩<sub>med</sub><sup>c</sup></th><th>⟨sSFR⟩<sub>med</sub><sup>d</sup></th><th>⟨E(B − V)⟩<sub>med</sub><sup>e</sup></th><th>γ (μm<sup>−1</sup>)</th><th>w<sub>λ</sub> (Å)<sup>f</sup></th><th>B<sub>k</sub><sup>g</sup></th><th>B<sub>A</sub> (mag)<sup>g</sup></th></tr></thead><tbody><tr><td>All</td><td>505</td><td>−1.020</td><td>10.177</td><td>−8.782</td><td>0.285</td><td>0.473 ± 0.007</td><td>224 ± 3</td><td>0.538 ± 0.008</td><td>0.149 ± 0.002</td></tr><tr><td>γ<sub>34</sub>-bin_1</td><td>84</td><td>−3.389</td><td>10.499</td><td>−8.838</td><td>0.361</td><td>0.436 ± 0.012</td><td>206 ± 5</td><td>0.802 ± 0.021</td><td>0.290 ± 0.007</td></tr><tr><td>γ<sub>34</sub>-bin_2</td><td>84</td><td>−2.234</td><td>10.300</td><td>−8.783</td><td>0.309</td><td>0.472 ± 0.014</td><td>224 ± 6</td><td>0.688 ± 0.019</td><td>0.218 ± 0.006</td></tr><tr><td>γ<sub>34</sub>-bin_3</td><td>84</td><td>−1.407</td><td>10.248</td><td>−8.791</td><td>0.286</td><td>0.465 ± 0.014</td><td>220 ± 6</td><td>0.590 ± 0.019</td><td>0.172 ± 0.005</td></tr><tr><td>γ<sub>34</sub>-bin_4</td><td>84</td><td>−0.678</td><td>10.121</td><td>−8.744</td><td>0.263</td><td>0.461 ± 0.018</td><td>218 ± 8</td><td>0.473 ± 0.019</td><td>0.124 ± 0.005</td></tr><tr><td>γ<sub>34</sub>-bin_5</td><td>84</td><td>0.047</td><td>9.984</td><td>−8.783</td><td>0.237</td><td>0.436 ± 0.023</td><td>206 ± 11</td><td>0.375 ± 0.021</td><td>0.088 ± 0.005</td></tr><tr><td>γ<sub>34</sub>-bin_6</td><td>85</td><td>0.918</td><td>10.019</td><td>−8.745</td><td>0.256</td><td>0.683 ± 0.068</td><td>325 ± 32</td><td>0.216 ± 0.021</td><td>0.050 ± 0.005</td></tr><tr><td>M<sub>*</sub>-bin_1</td><td>84</td><td>−0.118</td><td>9.754</td><td>−8.768</td><td>0.190</td><td>0.360 ± 0.023</td><td>170 ± 11</td><td>0.472 ± 0.034</td><td>0.088 ± 0.006</td></tr><tr><td>M<sub>*</sub>-bin_2</td><td>84</td><td>−0.285</td><td>9.964</td><td>−8.721</td><td>0.242</td><td>0.425 ± 0.023</td><td>201 ± 11</td><td>0.434 ± 0.023</td><td>0.103 ± 0.005</td></tr><tr><td>M<sub>*</sub>-bin_3</td><td>84</td><td>−0.805</td><td>10.100</td><td>−8.790</td><td>0.262</td><td>0.385 ± 0.017</td><td>182 ± 8</td><td>0.496 ± 0.023</td><td>0.129 ± 0.006</td></tr><tr><td>M<sub>*</sub>-bin_4</td><td>84</td><td>−1.302</td><td>10.255</td><td>−8.780</td><td>0.283</td><td>0.497 ± 0.015</td><td>235 ± 7</td><td>0.574 ± 0.018</td><td>0.161 ± 0.005</td></tr><tr><td>M<sub>*</sub>-bin_5</td><td>84</td><td>−1.806</td><td>10.450</td><td>−8.790</td><td>0.326</td><td>0.502 ± 0.014</td><td>238 ± 6</td><td>0.625 ± 0.017</td><td>0.206 ± 0.005</td></tr><tr><td>M<sub>*</sub>-bin_6</td><td>85</td><td>−2.013</td><td>10.800</td><td>−8.853</td><td>0.408</td><td>0.541 ± 0.016</td><td>257 ± 7</td><td>0.567 ± 0.015</td><td>0.226 ± 0.006</td></tr></tbody></table>

**Note.**

<sup>a</sup> The central wavelength ($\lambda_0 = 1/x_0$) was fixed to 2175 Å ($x_0 = 4.598 \ \mu\text{m}^{-1}$).

<sup>b</sup> Median values of the individual $\gamma_{34}$ measurements.

<sup>c</sup> Median values of the individual PDF-weighted median $\log(M_*/M_\odot)$ values derived from SED fitting.

<sup>d</sup> Median values of the individual log(sSFR<sub>UV,corr</sub>/yr<sup>-1</sup>).

<sup>e</sup> Median values of the individual PDF-weighted median $E(B-V)$ values derived from SED fitting.

<sup>f</sup> The FWHM in wavelength ($w_\lambda \approx \gamma_0 \lambda_0^2$).

<sup>g</sup> The errors include the effects of the uncertainties in the SED models (see Section 3.7).

spectra appear to be more bent at larger negative $\gamma_{34}$. The right panels show that the amplitude of the UV bump, $B_k$, decreases with increasing $\gamma_{34}$. The best-fit Drude parameters are given in Table 1. If $x_0 = 1/\lambda_0$ is treated as a free parameter, we found highly consistent values ($\Delta\lambda_0 = |\lambda_0 - 2175\AA| < 20\AA$; i.e., within 1%) when the bump is strong, as in the upper panels of Figure 12, while we measured $\lambda_0 = 2219 \pm 4\AA$ in the highest-$\gamma_{34}$ bin with the weakest bump. The amplitude $B_k$ and the width $\gamma$ hardly depend on whether or not the peak wavelength is fixed. Note that, in the highest $\gamma_{34}$ bin, we measured the width parameter $\gamma$ ($=0.68 \pm 0.07 \ \mu\text{m}^{-1}$; or $w_\lambda = 325 \pm 32\AA$) to be larger than the rest, with a relatively large error. This is presumably due to noise and the low $B_k$, rather than reflecting such a real wide bump. However, the bump amplitude $B_k$ itself is more or less accurate. Among these subsamples, the median $M_*$ decreases weakly with $\gamma_{34}$, and the median sSFR hardly varies (see Table 1). This means that the galaxies in a single $\gamma_{34}$ bin have a range of $M_*$ and sSFR. Even at a given $\gamma_{34}$, the variations in the intrinsic galaxy properties may bring a variation in the bump amplitude. We will see this later in Section 4.3.

Figure 13 indicates a tight negative correlation between $B_k$ and $\gamma_{34}$. The bump amplitude reaches to $B_k \approx 0.8$ in the lowest $\gamma_{34}$ bin. A concern with this approach, however, is that stacking those binned by the observed $\gamma_{34}$ may induce an artificial magnification of the bump strength, since we are constructing the stack on the basis of the quantity of interest itself. We thus repeated the analysis excluding those spectra with larger uncertainties in $\gamma_{34}$. We attempted two thresholds of $\sigma(\gamma_{34}) = 0.5$ and 1.0 while adopting the same binning grid as the original one; the numbers of galaxies in the bins are thus no longer equal. The results are shown together in Figure 13, demonstrating that the measurements of $B_k$ appear to converge in all but the lowest $\gamma_{34}$ bin, where the exclusion of lower-quality measurements results in a slightly lower $B_k$. The presence of a tight $\gamma_{34}-B_k$ correlation is robust. Our result thus confirms the result of Noll et al. (2009b, see their Figure 8),

whose $\gamma_{34}-B_k$ data points are plotted in the figure, and strengthens their statement that the more negative $\gamma_{34}$ could be associated with more prominent UV bump in the attenuation curves.

We next measure the bump strengths as a function of $\beta_b$ and $\gamma_{34}$ by stacking galaxies at grid points in the $\beta_b-\gamma_{34}$ plane in steps of 0.1 and 0.25, respectively, in $\Delta\beta_b$ and $\Delta\gamma_{34}$. Note that the size of the elliptic bin ($\Delta\beta_b = 0.2$ and $\Delta\gamma_{34} = 0.5$ in each radius) is larger than the grid separations to achieve a reasonable S/N in the stacked spectra. The left panel of Figure 14 shows a clear trend of the stacked $B_k$ that is similar to the simulations that are shown in the lower panel of Figure 8.

We also show in Figure 14 (middle panel) the observed amplitudes of the absolute excess absorption, $B_A$. The excess absorption reaches ~0.3 mag at 2175 Å in the population of lowest $\gamma_{34}$ and highest $\beta_b$. The observed $B_A$ increases with the distance from the locus that corresponds to $E(B-V) = 0$ in the upper left corner.

Lastly, the right panel of Figure 14 shows the median $E(B-V)$ in the corresponding bins. As expected, the $E(B-V)$ is tightly correlated with $\beta_b$ but almost independent of $\gamma_{34}$. Thus, we can conclude that the variations in $\gamma_{34}$ at a given $\beta_b$ reflect a real diversity of the bump strengths in the attenuation curves rather than a variation in the amount of overall reddening of the spectra.

### 4.3. UV Bump Strength versus Galaxy Stellar Mass and sSFR

In this subsection, we correlate the amplitude of $k_{\text{bump}}(\lambda)$ with galaxy properties. We first look at the global trend along the MS by stacking the individual spectra in six bins of the nominal stellar mass of the galaxies. Table 1 summarizes the bump parameters obtained from the fits, and Figure 15 shows the measured $B_k$ as a function of $M_*$. We found a tight positive correlation across $9.8 \lesssim \log M_*/M_\odot \lesssim 10.5$, though the correlation may not hold at the lowest- and highest-mass ends.

We next focus on both stellar mass and sSFR. In order to understand the dependence of the bump strength on $M_*$ and
<figure><img src="image_13.png" /><figcaption>Figure 12. Composite spectra in bins of the individual γ$_{34}$ measurements. Prominent narrow spectral features are marked as in Figure 9. Left panels: stacked VIMOS spectra (black lines) and coadded model spectra (blue lines). The median γ$_{34}$ is indicated with the number of galaxies in each panel. Right panels: coadded excess in the extinction curves, k$_{bump}$ (λ), in the same subset as the left panels. The best-fit Drude profile is overplotted in each panel (red curve).</figcaption></figure>

sSFR separately, we measure the k$_{bump}$ (λ) in stacked spectra of galaxies constructed at grid points in the log M$_{*}$–log sSFR plane in steps of 0.1 dex in M$_{*}$ and 0.05 dex in sSFR. At each grid point, objects within an elliptic bin of Δ log M$_{*}$ = 0.2 dex and Δlog sSFR = 0.1 dex in each radius were stacked. Therefore, the adjacent grid points partially share the same galaxies, and thus the measurements will be correlated. We limited bins to those containing ≥20 galaxies so as to achieve a reasonable S/N in the stacked spectra.

Figure 16 shows in the top row the stacked B$_{k}$ as a function of M$_{*}$ and sSFR. Across the accessible range, B$_{k}$ depends on both M$_{*}$ and sSFR, with a stronger dependence on sSFR. We find that log B$_{k}$ can be expressed by a linear function of log M$_{*}$

and log sSFR. The best fit is given by

$$
\begin{equation}
\begin{split}
\log B_k = {}& -0.212 + 0.191 \log(M_*/10^{10} M_{\odot}) \\
             & - 0.427 \log(\text{sSFR}/\text{Gyr}^{-1}),
\end{split}
\tag{16}
\end{equation}
$$

which is shown in the middle panel in the top row. The rms of the residuals in B_k from this fitted surface is 0.0399 (right panel). The result indicates that B_k decreases with increasing sSFR at fixed M_* and increases weakly with M_* at fixed sSFR. Note that this empirical functional form has been obtained from a limited region in the M_* versus sSFR space; thus, extrapolation of this fitted relation for galaxies lying outside the region may not be justified.
<figure><img src="image_14.png" /><figcaption>Figure 13. Upper panel: amplitude of the UV bump, $B_k$, vs. median $\gamma_{34}$ measured in the stacked spectra in six bins of the individual $\gamma_{34}$ measurements. The blue circles indicate the binned stacks using the full sample, while the open diamonds and squares are limited, respectively, to those with $\sigma(\gamma_{34}) < 1.0$ and $<0.5$. The horizontal error bars denote the 16th–84th percentiles in $\gamma_{34}$. The star indicates the stack of the entire sample. The gray triangles indicate the measurements from Noll et al. (2009b). Lower panel: median values of the individual PDF-weighted median $E(B-V)$ in the same bins as above. The vertical error bars present the 16th–84th percentiles in $E(B-V)$.</figcaption></figure>

Figure 16 also shows in the middle row the amplitude of the absolute excess absorption, $B_A$, in the same stacking bins. The amplitude $B_A$ ranges from $\approx 0.07$ to 0.25 mag, increases with increasing $M_*$, and decreases with sSFR. A linear fit to $\log B_A$ yields

$$ \log B_A = -0.857 + 0.457 \log(M_*/10^{10} M_{\odot}) - 0.303 \log(\text{sSFR}/\text{Gyr}^{-1}), \quad (17) $$

with the rms residual of $B_A$ of 0.0124 mag.

The bottom panels of Figure 16 show that the median $E(B-V)$ increases with both $M_*$ and sSFR. A linear fit yields

$$ \log E(B-V) = -0.602 + 0.251 \log(M_*/10^{10} M_{\odot}) + 0.150 \log(\text{sSFR}/\text{Gyr}^{-1}), \quad (18) $$

where the rms residual of $E(B-V)$ is 0.0092. Both positive coefficients on $M_*$ and sSFR indicate that the $E(B-V)$ is rather correlated with the absolute SFR. The presence of a tight positive correlation between the reddening and SFR (or the UV luminosity) is consistent with those found at $z \sim 2-3$ in previous studies (e.g., Meurer et al. 1999; Bouwens et al. 2009).

It may be noted that the dependence of $B_A$ (Equation (17)) can be reproduced trivially from a combination of those of $B_k$ and $E(B-V)$ (Equations (16) and (18)) and vice versa, as expected from the definition.

# 5. Discussion

## 5.1. Comparison with Local Starbursts

We have detected the broad excess absorption due to the UV bump in the attenuation curve of $1.3 \le z \le 1.8$ galaxies and mapped the dependence of this on $M_*$ and sSFR within the accessible region. The impact of the UV bump on the spectra is intermediate between that of the LMC2 supershell extinction law and that of a featureless attenuation law, like the Calzetti et al. (2000) curve or the SMC extinction law, without a bump. Our result is thus at face value different from what has been found for local starburst galaxies in the series of papers that led to the so-called Calzetti law (Calzetti et al. 1994; Calzetti 1997; Calzetti et al. 2000). We therefore here attempt to directly compare our high-$z$ sample of MS galaxies and the local starbursts from which the featureless attenuation curve has been derived. We refer back to Section 2.3 for the selection of the local starbursts.

The measurements of the absolute excess attenuation due to the UV bump are available for 23 starbursts (for 17 of them, the estimates of $M_*$ and sSFR are also available) in Calzetti et al. (1994) and thus can be compared with our $B_A$ measurements or $B_k$ by dividing by the reddening measurements of the sources, although there are some methodological differences in detail.<fnref n="18" /> An absorption excess, i.e., a positive bump amplitude, has been measured for 14/23 of the galaxies, or 9/17 if limited to those with an available $M_*$ and SFR. For the remainder, a negative amplitude has been obtained, presumably due to uncertainties in the spectral data.

In Figure 17, we show the positions of the local starburst galaxies in the $M_*$ versus sSFR diagram relative to the region in which our stacked spectra have yielded $B_k$ measurements. As noted earlier, the local “starburst” sources have sSFRs that are, on average, $\sim 0.5$ dex lower than our high-$z$ MS galaxy sample while spanning a wider range of stellar mass of $8 \lesssim \log M_*/M_{\odot} \lesssim 11$. The dependence of $B_k$ that we have observed within our high-redshift sample indicates that lower sSFR galaxies tend to have a stronger bump in the attenuation curves, and the extrapolation of this observed trend would clearly predict that the local sources should generally have much larger UV bumps. Since the majority of the local sample are located below the dashed line in Figure 17, we would definitely expect that the $B_k$ of the local sample should exceed the value observed in the stacked spectrum of our entire sample ($B_k = 0.538$). Considering the then-current data quality, it should have been easy to detect such a bump signature in the average attenuation law (see, e.g., Figure 20 of Calzetti et al. 1994).

The local starbursts, however, show on average much weaker bump strengths with the median amplitude $\langle B_k \rangle_{\text{med}} = 0.20$ ($\langle B_A \rangle_{\text{med}} = 0.028$ mag) for all 23 local sources and $\langle B_k \rangle_{\text{med}} = 0.06$ ($\langle B_A \rangle_{\text{med}} = 0.025$ mag) for the plotted 17 sources. These average excess absorptions are significantly smaller than the $B_k = 0.538$ ($B_A = 0.149$ mag) measured for our stack of the entire sample.

It is clear from the color-coding of Figure 17 that the $B_k$ values of the local starbursts are completely inconsistent with the trend established by the higher-redshift sample. We can thus conclude that local starbursts, or at least those used in the Calzetti et al. papers, have an attenuation curve in which the

<footnote n="18">¹⁸ Calzetti et al. (1994) parameterized the absorbed flux as $\eta \equiv \Delta \log f_{\lambda}$ at $\lambda_{\text{rest}} = 2175$ Å. This can be converted to our $B_A$ as $B_A = -2.5\eta$.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 14. Left panel: stacked measurements of $B_k$ as a function of $\beta_b$ and $\gamma_{34}$. The ellipse in the lower left corner indicates the bin size. The gray lines indicate the predictions for a representative model ($Z = 0.008$, delayed SFH with $\tau = 1$ Gyr, and age of 0.50 Gyr), and the purple dotted-dashed line indicates the linear fit to the individual measurements (Equation (15); same as in Figure 11). Middle panel: stacked measurements of $B_A$ in the same stacking bins. Right panel: median $E(B-V)$ in the same stacking bins.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Upper panel: $B_k$ vs. median $M_*$ in six bins of the individual $M_*$ estimates. The horizontal error bars present the 16th–84th percentiles in $M_*$. Lower panel: median values of the individual PDF-weighted median $E(B-V)$ in the same bins as above. The vertical error bars present the 16th–84th percentiles in $E(B-V)$.</figcaption></figure>

UV bump is, on average, substantially weaker than what is seen in $1.3 \le z \le 1.8$ star-forming galaxies, despite the latter's higher sSFR.

We here recall, however, that the local “starbursts” are, as their name implies, almost certainly undergoing a substantial, recent, and short-lived elevation of their SFR. In contrast, the galaxies in the higher-redshift sample are close to the MS and have probably been forming stars at a more or less steady rate.

To gain further insights into how the bump feature depends on the star formation activity, we now compare the high-z sample and the local starbursts by renormalizing their sSFRs to that of the MS at the appropriate epoch; i.e., we consider the $\Delta\log$sSFR = $\log(sSFR/\langle sSFR \rangle_{MS})$ relative to the appropriate MS. We utilize the MS relation at $z = 1.5-2.0$ derived by Whitaker et al. (2014) to compute $\Delta\log$sSFR for our high-z sample and a relation combining the Speagle et al. (2014) relation at $z = 0.01 (M_* \ge 10^{9.66} M_\odot)$ and the one from Renzini & Peng (2015; $<10^{9.66} M_\odot$) for the local sources (see Figure 5).

We then rebin our sample in the $\log M_*$ versus $\Delta\log$sSFR plane in steps of, respectively, 0.1 and 0.05 dex using elliptic apertures with radii of 0.2 and 0.1 dex in each axis.

Figure 18 summarizes the new bump measurements, $B_k$ and $B_A$, and the reddening $E(B-V)$ for our high-z MS galaxies, together with the local starbursts, in the $\log M_*$ versus $\Delta\log$sSFR plane. Our sample covers $\pm 0.3$ dex with respect to the MS. It may be worth providing the empirical fits to these measurements as a function of $M_*$ and $\Delta\log$sSFR. The best linear fits for our sample are given as follows:

$$ \log B_k = -0.319 + 0.245 \log(M_*/10^{10} M_\odot) \\ - 0.413 \Delta\log$sSFR, \qquad (19) $$

$$ \log B_A = -0.933 + 0.497 \log(M_*/10^{10} M_\odot) \\ - 0.268 \Delta\log$sSFR, \qquad (20) $$

$$ \log E(B-V) = -0.602 + 0.251 \log(M_*/10^{10} M_\odot) \\ + 0.150 \Delta\log$sSFR. \qquad (21) $$

The rms residuals in $B_k$, $B_A$, and $E(B-V)$ are, respectively, 0.0378, 0.0113, and 0.0092.

It is obvious in the top panel of Figure 18 that the majority of the local starbursts differ from the local MS, lying well above the commonly used threshold of $\Delta\log$sSFR > 4 (0.6 dex) for starbursts. Interestingly, the observed UV bump strengths, $B_k$, in the local starbursts are now in much better agreement, at least qualitatively, with the extrapolation of the trends in our high-z sample that predicts weaker bumps for higher $\Delta\log$sSFR. This suggests that the weak or absent UV bump in the local starbursts is linked to their large positive offset from the MS, i.e., their $\Delta\log$sSFR, rather than their absolute SFRs or sSFRs. In the following subsection, we will come back to this point in interpreting the observed behavior of the bump feature.

## 5.2. Interpretation of the Dependence of $B_k$

The variations in bump strength can be produced by the changes in the nature of dust grains, such as size distribution and/or composition, or the radiative transfer effects in different dust-to-star geometries. Here we attempt to interpret the observed trends of $B_k$, first in terms of the underlying grain properties and later in terms of geometrical effects in the subsequent subsection.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Top row: stacked measurements of $B_k$ as a function of $M_*$ and sSFR (left panel), the linear fit to $B_k$ (Equation (16); middle panel), and residuals from the fit (right panel). The large ellipse indicates the bin size with radii of 0.2 dex in $M_*$ and 0.1 dex in sSFR. Small circles indicate the individual galaxies. For the residuals, the color range matches that of the left panel, but the middle value is shifted to zero. Middle row: same as the top row but for the stacked measurements of $B_A$ in the same stacking bins. Bottom row: for the median $E(B-V)$ in the same stacking bins.</figcaption></figure>

The identification of carriers of the UV bump is a long-standing topic of controversy (see Draine 1989 for a review). Carbonaceous grains containing sp²-bonded structures are the most widely accepted materials because they exhibit a broad excess in the absorption cross section at $\lambda \sim 2200$ Å due to resonant absorption in the Rayleigh limit (i.e., grain size $a \lesssim 0.01$ µm $\ll \lambda$; e.g., Gilra 1971). Multiple possible forms of the carbonaceous bump carrier candidates have been considered, including (partially) graphitized particles, a random assembly of microscopic sp² carbon chips (Papoular & Papoular 2009), and polycyclic aromatic hydrocarbons (PAHs; e.g., Mathis 1994; Steglich et al. 2010; Hirashita & Murga 2020). In the following discussion, we assume that the UV bump is attributed to small ($a \lesssim 0.01$ µm) carbonaceous grains, whatever their form.

The strength of the bump feature will reflect the destruction and creation of the bump carriers. The destruction of small-sized ($a \lesssim 0.01$ µm) grains that may be able to give rise to the bump feature, whatever the chemical composition, is accelerated by an enhanced frequency of supernovae (SNe) and/or intense hard radiation fields (Gordon & Clayton 1998; Sloan et al. 2008).

Harder and more intense radiation fields may arise due to the lower metallicity in less massive and/or higher sSFR galaxies (e.g., Zahid et al. 2014; Sanders et al. 2018), and these may more efficiently destroy small grains (Sloan et al. 2008).

The production rate of small grains is also a factor. It is thought that dust grains of size $a \sim 0.1$ µm are produced by SNe and asymptotic giant branch (AGB) stars. Carbonaceous grains that could be the precursors of bump carriers would thus be attributed mainly to carbon-rich AGB stars (Gall et al. 2011). Because only moderate-mass stars ($\lesssim 3 M_\odot$) can be carbon stars (Renzini & Voli 1981), the dust production in the AGB phase achieves its peak efficiency $\sim 1$ Gyr (the lifetime of $\approx 2 M_\odot$ stars) after star formation (Zhukovska et al. 2008). Contrarily, the time delays in the appearance in Type II SNe and oxygen-rich AGB stars, whose progenitor masses are $> 8 M_\odot$ and $\sim (3-8) \times M_\odot$, respectively, are significantly shorter than 1 Gyr. The ejected carbonaceous grains then need to be shattered into small grains and, if amorphous, at least partially graphitized and/or aromatized to be bump carriers. The timescales of these processes, however, are both thought to be shorter ($\sim 100$ Myr) than the time delay in the appearance of the carbon-rich AGB stars (Sorrell 1990; Hirashita 2015).
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Bump amplitude, $B_k$, of the local starbursts from the Calzetti et al. papers (Section 2.3) in the $M_*$ vs. sSFR plane compared with our results (colored grid). The filled diamonds indicate the local starbursts and are color-coded by the bump amplitude, but the black filled diamonds correspond to the negative measurements of the bump. The white diamonds indicate the local sources for which no individual bump measurements are available. The color scale has been changed from Figure 16 to cover the measurements of the local sources down to zero. The dashed line indicates Equation (16) at $B_k = 0.538$, which is obtained for our stack of the entire sample. Note that a single red diamond has $B_k=3.96$, as labeled, which is largely exceeding the upper limit of the color range.</figcaption></figure>

We attempt to connect the observed dependence of $B_k$ to possible processes relevant to the evolution of dust grains. In the Local Group (i.e., among the MW, LMC, and SMC), more massive and/or metal-rich galaxies tend to have extinction curves that exhibit stronger bumps. At $z \sim 2$, Shivaei et al. (2020) reported that attenuation curves of higher-metallicity galaxies exhibit stronger bumps. To first order, the observed $M_*-B_k$ correlation (Figure 15) appears to be consistent with these trends, given the metallicity being tightly correlated with $M_*$ (e.g., Zahid et al. 2014; Kashino et al. 2017). The abundance of PAHs, estimated from the mid-IR luminosity, is known to depend strongly on the metallicity of the galaxies, which may be linked to the possible metallicity dependence of the bump strength given PAHs being a plausible candidate of the bump carrier (e.g., Engelbracht et al. 2005; Galliano et al. 2008; Shivaei et al. 2017). The PAH-to-dust abundance ratio increasing with metallicity is also predicted by theoretical models in terms of the changes in the efficiency of shattering and coagulation processes, which both play roles in producing small grains (Seok et al. 2014; Rau et al. 2019).

Moreover, our results, obtained from binning the sample onto the $M_*$-sSFR plane, indicate that the $M_*$ dependence is only moderate when fixing sSFR. This implies that an apparently stronger $M_*$ dependence may be derived if measuring $B_k$ only along the $M_*$ axis due to the anticorrelation between $M_*$ and sSFR, as well as between $B_k$ and sSFR within a representative sample of MS galaxies. From now on, we focus on the dependence of $B_k$ on sSFR at fixed $M_*$ and will propose a simple approach to simultaneously explain our findings and the absence of the bump in local starbursts.

Knowledge of the physical processes playing on dust grains, as noted above, suggests that a key quantity in the balance between production and destruction of the bump carriers may be the recent SFR as measured at $\sim 10^{7-8}$ yr compared with the SFR of order 1 Gyr ago (i.e., the time delay in the appearance of carbon-rich AGB stars). For simplicity, we may denote, following Wang & Lilly (2020), the ratio of the current instantaneous SFR (measured on $\sim 10^{7-8}$ yr timescales) to that

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Bump amplitude, $B_k$ (top) and $B_A$ (middle), and mean best-fit $E(B-V)$ (bottom) as a function of $M_*$ and $\Delta\log$sSFR = $\log(s\text{sSFR}/\langle\text{sSFR}\rangle_{\text{MS}})$ for our sample (colored grid) and the local starbursts (diamonds). The horizontal dotted lines indicate $\Delta\log$sSFR = 0. The dashed line in the top panel indicates Equation (20) at $B_k = 0.538$, which is obtained for our stack of the entire sample.</figcaption></figure>

1 Gyr ago as SF$_{79}$. It should be noted that this definition is actually slightly different from that in Wang & Lilly (2020), who defined SF$_{79}$ to be the ratio of the instantaneous SFR to that averaged over the previous 1 Gyr.

We may then look at the behavior of SF$_{79}$ in two idealized cases. First, for galaxies with a constant sSFR, both the stellar mass$^{19}$ and the SFR increase exponentially with an e-folding

<footnote>$^{19}$ Here we are referring to an sSFR computed using the stellar mass as the integral of the past SFH, so that the $sSFR^{-1}$ is the mass doubling timescale. This is different from the definition of the stellar masses estimated for the individual galaxies in our sample, which denotes the mass in stars that have survived to the time of observation.</footnote>
timescale given by the inverse sSFR⁻¹. The change in SFR will therefore depend on the sSFR multiplied by the time interval of interest, which in this case is of order 1 Gyr:

$$ \ln SF_{79} \sim sSFR \times 1 \text{ Gyr.} \qquad (22) $$

It is then easy to see that different galaxies with different constant sSFRs will have different SF₇₉. We can approximately write these in terms of some average fiducial value as

$$ \Delta \log SF_{79} \sim \Delta \log sSFR \times sSFR_{\text{Gyr}} \qquad (23) $$

if the offsets are small (≪1 dex). Here sSFR<sub>Gyr</sub> = sSFR × 1 Gyr, i.e., the sSFR in units of Gyr⁻¹. When the inverse sSFR⁻¹ timescale becomes comparable to the gigayear timescales of interest, i.e., sSFR<sub>Gyr</sub> ~ unity, as is certainly the case at high redshift z ~ 2, variations in the (steady) sSFR will cause significant variations in SF₇₉, with a consequent effect expected on the amplitude of the absorption bump. This effect alone could conceivably account for the observed trend of decreasing bump strength with sSFR that is seen within our high-redshift sample of MS galaxies. Note, however, that at low redshift, where sSFR ≪ 1 Gyr⁻¹, it can be seen that variations in the steady sSFR will have a much smaller effect on SF₇₉.

However, a potentially much larger effect on SF₇₉ will be produced by any rapid temporal changes in the sSFR. In particular, even at low redshift, when most galaxies have sSFR ≪ 1 Gyr⁻¹, a short sharp increase in the SFR, as in a “starburst,” could lead to a corresponding increase in SF₇₉. If the burst is of short duration, much less than 1 Gyr, then it simply follows (independent of the value of the fiducial sSFR) that

$$ \Delta \log SF_{79} \sim \Delta \log sSFR. \qquad (24) $$

This second effect is almost certainly the one that is relevant for the low-redshift sample of starburst galaxies. Remarkably, it can be seen that Equations (23) and (24) have exactly the same form if Equation (23) is evaluated in the regime where sSFR ~ 1 Gyr⁻¹, i.e., at z ~ 2.

This undoubtedly simplified picture therefore provides a natural explanation of the effects in Figure 18, in which high-redshift MS galaxies and local starburst galaxies are combined; high-redshift MS galaxies in a quasi-steady state and local starburst galaxies undergoing a short temporal elevation in their sSFR happen to display the same relation when the amplitude of the bump feature is plotted against the Δlog sSFR relative to the appropriate MS but quite different relations when plotted against the sSFR. The success of this simple approach in explaining the observations presented in this paper adds support to the idea that the abundance of the carrier grains of the 2175 Å excess absorption is strongly linked to variations in the recent SFH of the galaxies over the last billion years or so.

This picture, however, appears not to apply, at least immediately, to local MS galaxies; a prediction from Equation (23) is that we would expect to see a very much weaker or absent trend of bump strength with Δlog sSFR for MS galaxies at low redshift with sSFR<sub>Gyr</sub> ≪ 1. On the other hand, Salim et al. (2018) found that more massive galaxies tend to have attenuation curves that are steeper and exhibit weaker bumps, and as a secondary dependence, the slope and the bump strength become, respectively, steeper and stronger toward both the upper and lower envelopes of the MS at fixed M*. In parallel, no dependence on the gas-phase metallicity was found

in both the slope and the bump strength at fixed M*. These findings appear to conflict with the trend seen in the extinction curves in the Local Group galaxies and M* and metallicity trends seen at high redshifts (this paper and Shivaei et al. 2020). As discussed in Salim et al. (2018), their results suggest that the dust-to-star geometries may play an important role in regulating the attenuation curves in these local galaxies. We will discuss the effects of geometry in the next subsection.

We should also mention a caveat in the assumption that is a base of our interpretation. Here the bump carriers are assumed to be small carbonaceous grains that were initially supplied by carbon-rich AGB stars and then processed in the ISM to be able to give rise to the bump feature. Although many theoretical models employ this general picture (e.g., Zhukovska et al. 2008; Asano et al. 2013; Hirashita & Murga 2020), on the other hand, there is a known paradoxical problem that the bump strength appears to be anticorrelated with the abundance ratio of carbon-rich to oxygen-rich stars (the C/M ratio); the SMC (with no bump) presents the highest C/M, while the MW (with a strong bump) presents the lowest among the Local Group objects (e.g., Cook et al. 1986; Groenewegen 1999; Mouhcine & Lancon 2003). This established observational fact appears to disfavor carbonaceous grains being the bump carriers, or at least scenarios in which the supply of bump carriers is owed to carbon-rich AGB stars. There thus remains a tension between physical arguments favoring carbon grains and astrophysical arguments disfavoring them.

### 5.3. Effects of Dust-to-star Geometry

In the last subsection, we proposed an interpretation of the sSFR dependence of B<sub>k</sub> in terms of the production and destruction of the bump carriers, or, in other words, the actual variations in the underlying extinction curves. However, even if the extinction curve is constant, the radiative transfer effect in different dust-to-star geometries can produce a wide variation in the shape of the attenuation curves. Radiative transfer calculations have demonstrated that, in general, the attenuation curve tends to be shallower (or flatter) toward the FUV part and exhibit a weaker bump for a higher dust column density (or higher overall attenuation), a clumpier ISM, and a geometry in which the stars and dust are mixed (instead of a geometry with a screen or shell of dust; e.g., Witt & Gordon 1996, 2000; Seon & Draine 2016; Narayanan et al. 2018). The radiative transfer effects have also been indicated in observations, especially through detecting the correlation between the slope and the bump strength (shallower slope ⇔ weaker bump) and/or between the slope and the optical depth (higher A<sub>V</sub> ⇔ shallower slope; Kriek & Conroy 2013; Chevallard et al. 2013; Salmon et al. 2016; Reddy et al. 2018; Salim et al. 2018; Tress et al. 2018; Battisti et al. 2020).

It may be natural to expect that galaxies with higher sSFRs tend to have more complex, clumpy geometries of stars and dust. Indeed, Kriek & Conroy (2013) showed that galaxies having a higher Hα equivalent width tend to have attenuation curves that are flatter and exhibit weaker bumps. Battisti et al. (2020) found a similar correlation with sSFR for z ~ 0.1 galaxies. Furthermore, Shibuya et al. (2016) found that the fraction of clumpy galaxies increases with sSFR across 0 ≲ z ≲ 3. The effects of dust-to-star geometry are therefore to be considered as an alternative explanation of the observed sSFR–B<sub>k</sub> anticorrelation, which does not require any changes in the underlying extinction curve, i.e., the nature of dust grains.
Distinguishing the effects of different geometries and grain compositions needs more detailed investigation.

### 5.4. The Assumed Baseline Attenuation Curve

In this work, we adopted the Calzetti et al. (2000) law as the baseline of the attenuation curves for all galaxies in the sample, which is the same approach as employed by Noll et al. (2009b). The evolution of the dust population, however, implies that not only the bump profile but also the overall shape of the attenuation curve may be different at high redshift, as is the case locally. As already mentioned above, the global shape of the attenuation curve also depends on the geometrical configuration of dust and stars (e.g., Witt & Gordon 2000; Narayanan et al. 2018).

For a given UV continuum, the measurement of $B_k$ depends on the overall slope of the baseline attenuation curve. Application of a baseline curve that has a steeper rise toward the FUV results in a lower $E(B-V)$ for a given observed UV slope and thus a higher $B_k$ for a given level of absolute excess absorption due to the UV bump. Possible variations in the underlying attenuation curve also have impacts on the SFR derived from the UV luminosity. Applying a steeper attenuation curve will lead to a lower level of dust attenuation and thus a lower SFR for a given UV continuum.

Observationally, Scoville et al. (2015) found that star-forming galaxies at $2 < z < 6.5$ have an average attenuation curve that is very similar to the Calzetti et al. (2000) curve in the overall shape and present a moderate UV bump feature. On the other hand, there have been several claims that $z \sim 2$ star-forming galaxies have an attenuation curve that is steeper than the Calzetti et al. (2000) curve (Buat et al. 2012; Salmon et al. 2016; Reddy et al. 2018; Battisti et al. 2020). For example, Buat et al. (2012) found an average $\langle B_k \rangle = 1.6$ and $\langle \delta \rangle = -0.27$ for galaxies at $1 \lesssim z \lesssim 2$ applying a “modified” Calzetti et al. (2000) curve,

$$k_{\text{mod}}(\lambda) = \frac{R_V}{4.05} k_{\text{Cal}}(\lambda) \left( \frac{\lambda}{5500 \AA} \right)^{\delta} + k_{\text{bump}}(\lambda), \quad (25)$$

where $\delta$ modifies the slope. Their $\langle B_k \rangle = 1.6$ is higher than our measurements at face value. If we adopted the modified baseline curve using their average value of $\delta = -0.27$, however, we obtained $B_k \approx 1.5$ from the stack of the entire sample, which is in agreement with the average value of Buat et al. (2012).

If all galaxies in the sample can be assumed to follow a common baseline attenuation, whatever the slope, then the dependence of $B_k$ that we found would qualitatively still hold. On the other hand, there are some claims that lower sSFR galaxies tend to have a steeper attenuation curve (Kriek & Conroy 2013; Battisti et al. 2020). If this is true, then our analysis may have underestimated $B_k$ and overestimated SFR in lower sSFR galaxies. Subsequently, we may have artificially narrowed both the range in sSFR (or the width of the MS, i.e., the range of $\Delta \log$ sSFR) and the range in $B_k$; thus, the effects on the coefficients in Equations (16) and (20) would be small.

Though the conversion to the amplitude in $k(\lambda)$ depends on the assumption, we stress that the presence of the UV bump in the attenuation curves has been robustly confirmed through the pure observables such as $B_A$ and $\gamma_{34}$ for our sample. Attenuation curves without a UV bump cannot explain the shape of the stacked spectra that are bent at $\lambda \approx 2175 \AA$. The observed trend in $B_A$ (Figures 14, 16, and 17), therefore, is

essentially independent of the assumption of the shape of the baseline attenuation curve and is thus useful for predicting the absolute excess absorption for galaxies at similar redshifts.

The shape of the attenuation curve potentially has significant impacts on the SEDs of galaxies and thus derived fundamental quantities such as SFR, as mentioned above. In particular, the potential variations in the FUV slope of the attenuation curves also imply that the attenuation of the Lyman continuum photons may even largely vary from galaxy to galaxy. This, for instance, may have a significant impact in estimating SFR from the $H\alpha$ flux (or whatever the recombination line flux), since dust absorption in the Lyman continuum reduces the number of produced $H\alpha$ photons (Puglisi et al. 2016). Variations in the Lyman continuum absorption may also affect the measurements of the ionizing photon escape fraction of galaxies. Accurate determination of the overall shape of the attenuation curve is thus essential for better understanding the evolution of galaxies.

## 6. Summary

We have investigated the strength of the $2175 \AA$ UV bump feature in the attenuation curves of a sample of 505 star-forming galaxies at $1.3 \le z \le 1.8$ in the zCOSMOS-deep survey. Approximately 30% of the galaxies exhibit a robust signature of the UV bump ($\gamma_{34} < -2$) in their individual spectra (Section 4.1). Significant intrinsic scatter in the observed $\gamma_{34}$ at a given UV slope clearly indicates the presence of a real diversity in the bump strength in the attenuation curves across the sample (see Section 4.1).

To increase the S/N, we have also measured the UV bump profiles in stacked spectra representing the whole sample, in subsamples of galaxies selected in ($\gamma_{34}, \beta_b$) space, and in stellar mass and sSFR. The attenuation profiles are all well described by Drude profiles with a center wavelength of $2175 \AA$ but varying amplitudes. The derived bump amplitudes vary across the range $B_k \approx 0.2-0.8$ (Section 4.2) with an inverse correlation with $\gamma_{34}$. Using the stacks in $M_*$ bins, we found a tight positive $M_*-B_k$ correlation across $9.8 \lesssim \log M_*/M_\odot \lesssim 10.5$ along the MS, though the correlation may not hold at the low- and high-mass ends (Figure 15).

Binning the sample in the $M_*-\text{sSFR}$ plane, we found that there is a strong negative trend between $B_k$ and sSFR at fixed $M_*$ while $B_k$ increases moderately with $M_*$ at fixed sSFR (Section 4.3). This correlation with sSFR in the high-redshift sample is strikingly counter to the observed absence of the UV bump in the average attenuation curve of local starburst galaxies, since these local sources actually have lower sSFRs than our high-$z$ MS sample. We found, however, that the two samples empirically come into much better agreement if we plot the bump strength against the $\Delta \log$ sSFR relative to the evolving MS at the appropriate epoch, rather than against the sSFR itself (see Figure 18).

We have interpreted these findings in terms of the recent SFH of the galaxies, especially considering the changes in the underlying grain properties. The bump strength is determined by the balance between the destruction and production of the bump carriers. The former may be accelerated by a higher frequency of SNe and/or more intense radiation fields in galaxies with higher instantaneous sSFRs, whereas the latter only reflects the SFH 1 Gyr before, but not the current SFR, for the time delay in the appearance of carbon-rich AGB stars (which are here assumed to be the main suppliers of precursor
carbonaceous grains of the bump carriers) from the onset of star formation. This suggests that the bump strength should be largely determined by the ratio of the current SFR measured on order $10^{7-8}$ yr timescales to that of order 1 Gyr ago, which we denote as SF$_{79}$ (Section 5.2).

We therefore explored how the SF$_{79}$ would be expected to vary with the sSFR of a galaxy in two different regimes: (i) that of a quasi-constant sSFR and (ii) that of a short-lived rapid increase in the SFR. This reveals an interesting effect. High-redshift MS galaxies with quasi-constant sSFRs of order of sSFR ~ 1 Gyr$^{-1}$ (appropriate for the high-redshift sample at z ~ 2) and local starburst galaxies that are undergoing a short-lived sharp increase in SFR should both exhibit the same relations between SF$_{79}$ and ΔsSFR but quite different relations between SF$_{79}$ and sSFR.

The fact that the bump strength is observed to behave in this same way therefore adds weight to the idea that it is the variation in SFR over the last roughly 1 Gyr, as parameterized by SF$_{79}$, that is responsible for the observed variations in the bump strength through the creation and destruction of the carrier grains responsible for the bump.

A prediction of this no doubt oversimplified picture is that MS galaxies at low redshift should show a very much weaker trend of bump strength with Δ log sSFR than their high-redshift counterparts. This prediction, however, appears to be counter to recent results of Salim et al. (2018). This apparent conflict may indicate that the radiative transfer effects in different dust-to-star geometries and dust column densities play roles in shaping the attenuation curves.

In the future, next-generation multi-object spectrographs, such as VLT/MOONS and Subaru/Prime Focus Spectrograph, will provide us with rest-frame UV spectra of order $10^{5-6}$ galaxies at high redshifts. This will enable us to investigate the attenuation curves and the bump feature against various galaxy properties using both individual spectra and stacks and thus to understand much better the nature of dust through cosmic time.

We are grateful to Irene Shivaei for useful discussions and thank the anonymous referee for constructive comments and suggestions. This research is based on observations undertaken at the European Southern Observatory (ESO) Very Large Telescope (VLT) under Large Program 175.A-0839 and has been supported by the Swiss National Science Foundation (SNF). This work was supported in part by KAKENHI (JDS 26400221) through the Japan Society for the Promotion of Science (JSPS) and the World Premier International Research Center Initiative (WPI), MEXT, Japan. This work uses data collected at the Subaru telescope, which is operated by the National Astronomical Observatory of Japan.

## Appendix A
### Reanalysis with the UV+IR-based SFRs

In our main analysis, we used total SFRs that were estimated from the UV luminosity with appropriate correction for dust absorption (see Section 3.3). However, for a nonnegligible fraction of the galaxies, the dust correction becomes quite large (attenuation at rest frame 1600 Å, $A_{1600} \gtrsim 2.5$); thus, the dust-corrected UV luminosity may be uncertain. Another concern is that the rest-frame UV emission does not trace highly dust-obscured star formation (Puglisi et al. 2017); thus, dusty starburst galaxies could be included in our sample as normal

star-forming galaxies. Therefore, here we reestimate the total SFRs of our sample galaxies by incorporating the rest-frame far-IR-to-millimeter photometry and present the results from reanalysis using the new SFR estimates.

### A.1. SFR from IR Luminosity

We utilized a “super-deblended” far-IR-to-millimeter photometric catalog presented by Jin et al. (2018). This catalog contains point-spread function fitting photometry at fixed prior positions including 88,008 galaxies detected in VLA 1.4 GHz, 3 GHz, and/or MIPS 24 µm images. To derive the total IR luminosity, we use the available photometry in the five Herschel PACS (100 and 160 µm) and SPIRE (250, 350, and 500 µm) complemented with JCMT/SCUBA2 850 µm, ASTE/AzTEC 1.1 mm, and IRAM/ MAMBO 1.2 mm. We fit these photometric fluxes with a coupled modified blackbody plus mid-IR truncated power-law component using the prescription given in Casey (2012). The total IR luminosity, $L_{IR}$, is then taken from the rest frame 8–1000 µm and converted to the IR-based SFR by employing a relation in Madau & Dickinson (2014) converted to a Chabrier (2003) IMF:

$$SFR_{IR} (M_{\odot} \text{ yr}^{-1}) = 2.64 \times 10^{-44} L_{IR} (\text{erg s}^{-1}). \quad (A1)$$

The total SFR is then computed as $SFR_{UV+IR} = SFR_{UV} + SFR_{IR}$, where $SFR_{UV}$ is computed from the observed UV luminosity not corrected for dust absorption.

### A.2. The Bump Strength as a Function of $M_*$ and sSFR

Cross-matching our sample with the far-IR catalog, we found a counterpart for 429 among our sample of 505 galaxies. For 70 of these, we measured $L_{IR}$ at $S/N \ge 3.0$, ranging across $log L_{IR}/L_{\odot} \sim 11.5-12.5$. In Figure 19, we compare the new SFR$_{UV+IR}$ with the fiducial SFR$_{UV,corr}$ for these 70 galaxies. The data points are color-coded by the dust attenuation at rest frame 1600 Å, $A_{1600}$. Although these two SFRs are in broad agreement with each other, there is a substantial scatter. Particularly, there are some objects having relatively low $A_{1600}$ whose SFR$_{UV+IR}$ exceeds SFR$_{UV,corr}$ by $\gtrsim 0.4$ dex, suggesting that these galaxies contain star-forming regions that contribute largely to the total SFR but are heavily dust-obscured, and thus their total SFR is not fully recovered in SFR$_{UV,corr}$ even by applying dust absorption correction. In contrast, a couple of sources are located well below the one-to-one relation with relatively large $A_{1600}$, suggesting that their SFR$_{UV,corr}$ may be overestimated due to overcorrection for dust. This implies that the uncertainties in the total SFR$_{UV,corr}$ could, in general, be larger than the nominal error bars, especially in those with larger $A_{1600}$.

For the reanalysis, we replaced the fiducial SFR$_{UV,corr}$ with the total SFR$_{UV+IR}$ for these 70 objects. For the remainder, we adopted SFR$_{UV,corr}$ if they have $A_{1600} < 2.5$ mag (see Equation (5)) but exclude the other 74 sources with $A_{1600} \ge 2.5$ mag from the sample. The final sample here contains 431 galaxies. In the lower panel of Figure 19, we show the sample in the $M_*$ versus sSFR plane. After replacing SFR$_{UV,corr}$ with SFR$_{UV+IR}$, the sample remains largely consistent with the same MS, with only a handful of objects located well above the MS ($\Delta \log \text{sSFR} \gtrsim 0.5$ dex). These objects could be heavily obscured starburst galaxies. In this paper, we do not specifically treat this type of galaxy because of their minor contribution to the whole sample and thus to the conclusions.
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Upper panel: comparison between SFR<sub>UV+IR</sub> and SFR<sub>UV,corr</sub> for 70 galaxies in our sample for which L<sub>IR</sub> is measured at S/N ≥ 3.0. Each symbol is color-coded with A<sub>1600</sub>. The diagonal line indicates a one-to-one relation. Lower panel: sSFR as a function of M<sub>*</sub>. Small blue circles indicate galaxies with A<sub>1600</sub> ≤ 2.5 mag for which SFR<sub>UV,corr</sub> is adopted for reanalysis. Orange circles indicate those for which L<sub>IR</sub> is measured at S/N ≥ 3.0 and thus SFR is replaced with SFR<sub>UV+IR</sub>, while white circles indicate the original SFR<sub>UV,corr</sub> which is not used for reanalysis.</figcaption></figure>

Using this sample, we carried out the same analysis described in Section 4.3. An exception is that we stacked galaxies at each grid point in the M* versus sSFR plane within a radius of 0.2 dex in both axes, instead of 0.1 dex in the sSFR axis, because the number density of the data points in the M* versus sSFR plane is reduced (but this does not change the results anyway). Figure 20 shows the results. The linear fits are given as

$$ \log B_k = -0.242 + 0.167 \log(M_*/10^{10} M_{\odot}) \\ - 0.299 \log(\text{sSFR/Gyr}^{-1}), \qquad (A2) $$

$$ \log B_A = -0.888 + 0.399 \log(M_*/10^{10} M_{\odot}) \\ - 0.211 \log(\text{sSFR/Gyr}^{-1}), \qquad (A3) $$

$$ \log E(B-V) = -0.628 + 0.231 \log(M_*/10^{10} M_{\odot}) \\ + 0.067 \log(\text{sSFR/Gyr}^{-1}). \qquad (A4) $$

The rms residuals in B<sub>k</sub>, B<sub>A</sub>, and E(B-V) are, respectively, 0.040, 0.010, and 0.017. The dependence of the bump

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Same as the right panels in Figure 16 but showing the results from the reanalysis using SFR<sub>UV+IR</sub> as described in the text. The ranges of the color bars are the same as in Figure 16.</figcaption></figure>

amplitude B<sub>k</sub> appears to be very similar to what is seen in Figure 16; B<sub>k</sub> is in a tight negative correlation with sSFR while moderately increasing with M*. This is also the case for B<sub>A</sub> and
the median E(B-V), as compared with the corresponding panels in Figure 16. The consistent result from the reanalysis makes our statements in this paper further robust.

## ORCID iDs

Daichi Kashino https://orcid.org/0000-0001-9044-1747

Simon J. Lilly https://orcid.org/0000-0002-6423-3597

John D. Silverman https://orcid.org/0000-0002-0000-6977

Alvio Renzini https://orcid.org/0000-0002-7093-7355

Emanuele Daddi https://orcid.org/0000-0002-3331-9590

Sandro Bardelli https://orcid.org/0000-0002-8900-0298

Olga Cucciati https://orcid.org/0000-0002-9336-7551

Jeyhan S. Kartaltepe https://orcid.org/0000-0001-9187-3605

Vincenzo Mainieri https://orcid.org/0000-0002-1047-9583

Roser Pelló https://orcid.org/0000-0003-0858-6109

David B. Sanders https://orcid.org/0000-0002-1233-9998

Elena Zucca https://orcid.org/0000-0002-5845-8132

## References

Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540

Asano, R. S., Takeuchi, T. T., Hirashita, H., & Nozawa, T. 2013, MNRAS, 432, 637

Asano, R. S., Takeuchi, T. T., Hirashita, H., & Nozawa, T. 2014, MNRAS, 440, 134

Battisti, A. J., Calzetti, D., & Chary, R. R. 2017, ApJ, 851, 90

Battisti, A. J., Cunha, E. d., Shivaei, I., & Calzetti, D. 2020, ApJ, 888, 108

Bekki, K., Hirashita, H., & Tsujimoto, T. 2015, ApJ, 810, 39

Bianchi, L., Clayton, G. C., Bohlin, R. C., Hutchings, J. B., & Massey, P. 1996, ApJ, 471, 203

Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936

Bradley, J., Dai, Z. R., Erni, R., et al. 2005, Sci, 307, 244

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buat, V., Giovannoli, E., Heinis, S., et al. 2011, A&A, 533, A93

Buat, V., Noll, S., Burgarella, D., et al. 2012, A&A, 545, A141

Calzetti, D. 1997, AJ, 113, 162

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582

Casey, C. M. 2012, MNRAS, 425, 3094

Chabrier, G. 2003, PASP, 115, 763

Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061

Cook, K. H., Aaronson, M., & Norris, J. 1986, ApJ, 305, 634

Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746

Draine, B. 1989, in IAU Symp. 135, Interstellar Dust, ed. L. J. Allamandola & A. G. G. M. Tielens (Dordrecht: Kluwer), 313

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Dwek, E. 1998, ApJ, 501, 643

Dwek, E., & Scalo, J. M. 1980, ApJ, 239, 193

Engelbracht, C. W., Gordon, K. D., Rieke, G. H., et al. 2005, ApJL, 628, L29

Fitzpatrick, E. L. 1986, AJ, 92, 1068

Fitzpatrick, E. L., & Massa, D. 1986, ApJ, 307, 286

Fitzpatrick, E. L., & Massa, D. 1990, ApJS, 72, 163

Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320

Gall, C., Hjorth, J., & Andersen, A. C. 2011, A&ARv, 19, 43

Galliano, F., Dwek, E., & Chanial, P. 2008, ApJ, 672, 214

Gilra, D. P. 1971, Natur, 229, 237

Gordon, K. D., Calzetti, D., & Witt, A. N. 1997, ApJ, 487, 625

Gordon, K. D., & Clayton, G. C. 1998, ApJ, 500, 816

Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279

Groenewegen, M. A. T. 1999, in IAU Symp. 191, Asymptotic Giant Branch Stars, ed. T. Le Bertre, A. Lebre, & C. Waelkens (San Francisco, CA: ASP), 535

Hirashita, H. 2015, MNRAS, 447, 2937

Hirashita, H., & Murga, M. S. 2020, MNRAS, 492, 3779

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236

Inoue, A. K. 2011, EP&S, 63, 1027

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Kashino, D., Silverman, J. D., Sanders, D., et al. 2017, ApJ, 835, 88

Kashino, D., Silverman, J. D., Sanders, D., et al. 2019, ApJS, 241, 10

Kriek, M., & Conroy, C. 2013, ApJL, 775, L16

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Le Fèvre, O., Saisse, M., Mancini, D., et al. 2003, Proc. SPIE, 4841, 1670

Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24

Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maraston, C., Pforr, J., Renzini, A., et al. 2010, MNRAS, 407, 830

Mathis, J. S. 1994, ApJ, 422, 176

McKee, C. 1989, in IAU Symp. 135, Interstellar Dust, ed. L. J. Allamandola & A. G. G. M. Tielens (Dordrecht: Kluwer), 431

Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64

Mouhcine, M., & Lanccon, A. 2003, MNRAS, 338, 572

Narayanan, D., Conroy, C., Davé, R., Johnson, B. D., & Popping, G. 2018, ApJ, 869, 70

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009a, A&A, 507, 1793

Noll, S., & Pierini, D. 2005, A&A, 444, 137

Noll, S., Pierini, D., Cimatti, A., et al. 2009b, A&A, 499, 69

Noll, S., Pierini, D., Pannella, M., & Savaglio, S. 2007, A&A, 472, 455

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845

Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141

Papoular, R. J., & Papoular, R. 2009, MNRAS, 394, 2175

Pforr, J., Maraston, C., & Tonini, C. 2012, MNRAS, 422, 3285

Prevot, M. L., Lequeux, J., Prevot, L., Maurice, E., & Rocca-Volmerange, B. 1984, A&A, 132, 389

Puglisi, A., Daddi, E., Renzini, A., et al. 2017, ApJL, 838, L18

Puglisi, A., Rodighiero, G., Franceschini, A., et al. 2016, A&A, 586, A83

Rau, S.-J., Hirashita, H., & Murga, M. 2019, MNRAS, 489, 5218

Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56

Renzini, A., & Peng, Y.-j. 2015, ApJL, 801, L29

Renzini, A., & Voli, M. 1981, A&A, 94, 175

Rix, S. A., Pettini, M., Leitherer, C., et al. 2004, ApJ, 615, 98

Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11

Salim, S., & Narayanan, D. 2020, ARA&A, 58, 529

Salmon, B., Papovich, C., Long, J., et al. 2016, ApJ, 827, 20

Salpeter, E. E. 1955, ApJ, 121, 161

Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99

Scodeggio, M., Franzetti, P., Garilli, B., et al. 2005, PASP, 117, 1284

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Scoville, N., Faisst, A., Capak, P., et al. 2015, ApJ, 800, 108

Seok, J. Y., Hirashita, H., & Asano, R. S. 2014, MNRAS, 439, 2186

Seon, K.-I., & Draine, B. T. 2016, ApJ, 833, 201

Shibuya, T., Ouchi, M., Kubo, M., & Harikane, Y. 2016, ApJ, 821, 72

Shivaei, I., Reddy, N., Rieke, G., et al. 2020, ApJ, 899, 117

Shivaei, I., Reddy, N. A., Shapley, A. E., et al. 2017, ApJ, 837, 157

Silverman, J. D., Kashino, D., Sanders, D., et al. 2015, ApJS, 220, 12

Sloan, G. C., Kraemer, K. E., Wood, P. R., et al. 2008, ApJ, 686, 1056

Sorrell, W. H. 1990, MNRAS, 243, 570

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Stecher, T. P. 1965, ApJ, 142, 1683

Steglich, M., Jäger, C., Rouillé, G., et al. 2010, ApJL, 712, L16

Steidel, C. C., Shapley, A. E., Pettini, M., et al. 2004, ApJ, 604, 534

Tress, M., Mármol-Queraltó, E., Ferreras, I., et al. 2018, MNRAS, 475, 2363

Vijh, U. P., Witt, A. N., & Gordon, K. D. 2003, ApJ, 587, 533

Wang, E., & Lilly, S. J. 2020, ApJ, 892, 87

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Witt, A. N., & Gordon, K. D. 1996, ApJ, 463, 681

Witt, A. N., & Gordon, K. D. 2000, ApJ, 528, 799

Yan, H., Lazarian, A., & Draine, B. T. 2004, ApJ, 616, 895

Zahid, H. J., Dima, G. I., Kudritzki, R.-P., et al. 2014, ApJ, 791, 130

Zhukovska, S., Gail, H. P., & Trieloff, M. 2008, A&A, 479, 453