R Gobat
email: rgobat@kias.re.kr

Emanuele Daddi

G Magdis

F Bournaud

M Sargent

M Martig

S Jin

A Finoguenov

M Bethermin

H S S Hwang

M B Éthermin

A Renzini

G W Wilson

I Aretxaga

M Yun

V Strazzullo

F Valentino

The

Universe, the processes responsible for the cessation of star formation were already very efficient. The termination of star formation in ETGs is usually attributed to the removal of cool gas reservoirs (e.g., by stellar or quasar feedback [START_REF] Di Matteo | Energy input from quasars regulates the growth and activity of black holes and their host galaxies[END_REF] ) and/or by the suppression of gas infall and cooling (e.g., by virial shocks or AGN feedback [START_REF] Birnboim | Virial shocks in galactic haloes?[END_REF][START_REF] Croton | The many lives of active galactic nuclei: cooling flows, black holes and the luminosities and colours of galaxies[END_REF] ). Alternatively, the growth of bulges and stellar spheroids is thought to stabilise gas reservoirs, making star formation inefficient compared to disk galaxies [START_REF] Martig | Morphological Quenching of Star Formation: Making Early-Type Galaxies Red[END_REF][START_REF] Tacchella | Evidence for mature bulges and an inside-out quenching phase 3 billion years after the Big Bang[END_REF] . If the latter plays an important role in galaxy quenching, we might expect substantial reservoirs of untapped gas to exist in galaxies that have recently turned quiescent. Detecting this residual gas at high redshift, close to the epoch of quenching for massive quiescent galaxies, is however very challenging 2 and all attempts have so far been unsuccessful.

Deep survey fields allow us to circumnavigate this issue. By combining data from large numbers of undetected quiescent galaxies, their average far-infrared (FIR) emission can be observed, which traces the cold dust present in their interstellar medium (ISM) and, indirectly, their SFR. We select 977 isolated, high-mass ( M = 1.1 × 10 11 M ) quiescent galaxies in the 2 deg 2 COSMOS field by combining two different photometric criteria (BzK 10 and UVJ [START_REF] Wuyts | What Do We Learn from IRAC Observations of Galaxies at 2 < z < 3.5?[END_REF] ) and keeping only objects that are individually undetected at observed mid-infrared (MIR) wavelengths (Methods). These criteria ensure that the sample contains only the least star-forming galaxies at z = 1.4 -2.5, with clear early-type morphologies (as implied by a high median S érsic index n ∼ 3.5; Methods). We extract cutouts centred at the position of each galaxy from the 24µm (MIR), 100 -500µm (FIR), 0.85 -1.1 mm (sub-millimetre; sub-mm), 10 and 20 cm (radio) observations of COS-MOS and perform a median stacking analysis at each wavelength. After correcting for the contribution from satellite galaxies and unassociated neighbours in the line of sight (Methods), we obtain a clean measure of the FIR and radio emission of the quiescent galaxies, with detections in all but two bands (100 and 160µm) ranging in significance from 3.5σ to 11.9σ (Figure 1).

We model both the FIR emission from cold dust [START_REF] Draine | Infrared Emission from Interstellar Dust. IV. The Silicate-Graphite-PAH Model in the Post-Spitzer Era[END_REF] and non-thermal radio emission from star formation, based on the FIR flux (Methods). While the models underpredict the observed radio fluxes by a factor ∼3 (corresponding to an excess of ∼ 5 × 10 22 W/Hz), they reproduce the FIR SED well, implying a luminosity of L IR = 2.9±0.9×10 10 L . Furthermore, sampling the peak (∼150 µm rest-frame) and Rayleigh-Jeans tail (> 400µm rest-frame; Figure 1) of the FIR emission allows for an accurate determination of the dust temperature, T dust ∼ 21 -25 K, and mass, M dust = 1 +0.6 -0.4 × 10 8 M . The temperature of the cold dust component, as suggested by the lack of detection at ∼ 100 -160 µm but not at ≥ 250 µm, is similarly up to 10 K lower than the main sequence of star formation (MS) at the same L IR [START_REF] Hwang | Evolution of dust temperature of galaxies through cosmic time as seen by Herschel[END_REF] , among the lowest values for star forming (SF) galaxies at any redshift but consistent with that of local ETGs [START_REF] Smith | The Herschel Reference Survey: Dust in Early-type Galaxies and across the Hubble Sequence[END_REF] . This low temperature also confirms that the measured Figure 1: Mid-infrared to radio SED of 24 µm-undetected pBzK+UVJp galaxies, from stacked cutout images (filled circles, the empty circles representing <3σ nondetections), with best-fit model from Draine & Li [START_REF] Draine | Infrared Emission from Interstellar Dust. IV. The Silicate-Graphite-PAH Model in the Post-Spitzer Era[END_REF] (solid red curve, dotted when extended into the radio regime; Methods). For comparison, the "best-fit" standard MS template [START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF] , which has hotter dust, is shown in blue, and the median best-fit dust component to the FIR SED of local dust-rich ETGs as a dashed green curve. The latter model has been redshifted to z = 1.76 and rescaled to the same stellar mass. All models have been broadened using the redshift distribution of the sample as kernel. The error bars on the MIR, FIR, and radio fluxes show the bootstrap uncertainties derived from the stack. FIR emission does not originate from star forming satellites, which would have T dust values closer to the MS (T dust,MS ∼ 30 K; Methods). These distant ETGs appear however to be much more dust rich than their local counterparts, with (M dust /M ) z=1.76 ∼ 8 × 10 -4 compared to (M dust /M ) z=0 ∼ 10 -6 -10 -5 locally 23,1 , a difference of 2 orders of magnitude. We convert the measured dust mass into a gas mass assuming a gas-phase metallicity close to the solar value and a metallicity-dependent gas-to-dust ratio (G/D) [START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF] , which yields M gas = 9.5 +5.6 -3.5 × 10 9 M (Methods). In nearby quiescent galaxies for which such measures exist, the relative content of molecular (H 2 ) to atomic (Hi) hydrogen gas ranges from M(H 2 )/M(Hi) ∼ 0.1 to 100, with molecular gas nevertheless dominating by a factor ∼5 on average [START_REF] Young | The ATLAS 3D project -IV. The molecular gas content of early-type galaxies[END_REF][START_REF] Young | The ATLAS 3D project -XXVII. Cold gas and the colours and ages of early-type galaxies[END_REF] . We thus assume that distant quiescent galaxies are similar to their local counterparts in this regard (Methods) and that M gas hereafter refers to molecular hydrogen rather than the combined neutral and H 2 gas mass.

The large amount of dust we detect, and the inferred gas mass, implies that high redshift ETGs are remarkably gas rich, with up to 11% of their baryonic mass being in the form of gas. As with the dust fraction, this represents an increase of more than two orders of magnitude with respect to the local Universe 1 (Figure 2), even when accounting for differences in stellar mass between local and high-redshift samples, and a cosmic evolution that vastly exceeds the similar increase in the gas fraction of Figure 2: Evolution of the molecular gas fraction M mol /M as a function of redshift for both quiescent and MS galaxies.The latter include low-redshift [START_REF] Leroy | The Star Formation Efficiency in Nearby Galaxies: Measuring Where Gas Forms Stars Effectively[END_REF][START_REF] Saintonge | an IRAM legacy survey of molecular gas in massive galaxies -I. Relations between H 2 , Hi, stellar content and structural properties[END_REF] , intermediate-redshift 16,17 , and z > 1 MS galaxies [START_REF] Daddi | Different Star Formation Laws for Disks Versus Starbursts at Low and High Redshifts[END_REF][START_REF] Tacconi | Phibss: Molecular Gas Content and Scaling Relations in z ∼ 1 -3 Massive, Mainsequence Star-forming Galaxies[END_REF][START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF] .

The distant ETG sample is shown by an orange dot, with error bars from the SED modelling (Methods), and the upper limit constraints from CO(2-1) observations of a massive ETG at a similar redshift 2 by a red arrow. The blue arrow marks the median upper limit on this ratio for ETGs in the ATLAS 3D survey, with the light blue shaded region showing the range of values for the detected subsample (22%) of ATLAS 3D . For comparison, local ETGs in the Herschel Reference Survey (HRS) are shown by magenta ellipses: open for gas masses from CO measurements [START_REF] Boselli | Cold gas properties of the Herschel Reference Survey. II. Molecular and total gas scaling relations[END_REF] , including upper limits, and filled for dust masses 1 converted using the same G/D, with the error bar showing the median absolute deviation of the sample. The black curve and grey strip show the evolution of an average MS galaxy [START_REF] Sargent | Regularity Underlying Complexity: A Redshift-independent Description of the Continuous Variation of Galaxy-scale Molecular Gas Properties in the Mass-star Formation Rate Plane[END_REF] with stellar mass 5 × 10 10 M , and a 0.2 dex scatter approximately corresponding to the 1σ scatter of the MS. The dotted line shows the relation 21 offset by a factor ∼ 6 and for a stellar mass of log M (z) + ∆M, where ∆M is the offset between the median mass of the total pBzk+UVJp sample and the M of the stellar mass function of quiescent galaxies at the median redshift of the sample. Finally, the green-yellow region shows gas consumption in the closed box case (no gas inflow or outflow), with a timescale t gas between 2 and 3 Gyr (lower and upper edge, respectively).

SF galaxies. In "normal" SF galaxies, the FIR and submm emission is a relatively pure tracer of star formation.

If the observed cold dust emission does indeed arise from star formation activity, converting [START_REF] Kennicutt | Star Formation in Galaxies Along the Hubble Sequence[END_REF] L IR into a SFR then yields 4.8 +1.8 -1.3 M yr -1 , in excellent agreement with the SFR of z ∼ 1.5 ETGs derived from optical spectroscopic diagnostics [START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] , which probe different timescales (a few tens of Myr for spectroscopic line emission, several times that for FIR emission). Given the large stellar mass of these galaxies, it implies that the distant ETGs have a sSFR that is a factor ∼ 30 lower than that of active galaxies on the MS at the same epoch (Figure 3). The presence of both such large gas reservoirs and low SFR would imply that the star formation efficiency (SFE) of these quiescent galaxies is 2 -3 times lower Figure 3: A): correlation between the molecular gas content and SFR of galaxies, including the same MS samples as in Figure 2. Our median measurement, shown as a orange filled circle with error bars (Methods), is offset by ∼ 3σ from the average SFR -M mol relation for MS galaxies [START_REF] Sargent | Regularity Underlying Complexity: A Redshift-independent Description of the Continuous Variation of Galaxy-scale Molecular Gas Properties in the Mass-star Formation Rate Plane[END_REF] (black dashed line). For comparison, the upper limit constraints from CO(2-1) observations of a massive ETG at a similar redshift 2 is shown by red arrows, local ETGs from the HRS 20,1 as blue crosses and arrows (in the case of upper limits on M mol ), and independent measurements from the local ATLAS 3D sample [START_REF] Young | The ATLAS 3D project -IV. The molecular gas content of early-type galaxies[END_REF][START_REF] Davis | The ATLAS 3D project -XXVIII. Dynamically driven star formation suppression in earlytype galaxies[END_REF] , when more than one dex below the MS in SFR, as light blue squares. As in Figure 2, the average SFR and dust-derived gas mass of HRS ETGs is indicated by a filled magenta ellipse, with error bars showing the median absolute deviation of the sample. The green-yellow region shows the same closedbox evolution from z = 1.76 to z = 0, as in Figure 2. B): SFR of star forming and quiescent galaxies as a fraction of the MS one [START_REF] Sargent | Regularity Underlying Complexity: A Redshift-independent Description of the Continuous Variation of Galaxy-scale Molecular Gas Properties in the Mass-star Formation Rate Plane[END_REF] , for the same data as in A).

than that of disk galaxies and comparable to that of local ETGs 28 (Figure 3). This value is conservative and thus could be even lower, if part of the FIR emission of these distant ETGs arises from, e.g., heating of cirrus dust by an evolved stellar population consistent with the median ETG age at this redshift. Whether massive, main sequence stars or post-main-sequence stars (i.e., ongoing of past star formation, respectively) are the dominant heating source of the dust in these ETGs then depends on its spatial distribution (concentrated in and around molecular clouds or diffuse), which is still unknown. The good agreement between optical and IR-based SFRs [START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] , though, suggests that residual star formation is an important contributing mechanism to the FIR emission. A low SFE is also consistent with the low dust temperature and small value of the dust-mass weighted luminosity U estimated from the IR SED, since the intensity of the radiation field can be related to the SFE and metallicity by the relation U ∝ SFE/Z [START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF] .

It also implies a longer gas consumption timescale, t gas = 2 -3 Gyr compared to 1 Gyr for MS galaxies.

A low efficiency in ETGs is predicted by high-resolution hydrodynamical simulations 8,29 (Methods), in which the SFE of a forming ETG, which starts off as a high-redshift disk galaxy, decreases by a factor 5 -10 once the gas drops below a critical value (Figure 4). In our simulations, of gas disks embedded in spheroids, this corresponds to a gas fraction of 20 -25%. This effect is linked to the presence of the spheroidal stellar component which, unlike a stellar disk, does not destabilise the gas disk. Conversely, in galaxies containing more than ∼ 25% of gas, regardless of their stellar structure, the gas remains unstable and capable of forming stars efficiently.

Our observations suggest that this process does indeed occur in high-redshift ETGs once the gas fraction has been reduced enough, and that the suppression of star formation thus does not require the complete removal of cool gas reservoirs. In principle, a passive evolution with no further accretion and a low SFE could then account for the low gas fractions of local ETGs, compared with their high-redshift counterparts, while only increasing their stellar mass by a small fraction between z ∼ 1.76 and z = 0 (Figure 2). However, gas in local ETGs is likely to have had a more complex history. Indeed, several channels exist through which these galaxies could have been replenished with gas. Given ETGs ages of < 3 Gyr inferred from spectroscopy [START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] , we can estimate that ∼ 4.5 × 10 10 M of metal-enriched gas were already returned to the ISM by stellar evolution, at a rate of 2.9 M yr -1 at z = 1.76, with another ∼ 10 10 M produced over the next 10 Gyr assuming passive evolution (Methods). Gas is also expected to be reaccreted through mergers with gas-rich satellite galaxies and through the cosmic web. This suggests that there must exist other processes to keep the cold gas fraction low in these ETGs for the rest of their cosmic evolution, either by efficiently removing the surplus gas or preventing it from cooling. There are three likely such "maintenance" mechanisms: first, the host halos of these distant ETGs are, at 10 13 M , already massive enough to support a hot atmosphere, as evidenced by X-ray stacking (Methods). This hot gas can efficiently prevent the infall of cold gas onto the central galaxy [START_REF] Birnboim | Virial shocks in galactic haloes?[END_REF] . The X-rays emitted by the halo can also remove, on short timescales and through sputtering, the dust produced by stellar evolution, thus contributing to the evolution of the dust content of ETGs between z > 1.4 and z = 0 (Methods). Figure 4: SFE as a function of gas fraction in highresolution hydrodynamic simulations. The black line shows the evolutionary track of an ETG with an embedded gas disk in a cosmological zoom-in simulation, with initial peaks of SFE corresponding to galaxy mergers, followed by gas consumption and a strong decrease of the SFE once the gas fraction gets lower than 20-25%. This effect, which has been dubbed "morphological quenching", does not occur in galaxies whose stellar structure remains disk-like [START_REF] Martig | Morphological Quenching of Star Formation: Making Early-Type Galaxies Red[END_REF] . The galaxy already has an ETG morphology at the beginning of the track (Methods) but only quenches ∼ 2 Gyr later when enough gas has been consumed. The red and blue symbols are from idealised simulations of galaxies of ETG and disk-dominated galaxies, respectively, where both gas structure and star formation are resolved at the parsec scale [START_REF] Martig | The ATLAS 3D project -XXII. Lowefficiency star formation in early-type galaxies: hydrodynamic models and observations[END_REF] (Methods). The error bars correspond to the r.m.s. fluctuation of the SFE between various snapshots of a given simulation. These confirm that the SFE drops below gas fractions of 20-25% only for ETGs, at a level consistent with our observations (filled green square with error bars). Local ETGs from ATLAS 3D and HRS have SFEs similar to the idealised simulations and high redshift sample, but with much lower gas fractions (3 × 10 -3 on average).

Second, the "radio-mode" feedback from AGN-driven jets (as tantalisingly suggested by the systematic radio excess, which would imply a duty cycle of nearly 1, as opposed to 1% locally; Figure 1 and Methods) can not only heat the halo gas but also regularly drive away a fraction of the cold gas reservoir [START_REF] Croton | The many lives of active galactic nuclei: cooling flows, black holes and the luminosities and colours of galaxies[END_REF] . Finally, gravitational interaction between satellites and the halo can also provide a source of heating for the gas [START_REF] Johansson | Gravitational Heating Helps Make Massive Galaxies Red and Dead[END_REF] .

The large amounts of dust we detect in our z > 1.4 sample, through its FIR emission, implies similarly large quantities of molecular hydrogen. This H 2 gas would thus appear to be abundant enough in these galaxies to be detectable with the Atacama Large Millimeter Array (ALMA) through sub-mm line tracers, although the confirmation of even a single galaxy can be expected to require several hours of observation. Nevertheless, our result indicates that semi-direct studies of gas in high-redshift quiescent galaxies are already feasible, although the large time expenditure required and stillunknown dependencies of the dust and gas fractions on stellar mass or environment could make their design challenging.

Methods 1 Conventions

Stellar masses and star formation rates (SFRs) derived from photometry assume a Salpeter [START_REF] Salpeter | The Luminosity Function and Stellar Evolution[END_REF] stellar initial mass function, in accordance with the literature on IMF variations in ETGs [START_REF] Grillo | On the initial mass function and tilt of the fundamental plane of massive early-type galaxies[END_REF][START_REF] Cappellari | Systematic variation of the stellar initial mass function in early-type galaxies[END_REF][START_REF] Conroy | The Stellar Initial Mass Function in Early-type Galaxies From Absorption Line Spectroscopy. II. Results[END_REF] . We assume for consistency the same IMF for galaxies on the main sequence of star formation (MS). For distances, cosmic times, and cosmological simulations, we assume a ΛCDM cosmology with H 0 = 70 km s -1 Mpc -1 , Ω M = 0.3, and Ω Λ = 0.7 [START_REF] Larson | Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Power Spectra and WMAP-derived Parameters[END_REF][START_REF]Planck 2015 results. XIII. Cosmological parameters[END_REF] .

Sample selection

We select galaxies from the COSMOS field [START_REF] Scoville | The Cosmic Evolution Survey (COSMOS): Overview[END_REF] using public UV-to-near infrared (NIR) photometry [START_REF] Mccracken | The COSMOS-WIRCam Near-Infrared Imaging Survey. I. BzK-Selected Passive and Star-Forming Galaxy Candidates at z 1.4[END_REF][START_REF] Muzzin | A Public K s -selected Catalog in the COSMOS/ULTRAVISTA Field: Photometry, Photometric Redshifts, and Stellar Population Parameters[END_REF] and a new point-source-matched catalogue (Jin et al., in prep.; 8 µJy r.m.s.) based on public Spitzer/MIPS 24 µm data [START_REF] Floc'h | Deep Spitzer 24 µm COSMOS Imaging. I. The Evolution of Luminous Dusty Galaxies -Confronting the Models[END_REF] . We use spectroscopic redshifts from zCOS-MOS [START_REF] Lilly | zCOSMOS: A Large VLT/VIMOS Redshift Survey Covering 0 < z < 3 in the COS-MOS Field[END_REF] when available and photometric redshifts [START_REF] Muzzin | A Public K s -selected Catalog in the COSMOS/ULTRAVISTA Field: Photometry, Photometric Redshifts, and Stellar Population Parameters[END_REF][START_REF] Strazzullo | Passive galaxies as tracers of cluster environments at z ∼ 2[END_REF] otherwise. Quiescent galaxies at 1.4 < z < 2.5 are first selected using the BzK criterion 10 , using photometry from the original 2010 catalogue [START_REF] Mccracken | The COSMOS-WIRCam Near-Infrared Imaging Survey. I. BzK-Selected Passive and Star-Forming Galaxy Candidates at z 1.4[END_REF] , for which colour selections have been extensively tested. The rest of the analysis makes use of more recent UltraVISTA photometry [START_REF] Muzzin | A Public K s -selected Catalog in the COSMOS/ULTRAVISTA Field: Photometry, Photometric Redshifts, and Stellar Population Parameters[END_REF] , however. In this first step, galaxies undetected at <3σ in the B band, implying a Bz colour compatible with the locus of quiescent galaxies, are also considered. Photometric redshifts calibrated for passive stellar populations [START_REF] Strazzullo | Passive galaxies as tracers of cluster environments at z ∼ 2[END_REF] are then used to compute rest-frame U-, V-and J-band magnitudes for the initial sample, so as to discard from the sample galaxies not also selected as quiescent by the high-redshift UVJ criterion [START_REF] Williams | Detection of Quiescent Galaxies in a Bicolor Sequence from z = 0 -2[END_REF] . As illustrated in Supplementary Figure 1, this double selection minimises the number of interlopers due to photometric uncertainties and highly-reddened star-forming galaxies. We then discard from the sample sources that are detected at ≥ 3σ significance in the 24 µm catalogue, about 14% of the sample, as it implies either considerable star formation (SF) or AGN activity. We then keep only galaxies with stellar masses higher than log M = 10.8 M , within uncertainties, so as to ensure that the stellar mass distribution of this sample is similar to that of a spectroscopically observed subsample [START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] (Supplementary Figure 2) for which accurate redshifts and stellar population parameters (age, SFR, metallicity) were estimated. This mass limit also corresponds to twice the completeness limit of the catalogue for z > 1.5 quiescent galaxies [START_REF] Gobat | Satellite content and quenching of star formation in galaxy groups at z ∼ 1.8[END_REF] . We also ignore possible galaxy pairs, i.e., with concordant redshifts (within ∆z phot ≤ 0.2) and separated by less than 60" ( 30% of the sample), to simplify the modelling of the far-infrared (FIR) emission (Supplementary Section 3). The final sample contains 977 passive galaxies selected from a combined BzK+UVJ diagram, with a median redshift of z = 1.76 and a median stellar mass of M = 1.1 × 10 11 M .

Far-infrared stack

The COSMOS survey includes imaging data at 24 µm 42 (Spitzer/MIPS), 100 and 160 µm (Herschel/PACS, from the PEP survey [START_REF] Lutz | PACS Evolutionary Probe (PEP) -A Herschel key program[END_REF] ), 250, 350 and 500 µm (Herschel/SPIRE, from the HerMES survey [START_REF] Oliver | The Herschel Multi-tiered Extragalactic Survey: HerMES[END_REF] ), 850 µm [START_REF] Geach | The SCUBA-2 Cosmology Legacy Survey: 850 µm maps, catalogues and number counts[END_REF] (JCMT/SCUBA2), 1.1 mm 50 (ASTE/AzTEC), 3 GHz 52 and 1.4 GHz 51 (VLA). From these data we extract cutouts centred on each source in the sample with sizes 5 times larger than the beam FWHM. These are then combined to create a median 2D image in each band (Supplementary Figure 3), using a thousand bootstrap resamplings of the data with sizes of half the initial sample to estimate the variance of the sample. At each bootstrap iteration, we estimate the flux of the central ETG by assuming that the resulting median image is produced by a combination of point-like emission, contribution from star-forming satellites, and a background term. The FIR emission from the ETGs' satellite halos is estimated in the following way [START_REF] Gobat | Satellite content and quenching of star formation in galaxy groups at z ∼ 1.8[END_REF] : we select, within 120" of each ETG (hereafter, "central"), star-forming galaxies that satisfy the high-redshift UVJ criterion [START_REF] Williams | Detection of Quiescent Galaxies in a Bicolor Sequence from z = 0 -2[END_REF] and with z L68 < z cen < z H68 , where z cen is the redshift of the central and z L68 (z H68 ) the lower (upper) 68% confidence limit to the photometric redshift of the putative satellites. Subtracting the background distribution [START_REF] Gobat | Satellite content and quenching of star formation in galaxy groups at z ∼ 1.8[END_REF] yields an excess of 1234 star-forming satellites, with an additional ∼ 400 quiescent ones, corresponding to an average of 1 -2 satellites per central and a stellar mass distribution consistent with field BzK -selected galaxies. As shown in Supplementary Figure 4, the resulting overdensity of star-forming satellite galaxies (hereafter, "satellite halo") becomes clear when averaging over the whole sample and reaches the level of the background at 30" from the central. We estimate the SFR of each satellite thus selected by modelling its rest-frame UV SED with constant star formation stellar population models [START_REF] Bruzual | Stellar population synthesis at the resolution of 2003[END_REF] , with and without dust extinction [START_REF] Calzetti | The Dust Content and Opacity of Actively Star-forming Galaxies[END_REF] . The obscured SFR of satellites is then given by the difference between the dustcorrected and uncorrected UV SFRs. We then divide the field around the centrals in annuli of width ∆ log r = 0.1", starting at r = 1 from the central (i.e., twice the seeing and more than 6 times the galaxies' effective radius; see Supplementary Section 4), and sum in each annulus the obscured SFR of satellites to estimate the radial SFR surface density profile of the satellite halo. We model this profile using a combination of a β-function and a multiplicative, two-parameter "quenching" term of the form min(r/r q , 1) αq (Supplementary Figure 4,B), where r is the distance to the central and r q and α q are free parameters, and extrapolate it to the central 2". As the formal error on these parameters contributes negligibly to the uncertainty budget, we fix them to their best-fit values. As the obscured SFR is linearly proportional to the infrared luminosity, the SFR surface density profile is then used as a model of the FIR emission of satellites. Quiescent satellites are not included in this estimate since their average number and stellar mass densities are low, and amount to (6 ± 0.7) × 10 9 M when integrated up to the virial radius. If the mass fraction and temperature of the dust in quiescent satellites are the same than for the quiescent centrals (see below), their contribution to the FIR emission is at most 5%, i.e., somewhat smaller than the uncertainties. This value could be even lower if, on the other hand, the quenching of SF in these satellites was mediated by mechanisms that removed the interstellar medium (ISM) at a greater rate than for the centrals. When the FWHM of the FIR instrumental beam is smaller than the extension of the satellite halo (∼20"; this applies to the MIPS, PACS, 250 µm SPIRE, SCUBA-2, AzTEC, and VLA data), the FIR flux of the central in each band is estimated from the decomposition of the observed signal into a point-source and a 2D extended component based on the SFR profile of the satellite halo, convolved with the instrumental beam. On the other hand, the beam size of SPIRE at 350 and 500 µm is larger than the extent of the satellite halo, precluding a decomposition of the stacked cutout image into two components. In this case, we only fit the central part of the cutout (within one beam FWHM) with a single point-source and correct the estimate for the fraction of flux in the beam arising from the satellite halo assuming a main-sequence FIR SED [START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF] for the latter. We note that, in the bands where the two-component decomposition is possible, the two methods yield consistent flux estimates. In both cases the flux of the central galaxy is essentially measured through a point source fit to the unconvolved image. We then take the variance of bootstrap iterations as uncertainty on the flux. Finally, the measured fluxes are corrected for clustering of randomly distributed sources, based on simulations performed on the FIR COSMOS data [START_REF] B Éthermin | Evolution of the dust emission of massive galaxies up to z = 4 and constraints on their dominant mode of star formation[END_REF] . The clustering-corrected fluxes and satellite corrections are given in Supplementary Table 1. These, and the total infrared luminosity L IR we find, are a factor ∼ 5 lower than those inferred recently for quiescent galaxies at similar redshift and stellar mass [START_REF] Man | Confirming the Existence of a Quiescent Galaxy Population out to z = 3: A Stacking Analysis of Mid-, Far-Infrared and Radio Data[END_REF] , but without probing the Rayleigh-Jeans tail of the cold dust emission. Likewise the recovered 1.1 mm is a factor ∼4 lower than the contribution of quiescent galaxies to the cosmic infrared background [START_REF] Viero | HerMES: The Contribution to the Cosmic Infrared Background from Galaxies Selected by Mass and Redshift[END_REF] . However, both analyses use a different selection criterion, which, along with our strict treatment of satellite emission, could account for this difference.

The recovered MIR and FIR SED, including the <3σ (non-) detections at 100 and 160 µm, is then modelled with a set of templates of dust grain emission 12 (DL07), assuming q PAH = 3.19% and U max = 10 6 (respectively, the fraction of dust mass in polycyclic aromatic hydrocarbon grains and maximum starlight intensity seen by the dust, relative to the local interstellar radiation field) but with varying γ and U min (respectively, fraction of dust enclosed in photodissociation regions and minimum starlight intensity), to estimate the total infrared luminosity L IR and dust mass M dust . On the other hand, the dust temperature T dust is determined using a simple modified blackbody model (MBB) with a varying effective emissivity β, which yields β = 1.8 ± 0.3. The derived pa-rameter estimates are given in Supplementary Table 2. We note that M dust estimates could vary substantially between a simple single-temperature MBB fit to the data and that of that of a full dust model like the one provided by DL07. However, when the adopted absorption cross section of the dust distribution (κ abs ) used in the MBB model is consistent both in the normalisation and in the spectral index with that of the full dust model (DL07), the two techniques yield comparable results [START_REF] Bianchi | Vindicating single-T modified blackbody fits to Herschel SEDs[END_REF][START_REF] Berta | Measures of galaxy dust and gas mass with Herschel photometry and prospects for ALMA[END_REF] .

The major source of uncertainty in the derivation of M dust is the appropriate value of κ abs , which remains unknown. However, the prescription used in the DL07 model provides M dust (and therefore M gas ) estimates that are consistent with independent gas measurements of CO excitation in star-forming galaxies at various redshifts [START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF][START_REF] Genzel | Combined CO and Dust Scaling Relations of Depletion Time and Molecular Gas Fractions with Cosmic Time, Specific Star-formation Rate, and Stellar Mass[END_REF] . Furthermore, the DL07 model has been used to derive M gas estimates for the comparison sample of local quiescent galaxies. Thus, our choice of the DL07 model for our analysis, even if the derived absolute value of M dust remains uncertain, facilitates the comparison between various local and high-z samples and avoids any model dependent systematics.

Finally, we also add to the best-fit DL07 model a power-law radio slope with spectral index α = 0.8 and a normalisation given by the FIR-radio correlation [START_REF] Magnelli | The far-infrared/radio correlation and radio spectral index of galaxies in the SFR-M plane up to z ∼ 2[END_REF] . The radio fluxes are however not included in the fit and both the data point and model are shown for reference only. We note that there appears to be a significant excess in the median flux at 1.4 GHz and 3 GHz, compared to the model, corresponding to ∼ 5 × 10 22 W/Hz. On the one hand, this is suggestive of the presence in our sample of low-excitation AGNs [START_REF] Best | On the fundamental dichotomy in the local radio-AGN population: accretion, evolution and host galaxy properties[END_REF] and might be a sign of widespread, persistent AGN feedback. In this case the duty cycle should be f BH 0.5, since the excess is detected in a median stack. A more accurate value can be estimated under the simple assumption that the AGN activity is randomly distributed, over the time interval (t z=1.4 -t z=2.5 ) probed by our sample, in single bursts of varying duration ∆t b and constant luminosity set by the low-redshift radio luminosity functions [START_REF] Best | On the fundamental dichotomy in the local radio-AGN population: accretion, evolution and host galaxy properties[END_REF] . For each ∆t b the radio luminosity is estimated through Monte Carlo simulations by sampling, for each galaxy, the single burst model at its observed redshift. Varying ∆t b to match the observed excess luminosity then yields f BH = ∆t b /(t z=1.4 -t z=2.5 ) ∼ 0.66. This would imply that the radio AGN is almost always "on" in these galaxies, as opposed to the local Universe where f BH < 0.01 [START_REF] Ciotti | Cooling Flows and Quasars. II. Detailed Models of Feedback-modulated Accretion Flows[END_REF] . The radio spectral index is ∼ -0.82, which is very close to the typical value for SF emission but can also be produced by AGN-related processes [START_REF] O'dea | The Compact Steep-Spectrum and Gigahertz Peaked-Spectrum Radio Sources, Pub. Astron[END_REF][START_REF] Richards | Spectral Energy Distributions and Multiwavelength Selection of Type 1 Quasars[END_REF] . On the other hand, the normalisation of the FIR-radio correlation (FRC) has a dispersion of a factor ∼ 2 [START_REF] Yun | Radio Properties of Infrared-selected Galaxies in the IRAS 2 Jy Sample[END_REF] , potentially lowering the significance of this excess, and in local ETGs the 1.4 GHz flux appears to be weakly correlated with FIR luminosity [START_REF] Nyland | Star formation in nearby early-type galaxies: the radio continuum perspective[END_REF] . However, for an entire galaxy population to deviate from the median relation would imply that the slope and normalisation of the FRC vary substantially with, e.g., star formation efficiency (SFE) or specific SFR (sSFR), which is not substantiated by recent high-redshift studies [START_REF] Ibar | Deep multi-frequency radio imaging in the Lockman Hole using the GMRT and VLA -I. The nature of the sub-mJy radio population[END_REF][START_REF] Thomson | An ALMA survey of submillimetre galaxies in the Extended Chandra Deep Field South: radio properties and the farinfrared/radio correlation[END_REF][START_REF] Magnelli | The far-infrared/radio correlation and radio spectral index of galaxies in the SFR-M plane up to z ∼ 2[END_REF] .

Converting M dust into M gas

While M dust can be obtained directly from the FIR emission, the molecular gas mass M gas is a more important quantity in the context of galaxy evolution. Since molecular hydrogen is essentially invisible, M gas is usually extrapolated from other observables, either from sub-mm molecular (e.g., CO) line emission or M dust , using conversion factors that are somewhat uncertain. In the case of M dust , a gas-to-dust ratio (G/D) is used, which is related to the gas-phase metallicity and to processes of dust creation and destruction. We use here G/D ∼ 95, based on the parameterisation of the G/D as a function of metallicity from Magdis et al. 2012 [START_REF] Magdis | The Evolving Interstellar Medium of Star-forming Galaxies since z = 2 as Probed by Their Infrared Spectral Energy Distributions[END_REF] and the assumption of solar metallicity for the gas in the ETGs, which is supported by two lines of indirect evidence: first, the average gas metallicity of star-forming galaxies with the same sSFR is close to (or slightly above) the solar value, up to z ∼ 1.7 [START_REF] Mannucci | A fundamental relation between mass, star formation rate and metallicity in local and highredshift galaxies[END_REF][START_REF] Kashino | The FMOS-COSMOS Survey of Star-forming Galaxies at z ∼ 1.6. IV. Excitation State and Chemical Enrichment of the Interstellar Medium[END_REF] ; second, a spectroscopic, rest-frame UV measurement of the metallicity of the last batch of stars in a subset of our ETG sample also yields ∼ Z 27 , the stellar metallicity being expected to be slightly lower, but otherwise comparable, to the gas one [START_REF] Halliday | GMASS ultradeep spectroscopy of galaxies at z ∼ 2. I. The stellar metallicity[END_REF] .

The picture is however complicated by processes which can produce dust, such as supernovae, asymptotic giant-branch stars, or accretion of gas-phase metals in the ISM, and those that can destroy it (e.g., X-ray sputtering [START_REF] Arimoto | The Iron Discrepancy in Elliptical Galaxies after ASCA[END_REF] ). While the bulk of cold dust in local ETGs does not appear to originate from stellar evolution [START_REF] Smith | The Herschel Reference Survey: Dust in Early-type Galaxies and across the Hubble Sequence[END_REF][START_REF] Agius | GAMA/H-ATLAS: linking the properties of submm detected and undetected early-type galaxies -I. z ≤ 0.06 sample[END_REF] , the G/D of ETGs might nonetheless not follow the same relations as MS galaxies. On the other hand, in local ETGs which have detected H 2 gas (i.e., the gas-rich population) the median G/D is somewhat higher (∼ 400), while in the undetected ones the median upper-limit on the G/D is slightly lower (∼ 80) than our adopted value [START_REF] Lianou | The dustier early-type galaxies deviate from latetype galaxies' scaling relations[END_REF] . Overall, the global G/D of local ETGs appears to follow a similar relation to that of star-forming galaxies [START_REF] Leeuw | Spatially Resolved Imaging at 350 µm of Cold Dust in Nearby Elliptical Galaxies[END_REF][START_REF] Smith | The Herschel Reference Survey: Dust in Early-type Galaxies and across the Hubble Sequence[END_REF][START_REF] Lianou | The dustier early-type galaxies deviate from latetype galaxies' scaling relations[END_REF] .

Finally, although knowing the exact state of the gas in our distant ETGs is not necessary for the purposes of this work, we make the simplification that all the hydrogen gas is in molecular (H 2 ) form, rather than atomic (Hi) or a mix of both. We base this assumption on the high average H 2 /Hi ∼ 2 -5 ratio seen in local ETGs 76,24,25,77 . Moreover, the H 2 /Hi ratio is expected to increase with redshift [START_REF] Lagos | Cosmic evolution of the atomic and molecular gas contents of galaxies[END_REF][START_REF] Lagos | Which galaxies dominate the neutral gas content of the Universe?[END_REF] and gas column density [START_REF] Sternberg | Hi-to-H 2 Transitions and Hi Column Densities in Galaxy Star-forming Regions[END_REF] . This ratio appears nevertheless to vary substantially in the local population [START_REF] Welch | The cool interstellar medium in elliptical galaxies. II. Gas content in the volume-limited sample and results from the combined elliptical and lenticular surveys[END_REF] , possibly as a result of differential evolution, with lenticular (S0) galaxies containing higher amounts of H 2 [START_REF] Sage | The cool ISM in elliptical galaxies. I. A survey of molecular gas[END_REF] . We note that, while high-redshift quiescent galaxies tend to have a more disky morphology (i.e., closer to S0s) than local ETGs [START_REF] Chevance | On the Shapes and Structures of High-redshift Compact Galaxies[END_REF][START_REF] Krogager | A Spectroscopic Sample of Massive, Quiescent z ∼ 2 Galaxies: Implications for the Evolution of the Mass-Size Relation[END_REF][START_REF] Bruce | The bulge-disc decomposed evolution of massive galaxies at 1 < z < 3 in CANDELS[END_REF] (see also Supplementary Section 4), the cold gas in at least some local S0s might have a different origin than that of high-redshift ETGs [START_REF] Davis | The ATLAS 3D project -X. On the origin of the molecular and ionized gas in early-type galaxies[END_REF][START_REF] Katkov | Decoupled gas kinematics in isolated S0 galaxies[END_REF] .

Morphologies

There exist several publicly available high-resolution imaging datasets on the COSMOS field, taken with the Hubble Space Telescope (HST). The most extensive was obtained with the Advanced Camera for Surveys and covers 1.7 degrees at 814 nm [START_REF] Koekemoer | The COSMOS Survey: Hubble Space Telescope Advanced Camera for Survey Observations and Data Processing[END_REF] , corresponding to rest-frame UV wavelengths at z > 1.4. This implies that HST /ACS imaging might, for this sample, trace the younger stars produced by the last significant SF episode more than the total stellar mass distribution. Moreover, the existence of colour gradients in high-redshift ETGs is well attested [START_REF] Guo | Color and Stellar Population Gradients of Passively Evolving Galaxies at z ∼ 2 from HST/WFC3 Deep Imaging in the Hubble Ultra Deep Field[END_REF][START_REF] Gargiulo | Spatially resolved colours and stellar population properties in early-type galaxies at z ∼ 1.5[END_REF][START_REF] Chan | Sizes, colours gradients and resolved stellar mass distributions for the massive cluster galaxies in XMMUJ2235-2557 at z = 1.39[END_REF] , which implies that their morphological parameters depend on the wavelength at which they are measured. Here we therefore use the much smaller CANDELS survey [START_REF] Koekemoer | CANDELS: The Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey -The Hubble Space Telescope Observations[END_REF] , which benefits from NIR (i.e., restframe optical) data at 1.2 and 1.6 µm. We select the 46 ETGs in our sample whose position overlaps with the HST /WFC3 mosaics of the CANDELS survey (a random selection of which is shown in Supplementary Figure 5), as well as 12 medium luminosity stars (14 < V < 16) in the same region. We extract cutouts around the targets from both the F125W and F160W images and combine them using the same method as in Supplementary Section 3, to generate median NIR images of the galaxies and point sources. These are then used in conjunction with the galfit code [START_REF] Peng | Detailed Decomposition of Galaxy Images. II. Beyond Axisymmetric Models[END_REF] and fitted with a S érsic profile. This yields a S érsic index of n = 3.5±0.1 and effective radius of r e = 1.7 ± 0.1 kpc, consistent with spheroidal morphology and on the stellar mass-size relation for earlytype galaxies at this redshift [START_REF] Mancini | High-redshift elliptical galaxies: are they (all) really compact ?[END_REF][START_REF] Van Der Wel | 3d-HST+CANDELS: The Evolution of the Galaxy Size-Mass Distribution since z = 3[END_REF] . The results of this morphological modelling are shown in Supplementary Figure 6. Furthermore, 40 of these galaxies have individually measured morphologies in the CANDELS catalogue 95 , with median values of n ∼ 3 and r e ∼ 1.3 kpc, respectively. Similar results and high S érsic indices are also found for the full sample; more details will be presented in a future publication. On the other hand, we note that the in 3 GHz median image the ETG stack appears less concentrated, with r e ∼ 3.4 kpc for a fit with a fixed n ∼ 1 (exponential) profile. The signal-to-noise ratio of the stacked image is however not sufficient for a two-parameter fit. While n ∼ 1 yields a good fit, it is not clear whether this reflects the structure of the galaxies (or parts thereof) or stems from, e.g., the averaging over the whole sample of randomly oriented jets.

Hydrodynamical simulations

The cosmological simulation used in this work (Figure 4) is performed in a zoom box focusing on a massive galaxy at the centre of its dark matter halo. The size of the box corresponds to twice the virial radius of the dark halo, without boundary conditions for mass infall. The latter is based on a large-scale cosmological simulation. The simulation of a forming ETG analysed here is from Martig et al. 2009 [START_REF] Martig | Morphological Quenching of Star Formation: Making Early-Type Galaxies Red[END_REF] , with a complete sample presented in Martig et al. 2012 [START_REF] Martig | A Diversity of Progenitors and Histories for Isolated Spiral Galaxies[END_REF] . The evolution of the morphology of both its gas and star components is shown in Supplementary Figure 7. The galaxy starts as disk dominated at high redshift (z > 1.5), then undergoes a series of mergers with mass ratios ranging from 10:1 to 3:1. These mergers have two effects: they convert the stellar structure into a spheroid and compress the gas, increasing the SF and gas consumption rates. Over the 3 Gyr long track used in Figure 4, the stellar morphology is consistently spheroidal, although sometimes disturbed by mergers in their early stages. Over this period, the stellar mass increases from ∼ 6 × 10 10 M to 1.2 × 10 11 M , consistent with our observed sample. About 40% of this increase comes from SF internal to the studied galaxies, while the rest is accreted during mergers with other, smaller galaxies.

The idealised simulations are performed with the Adaptive Mesh Refinement code RAMSES [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF] . They are similar to those presented in Martig et al. 2013 [START_REF] Martig | The ATLAS 3D project -XXII. Lowefficiency star formation in early-type galaxies: hydrodynamic models and observations[END_REF] , but have double the resolution, reaching 3 pc in the densest gas regions 2 . The gas fraction coverage has been extended with new simulations, fully identical to the previous ones for all other parameters. These idealised simulations do not follow the formation of a given galaxy in a cosmological context. Instead, the stellar and gas mass, as well as their spatial distributions, are arbitrarily chosen to be representative of a certain type of galaxy. The simulations reach a very high resolution in order to accurately describe the fragmentation of gas in the dense clouds where stars form [START_REF] Martig | The ATLAS 3D project -XXII. Lowefficiency star formation in early-type galaxies: hydrodynamic models and observations[END_REF] . They also allow for direct comparisons of the SFE between disk-dominated galaxies and ETGs, with all other parameters being identical. The SFE is determined by the gas surface density, principally, and the rotation curve of the galaxy where the gas is located. In our simulated high-redshift galaxies the critical value is 70 -100 M pc -2 at the half-mass radius of the gas. This then corresponds to f gas ∼ 25% at high redshift, compared to f gas = 2 -3% at z ∼ 0, due to the expected ∼ 3 -4-fold increase in the size of gas disks (i.e., a factor 9 -16 in surface) in high-redshift ETGs compared to the compact ones seen in local ETGs [START_REF] Crocker | Molecular gas and star formation in early-type galaxies[END_REF][START_REF] Davis | The ATLAS 3D project -X. On the origin of the molecular and ionized gas in early-type galaxies[END_REF][START_REF] Van Der Wel | 3d-HST+CANDELS: The Evolution of the Galaxy Size-Mass Distribution since z = 3[END_REF] . The error bars in Figure 4 correspond to the r.m.s. fluctuation of the SFE between various snapshots of a given simulation. The agreement between these idealised simulations and the cosmological simulation confirms that the process highlighted in the cosmological simulation is not redshift-dependent: while the reduction of the SFE in a spheroid-dominated galaxy (once enough gas has been consumed) happens at about redshift z ∼ 1 in this simulation, it can just as well occur at somewhat higher redshift, as in our observations. Finally, we note that the idealised simulations do not include diffuse gas reservoirs in the galaxies' outskirts, yielding slightly lower gas masses and thus slightly higher values of the SFE than the cosmological simulation.

Gas and dust return rate

The amount of gas returned to the ISM by a stellar population depends on both its IMF and the star formation history (SFH). It is generally among the parameters estimated by stellar population synthesis models for single age populations [START_REF] Bruzual | Stellar population synthesis at the resolution of 2003[END_REF][START_REF] Maraston | Evolutionary population synthesis: models, analysis of the ingredients and application to high-z galaxies[END_REF] (SSP). The typical mass-weighted age of BzK -selected quiescent galaxies inferred from spectroscopy is t SF = 2 -3 Gyr [START_REF] Onodera | The Ages, Metallicities, and Element Abundance Ratios of Massive Quenched Galaxies at z ≥ 1.6[END_REF][START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] . We assume a delayed, exponentially declining SFH, arbitrarily starting at z = 10, as the simple parametric form that is closest to the probable "true" SFH of these high-redshift ETGs (rising as the galaxy grows on or above the MS, then decreasing during quenching). In this case, the SFH that would reproduce the SFR inferred from the cold dust emission has t SF ∼ 2.24 Gyr. The amount of gas returned to the ISM over a certain time span can then be computed easily as the convolution of the SFH and mass loss function due to stellar death. For our chosen SFH and IMF, this would amount to ∼ 4.5 × 10 10 M , or 37 +2 -1 % of the observed stellar mass at z = 1.76, and, given the median stellar mass of ETGs in this sample, an average return rate of 1.3 +1. [START_REF] Birnboim | Virial shocks in galactic haloes?[END_REF] -1.3 M yr -1 over the last 500 Myr. Similarly, assuming dust yields for solar-metallicity asymptotic giant branch stars [START_REF] Ventura | Dust from asymptotic giant branch stars: relevant factors and modelling uncertainties[END_REF] , ∼ 5 × 10 8 M would have been returned to the ISM, corresponding to ∼ 1% of the returned gas mass. We note that single burst or truncated SFHs, also commonly used in spectral and photometric modelling, would yield similar total gas fractions but negligible return rates at z = 1.76.

The typical best-fit SFH we infer for these ETGs [START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] would produce a large amount of post-main-sequence stars at z ∼ 1.76, which could then contribute significantly to the radiation field to which the dust is exposed. Whether the FIR emission we see is reprocessed radiation from massive, young stars (i.e., SF) or evolved stars then depends on the geometric distribution of the dust, which is unknown. We can, however, estimate the amount of FIR emission produced by non-SF sources in a simple case using the class of delayed exponential SFHs mentioned above, however truncated for the last 500 Myr so as to have no ongoing SF. We then use these SFHs to produce composite stellar population models [START_REF] Bruzual | Stellar population synthesis at the resolution of 2003[END_REF] with which we model the UV-optical SEDs of our ETGs. We assume a standard dust extinction law 100 with a fixed steep slope of δ = -0.4, typical of low sSFR galaxies, and a UV bump amplitude linearly dependent on δ 101 .

In this case, we find that the median difference between the integrated flux of models before and after applying extinction, i.e., the maximum luminosity that can arise from evolved stellar populations absorbed by the dust, is L abs = (9 ± 5) × 10 9 L , or less than half of the measured L IR .

Finally, assuming a closed-box evolution to z = 0, an additional ∼ 10 10 M of gas and ∼ 6 × 10 8 M of dust should be produced by stellar evolution, vastly more than currently observed in local ETGs [START_REF] Smith | The Herschel Reference Survey: Dust in Early-type Galaxies and across the Hubble Sequence[END_REF][START_REF] Lianou | The dustier early-type galaxies deviate from latetype galaxies' scaling relations[END_REF] and highlighting the importance of dust and gas removal mechanisms [START_REF] Ciotti | Winds, outflows, and inflows in X-ray elliptical galaxies[END_REF][START_REF] Arimoto | The Iron Discrepancy in Elliptical Galaxies after ASCA[END_REF][START_REF] Ciotti | Cooling Flows and Quasars. II. Detailed Models of Feedback-modulated Accretion Flows[END_REF][START_REF] Croton | The many lives of active galactic nuclei: cooling flows, black holes and the luminosities and colours of galaxies[END_REF] in the evolution of the ISM of ETGs.

X-ray stacking

We use deep X-ray observations of the COSMOS field by the Chandra and XMM-Newton observatories [START_REF] Finoguenov | The XMM-Newton Wide-Field Survey in the COSMOS Field: Statistical Properties of Clusters of Galaxies[END_REF][START_REF] Elvis | The Chandra COSMOS Survey. I. Overview and Point Source Catalog[END_REF] to constrain the total mass of the halos hosting our sample ETGs. This procedure is similar to that described in B éthermin et al. 2014 [START_REF] B Éthermin | Clustering, host halos, and environment of z ∼ 2 galaxies as a function of their physical properties[END_REF] and Gobat et al. 2015 [START_REF] Gobat | Satellite content and quenching of star formation in galaxy groups at z ∼ 1.8[END_REF] , and we refer to these articles for additional information. There are 585 sample ETGs in the Chandra COSMOS Survey field and located in zones free from emission. After subtracting both background and point sources, we find a residual flux in the 0.5-2 keV band of (2.3 ± 0.6) × 10 -17 erg cm -2 s -1 , corresponding to a luminosity of L X = (5.1 ± 1.3) × 10 42 erg cm -2 s -1 and a mass of M 200c = (1.2 ± 0.2) × 10 13 M .

2 (

 2 time since last major merger tM in Gyr)

Acknowledgements

The authors thank S. Lianou for providing models of dust emission in local ETGs and V. Smolcić for the 3 GHz radio data. S.J. acknowledges China Scholarship Council funding. The new simulations presented in this work were performed on GENCI resources (allocations 2016-04-2019 and 2017-04-2192).

Data Availability

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Author Contributions

R.G. and E.D. devised the project. R.G. analysed the data and wrote the manuscript. G.M. modelled the FIR emission. F.B. and M.M. carried out and analysed the simulations. M.S. and M.B. provided some of the theoretical framework. S.J. provided the MIR catalogue. A.F. analysed the X-ray observations. G.W.W., I.A., and M.Y. provided sub-mm data. H.S.H., A.R., V.S., and F.V. provided critical feedback that helped shape the manuscript. Supplementary Figure 1: Colour-colour diagrams showing the z > 1.4 ETG sample selected using the observed-frame BzK (A) and rest-frame UVJ (B) criteria. The grey distributions shows the position of sources selected using only the UVJ (A) or BzK (B) criteria, including, in the case of the latter, sources undetected at 5σ in the B or z band. Supplementary Figure 2: Distribution of photometric redshifts (A) and stellar masses of centrals (B) and satellites (C) derived from UltraVISTA photometry. In (A) and (B), the solid histogram shows the distribution of the spectroscopic subsample [START_REF] Gobat | In and out star formation in z ∼ 1.5 quiescent galaxies from rest-frame UV spectroscopy and the far-infrared[END_REF] , while the dotted one represents the distribution of sample ETGs in the CANDELS sub-field of COSMOS (see Supplementary Section 4). In (C), the solid histogram shows the distribution of UVJ-selected star-forming satellites, while the dotted one shows, for comparison, the distribution of all BzK-selected galaxies in the COSMOS field, rescaled to match the number of satellites. Supplementary

Supplementary