Early Growth of the Star Formation Rate Function in the Epoch of Reionization: An Approach with Rest-frame Optical Emissions

We present a star formation rate function (SFRF) at z ∼ 6 based on star formation rates (SFRs) derived by spectral energy distribution fitting on data from rest-frame UV to optical wavelengths of galaxies in the CANDELS GOODS-South and North fields. The resulting SFRF shows an excess compared to the previous estimations by using rest-frame UV luminosity functions (LFs) corrected for the dust attenuation and is comparable to that estimated from a far-infrared LF. This suggests that the number density of dust-obscured intensively star-forming galaxies at z ∼ 6 has been underestimated in the previous approach based only on rest-frame UV observations. We parameterize the SFRF using the Schechter function and obtain the best-fit parameter of the characteristic SFR (SFR*) when the faint-end slope and characteristic number density are fixed. The best-fit SFR* at z ∼ 6 is comparable to that at z ∼ 2, when the cosmic star formation activity reaches its peak. Together with SFRF estimations with a similar approach using rest-frame UV to optical data, the SFR* is roughly constant from z ∼ 2 to ∼6 and may decrease above z ∼ 6. Since the SFR* is sensitive to the high-SFR end of the SFRF, this evolution of SFR* suggests that the high-SFR end of the SFRF grows rapidly during the epoch of reionization and reaches a similar level observed at z ∼ 2.


INTRODUCTION
Galaxies evolve through converting their gas into stars, or star formation activities.The investigation on the star formation history across the cosmic time is thus essential to understand the galaxy formation and evolution.To quantitatively evaluate the cosmological evolution of star formation activities in galaxies, it is crucial to examine the star formation rate function (SFRF) of galaxies at a wide range of redshift.
The SFRF is defined as the number density of galaxies written as a function of the star formation rate (SFR) at a given redshift.Therefore, the SFRF describes the ongoing evolution of galaxies at a period of the universe, and the redshift evolution of the SFRF directly gives a picture of cosmological evolution history of galaxies.
Moreover, integrating the SFRF over the SFR provide the star formation rate density (SFRD) at the redshift, which is also useful and important to quantify the stellar mass assembly of galaxies across the cosmic time.
In the high-z (z ≳ 4) universe, the SFRF and/or SFRD is mainly investigated through the rest-frame UV luminosity functions (LFs) (see e.g., a review by Madau & Dickinson 2014).Since the rest-frame UV luminosity is one of the tracers of the SFR of galaxies, the restframe UV LF can be converted into the SFRF.However, the rest-frame UV light can be easily attenuated by dust, since it is necessary to correct for the loss of light due to dust attenuation.In this correction, the rest-frame UV spectral index β (f λ ∝ λ β ) is commonly used to estimate the amount of dust attenuation with the IRX-β relation (or more directly with A UV -β relation), which links the spectral slope β and the infrared excess defined as IRX = L TIR /L UV (see e.g., Meurer et al. 1999).
On the other hand, recent far-infrared (FIR) observations have been available even at z ≳ 4 (e.g., Rowan-Robinson et al. 2016;Koprowski et al. 2017).FIR observations can probe the dust-obscured star formation activities, and investigations of the SFRF and/or SFRD at z ≳ 4 with FIR observations have been carried out.However, part of previous studies on FIR observations at z ≳ 4 suggest that the contribution of dust-obscured star formation to the total is negligible and the estimation with rest-frame UV-based approach is well corrected for the dust attenuation (e.g., Koprowski et al. 2017), while others suggest the dust-obscured star formation dominates over the dust-unobscured star formation and the rest-frame UV-based approach underestimates the total star formation activities (e.g., Rowan-Robinson et al. 2016).Thus, on the contribution of the dust-obscured star formation at z ≳ 4, no consensus has been reached yet, and an independent examination is desired.
As an independent examination, very recently, investigations on the SFRF and/or SFRD using not only rest-frame UV but also optical data have been reported (Asada et al. 2021;Asada & Ohta 2022;Rinaldi et al. 2022).The rest-frame optical light suffers much less from dust attenuation than rest-frame UV light, thus utilizing rest-frame optical data is expected to derive the properties of galaxies more properly.Additionally, in the rest-frame optical wavelength, a tracer of the SFR, the Hα emission line, is available.At z ≳ 4, rest-frame optical light including the Hα line is redshifted to midinfreared (MIR) wavelength, and observations with In-fraRed Array Camera (IRAC) on Spitzer are the most suitable among the current facilities.In particular, the Hα emission lines from high-z galaxies are strong enough to boost the IRAC broadband photometry (e.g., Yabe et al. 2009;Stark et al. 2013), and can be recognized through the IRAC color excess.Therefore, using not only rest-frame UV but also optical data including the Hα emission line with Spitzer/IRAC will probe the star formation activities at z ≳ 4 in an independent way of UV-based studies.
As discussed by Asada & Ohta (2022, A22 hereafter), to measure the flux of the Hα emission line from the IRAC broadband photometry, it is necessary to estimate the contribution from the continuum accurately, and photometry at a band free from any other emission lines is essential.Particularly, the observation of the continuum at a longer wavelength band than where the Hα emission line falls is crucial.Without an observation at a longer wavelength, the boosted flux and the IRAC color excess due to the emission line can also be interpreted as a red color of stellar continuum or dust reddening.
There are four broadband filters on IRAC,3.6,4.5,5.8,and 8.0 µm band.To observe both the Hα emission line and the continuum emission at longer wavelength with IRAC, the target redshift must be z ∼ 4.5 (3.9 ≲ z ≲ 4.9), 5.8 (5.1 ≲ z ≲ 6.6), or 7.8 (6.9 ≲ z ≲ 8.6), where the Hα emission line falls into the 3.6, 4.5, and 5.8 µm band, respectively.At z ∼ 4.5, an SFRF is obtained using not only rest-frame UV but also optical data by the spectral energy distribution (SED) fitting method taking the Hα emission into account by Asada et al. (2021, A21 hereafter).The SFRF at z ∼ 4.5 obtained with this approach shows an excess compared to that estimated from the rest-frame UV-based approach, and suggests the (dust-obscured) intensive star formation may play a more important role in the evolution of galaxies at z ∼ 4.5 than previously expected.At z ∼ 7.8, A22 obtained a constraint on the SFRF and SFRD through the Hα LF derived from the nondetection of Hα emission line in the 5.8 µm band.However, the SFRF at z ∼ 5.8 has not been investigated with this approach yet, and the redshift evolution of the SFRF obtained with the approach at z ≳ 4 is not known.
In this work, we aim at investigating the SFRF at z ∼ 5.8 using rest-frame UV and optical data including the Hα emission line.At this redshift, the Hα emission line is observed with IRAC 4.5 µm band, and the continuum emission is observed with the 5.8 and 8.0 µm band.Thus, we utilize the photometric data with all the four broadbands on IRAC, and perform SED fitting to obtain the SFRF.In SED fitting, as similar to A21, we extensively examine the assumptions on the model SED including various SFHs, dust attenuation laws, and a two-component model that consists of a young starforming population and old quenched population, and we take the uncertainty of these model assumptions into account in the resulting SFRF.In addition, we exploit the resulting SFRF at z ∼ 5.8 and other recent estimations of the SFRF at z ≳ 4 with rest-frame optical emissions, and examine the redshift evolution of the SFRF and SFRD at high-z universe with this independent approach.
This paper is structured as follows.In Section 2, we describe the data and the sample selection used in this work.In Section 3, we derive the physical parameters of the sample galaxies through SED fitting, and we obtain the SFRF using the result of SED fitting in Section 4. In Section 5, we discuss the implications from the resulting SFRF including the redshift evolution of the SFRF and SFRD at z ≳ 4. In Section 6, we give further discussions on results in this work, and we summarize this paper in Section 7. Throughout this paper, all magnitudes are quoted in the AB system (Oke & Gunn 1983), and we assume the flat cosmological parameters of H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3.and Ω Λ = 0.7.

DATA AND SAMPLE SELECTION
Among surveys carried out with Spitzer, the Great Observatories Origins Deep Survey (GOODS) fields are one of the deepest and widest fields.Additionally, part of these observations were conducted during the cryogenic mission and thus observations with IRAC 5.8µm and 8.0µm band are available, which are essential to achieve our goals.We focus on both the GOODS-South and -North fields to maximize the survey volume, and use the photometric catalog in the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CAN-DELS) program.

GOODS-South field
In GOODS-South field, we use a photometric catalog given by Guo et al. (2013), which covers a wavelength range from UV to MIR.Here we briefly summarize the method for the source extraction and photometry in the catalog, but we recommend readers to refer Guo et al. (2013) for details.The catalog includes multiwavelength band photometry consisting of observations with the Advanced Camera for Surveys (ACS) and Wide Field Camera 3 (WFC3) on Hubble Space Telescope (HST), IRAC on Spitzer, and several ground-based observatories.The optical data are composed of observations with F435W, F606W, F775W, F814W and F850LP bands on HST/ACS.The near-infrared (NIR) data are composed of observations with HST/WFC3 F098M, F105W, F125W and F160W bands.The MIR data contain observations with Spitzer/IRAC 3.6µm, 4.5µm, 5.8µm, and 8.0µm bands.From ground-based facilities, the catalog also contains VLT/VIMOS and CTIO/MOSAIC U -band data and VLT/ISAAC and VLT/HAWK-I K sband data.
Photometry on HST images was performed using SExtractor's dual-image mode.The combined maxdepth mosaic in the F160W band was used for source detection, and photometry was performed for pointspread-function (PSF)-matched images.On groundbased and Spitzer images, photometry was made through TFIT (Laidler et al. 2007).In this work, we use all the photometry in the 17 bands except for the CTIO U -band data from this photometry catalog since this band was revealed to have a red leak (Guo et al. 2013).
We utilize the CANDELS photometric redshift catalog (Dahlen et al. 2013).For the entire sources in the photometric catalog, photometric redshift is derived using a hierarchical Bayesian approach that combines the probability distribution functions (PDFs) estimated from several manners.

GOODS-North field
In GOODS-North field, we use a catalog given by Barro et al. (2019), which contains photometry from UV to far-infrared.The catalog contains not only the photometry but also photometric redshifts and other stellar parameters, but, in this work, we only use the photometry from UV to MIR and the photometric redshift.The catalog includes mutiwavelength photometry on observations with ACS and WFC3 on HST, IRAC on Spitzer, and several ground-based observatories.For optical data, observations with HST/ACS F435W, F606W, F775W, F814W, and F850LP bands are available.The NIR data consist of observations with HST/WFC3 F105W, F125W, F140W, and F160W bands.The MIR data are composed of observations with Spitzer/IRAC 3.6µm, 4.5µm, 5.8µm, and 8.0µm bands.From groundbased facilities, KPNO/Mosaic and LBT/LBC U -band, Subaru/MOIRCS K s -band, and CFHT/Megacam Kband observations are included.This catalog also includes optical medium-band data that covers λ = 2500 to 9500 Å from the GTC SHARDS survey and WFC3 IR spectroscopic observations with the G102 and G141 grisms, but none of our sample galaxies are detected in these observations (c.f., sample selection in Section 2.3) and we do not use data from these observations in the following analysis.
Photometry in this catalog was carried out in a similar way as in the previous CANDELS catalogs including that by Guo et al. (2013) in GOODS-South field, which is introduced in the previous subsection.The source detection was performed in the F160W-band image, and photometry on all the HST band images was made for PSF-matched images with SExtractor's dual-image mode.For the ground-based and Spitzer/IRAC observations, photometry was conducted through TFIT.
The catalog also contains photometric redshifts for the entire sources in the catalog.The photometric redshifts are also derived in a similar approach as in the previous CANDELS catalogs: the photometric redshift is derived by combining the PDFs estimated from several manners.

Sample Selection
To make a sample of galaxies at 5.09 < z < 6.62 whose Hα emissions are redshifted into the IRAC 4.5µm band, we set the sample selection criteria as follows.
First, we select galaxies based on the photometric redshift.Specifically, we extract galaxies whose bestestimated photometric redshift z best is in our target range and whose relative uncertainty on the photo-z es-timation is less than 5%: 5.09 < z best < 6.62, (1) and ∆z 1 + z < 0.05. (2) Next, we exclude X-ray sources to remove active galactic nuclei (AGNs) from our sample.In GOODS-S field, we eliminate AGNs identified by Hsu et al. (2014).The X-ray sensitivity limits of this AGN catalog are typically 3.2 × 10 −17 , 9.1 × 10 −18 , and 5.5 × 10 −17 erg s −1 cm −2 for the full (0.5-8 keV), soft (0.5-2keV), and hard (2-8keV) bands, respectively.These correspond to luminosity limits of 1.2 × 10 43 , 3.4 × 10 42 , and 2.0 × 10 43 erg s −1 with assuming the redshift z = 5.8, respectively.In GOODS-N field, we eliminate X-ray sources identified by Alexander et al. (2003).The sensitivity limits are about 2.5 × 10 −17 and 1.4 × 10 −16 erg s −1 cm −2 for soft (0.5-2.0 keV) and hard (2-8 keV) bands, which correspond to luminosity limits of 9.2 × 10 42 and 5.2 × 10 43 erg s −1 by the assumption of redshift, respectively.Among the 34,930 (35,445) objects in the CANDELS GOODS-S (GOODS-N) catalog, 300 (254) galaxies pass the criteria above (we refer them as "phot-z sample").We apply additional cuts to secure the reliability on SEDs of galaxies in our sample and to make a robust estimation of their physical properties by SED fitting.As introduced in Section 1, the detection of the restframe optical continuum is important to estimate the contribution of the continuum in 4.5µm band.Thus, we require the signal-to-noise ratio (S/N) both in IRAC 5.8µm and 8.0µm bands larger than 3.Among the 554 (= 300 + 254) galaxies in the phot-z sample, 23 galaxies have S/N > 3 both in IRAC 5.8µm and 8.0µm bands.We also remove galaxies that are heavily blended with the neighboring objects in the IRAC images.Although the photometry in IRAC bands was performed with a deblending tool, TFIT, heavy contamination can lead to a systematic uncertainty in the photometry.To this end, we eliminate all the galaxies whose separation from the nearest object in the photometric catalog is < 2 ′′ , considering the FWHMs of IRAC images are ∼ 2 ′′ .Finally, we conduct visual inspections on the galaxies that meet criteria above, and obtain a sample of galaxies that contains nine galaxies (six from GOODS-S and three from GOODS-N).A list of these nine galaxies is shown in Table 1, and the distribution of apparent magnitudes at H 160 and 5.8µm band is shown in the lower left panel of Figure 1.In Figure 2, postage stamps of galaxies in the final sample are shown.In Appendix A, the SEDs of the nine galaxies are also shown.

SED fitting
To estimate the physical properties including the SFR, we first carry out SED fitting on the galaxies in our sample.We use the population synthesis code Pégase.3(Fioc & Rocca-Volmerange 2019) to make model SEDs.This code includes the nebular emissions in the model spectrum, and follows the chemical evolution in a selfconsistent way, which is used to calculate the metallicity of the interstellar medium and the next generation stars 1 .We adopt Chabrier03 IMF (Chabrier 2003) with a mass range of 0.08-120 M ⊙ , and the attenuation by intergalactic neutral hydrogen is modeled following the prescription given by Madau (1995).The internal dust attenuation is modeled with Calzetti law (Calzetti et al. 2000) and Small Magellanic Cloud (SMC; Pei 1992) law.We assume R V = 4.05 and R V = 2.93 for Calzetti and SMC law, respectively.As for the star formation history (SFH), we consider seven models; constant star forma- tion (CSF), exponentially declining (∝ e −t/τ ), and delayed exponential (∝ te −t/τ ) models with the e-folding time of τ = 10, 100, 1000 Myr.
Under these assumptions, we search for the bestfitting model SED with the free parameters of age and color excess E(B −V ) by χ 2 minimization.We allow the age of galaxies to vary with ∼ 60 steps from 1 Myr to the age of the universe at the redshift where each galaxy is located.The color excess is taken from 0.0 to 0.8 mag at an interval of 0.01 mag.The ratio between stellar and nebular attenuation is assumed to be 1, since this ratio is suggested to be about unity in the high-z universe (e.g., Cullen et al. 2014;Reddy et al. 2015;Theios et al. 2019).Note that, we find that the best-fitting E(B − V ) is larger than 0.8 for only one galaxy (29111-GN), and thus we allow E(B −V ) to vary up to 1.5 for this galaxy.We do not use the photometry at the wavelength shortward to Lyα in calculating the χ 2 .Consequently, we obtain 14 best-fitting model SEDs for each galaxy in the final sample; two types of assumptions on dust attenuation law and seven on SFHs, and each of them gives a respective set of best-fitting physical parameters including age, E(B − V ), stellar mass M ⋆ , and SFR.

Two-component model SED fitting
a Name represents the ID given in the CANDELS catalogs and the field where the galaxy is located.
b The rest-frame UV spectral slope β is directly measured from the rest-frame UV photometry.
c The spectroscopic redshift is available only for 3073-GS, zspec = 5.563.
d Some sources have a relatively large magnitude errors on 5.8-and 8.0-µm band photometry although they pass the S/N cut in the sample selection.This is because we used the limiting magnitudes as defined by averaged rms values given in the catalog for the S/N calculation while the magnitude errors in this table are nominally taken from the photometric errors in the catalog.
As discussed by A21, the SED fitting on high-z galaxies with a red SED from rest-frame UV to optical wavelength range can result in a significantly dusty and highly star-forming solution, but the red color in SEDs can be explained not only by the dust attenuation but also by the aging of stellar population.A21 showed that considering a two-component model that consists of old quenched stellar population and young star-forming population in SED fitting can reduce the best-fitting E(B − V ) and SFR particularly for highly star-forming solutions, leading to an effect on the resulting SFRF.We thus also consider such a two-component model.Here, we adopt a similar approach as that by A21.For simplicity, we assume there is no dust attenuation in the old stellar population, and fix the spectrum of the old stellar component to that is expected for the reddest in the wavelength range of rest-frame UV to optical within the age of the universe.Since the age of the universe at z ∼ 5.8 is ≲ 1 Gyr, the age of the old component is fixed to 500 Myr2 .We compare the continuum spectra at the age of 500 Myr with four SFHs; instantaneous burst and CSF with a duration of ∆τ = 10, 50, 100 Myr, and find that the CSF with ∆τ = 100 Myr gives the reddest spectrum, which is consistent with the result shown by A213 .Thus, we adopt the spectrum from the CSF model with ∆τ = 100 Myr aged 500 Myr as the old stellar population spectrum.
Using this spectrum, we subtract the old population SED with a fraction of f old described below from the observed SED, and perform SED fitting as described in Section 3.1 to the residual, which gives the best-fitting SED of "young star-forming population".Here, we allow the contribution of the old stellar population to the whole SED to vary; specifically, we allow the fraction of flux in the 5.8µm band by the old stellar population to the total (observed) flux (f old ) to be changed from 0 to 0.95.Thus, in this two-component model, free parameters and the assumptions for the young starforming population are the same as those in the normal, one-component model that is described in Section 3.1, and one additional free parameter f old is considered.We obtain 14 best fit SEDs and the number of sets of best-fitting physical parameters for each galaxy with this two-component SED fitting as well.The M ⋆ obtained in this two-component model is the estimation for the whole system, and the SFR, age, and E(B − V ) are those for the young component.Note that, since the SFR of the old component is assumed to be zero, the SFR obtained here can be regarded as the estimation for the whole system.
Among the 14 sets of best-fitting physical parameters from this two-component model, in the following analysis, we only use results that meet following criteria: I.The best-fitting age for the young population is less than 400 Myr.
II. f old is larger than 0.5.
In this model, the old population is assumed to be quenched 400 Myrs ago, so the young component must be younger than 400 Myr (I).If the best-fit f old is small, that means old stellar population is almost negligible and/or observed SED is well fitted without such an old population (II).

Physical parameter estimation
From the results in Section 3.1 and 3.2, we obtain 14 or less sets of estimations4 of physical parameters (i.e., age, E(B − V ), SFR, and M ⋆ ) for each attenuation law for each galaxy in the final sample.
Using these sets of estimations, we aim at deriving the best estimation of each physical parameter taking into account the goodness-of-fit of the model SEDs under the assumption of the attenuation law for each galaxy.
First, we evaluate the goodness-of-fit for each of the 14 or less estimations using the Bayesian Information Criterion (BIC) that is defined as where χ 2 min is the minimized χ 2 given by the best fit model SED, q is the number of free parameters (q = 2 for results from Section 3.1 and q = 3 for results from Section 3.2), and m is the number of independent observations.The smaller value of BIC means the model SED fits better to the observed SED, and the likelihood L of the model SED given the observed data is L ∝ exp(−BIC/2).Naively, the estimation that gives the smallest BIC is the most likely, and thus one may think adopting the physical parameter estimation with the smallest BIC (which gives the maximum likelihood) is reasonable.However, the maximum likelihood value is not distinctively larger than those in other cases, and the difference in physical parameter estimations is not negligible between the maximum likelihood model and other models with comparable likelihood (c.f., Figure 13 in Appendix B).Thus adopting only the estimation with the maximum likelihood and discarding all the other estimations seem to be risky.To take this into account, for each physical parameter, we calculate the weighted mean and weighted standard deviation using exp(−BIC/2) as the weight, and use them as the best estimations and its uncertainties for each physical parameters.We derive weighted means and weighted standard deviations for each of attenuation law assumptions (Calzetti law or SMC law).This gives the physical parameter estimations under the assumption of attenuation law.We also derive these values using all the set of best-fit SEDs, including both of Calzetti and SMC assumption, to obtain a best estimation taking the uncertainty of dust attenuation law into account (we call this estimation as "Calzetti+SMC estimation").
As a result, for each physical parameter, we derive three estimations and the associated uncertainties: Calzetti, SMC, and Calzetti+SMC estimations.In Appendix B, we show an example of this physical parameter estimation described in this subsection for one galaxy in the final sample (1411-GS).
In Figure 3, we show the resulting stellar mass against SFR obtained with Calzetti+SMC, Calzetti, and SMC estimation in the left, middle, and right panel, respectively.For comparison, we show a main sequence (MS) relation at z ∼ 5.8 predicted by the redshift evolution of the MS relation by Schreiber et al. (2015).Schreiber et al. (2015) investigated the MS relation from z = 0 to z = 4 and derived the redshift evolution.Here we extrapolate the evolution up to z ∼ 5.8 to plot the MS relation in the figure.We can see that most galaxies in the final sample are aligned with the MS relation without a few exceptions; one starburst and one or two quiescent galaxies.Thus, the final sample mainly contains (massive) MS galaxies at this redshift and a few starburst or quiescent galaxies.Note that the SED of 30423-GS is well fitted with quiescent solutions, and its SFR value is so low that the galaxy is missing in Figure 3 and 4. The MS relation at z ∼ 5.8 does not differ significantly between literature, and the consequence here persists even if we compare with MS relations obtained in other studies (e.g., Speagle et al. 2014;Santini et al. 2017;Popesso et al. 2022).In Appendix A, we show the results of SED fittings and the resulting estimations of physical parameters for all the galaxies in the final sample.
The methodology described in this subsection is essentially for estimating physical parameters taking various SFHs into account.This method is simple and easy to be applied even with results from SED fitting code assuming parameterized SFH, but the statistical validity may be unclear.Thus, it is worth comparing with other non-parametric SFH SED fitting approaches.We use a SED fitting code dense basis (Iyer & Gawiser 2017;Iyer et al. 2019), which assumes non-parametric SFH, and examine how it affects the results in this paper.Results in this paper persist when we use dense basis estimations instead of the weighted-mean estimations (see Appendix C).

STAR FORMATION RATE FUNCTION
In this section, we describe the method for deriving the SFRF using the result in Section 3. We use a similar approach to derive the SFRF as by A21, but slightly modified to apply to the data set in this work.

Method for deriving SFRF
It is necessary to correct for the incompletness of the sample to properly derive the SFRF.As similar to A21, we consider three factors causing the incompleteness: (i) detection in the H 160 band, (ii) S/N cut in the 5.8 µm and 8.0µm band, and (iii) elimination of the blended objects in the IRAC images.
To correct for (i), the detection rate in the H 160 -band image for each field is necessary.In GOODS-S, Duncan et al. (2014) derived the detection rate of the same CANDELS GOODS-S catalog as we use in this work by performing mock observations, and we use the result.The solid lines in the top panel of Figure 1 show the detection rate as a function of H 160 magnitude.Duncan et al. (2014) obtained the detection rate by dividing the GOODS-S field into four subregions according to the exposure time; HUDF, ERS, DEEP, and WIDE.In GOODS-N, we estimate the H 160 -band detection rate from the differential number counts for objects in the catalog following the prescription by Guo et al. (2013).Barro et al. (2019) derived the differential number density and the best-fit power law function dividing the field into two subregions (DEEP and WIDE).We use the result and estimate the detection rate in the faint region by calculating the ratio of the observed number density to that predicted by the best-fit power law.The dashed lines in the top panel of Figure 1 show the resulting detection rate.
As for the incompleteness due to (ii), since the exposure times in IRAC 5.8µm and 8.0µm band are also inhomogeneous, we correct for the effect of S/N cut considering this inhomogeneous depth of these images.To this end, we calculate the 3σ limiting magnitudes in the 5.8µm band at the locations of galaxies in the phot-z sample (see Section 2.3 for the definition of the phot-z sample), and make the normalized cumulative histogram of them in each GOODS-S and -N fields.In the CANDELS GOODS-S catalog, 1σ limiting magnitudes at the locations of the objects are available, and we use them to calculate the 3σ limiting magnitudes to make the normalized cumulative histogram.In the CANDELS GOODS-N catalog, on the other hand, they are not available, thus we use the flux uncertainties of galaxies in the phot-z sample instead.Since the photz sample are composed of galaxies at z ∼ 5.8 without AGN activities, their photometric errors are expected to be background limited5 .As discussed by A21, we can use the normalized cumulative histogram of the 3σ limiting magnitudes in the 5.8 µm band of galaxies in the phot-z sample to correct for the incompleteness due to S/N cut in the 5.8µm band.In the sample selection, we applied S/N cut not only in 5.8µm band but also in 8.0µm band, so one might think an additional correction is required.However, firstly, the limiting magnitudes at 5.8µm and 8.0µm band are almost identical.Secondly, the observed color of [5.8] − [8.0] is typically 0 or more (c.f. Figure 12 in Appendix A).Finally, the exposure maps in 5.8µm and 8.0µm band are similar to each other, and regions with a deep observation in 5.8µm band is also expected to have a deep observation in 8.0µm band.Considering these factors, galaxies that are detected in the 5.8µm band with S/N> 3 should be detected also in the 8.0µm band (c.f., Figure 2), and thus an additional correction should not be necessary.In the right panel of Figure 1, the resulting normalized cumulative histograms in the GOODS-S and -N fields are shown.
We correct for the incomleteness due to (iii) using the fraction (f sel ) of isolated objects both in 4.5-and 5.8-µm band IRAC images at z ∼ 5.8 as similar to A21.To evaluate f sel , we first randomly pick out ∼ 200 galaxies from the phot-z sample, and categorize them into the following three groups by visual inspections; galaxies that are detected in the 5.8µm band image and not blended with the neighbors (group A), galaxies that are detected in the 5.8µm band image but heavily blended with the neighbors (group B), and objects that are not detected in the 5.8µm band image.We then calculate f sel by #(A)/(#(A) + #(B)), and obtain result in f sel ∼ 0.52.
Taking these incomleteness factors into account, we can estimate the SFRF ϕ( Ṁ⋆ ) [Mpc −3 dex −1 ] in the bin from Ṁ⋆,1 to Ṁ⋆,2 with a bin width of d(log Ṁ⋆ ) as follows; where the subscript i represents each galaxy in the final sample, N gal is the number of galaxies (= 9), P i ( Ṁ⋆ ) is the PDF of the SFR for the galaxy i, and the V eff i is the effective volume of this survey for the galaxy i.We assume Gaussian profiles for P i ( Ṁ⋆ ), with the best estimations and their uncertainties obtained in Section 3.3 being used as the mean and sigma of the profile.This effective volume V eff i can be calculated as where subscript k represents the regions in the fields, N region is the number of regions (= 6), f 160 k is the detection rate in the H 160 band as a function of the H 160 -band magnitude (i), f ch3 k is the value in the cumulative histogram of 3σ limiting magnitudes as a function of the 5.8µm-band magnitude (ii), f sel is the fraction of isolated objects in IRAC images (iii), m 160,i and m ch3,i are the observed magnitudes for the galaxy i in the H 160 and 5.8µm band, respectively, Ω k is the solid angle of the region k, and r z=a is the comoving distance from z = 0 to z = a.
We also estimate the completeness limit of the SFR, above which the SFRF obtained with Equation ( 4) is well-corrected for the incompleteness.To estimate the completeness limit of the SFR, we show the diagram of 5.8µm-band magnitude versus the SFR (Figure 4).We can see a rough anti-correlation between the [5.8] and SFR regardless of the assumption of the dust attenuation.However, due to the small number of the sample galaxies, it is difficult to estimate the completeness limit only from this distribution.In the figure, we also plot a guideline that shows the upper limit on the SFR by the black dotted line.The rest-frame Hα equivalent width (EW) cannot be larger than ∼ 4000 Å even if most extreme stellar populations (e.g., PopIII stars) are considered (e.g., Inoue 2011).If we assume the 5.8µm flux density to be the continuum around Hα at z ∼ 5.8, we then can obtain an upper limit on the Hα luminosity for a given magnitude in the 5.8 µm band, which leads to an upper limit on the SFR6 .In spite of a rough estimation of the upper limit, this constraint is consistent to the SFRs derived in this work.Considering this upper limit and the overall distribution of galaxies in the [5.8]-SFR diagram, galaxies with SFR ≳ 10 2.5 M ⊙ yr −1 should be brighter than ∼ 24.5 mag in the 5.8µm band, but those with SFR ≲ 10 2.5 M ⊙ yr −1 can be fainter than ∼ 24.5 mag in the band.Thus, the limiting magnitude of ∼ 24.5 mag in the 5.8µm band (c.f.lower right panel in Figure 1) corresponds to the completeness limit of SFR ∼ 10 2.5 M ⊙ yr −1 .

Star Formation Rate Function
In Figure 5, we show the resulting SFRFs obtained under the each assumption of dust attenuation law (Calzetti+SMC, Calzetti, and SMC estimation).The uncertainty of the SFRF is calculated assuming the Poisson uncertainty by Gehrels (1986).As discussed in Section 4.1, the completeness limit of the sample is ∼ 10 2.5 M ⊙ yr −1 , thus the data point below that value is treated as a lower limit.In the figure, we also show several SFRFs from literature for comparison converted to the same IMF as we use.Smit et al. (2016, S16 hereafter) derived an SFRF based on the rest-frame UV LF at z ∼ 6 after corrected for the dust attenuation with assuming the A UV -β relation.S16 assumed Meurer et al. (1999)-and SMC-type A UV -β relation and derived the respective SFRFs.For a fair comparison, we compare our result to SFRFs in the literature under the corresponding assumption of dust attenuation law.For result using Calzetti+SMC estimation (left panel in Figure 5), we simply compare to the average of the SFRFs with Calzetti and SMC law from the literature.In addition, we show the SFRF obtained from FIR observations at z ∼ 6 estimated by Rowan-Robinson et al. (2016) in the figure .We have to keep it in mind that this SFRF from FIR observations was measured only in the brightest region (SFR ∼ 10 4 M ⊙ yr −1 ) and obtained by extrapolating down to the fainter region.
The SFRF obtained in this work (red circles) clearly shows an excess compared to that estimated from dustcorrected rest-frame UV LF (blue solid line).Although the data points of the SFRF above SFR = 10 2.5 M ⊙ yr −1 contain a small number of galaxies (2 galaxies) and the uncertainties on these data points are large, the lower limit at the bin of SFR = 10 1.5−2.5 M ⊙ yr −1 also suggests that the SFRF obtained in this work is in an excess compared to the rest-frame UV-based estimation, and thus the excess seems to be rather robust.This result is similar to that by A21, which showed the SFRF at z ∼ 4.5 derived from SED fitting shows an excess compared to that estimated from the dust-corrected rest-frame UV LF.On the contrary, the resulting SFRF is roughly consistent with the FIR estimation given by Rowan-Robinson et al. ( 2016) in spite of their large extrapolation.
As discussed by A21, the difference in the SFRFs between that obtained with SED fitting on data from rest UV to optical and that estimated from the rest-frame UV LF stems from the difference in the estimation of the dust attenuation.In Figure 6, we compare the estimated amount of dust attenuation A 1600 in this work and that from the A UV -β relation.In the figure, we show the comparison with assuming Calzetti attenuation law, and use the A UV -β relation given by Meurer et al. (1999): The rest-frame UV spectral slope β is measured by fitting a power law to the rest-frame UV photometry for each object in the final sample.The figure shows galaxies with larger SFRs tend to suffer from heavier dust attenuation, and the amount of dust attenuation on such galaxies is likely to be underestimated with the A UV -β relation.These results suggest that (1) previous approach to estimate the SFRF based only on the rest-frame UV observations underestimate the contribution from dust-obscured SF particularly in the brighter (larger SFR) region and (2) the contribution from dustobscured SF can be well corrected for by using not only rest-frame UV but also optical data.
One may think the difference stems from a difference in sample selection.Although we use catalogs selected in the rest UV wavelength (i.e., F160W detected), we make a sample based on photometric redshifts which is not exactly the same way as in part of rest-frame UV based works that uses Lyman break selection.Practically, the Lyman break selection requests relatively blue in the rest-frame UV slope, while photometric redshift selection does not necessarily require.In particular, two galaxies are mainly contributing to the highest-SFR bin where the excess is most significant, and these two galaxies can be the key origin of the excess.However, both of the two galaxies meet the rest-frame UV slope crite- ria in i 814 -dropout selection (e.g., Atek et al. 2015) and their color can be classified as LBGs at this redshift.Thus, the difference in sample selection is not expected to be the main reason for the difference in the estimated SFRFs.

Redshift evolution of the SFRF
To see the redshift evolution of the SFRF, we also compare the result with the SFRF at z ∼ 4.5 derived from SED fitting by A21 (black dashed lines in Figure 5).The number density derived in this work at the high-SFR end region is not largely different from that at z ∼ 4.5.This suggests that the SFRF obtained with an approach using data from rest-frame UV to optical does not significantly evolve from z ∼ 4.5 to z ∼ 5.9 at the high-SFR end region.In this subsection, we will further investigate the redshift evolution of the SFRF obtained with this approach.
One way to quantify the evolution of the SFRF is to parameterize the SFRF with an analytical form such as Schechter function.Thus, we first fit the Schechter function to the SFRF derived in this work.The Schechter function can be written as where SFR * , ϕ * , and α is the characteristic SFR, the characteristic number density, and the low-SFR end slope, respectively.Since the SFRF obtained in this work puts a constraint only in the high-SFR end region and the number of data points is small, we fix the α and ϕ * and allow only SFR * to vary.We fix α to the same as that of the SFRF from rest-frame UV LF at z ∼ 6 (S16) since only rest-frame UV observations give constraints on α at this redshift: α = −1.63 and −1.72 for Calzetti and SMC, respectively.We also fix ϕ * to the same as that estimated in a similar approach at z ∼ 4.5 (A21): ϕ * = 7.24 × 10 −5 and 3.89 × 10 −5 for Calzetti and SMC, respectively.We search for the best-fit parameter of SFR * by minimizing χ 2 , and the lower limit on the SFRF at the bin of SFR = 10 2.0 M ⊙ yr −1 is also taken into account in calculating χ 2 by reference to the approach given by Sawicki (2012): Here, subscript i and j denotes bins where the observed SFRF is complete and not complete (lower limit), respectively.ϕ d,i is the observed SFRF at the ith bin, σ i is its uncertainty, ϕ lim,j is the lower limit on the SFRF at the jth bin, and ϕ m,i or ϕ m,j is the SFRF calculated from the Schechter function at the corresponding bin.In this calculation using Equation (8), we treat the lower limit at the bin of SFR = 10 2.0 M ⊙ yr −1 as 1σ lower limit.
The resulting best-fit SFR * is shown in Table 2. Regardless of the assumption on the dust attenuation law, the best-fit SFR * does not change significantly and is consistent with each other within the uncertainties.Since recent ALMA observations in the high-z universe suggest SMC-like attenuation law can be more appropriate in modeling high-z SF galaxies (e.g., Fudamoto et al. 2020;Schouws et al. 2021), we use the result with SMC estimation to discuss the evolution of the SFRF hereafter in this subsection.
To see the redshift evolution, we next compile measurements of the SFRF at a wide range of redshift using various probes of the SFR, and compare the best-fit Schechter parameters.At z ≲ 3, investigations with various probes are available and we compile results using Hα (Sobral et al. 2013), infrared (Magnelli et al. 2011), and bolometric (Bell et al. 2007;Reddy et al. 2008) observations.At z ≳ 3, on the other hand, we compare results using (dust-corrected) rest-frame UV LFs (S16) and those with not only rest-frame UV but also optical data including this work (A21; A22).For a fair compari-  son, from S16, we use results using an SMC-type A UV -β relation in the dust correction of the UV LFs: Results in the literature using other IMF are converted to the same IMF as we use in this work.When bestfit Schechter parameters for the SFRF are not quoted in literature, we derive the best-fit parameters by ourselves using the result in the literature.The detailed procedure for this fit is given in Appendix D.
Figure 7 shows the resulting evolution of the SFR * .At z ≲ 3, results with various probes are roughly consistent with each other, although the estimations from Hα LF give slightly lower values.As discussed by Smit et al. (2012), this may be originated from an incomplete sampling or inappropriate correction for the dust attenuation of the Hα LFs.The SFR * decreases from z ∼ 2 to the local by ∼ 1 dex, which means the star formation activities in galaxies are typically suppressed from z ∼ 2 to the present universe.Recent hydrodynamical cosmological simulation claimed that the main mechanism of this quenching is AGN feedback (e.g., Katsianis et al. 2017).Other observational studies (Harikane et al. 2022) claimed this quenching is attributed to a decrease of the dark-matter accretion rate onto halos due to the cosmic expansion.
At z ≳ 3, the SFR * s estimated with the approach using rest-frame UV to optical data (red symbols in Figure 7) show as high values as that at z ∼ 2, while those from rest-frame UV LFs (blue symbols in Figure 7) show lower values than at z ∼ 2 by ∼ 1 dex throughout z ≳ 3.With the previous measurements using rest-frame UV LFs, SFR * is thought to decrease from z ∼ 2 to z ∼ 8 (e.g., Smit et al. 2012).However, the SFRF measurements using rest-frame optical data indicate the SFR * is roughly constant from z ∼ 2 to z ∼ 6 and may decrease above z ∼ 6.Since the parameter SFR ⋆ is sensitive to the high-SFR end of the SFRF, this evolution of the SFR ⋆ suggests that the high-SFR end of the SFRF grows as early as by z ∼ 6, or in the epoch of reionization, to a similar value as at z ∼ 2 when the star formation activity in the universe reaches its peak.

Star Formation Rate Density
The star formation rate density (SFRD) is usually measured by integrating the SFRF down to a lower bound.The rest-frame UV absolute magnitude of M UV = −17.0mag is commonly adopted as the lower bound, which corresponds to a SFR of ∼ 0.22 M ⊙ yr −1 .The SFRF measurement in this work is only in much larger SFR region than this lower bound of integration.We here calculate the SFRD in this work as follows: where the notations are the same as those in Equation (4).Since the summation is taken for galaxies in the final sample which is complete only above the completeness limit (SFR = 10 2.5 M ⊙ yr −1 ), this values gives a lower limit on the SFRD.We compare this limit with the previous measurements using the lower bound of integration, M UV = −17.0mag.The results are shown in Figure 8 and Table 2.In the figure, for a comparison, we also show recent measurements of SFRD from literature using FIR (Rowan-Robinson et al. 2016;Gruppioni et al. 2020), rest-frame UV (Bouwens et al. 2020), both rest-frame UV and optical (A21), and Hα (A22) observations.At z ≳ 5, part of FIR observations (e.g., Rowan-Robinson et al. 2016;Gruppioni et al. 2020) indicate a more intensive star formation activities than dust-corrected rest-frame UV observations (e.g., Madau & Dickinson 2014;Bouwens et al. 2020, and see Figure 8).With Calzetti+SMC (the red filled circle) or Calzetti (the red filled square) estimation, the lower limit obtained in this work suggests that the estimations from rest-frame UV observation are marginally acceptable or underestimate the SFRD, and the measurements with FIR observations are more likely, even though the lower limit does not include the contribution from fainter (lower SFR) galaxies.On the contrary, the constraint obtained with SMC estimation (the red filled diamond) is consistent with both the measurements with FIR and rest-frame UV observations and not tight enough to determine which measurements are more likely.
It is worth noting that the assumption on the dust attenuation law can affect the SFRD estimations.Particularly, if SMC law is assumed instead of Calzetti law, the estimated SFRD can decrease.Madau & Dickinson (2014) (black solid line in the figure) calculated the correction factor for dust attenuation using Equation ( 6), which is almost the same as a correction using Calzetti law.From A21, we show the estimation using Calzetti law in Figure 8 (the open red hexagon), but A21 discuss that the SFRD can be decreased by 0.2-0.6 dex when SMC law is assumed instead.Bouwens et al. (2020) (open blue circles) assumed a Calzetti-like and SMC-like A UV -β relation for high-mass and lowmass galaxies, respectively.Thus, the comparison of our result especially with SMC law to the literature may not be fair in that different assumptions are adopted.Nevertheless, the consequences from comparisons with previous measurements using several probes above do not change.At z ∼ 6, the dust correction calculated from rest-frame UV observations (i.e., A UV -β relation) is ∼ 0.20 dex (Bouwens et al. 2020) and thus the dustuncorrected SFRD is smaller than the dust-corrected value shown in the figure by ∼ 0.20 dex.Even with the dust-uncorrected value, the lower limit obtained in this work with the SMC estimations is still consistent, which is the same consequence as above.For comparisons with Calzetti+SMC and Calzetti estimation, the results are strengthened if the estimations with rest-frame UV observations are decreased.

Comparison of other physical properties to the previous studies
In Section 4, we derived an SFRF based on SED fitting, and showed the SFRF is in an excess compared to that estimated only from rest-frame UV observations.This difference may be originated in a difference of the sample (e.g., field-to-field variations) or a systematic difference in the SED modelings in this work.To examine these possibilities, first, we derive the rest-frame UV LFs using the sample of galaxies in this work and compare them with those in the previous works.We derive the rest-frame UV LFs with both the phot-z sample and final sample.The (dust-uncorrected) rest-frame UV absolute magnitude for each galaxy is calculated using the  photometric redshift z best and the observed flux in the rest-frame UV wavelength.With the final sample, we derive the UV LF using the similar approach as we used in deriving the SFRF (Equation ( 4)).With the photz sample, since this sample is not subject to selections in IRAC photometry (i.e., the S/N cut in 5.8-and 8.0µm band and the elimination of blended source in IRAC images), the UV LF is derived without the corresponding incompleteness correction factors: f sel and f ch3 k in Equation ( 5).
The resulting UV LFs with the final and phot-z sample are shown in the top and bottom panel of Figure 9, respectively.For comparison, best-fit UV LFs in literature are also shown (Bouwens et al. 2015;Harikane et al. 2021).The SFRF derived by S16 (blue solid line in Figure 5) is based on the UV LF by Bouwens et al. (2015), and thus we show it for comparison (blue dotted line in Figure 9).Harikane et al. (2021) found a significant excess of the rest-frame UV LFs compared to the previous best fit UV LFs written in the analytical form of Schechter function in the bright-end region, and suggested that the rest UV LF is well fitted by a double power law (DPL) rather than Schechter function.To see if the excess in the SFRF found in this work stems from this excess in the UV LFs, we also show the UV LF by Harikane et al. (2021) in the figure (blue solid line).
The UV LFs obtained with both the phot-z sample and final sample agree well with those in the literature.In particular, the difference between DPL and Schechter function fit (blue solid and dotted lines) emerges at the brighter population of galaxies (M UV < −23 mag) than that is used in this work.This indicates that the excess of the SFRF is not caused by the difference of the (dust-  2021) and Bouwens et al. (2015), respectively.In the bottom panel, we also show the UV LFs with galaxies in the phot-z sample in GOODS-South (North) field by red (green) squares.
uncorrected) rest UV LF itself but by the amount of dust correction.
In addition, we derive the UV LFs with the photz sample in GOODS-S (red squares) and GOODS-N (green squares) fields.The resulting UV LFs are in a good agreement with each other and with that of the phot-z sample as a whole (black circles).This indicates that there is no field-to-field variance between GOODS-S and -N, and our approach to derive the statistical properties such as SFRF and UV LF with combining sample of galaxies in GOODS-S and -N works well.Although the UV LF with GOODS-N phot-z sample seems to be systematically underestimated at the fain-end region (M UV ∼ −19 mag), this is simply because the GOODS-N sample is incomplete at this region while the GOODS-S sample (and thus the phot-z sample as a whole) is still complete.The max-depth of the rest-frame UV observation by HST is different in the GOODS-S and -N due to the presence of Hubble Ultra Deep Field region in the GOODS-S field (see the top panel in Figure 1), and thus the completeness limit of the UV LF for GOODS-N sample is shallower than those for GOODS-S and phot-z sample.
Next, to examine if there is a systematic difference of SED modelings in this work and those in the previous works, we derive the galaxy stellar mass function (GSMF) using the result of SED fitting obtained in Section 3.3, and compare it with other studies.Again, we use the same approach to derive the GSMF as we used to derive the SFRF (Equation ( 4)), and derive the GSMF under the each assumption of dust attenuation law (Calzetti+SMC, Calzetti, and SMC).However, the difference in the resulting GSMFs with the different attenuation laws is much smaller than that in the SFRFs and the three GSMFs are almost identical.Thus, we only use the resulting GSMF with the Calzetti estimation hereafter.
The result is shown in Figure 10.In the figure, the previous measurements of GSMF at z ∼ 6 are also shown for comparison (Duncan et al. 2014;Song et al. 2016;Deshmukh et al. 2018).The resulting GSMF is in a good agreement with the previous estimations of the GSMF from the literature except for the lowest stellar mass bin (M ⋆ = 10 9.5 M ⊙ ), where our sample is incomplete since the detection limit in 5.8µm (and 8.0µm) band in this work is shallow.Since the GSMFs in the previous works are also derived based on SED fitting, significant difference in the SED modeling between this work and literature is not indicated from this comparison.

Consistency with observations in other wavelengths range
We have used photometry from rest-frame UV to optical and derived the SFRF and other physical properties.On the contrary, part of galaxies in the final sample are also observed in other wavelengths range such as farinfrared (FIR) and radio.In this subsection, we compare our result of SED fittings with the FIR and/or radio observation to see the consistency between the estimations with SED fitting and observations.
To compare with FIR observations, we calculate the total FIR luminosity (L TIR ) from the result of SED fitting.Namely, we calculate the total energy attenuated by dust for each of the best-fit model SEDs obtained in Section 3.1 and 3.2, and regard this energy to be re-emitted in FIR through the dust thermal emission assuming the energy conservation.Using the bestfitting L TIR , we derive three set of weighted means and weighted standard deviations of log 10 (L TIR ) for each of galaxies in the final sample with the Calzetti, SMC, and Calzetti+SMC estimation as we did for other physical parameters such as the SFR and stellar mass in Section 3.3.Here, we calculate the weighted mean of log 10 (L TIR ) rather than L TIR itself since the estimated total energy from SED fitting largely changes by several order of magnitudes with the model assumptions.When the best-fitting L TIR is 0 (i.e., a dust-free solution E(B − V ) = 0), its logarithm is treated as a sufficiently small number (log(L TIR /L ⊙ ) = 9.3, which corresponds to SFR IR ∼ 0.2 M ⊙ yr −1 ).The resulting estimations are shown in Table 3.

FIR and/or radio observations
In this subsection, we introduce the observations in FIR and/or radio wavelengths available from literature.
• 3073-GS The galaxy 3073-GS is also selected as the target in the ALPINE-ALMA survey 7 (Le Fèvre et al. 2020).The galaxy is observed in ALMA Band 7 and its dust continuum emission is not detected with an upper limit of log 10 (L TIR /L ⊙ ) < 11.31 (Béthermin et al. 2020).
• 17541-GS and 30423-GS These two galaxies of 17541-GS and 30423-GS are located in the region covered with an unbiased survey in ALMA Band 6 by GOODS-ALMA survey (Franco et al. 2018(Franco et al. , 2020)).There are no counterparts of 17541-GS or 30423-GS in the GOODS-ALMA catalog, and thus these two galaxies are not detected in ALMA Band 6.The limiting flux of the catalog is ∼ 640 µJy.Using an FIR SED template given by Schreiber et al. (2018), this corresponds to an upper limit of log 10 (L TIR /L ⊙ ) < 12.59.
• 29111-GN The galaxy 29111-GN is detected in a wide wavelength range from MIR to radio.The radio emission from the galaxy is detected at 1.4 GHz with the Very Large Array (VLA) with a flux of f ν = 25.4 ± 3.2 µJy (Owen 2018).Liu et al. (2018, L18 hereafter) used the photometry at Spitzer/MIPS 24 µm and VLA 1.4 GHz in GOODS-N field as the prior information to deblend the highly confused images by Herschel and ground-based FIR to (sub-)millimeter observations.For this galaxy 8 , L18 obtained an SED from MIR to radio and estimate the total FIR luminosity of log 10 (L TIR /L ⊙ ) = 13.24 ± 0.03.
In addition to the total FIR luminosity, for the galaxy, we compare the SFR derived from SED fitting with that directly estimated from the radio observation.Given the photometric redshift z ph = 5.67, the radio emission observed at 1.4 GHz corresponds to the emission at a rest-frame frequency of 9.3 GHz, where the non-thermal (synchrotron) emissions from supernova remnants are dominant in bright starburst galaxies.Condon (1992) obtained a conversion factor from (non-thermal) radio luminosity L NT to the SFR for massive (> 5 M ⊙ ) stars assuming the extended Miller-Scalo IMF.For a fair comparison, starting with the formulation given by Condon (1992) that links the supernova rate and (non-thermal) radio luminosity L NT , we calculate the conversion factor between L NT and SFR for the mass range of 0.08-120 M ⊙ under the same assumption of IMF (i.e., Chabrier03 IMF): 7 The ID of this galaxy in ALPINE survey is CAN-DELS GOODSS 14. 8 The ID of this galaxy in L18 is 14914 Figure 11.Comparisons of the SFR (left panel) and LTIR (right panel) derived from SED fitting using data from restframe UV to optical (red circles) with those derived from independent observations in other wavelengths range (black squares).We show the comparison with Calzetti, SMC, and Calzetti+SMC estimations in the top, middle, and bottom panel, respectively.As for the SFR, we compare the SFR estimated from SED fittings and that from the radio flux density using Equation (11).See text for details.
where ν rest is the rest-frame frequency of the observation and α is the spectral index (see Appendix E for the detail derivation of this conversion factor).Here we use α = 0.8 (Condon 1992) and obtain SFR rad = 5286 ± 666 M ⊙ yr −1 from the VLA observation.

Comparisons to FIR and/or radio observations
We compare the SFR and L TIR derived from SED fitting using data from rest-frame UV to optical with those indicated from independent observations in other wavelength described in Section 6.2.1.The result of comparisons is shown in Figure 11.In the figure, estimations from SED fitting and constraints from independent observations are shown by red and black symbols, respectively.
When Calzetti law is assumed as the dust attenuation law in SED fitting (top panels), L TIR seems to be overestimated.This indicates that using Calzetti law with SED fitting overestimates the amount of dust attenuation, which results in a larger value of the intrinsic rest-frame UV luminosity, possibly overestimation of the SFRF as well.
On the contrary, when SMC law is assumed (middle panels), the total FIR luminosity predicted from SED fitting is in a good agreement with the observations (middle right panel), while the SFR predicted from SED fitting using SMC law is in tension with that indicated from the radio observation (middle left panel).This agreement in the total FIR luminosity is consistent with the previous results from recent ALMA observations in the high-z universe based on analyses of the IRX-β relation (e.g., Fudamoto et al. 2020;Schouws et al. 2021).It is worth noting that these previous approaches based on the IRX-β relation use only rest-frame UV and FIR observations and assume a simplified spectrum of galaxies by a power-law function in the rest-frame UV wavelength range, while the comparison in this work is a consequence from rest-frame UV, optical, and FIR observations by assuming more realistic spectra of galaxies.Thus, this comparison further supports these recent studies that claim the SMC-like law is better representation for high-z galaxies.
Lastly, with the Calzetti+SMC estimation (bottom panels), the SFR estimated from SED fitting and radio observation agree with each other (bottom left panel), while the total FIR luminosity is overestimated with SED fitting (bottom right panel).
We have to keep it in mind that, in calculating the SFR radio using Equation (11), we do not consider the contribution from thermal radio emissions but regard all the radio flux observed with VLA as the non-thermal emission.The contribution of thermal emission to the total radio luminosity at ν rest = 9.3 GHz is roughly ∼ 25% in the local starburst galaxy M82 (Condon 1992).Inferring from this example, the non-thermal radio luminosity L NT may get smaller than the total radio luminosity at ν rest = 9.3 GHz observed with VLA by a factor of ∼ 0.75, which results in decreasing SFR radio (black squares in left panels) by the same factor.However, this displacement (∼ 0.1 dex) affects the SFR radio only slightly and changes the SFR radio within its uncertainty, thus this effect is negligible.
Considering these results of comparisons, estimations using SMC law as the dust attenuation law in SED fitting is consistent with independent observations.However, estimations using Calzetti law or Calzetti+SMC seem to be in tension, and their resulting SFRF (and SFRD) may be overestimated.Since the comparison is only for part of the final sample of this work, a further investigation is needed to examine more robustly including much deeper FIR observations and enlarging the sample size.

AGN contamination
In this work, we eliminate AGNs by removing X-ray sources brighter than ≳ 1 × 10 43 erg s −1 (c.f., sample selection in §2.3), thus the lower-luminousity AGNs may contaminate our sample.It should be useful to roughly estimate the effect of these AGNs that may be remained in the sample.Considering typical SEDs of AGNs (e.g., Elvis et al. 1994;Koratkar & Blaes 1999;Nemmen et al. 2014), the H 160 -band magnitude of an AGN at z ∼ 6 with X-ray luminosity of ∼ 1×10 43 erg s −1 is ∼ 27 mag, thus AGNs fainter than this can be detected if they locate in the deepest survey region (c.f., see the upper panel in Figure 1).The X-ray luminosity limit corresponds to a UV absolute magnitude limit of M UV ∼ −19 mag.This means that our sample after the removal of Xray sources (phot-z sample and final sample) is expected not to contain AGNs brighter than M UV ∼ −19 mag, and therefore only the faintest bin in the bottom panel of Figure 9 may be contaminated by low-luminousity AGNs.
Nevertheless, it is also meaningful to estimate the effect of AGNs using AGN LFs at z ∼ 6 without assuming the AGN SEDs, because the SED templates are derived from low-z AGNs and the SED may be different for AGNs at z ∼ 6. Matsuoka et al. (2018) derived an AGN LF at z ∼ 6 over the range of M UV from −30 to −22 mag.The number density of AGNs was estimated to be ∼ 2 × 10 −8 Mpc −3 mag −1 even at their most abundant and faintest bins of M UV = −22 or −23 mag.This number density is less than 1% of that obtained for the sample in this work at the magnitude bin (Figure 9).This magnitude bin of M UV ∼ −22-−23 mag is the brightest bin for our sample, and the fractions of AGNs in the fainter bins where most galaxies in our sample are included are expected to be several orders of magnitude smaller than 1%, if we suppose the extrapolated fainter part of the AGN LF (Figure 12 by Matsuoka et al. (2018)).Thus, the effect of AGN contamination is considered to be negligible.
Indeed, we can estimate the expected number of AGNs in each magnitude bin in Figure 9 by extrapolating the AGN LF by Matsuoka et al. (2018).We calculate the maximum effective volume for the search of contaminating AGNs by assuming incompleteness correction factors in Equation ( 5) to be unity.With this effective volume and the extrapolated AGN LF, the expected number of AGNs is ∼ 0.02 at the bin of M UV = −22.5, and ∼ 0.07 at M UV = −17.5.Taking the summation of these expected numbers suggests the maximum expected number of AGNs contaminating the phot-z sample is ∼ 0.25.We have to note that this estimation does not consider the effect of X-ray luminosity cut.The low-luminousity AGNs which may remain in the sample even after the Xray luminosity cut should be included in the faint magnitude bins, and the effective volume for the search of these faint sources is smaller than that we used here.Therefore, the expected number of 0.25 obtained here is an upper limit, and the actual number must be much smaller than this.By considering the sample size of the phot-z sample (∼ 550), the effect of AGNs on the resulting SFRF is expected to be negligible.

Incompleteness of the lowest-SFR bin
As discussed in Section 4.2 and in Figure 5, the SFRF obtained in this work shows an excess as compared to that from dust-corrected rest-frame UV LFs.Because not only the highest SFR bin(s) whose uncertainty is large but also the lower limit at the bin of SFR = 10 1.5−2.5 M ⊙ yr −1 suggests the excess, we conclude the excess seems to be rather robust.However, given that the incompleteness at the lowest-SFR bin of SFR = 10 1.5−2.5 M ⊙ yr −1 is based on the upper limit on rest-frame Hα EW < 4000 Å (c.f., Section 4.1) and this maximum value is achievable only with an extreme stellar population such as Pop III stars, the incompleteness at this bin may be less reliable.
To see if the lowest-SFR bin is robustly incomplete, we examine how the SFR completeness limit changes with different assumptions of the maximum EW(Hα) value.As shown in Section 4.1, we use the relation between the 5.8µm band magnitude and SFR that can be obtained when an EW(Hα) value is assumed, and estimate the SFR completeness limit that corresponds to the 5.8µm band detection limit (∼ 24.5 mag).In Section 4.1, we used the most extreme case of EW(Hα)= 4000 Å as the upper limit to obtain the completeness limit of SFR ∼ 10 2.5 M ⊙ yr −1 .Here, we examine how small the maximum EW(Hα) should be to make the lowest-SFR bin complete.We obtain that the SFR completeness limit can be as small as SFR = 10 1.5 M ⊙ yr −1 and the lowest-SFR bin can be complete only when we assume the maximum EW(Hα) is ∼ 500 Å.However, star-formging galaxies (SFGs) at this redshift are expected to have EW(Hα)∼ 500 Å on average (e.g., Faisst et al. 2016), and non-negligible fraction of SFGs can have EW(Hα)> 500 Å.This indicates that the SFR bin cannot be complete, and thus the arguments in Section 4.2 should be rather robust.

SUMMARY
In the high-z (z ≳ 4) universe, the contribution of dust-obscured star formation to the total star formation activity is still unclear, and an investigation of the SFRF is desired in an independent way of the previous approaches based on rest-frame UV observations or FIR observations.Very recently, investigations at z ≳ 4 with an approach using not only rest-frame UV but also optical data including Hα emission line have been reported.In this work, we derived an SFRF at z ∼ 5.8 with SED fitting using data from rest-frame UV to optical wavelength of galaxies in the CANDELS GOODS-South and -North fields.In deriving the SFRF, we extensively examined the assumptions on the SED modeling including various SFHs, dust attenuation law, and a two-component model.We obtained three resulting SFRFs by assuming the Calzetti attenuation law, the SMC law, and a combination of Calzetti and SMC laws.We also examined the redshift evolution of the SFRF and SFRD at z ≳ 4 obtained with the independent approach using rest-frame optical emissions.Our main results are as follows: 1.The resulting SFRFs at z ∼ 5.8 show an excess compared to those estimated from rest-frame UV observations but are roughly consistent with those estimated from FIR observations, although the SFRF from FIR observations is actually measured only at the brightest region and obtained with a large extrapolation down to the fainter region (Figure 5).The excess in the resulting SFRFs is originated from the difference of the estimated amount of the dust attenuation (Figure 6).This suggests that the contribution from the dustobscured intensively star-forming galaxies to the total star formation activities at z ∼ 6 is underestimated with the previous approach based only on the rest-frame UV observations (Section 4.2).
2. The SFRF is parameterized by assuming the Schechter functional form.The low-SFR end slope and characteristic number density are fixed and best-fit SFR * is estimated (Table 2).Together with SFRFs at z ≳ 4 in literature estimated with the approach using rest-frame optical emissions, the characteristic SFR (SFR ⋆ ) is roughly constant at z ∼ 2 to z ∼ 6 and may decrease above z ∼ 6.Since the parameter SFR ⋆ is sensitive to the high-SFR end of the SFRF, this suggests an early growth of the high-SFR end of the SFRF in the epoch of reionization, and the high-SFR end at z ∼ 4-6 is almost comparable to that at z ∼ 2 when the star formation activities reaches its peak (Section 5.1 and Figure 7).
3. Although the SFRFs obtained in this work are complete only in the high SFR region (SFR > 10 2.5 M ⊙ yr −1 ), a lower limit on the SFRD is obtained simply by taking the summation of SFRs.
From the estimation by assuming Calzetti law or Calzetti+SMC estimation, the previous measurements of the SFRD at z ∼ 6 based on rest-frame UV observations are marginally consistent or violate the lower limit obtained in this work, while the measurements based on FIR observations are acceptable and more likely.On the contrary, from the estimation by assuming SMC law, the lower limit on the SFRD is not tight enough to determine if the SFRD is underestimated or not with the previous measurements using rest-frame UV observations (Section 5.2 and Figure 8).
We thank the anonymous referee for useful comments and suggestions to improve this paper.In this appendix, we show the physical properties of all the galaxies in the final sample including their SEDs and the best estimations obtained in Section 3.3.
Figure 12 shows the SEDs of galaxies in the final sample.We also show one of the best fit SEDs derived in Section 3.1 and 3.2.Among the best-fitting SEDs, we only overplot the model SED that gives the minimum BIC.
In Table 3, all the physical parameters obtained in this work are presented for each of the galaxies in the final sample with Calzetti, SMC, and Calzetti+SMC estimation.In the physical parameter estimation, we also utilize the result of the two-component SED fitting (Section 3.2).As noted in Section 3.2, the best-fit parameters of SFR, M ⋆ , and L TIR in the two-component model are estimations for the whole system that includes both the old stellar population and the young star-forming population, but age and A V are estimations for the young star-forming population9 .Therefore, the best estimations for age and A V , which are obtained using this two-component modeling, are not estimations for the whole system.In particular, the best estimations of age obtained in this work do not denote time intervals since the onset of the first star formation episode that the galaxy experienced, but should be regarded as the time intervals since the onset of the last star formation episode.

B. AN EXAMPLE OF PHYSICAL PARAMETER ESTIMATIONS
In this appendix, we show an example of physical parameter estimation described in Section 3.3 using a galaxy in the final sample (1411-GS).Figure 13 shows the distribution of the best-fitting physical parameter versus the weights (exp(−BIC/2)) derived by SED fittings in Section 3.1 and 3.2 under each assumption of dust attenuation laws.Panels in the top row contain both results of SED fitting with assuming Calzetti and SMC law.Those in the middle (bottom) row contain only the results of SED fitting with assuming Calzetti (SMC) law.Note that, in these panels, the results of two-component model SED fitting that do not meet the criteria I or II (Section 3.2) are not shown, and are not used in the following analysis.The red solid line in the figure shows the weighted mean value, and the shaded region presents their uncertainties obtained from the weighted standard deviations.
Figure SEDs of the galaxies in the final sample and best-fitting model SEDs.In each panel, the black points with error bars represent the observed SED, the blue solid line shows the best fit model spectrum, and red crosses show the best fit model SED.Among the 28 best-fitting models for each galaxy, in this figure, we only show the best-fitting model that gives minimum BIC (see also Section 3.3).For galaxies that the two-component model give the minimum BIC, we also show the best-fitting spectrum of the old stellar population and young star-forming population with red and blue dotted lines, respectively.At the top left in each panel, we also show the name of the galaxy and the SFR with Calzetti+SMC estimation.

C. COMPARISONS OF RESULTS WITH NON-PARAMETRIC SFH APPROACH
In this paper, we calculated the weighted means of best-fit SED models to obtain the best estimations for physical parameters, taking various SFHs into account.In this appendix, we utilize a public SED fitting code using a non-parametric SFH approach, dense basis (Iyer & Gawiser 2017;Iyer et al. 2019), and examine if the main results in this paper persist when dense basis method is adopted.We perform SED fittings using dense basis on galaxies in the final sample to obtain the PDFs of each physical parameter including the SFR, and followed the same procedure to calculate the SFRF, best-fitting SFR ⋆ , and lower limit of SFRD from the PDFs obtained with dense basis.For a fair comparison, we assumed flat prior distribution in dense basis configuration and used Calzetti attenuation law since the Calzetti law is implemented in dense basis.
Figure 14 shows the result of comparisons.The resulting SFRF at z ∼ 5.8 (left panel) agrees with each other within their uncertainties, although the number density in the highest SFR bin may be overestimated with the weighted-mean estimations.This agreement results in a persistence of estimations of the bestfitting SFR ⋆ value (middle panel) and lower limit on SFRD (right panel).The best-fitting SFR ⋆ value is log 10 (SFR ⋆ ) = 2.36 +0.12 −0.90 and the lower limit on the SFRD is log 10 (ρ lim SFR ) > −1.96 with weighted-mean method, while they are log 10 (SFR ⋆ ) = 2.63 +0.21  −0.69 and log 10 (ρ lim SFR ) > −1.91 with dense basis method.D. SFR ⋆ CALCULATION IN LITERATURE The SFR ⋆ by Bell et al. (2007), S16, and A21 in Figure 7 are given in the literature, but other previous studies do not give SFR ⋆ .Thus we derive SFR ⋆ by ourselves as follows.Sobral et al. (2013) gave the best-fit Schechter parameters for (dust-corrected) Hα LFs.Since the Hα luminosity can be converted into the SFR, we simply convert the characteristic luminosity (L ⋆ Hα ) to the characteristic SFR (SFR ⋆ ) using the following equation by assuming Chabrier03 IMF: L Hα = 2.1 × 10 41 SFR M ⊙ yr −1 erg s −1 (D1) Magnelli et al. (2011) measured a total FIR LF, but SFR ⋆ for the SFRF is not derived.We thus convert the total FIR LF into an SFRF using the following equation and obtain the best-fit Schechter parameters including SFR ⋆ : L TIR = 1.0 × 10 10 SFR M ⊙ yr −1 L ⊙ (D2) Reddy et al. (2008) measured rest-frame UV, Hα, and IR LFs, but the total SFRF is not derived.We use the total IR LF and obtain the best-fit Schechter parameters as well.
A22 derived a constraint on the Hα LF via Schechter function.Specifically, they obtained a forbidden region in the parameter space of the Hα luminosity versus the number density, and derived a constraint on the Schechter parameters of the Hα LF as follows.They first assumed a relation between dust-uncorrected (observed) rest-frame UV luminosity (L UV,obs ) and dust-corrected (intrinsic) rest-frame UV luminosity (L UV,int ) as and this correction factor η depends on the rest-frame UV luminosity: With this assumption, the characteristic luminosity of the dust-corrected rest-frame UV LF (L ⋆ int ) can be expressed using that of dust-uncorrected rest-frame UV LF (L ⋆ ): A22 derived an upper limit on the Hα luminosity function by converting this dust-corrected rest-frame UV LF into a Hα LF.In this work, we use the same L ⋆ as by A22 and their resulting parameters (a, b) to calculate  Figure 8) except for filled black circles or squares, which show the result using dense basis method.In the middle and right panel, results in this paper at z ∼ 5.8 is displaced by ∆z = ±0.1 for clarity.
the upper limit on L ⋆ int by Equation (D5), and convert it into an upper limit on the SFR ⋆ using the following equation: L UV = 1.3 × 10 28 SFR M ⊙ yr −1 erg s −1 Hz −1 (D6)

E. CONVERSION FACTOR BETWEEN THE RADIO LUMINOSITY AND SFR
In this appendix, we derive the conversion factor between the radio luminosity and the SFR (Equation 11).Condon (1992) derived a conversion factor from (nonthermal) radio luminosity L NT to the SFR for massive stars (M ⋆ > 5M ⊙ ) assuming extended Miller-Scalo IMF.

Figure 1 .
Figure 1.Top: detection rates at the H160 band in GOODS-S and -N fields.The different linestyle shows different fields, and different color shows different regions in a field.Lower left: the distribution of apparent magnitudes at H160 versus 5.8µm band of galaxies in the final sample.The circles (squares) represent galaxies in the GOODS-S (-N) field.Lower right: normalized cumulative histogram of the 3σ limiting magnitudes in the 5.8µm (Ch3) band, which effectively shows the detection rate in the band.(see text for details).

Figure 2 .
Figure 2. Postage stamps (16 ′′ × 16 ′′ ) of galaxies in the final sample.From left to right, images in ACS F435W, F606W, F814W, WFC3 F105W, 125W, F160W, IRAC 3.6 µm, 4.5 µm, 5.8 µm, and 8.0 µm bands are shown for each galaxy (top to bottom).The HST and Spitzer images are taken from publicly available science images created by the Hubble Legacy Fields project and by the GOODS Reionization Era wide-Area Treasury from Spitzer project (Stefanon et al. 2021), respectively.

Figure 3 .
Figure3.The stellar mass versus SFR diagram for galaxies in the final sample.The best estimations of the stellar mass and SFR obtained in Section 3.3 using Calzetti+SMC, Calzetti, SMC estimation are used to plot in the left, middle, and right panel, respectively.The black solid line denotes the MS relation at z ∼ 5.8 predicted from the redshift evolution of the MS obtained bySchreiber et al. (2015).The thick (thin) shaded region show a typical scatter of ±0.3 (±0.5) dex from the MS.

Figure 4 .
Figure4.The SFRs versus 5.8µm band magnitudes for galaxies in the final sample.Left, middle, and right panel shows the diagram with Calzetti+SMC, Calzetti, and SMC estimations, respectively.The black dotted line represents a guideline that corresponds to the upper limit on the SFR obtained from a maximum value of rest-frame Hα EW (EW(Hα) = 4000 Å), and the region above this line is expected to be forbidden (see the text for details).

Figure 5 .Figure 6 .
Figure 5.The SFRF at z ∼ 5.8 obtained by SED fitting in this work (red circles).Left, middle, and right panel shows the SFRF using the results with Calzetti+SMC, Calzetti, and SMC estimation, respectively.The vertical error bars are 1σ Poisson uncertainties by Gehrels (1986).The blue and green solid line show SFRFs estimated with rest-frame UV and FIR observations at z ∼ 6, respectively.The black dashed line shows the SFRF at z ∼ 4.5 derived in a similar way to this work.

Figure 7 .
Figure 7.The redshift evolution of the characteristic SFR SFR * .At z ≳ 3, investigations based only on rest-frame UV observations are shown by blue symbols, and those using rest-frame UV to optical observations are shown by red symbols including the result in this work.

Figure 8 .
Figure 8.The SFRD estimation and comparisons with the literature.The lower limits obtained in this work using results by Calzetti+SMC (red filled circle), Calzetti (red filled square), and SMC (red filled diamond) estimation are shown.Recent SFRD estimations using various probes are also shown for comparison.For clarity, our results with Calzetti (SMC) law is shifted by ∆z = 0.2 (−0.2), and the result at z = 5.25 by Gruppioni et al. (2020) is shifted by ∆z = −0.15.

Figure 9 .
Figure 9.The (dust-uncorrected) rest-frame UV LFs derived with the final (top) and phot-z (bottom) sample (black circles).The solid and dotted blue lines show the best-fit UV LF by Harikane et al. (2021) andBouwens et al. (2015), respectively.In the bottom panel, we also show the UV LFs with galaxies in the phot-z sample in GOODS-South (North) field by red (green) squares.

Figure 10 .
Figure10.The GSMF measurement at z ∼ 6 and comparison to the previous works.In this figure, only the GSMF obtained using the result of Calzetti estimations is shown, but the difference of dust attenuation law has negligible effects on the resulting GSMF.

Figure 13 .
Figure 13.Distribution of best-fitting physical parameters and the weights (exp(−BIC/2)) derived from SED fitting in Section 3.1 and 3.2 for a galaxy in the final sample (1411-GS).The red (blue) symbols represent the results by assuming Calzetti (SMC) attenuation law for the dust attenuation, and circles (squares) represent the results with one-(two-)component SED fitting.For one-component SED fitting given in Section 3.1, all the results with seven SFHs are shown, while for two-component SED fitting, results that do not meet either of the criteria I or II in Section 3.2 are not shown.Top row contains both results by assuming Calzetti and SMC law, and middle and bottom row show results with Calzetti and SMC law, respectively.From left to right, physical parameters of age, AV , stellar mass, SFR, and the total infrared luminosity are shown.The vertical red solid lines and shaded regions show the best estimation (weighted mean) of each physical parameter and their uncertainties (weighted standard deviation), respectively.

Figure 14 .
Figure14.Comparisons of results in this paper when the weighted-mean estimation and a non-parametric SFH SED fitting code (dense basis) is used.Left, middle, and right panel show the resulting SFRF, SFR * , and SFRD estimation at z ∼ 5.8 in this paper, respectively.Symbols and lines are the same as in corresponding figures (left: Figure5, middle: Figure7, and right: Figure8) except for filled black circles or squares, which show the result using dense basis method.In the middle and right panel, results in this paper at z ∼ 5.8 is displaced by ∆z = ±0.1 for clarity.

Table 1 .
Galaxies in the final sample

Table 2 .
The Best-fit Characteristic SFR and the Lower Limit on the SFRD Y.A. is supported by a Research Fellowship for Young Scientists from the Japan Society of the Promotion of Science (JSPS).K.O. is supported by JSPS KAKENHI Grant Number JP19K03928.This work is based on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI).STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.