A Large Catalogue of Multi-wavelength GRB Afterglows I: Color Evolution And Its Physical Implication

The spectrum of gamma-ray burst (GRB) afterglows can be studied with color indices. Here we present a large comprehensive catalogue of 70 GRBs with multi-wavelength optical transient data on which we perform a systematic study to find the temporal evolution of color indices. We categorize them into two samples based on how well the color indices are evaluated. The Golden sample includes 25 bursts mostly observed by GROND, and the Silver sample includes 45 bursts observed by other telescopes. For the Golden sample, we find that 95\% of the color indices do not vary over time. However, the color indices do vary during short periods in most bursts. The observed variations are consistent with effects of (i) the cooling frequency crossing the studied energy bands in a wind medium (43\%) and in a constant density medium (30\%), (ii) early dust extinction (12\%), (iii) transition from reverse shock to forward shock emission (5\%), or (iv) an emergent supernova emission (10\%). We also study the evolutionary properties of the mean color indices for different emission episodes. We find that 86\% of the color indices in the 70 bursts show constancy between consecutive ones. The color index variations occur mainly during the late GRB-SN bump, the flare and early reversed-shock emission components. We further perform a statistical analysis of various observational properties and model parameters (spectral index $\beta_{o}^{\rm CI}$, electron spectral indices $p^{\rm CI}$, etc.) using color indices. Overall, we conclude that $\sim$ 90\% of colors are constant in time and can be accounted for by the simplest external forward shock model, while the varying color indices call for more detailed modeling.


INTRODUCTION
The first observation of gamma-ray burst (GRB) optical afterglow was that of the BeppoSAX-detected GRB 970228 (van Paradijs et al. 1997). Since the launch of the Swift satellite more than 10 yr ago (Gehrels et al. 2004), many ground-based optical telescopes with increasing sensitivity have accumulated a rich collection of optical afterglows. In recent years, optical afterglows have been studied extensively either in terms of either their multiband light curves (e.g., Panaitescu & Kumar 2002;Huang et al. 2006;Zeh et al. 2006;Kann et al. 2010Kann et al. , 2011Li et al. 2012;Liang et al. 2013;Wang et al. 2015;Roming et al. 2017, and references therein), or their spectral energy distributions (SEDs) derived from simultaneous multiband photometry (e.g., Stratta et al. 2004;Starling et al. 2007;Kann et al. 2006Kann et al. , 2010Kann et al. , 2011Roming et al. 2017). The latter, if studied in a broad bandpass, provides tight constraints on the radiation mechanism and circumburst environment of GRBs. Given the rapid progress made in the detection of afterglows in the Swift era, unforeseen temporal and spectral features were discovered, which challenge our understanding of GRB afterglow physics (Zhang 2007;Gehrels et al. 2009, and references therein).
The temporal and spectral properties of optical transients can be studied by color indices (hereinafter CIs), defined as the magnitude difference between two filters. They can be used to study the SED with a good temporal resolution, even when high-resolution spectra are not available. For instance, CIs can resolve small variations in the spectral profiles of optical transients (Šimon et al. 2001;Šimon et al. 2013). They are also helpful in identifying different radiation mechanisms or progenitors in a single burst. As an example, CIs were used as an indicator of the underlying Type Ib/c supernovae (SNe; Simon et al. 2004).
Theoretically, the conventional afterglow model has been discussed by many authors (e.g., Mészáros & Rees 1997;Sari et al. 1998;Huang & Cheng 2003), and the relationship between temporal and spectral indices (the so-called closure relation) has been extensively reviewed in Zhang & Mészáros (2004) and Gao et al. (2013). According to this model, which assumes that synchrotron radiation dominates the afterglow emission, the color index should not change with time. Instead, if a change is observed, there may be several possible reasons (Melandri et al. 2017), for instance: (i) the cooling frequency crosses the studied energy bands (e.g., Filgas et al. 2011a). (ii) an additional emission component emerges, e.g., the SN counterpart observed in some GRBs, or the mixture of different emission stages, such as forward and reverse shock (Kobayashi & Zhang 2003) in the context of the fireball model (Piran 2004). (iii) after being destroyed by the initial intense radiation, host extinction could be changed owing to dust photodestruction (Perna & Lazzati 2002;Morgan et al. 2014). Some other possible reasons, such as a structured jet, could also produce the change in the CIs.
In this paper, we have made an extensive effort to collect publicly available optical and near-infrared photometric data. We restricted the study to 70 bursts with high-quality multiband observations, which were detected during both the pre-Swift and Swift eras. Our goal is to study the temporal variability of CIs and explore its physical implications. In a following paper, we will present a detailed time-resolved analysis of afterglow spectra using X-ray and optical data, which will be used to study the extinction curve and explore the physics in the circumburst medium in the host galaxy (Li et al 2017, in preparation).
This paper is organized as follows. The sample selection and data analysis are presented in Section 2. The statistical properties of CIs are summarized in Section 3. The evolution of the CIs is presented in Section 4. The theoretical implications are discussed in Section 5. Finally, conclusions are provided in Section 6.

SAMPLE SELECTION, CORRECTIONS, AND SCALING
We extensively compiled as many optical/near-IR (NIR) photometric data as possible, detected by groundbased telescopes from both pre-Swift and Swift-era afterglows. The bursts in our sample were observed before 2016 1 . We collected observational data from published papers and the Gamma-ray Burst Coordinates Networks Circular Service (GCN) 2 . We initially selected GRBs with known redshift that have high-quality data in at least three optical bands. We then restricted the study to the GRBs that were observed in optical wavelength for at least 1 hr, facilitating the study of CIs temporal evolution. A total of 70 bursts satisfy these criteria (from 1997 February 28 to 2015 December). To obtain a trustworthy data set, the data sources should be well calibrated. We made corrections to the data as follows.
We studied the temporal properties of the CIs in the rest frame, t rest = t obs /(1 + z). For some specific cases (GRB 970508, GRB 990510, GRB 990712, GRB 060614, GRB 060908, GRB 090426, GRB 091127, GRB 100621A, GRB 100814A, GRB 110918A, GRB 120422A, GRB 120729A, GRB 130702A, GRB 130925A), we subtracted the flux contribution at very late epochs that we interpret as coming from the host galaxy. This contribution is determined by identifying constant flux values at a late time, ∼ 10 6 s after the GRB trigger.
We calculated the Galactic extinction correction from the reddening map presented in Schlafly & Finkbeiner (2011) for optical and NIR magnitudes, and we assume 3 R v =3.1.
The observed spectra are k-corrected, where k is defined as m T = m O −k (e.g., Oke & Sandage 1968;Peterson 1997), with k = 2.5(β o − 1) log(1 + z). Here m O and m T are the observed and true magnitudes, respectively. β o is the spectral index 4 of the optical afterglow, and z is the redshift.

Corrections for Host Dust Extinction
The extinction due to dust from the host galaxy can significantly affect the observations. For wavelength λ i , the dust extinction is defined as (1) The observed SED of the optical afterglow is modeled by an absorbed power law in frequency (Kann et al. 2006) where τ (ν host ) = 1 1.086 in which β o is the intrinsic power-law slope of the SED, F 0 is the normalization constant, η(ν host ) = A host λ /A host V is given by the extinction model assumed for the GRB host galaxy, and A host λ is the host galaxy extinction correction (known only with large uncertainties). We do not carry out spectral fitting in this paper. Instead, the host galaxy extinction A host V , the spectral index β o , and the redshift z are acquired from the literature, and their values are given in Table 1.
Specifically, A host V is chosen by the following three steps. First, we selected the values obtained by Kann et al. (2006) and Kann et al. (2010) corresponding to the Small Magellanic Cloud (SMC) dust model. Second, when they are not available in Kann et al. (2006) and Kann et al. (2010), the extinction parameters are retrieved from other papers (given as a reference in Table 1) from the SMC dust model. This model is our first choice since it is shown to best describe dust in GRB environments (Kann et al. 2010). When the information is not available for this model, the value of A host V is obtained from either from the Milky Way (MW) or the large Magellanic Cloud (LMC) models. Finally, for a few cases, when the value of A host V is unavailable for any model, we use the value obtained by fitting a Gaussian to the extinction distribution (see Figure 1), log 10 (A host V )=-0.82±0.41; here -0.82 and 0.41 are the mean value and standard deviation of the Gaussian fit 5 , respectively. After obtaining the value of A host V , we derived the extinction values in any optical band, A host λ , in the GRB host galaxies following Pei (1992). Figure 2 illustrates the spectral behavior of A host λ /A host V for the three dust models and displays the position of the filters considered in this paper.

Light-curve Fitting
In this section, we describe the fits to the GRB light curves in our 70 bursts. The purpose of this procedure is twofold. First, it allows us to identify different emission phases. Second, the results are used to interpolate the magnitude, in order to compute the CIs, when there are no simultaneous multi-wavelength observations available.
The observed magnitude m is obtained from the flux F as where F 0 is the flux of an object with magnitude zero. Note that in this paper we adopt the AB magnitude system. After taking into account all the correction factors (our Galaxy and the host galaxy extinction corrections, and the spectral k-correction), we fit the light curves in each band separately with a model of multiple components (Li et al. 2012). The basic model is either a single power law (SPL) or a smoothly broken power law (BKPL) for the flux, resulting in 6 : m T (t) = −2.5log 10 (10 −0.4mc t −α + 10 −0.4mG + 10 −0.4mK +10 −0.4m host + 10 −0.4mGF ), m T (t) = −2.5log 10 {10 −0.4mc [(t/t b,1 ) α1ω1 + (t/t b,1 ) α2ω1 ] −1/ω1 +10 −0.4mG + 10 −0.4mK + 10 −0.4m host + 10 −0.4mGF }, here α, α 1 , and α 2 are the temporal slopes, t b1 is the break time, and ω 1 describes the sharpness of the break (the smaller the value, the smoother the break). For most GRBs, we fix this last parameter to ω 1 = 3. For a few cases (GRB 030226, GRB 050922C, GRB 061007, GRB 080603A) that have a smoother break, we fix it to ω 1 = 1, which improves the fit. Equation (2.2) can be used for both for the light curves with a break and those who exhibit an increasing behavior (e.g., afterglow onset indicated by the dashed line in Fig.3b). In addition, m c is the magnitude constant value, m G is the extinction magnitude for our Galaxy, m K is the spectral k-correction, m host is the extinction magnitude for the host galaxy, and m GF is the possible contribution from the host galaxy at late time (if identified). Note that m G , m K , and m host are known before fitting. We notice that for two bursts (GRB 061126 and GRB 080319B), the optical-afterglow light curves have a steeper decay slope α 1 in the pre-break segment (which might originate from the reverse-shock emission) than the decay slope α 2 in the post-break segment (possibly from the normal decay dominated by the external forward shock). However, due to the limitation of the BKPL model, it is impossible to fit the light curves whose pre-break slope is steeper than the postbreak slope. In this case, we adopted the BKPL model function with a negative sharpness of the break ω=-3 (Fig.3c).
A double-BKPL light curve (Fig.3d) is also expected in some afterglow models. For example, it is theoretically expected that the afterglow light curves may have a shallow segment owing to energy injection at an early time, which then changes into a normal-decay segment when the energy injection is over, and finally steepens owing to a jet break (Zhang et al. 2006). We therefore consider a smooth triple-power-law function (TPL) to fit the light curves (e.g., Liang et al. 2008;Li et al. 2012): +10 −0.4mG + 10 −0.4mK + 10 −0.4m host + 10 −0.4mGF }.
where ω 2 is the sharpness parameter at the second break time t b,2 .
The fits are performed with the IDL routine "mpfitfun.pro" 7 (Markwardt 2009), which uses the Levenberg-7 http://purl.com/net/mpfit Marquardt algorithm to achieve minimization. To select the best model, we compare the reduced χ 2 values and choose the one that has a more reasonable value (close to 1). For example, the χ 2 r of GRB 990510 (R band) is 31/39 (a little less than 1) for the BKPL model and 169/36 (much larger than 1) for the SPL model. Therefore, we adopt the BKPL as the best model for GRB 990510. Second, if both models have a χ 2 r close to 1, our principle is to choose the simpler one (fewer parameters). For instance, for GRB 060908, a χ 2 r is equal to 38/49 (∼ 1) for the BKPL model and 52/46 (∼ 1) for the SPL model. In this case, we choose the SPL model.
For multicomponent light curves (see §2.3), we introduce the minimum number of components by eye inspection of the temporal features. If the χ 2 r is still much larger than 1, we continue to add more components and redo the fit, until the χ 2 r becomes close to 1 (Li et al. 2012). For instance, the value of the χ 2 r of GRB 071025 (J band) for the BKPL is 474/39 (much larger than 1), while it is 58/34 (close to 1) for the double-BKPL model. As a result, the double-BKPL model is selected as the best model for this burst. Examples of light-curve fitting with various models or their composition are shown in Figure 4.

Component Identification
Optical afterglows have complicated light curves with up to eight emission episodes (Li et al. 2012;Li et al. 2015;Liang et al. 2006;Liang et al. 2013;Nardini et al. 2006;Kann et al. 2006Kann et al. , 2010Kann et al. , 2011Panaitescu & Vestrand 2008, 2011. A synthetic optical-afterglow light curve is presented in Li et al. (2012). To study the spectral behaviors in these different emission phases, we investigate the temporal evolution of their mean CIs. Here we explain which components are considered and how they are identified based on their temporal and spectral features.
Ia: the prompt optical flares (Prompt-Optical); in the very early time of some bursts when the prompt GRB emission is still going on, a highly variable optical emission component may be observed as in GRB 080319B 8 . They are selected such that the observation time is smaller than T 90 and the temporal slope 9 α >2.0.
Ib: the early optical flares (Reverse-Shock); in a few cases, the early light curves have steep slopes, which likely indicate an early reverse-shock emission component, such as the early phase of GRB 061126. They were selected such that their typical temporal index α ∼ 1.7 and the peak time t p is a few hundreds seconds.
II: an early shallow-decay component (Energy-Injection); the afterglow light curve might show an initial shallow-decay segment followed by a normaldecay/post-jet-break segment, which is likely due to energy injection from a long-lasting spinning-down central engine or piling up of flare materials into the blast wave (Liang et al. 2007;Li et al. 2012;Li et al. 2015). This component is described as the energy injection phase. They are identified by their temporal index α Shallow <α Normal with a typical value α Shallow ∼ 0.5 and with a typical value of break time t b ∼ 10 4 s. Here α Shallow is the temporal index during the shallow-decay segment and α Normal is the temporal index during the normal-decay phase.
III: the standard afterglow component (Onset/Normal-Decay); the light curves sometimes have an early onset rising segment followed by a normal decay. In most of the cases, a lack of observations in the early time lead to only a single normal decay. Afterglow onsets are identified with an early smooth bump with a typical value of peak time t p of several hundreds of seconds, coupled to the following normal decay with α ∼ 1.2.
IV: the jet-break component (Jet-Break); the light curves break into a steeper decay. They are identified with α ∼ p ∼ 2.5 (Zhang et al. 2006) and a break time ∼ 10 5 s. Here p is the electron index.
V: the late optical flares (Flare); the light curves have prompt-like flares during the afterglow phase when the prompt emission is turned off. They indicate the latetime activities of the central engine. The late optical flares are characterized by a very sharp temporal index α > 2.0.
VI: the late re-brightening bumps (Re-Bump); the late bumps would emerge at late times, which is distinguished from the onset bump of the early afterglow. Both are likely involved in the jet component that produces the re-brightening bump, which seems to be onaxis and independent of the prompt emission jet component (Liang et al. 2013). They are described by a smooth bump around t p ∼ 10 5 seconds.
VII: the late SN bumps (SN-Bump); in some cases, the optical transient light curves at late time show an SN bump. They form a late smooth bump at t p ∼ 10 6 s. In addition, we checked that the GRB-SN associations were confirmed in the literature (see §5.2.3).
All the assigned emission components in our 70 bursts are presented in Table 2 and 3, and are marked with arrows in Figures 5 and 6.

Definition of Color Index Variation
Once the fits are obtained, one can derive the CIs. They are defined as the magnitude difference between two filters, CI = m(λ 2 ) − m(λ 1 ), where λ 2 < λ 1 is required so that the smaller the CI, the bluer (or hotter) the optical afterglow. In this paper, we uniquely focus on the CI between two adjacent bands 10 . CI can be affected by optical extinction, which can be characterised by the color excess (CE), defined by the difference between the observed color index (CI O ) and the intrinsic color index (CI T )

Definition of the Golden and Silver Samples
To compute the CIs, simultaneous observations or interpolated values are required. We also note that different telescopes have different photometric systems (i.e., different light filters and spectral response functions), which are either the common Johnson-Cousins UBVRI photometry system 11 or the Sloan Digital Sky Survey (SDSS) ugriz photometry system 12 . Instead of converting one system into another with empirical correlations (e.g., Jordi et al. 2006), we group the bursts by photometric systems, which are studied separately.
Therefore, we sort all GRBs into two samples based on the photometric system used for the observation. This also naturally separates bursts for which interpolation is required from those for which it is not.
• Golden sample: The 25 bursts have simultaneous observations in multiple filters in the optical/IR bands (g, r, i, z, J, H, Ks) for a wide time window. They are mostly observed by the Gamma-Ray Burst Optical/NIR Detector (GROND; Greiner et al. 2008), with some events contributed from observations of the RAPid Telescopes for Optical Response (RAP-TOR;Wren et al. 2010) and the Peters Automated IR Imaging Telescope (PAIRITEL; Bloom et al. 2006). These define our Golden sample. Obtaining the CI 13 for these bursts is straightforward. The afterglow multicolor light curves and the time evolution of CIs are displayed in Figure 5.
• Silver sample: The 45 remaining bursts in our sample do not have simultaneous observations in multiple filters but have instead multiple observations in different filters from different optical instruments at different times. For those bursts, we use the lightcurve fits ( §2.3) and interpolate the missing data points to obtain concurrent values in different filters (UVW2, UVM2, UVW1, U, B, V, R, I, J, H and K). This group of bursts defines the Silver sample. It has wider observed energy bands than the Golden sample. We try to include a maximum number of possible multiple data to represent the whole light curve. The multicolor light curves 14 and the time evolution of CIs for each burst are shown in Figure 6.
The most accurate CIs are obtained in cases where there are simultaneous observations in multiple filters. This is the case for the bursts in the Golden sample. For many of these bursts there are additional nonsimultaneous data as well, which are discarded in order to keep the sample pure. Three criteria are applied to consider whether to add the additional data sources into the Golden sample: (i) there is at least an order-ofmagnitude extension in the duration of observation; (ii) the number of data points should be greater than (or at least equal to) that in the original data source; (iii) at least three new bands have good observational data. Most cases do not satisfy these criteria. However, in two particular cases, GRB 090426 and GRB 130427A, there is a significant amount of additional data. We therefore use these cases to compare the results from the Golden and Silver samples.
For example, we add the R-band data to the multicolor light curves in the Golden sample for GRB 090426 ( Figure 5). It provides a much earlier observation that can be compared to the GROND griz observation.

Methods to Determine the Temporal Behavior of Color Indices
We further define sets of CIs. The first one contains time-resolved data (Data Set I). The second set is made of component-wise averaged CIs (Data Set II). Finally, 14 We uniformly added different numbers to make a fixed set of offsets to different bands. For the Golden sample: g: 6, r: 5, i: 4, z: 3, J: 2, H: 1, K: 0; For the Silver sample: UVW2: 12, UVM2: 9, UVW1: 8, UVM1: 7, U: 7, B: 5, V: 4, R: 3, I: 2, J: 1.5, H: 1, K: 0, u: 14, g: 11, r: 10, i: 9, z: 8. the mean CIs over the entire burst duration give Data Set III. Hereinafter all the analysis (such as the figures and the tables) will be marked with which sets of data are used.
According toŠimon et al. (2001), CIs should not be derived from fits to the light curve. They claimed that any fits to the data would distort the CIs. The Golden sample gives directly CIs directly while the Silver sample CIs are obtained with a fit to the light curve. Studying the difference between Golden and Silver samples allows us to check the consistency between our results and that ofŠimon et al. (2001).
To account for the missing data points, two methods can be employed: (i) interpolation of the light curve from the fits (applied to the Silver sample), and (ii) averaging of existing data points as inŠimon et al. (2001). Interpolating missing data points results in large uncertainties but has the advantage of preserving a maximum number of simultaneous CIs. If the afterglow emission is smooth enough, the uncertainties should not be too large. On the other hand, averaging data results in more accurate estimation of the CIs, but strongly reduces their number. As pointed out byŠimon et al. (2001), averaging is also not reliable when light curves have rapid changes. Interestingly, after comparing both methods, we find that the results are consistent. We note that our Set II, which is obtained by averaging the CIs on each afterglow component, is similar to averages in arbitrary time intervals, as inŠimon et al. (2001), except that, in our paper, the time intervals are imposed by expected physical variations of the light curves.
We investigate the correlation of the CIs from GRB 130427A using these two methods. This burst has a well-observed afterglow with a great number of data points. In Figure 7, we compare the CIs in the same observation duration between the Golden and the Silver samples. The CIs from the Golden sample are derived directly from RAPTOR (Vestrand et al. 2014) with simultaneous observations in multiple filters, while the CIs from the Silver sample are obtained by an interpolation method from the RATIR (Becerra et al. 2017), which in general do not have simultaneous observation, compared with the GROND observation in the corresponding period. We find that the data cluster around the equal line, which indicates that the results of both methods are consistent with each other.

Determination of Time Evolution
To determine whether the CI changes with time, we consider the following quantitative analysis. First, we calculate the mean color index m and the standard deviation σ s for each CI. We do this for Data Set I and II defined in Section §2.4.2. If a significant variation exists for the ith CIs, the criteria | m i − m |> σ mi + σ s should be met; here σ mi is the error on the ith CI. Note that both σ mi and σ s are positive. Therefore, we define two grades for the CI: • Constant: the CI is consistent with the overall mean value.
• Variable: the CI presents a significant variation. Table 4 summarizes the statistical properties of the CIs in our sample, which include the total number of data points, the distributed ranges, and their mean values. In Figure 8, we present all CIs as functions of time in one panel to show the global behavior. We find that more than 90% of the CIs distribute within [-1.0, 1.0] in the observed time intervals.

STATISTICAL PROPERTIES OF COLOR INDICES
We analyze the distributions of the CIs (Data Set I) for both the Golden and Silver samples. They are all well fitted with Gaussian functions. Figure 9 displays the distributions together with the best Gaussian fits and with the mean CI values, as well as their standard deviations. For the Golden sample (Fig.9a), we have g-r=0.15±0.14, r-i=0.09±0.10, i-z=0.13±0.11, J-H=0.23±0.10 and H-Ks=0.25±0.17, respectively. We find that the distributions 15 are consistent with each other with a mean value around 0.2. We note the presence of a tail at large values for the CI g-r. For the Silver sample (Fig.9b), we obtained UVM2-UVW1=0.16±0.27, UVM2-UVW1=-0.19±0.24, UVW1-U=0.93±0.39, U-B=0.09±0.23, B-V=0.23±0.23, V-R=0.15±0.16, R-I=0.22±0.07, I-J=0.17±0.08, J-H=0.28±0.14 and H-K=0.14±0.22, respectively. We find that the distributions are still consistent with each other with a mean value also around 0.2 for most of the CIs. This indicates that most of the CIs could be explained by the standard afterglow model with a single power-law spectral slope.
Note, however, that there are two inconsistent peaks in the Ultraviolet (UV) bands, with the UVM2-UVW1 CIs at -0.19 and UVW1-U CIs at 0.93. The difference between the distributions could originate from the method that the CIs are obtained for the Silver sample. Because of the fitting and interpolation procedures, they carry large uncertainties, not taken into account by the histograms. A second reason for this difference could be underestimation of the extinction corrections. This is because the extinction curve in the host galaxy 15 We note that these results are similar to the ones if only CIs with simultaneous observations for all five CIs are selected. is largely uncertain 16 and therefore the corrected data could deviate from the actual values. This implies that after making the correction to the data, there could also be an additional reddening. Further investigation of the extinction property for each individual burst will be presented in a following paper (Li et al 2017, in preparation). Stronger deviations are expected in high-energy bands in the Silver sample, tentatively explaining the possible reddening in those bands (Fig.9b). Also, very few GRB host galaxies show MW/LMC-like 2175 • A features. However, the redshifts are distributed over a wide range of values, so even if such features exist, they would appear in different observer-frame bands. Therefore, the effect would be completely diluted.
To test whether there exists an additional reddening, we analyze the SED, using the average CIs (Data Set III, regardless GRBs) and assuming a typical magnitude for the r-band (Golden sample)/R band (Silver sample). In Figure 10, we show the SED of the afterglow from a wide-band observations. We note that the difference in the distributions between the Golden and Silver samples are negligible in the SED; thus, it is reasonable to analyze the SED of the Silver sample. However, we find that the SED is still not completely consistent with the theoretical expectations (a single power-law spectral slope). This implies that the traditional extinction models (LMC, SMC, MK, etc.) are not precise in fully describing the extinction character for some GRBs.
We conclude by estimating the influence of redshift on the CI. The redshifts of our bursts lie between z =0.007 and 5.2. In Figure 11, we plot the CIs against the redshifts. We obtain a Spearman correlation coefficient r =0.21, corresponding to a chance probabilityp= ≤ 10 −4 . This indicates that the possible dependence of the CIs on the redshift is very weak. Following the definition of time evolution outlined in Section §2.4.3, we find that there are 20 out of 25 GRBs in the Golden sample that have at least more than two data points that are variable CIs. This indicates that the evolution of the CIs in the afterglow is very common. GRB 080413B is taken as an example, and its CIs as functions of time are displayed in Figure 12(a). 16 The uncertainty of the value stems from the method of obtaining the host extinction. If the extinction value of each energy band is constant, then the derived CIs will also be constant.
The variable CIs are indicated by arrows. We note that the evolution of CIs tends to appear at a late time. In Figure 12(b), we compare the distributions of constant and variable CIs. We see that the nonvarying CIs are distributed similarly from one color to the others. Fitting the histograms by a Gaussian function, one has g-r=0.15±0.14 (mean ± standard deviations), r-i=0.08±0.09, i-z=0.13±0.10, J-H=0. respectively. Instead, the group of varying CIs shows a broad distribution for different colors, possibly indicating different underlying phenomena. We investigate the connection between the values of CIs and the afterglow light-curve components to determine their temporal evolution from one component to another. The components were identified in §2.4. We calculate the mean value of CIs for each emission component 17 . The results are presented in Table 2 for the Golden sample and in Table 3 for the Silver sample.
In Figure 13(a), we show the distribution of CIs among various optical emission components (Data Set II, regardless of the bursts and CIs) and the distribution of ∆CI (results from Gaussian fits are summarized in Table  6) between one component and the following one (Figs. 13b). ∆CI is defined by where t 2 > t 1 . This means that an increasing CI, with a red-to-blue color change, corresponds to ∆CI > 0; in contrast, ∆CI < 0 represents a decrease in CI with a blue-to-red color change.
We defined the global sample as the set of all CIs from both the Golden and Silver samples. We find that various emission components generally have similar distributions, with typical values ∼ 0.2, except the SN component, which presents a different ∆CI distribution. Comparing the CIs, we find that the distributions are consis-17 Note that if several similar emission phases exist at different times in a single burst, we merge them together and calculate the mean value. This is applied for GRB 021004 and GRB 071031 (with two late flares). tent, though the number of data points is small. Note, however, that the distribution is very broad.
We present in Figure 14 the CI-CI correlations between consecutive components, regardless of the CIs. We compute the correlation index r in all cases and find that it is systematically close to 1. The two smallest values are the transition to an SN component with r=0.57 and a flare component with r =0.60. In addition, a linear fit to the slope for most correlations is also close to 1. This indicates that the CIs do not have a significant change between consecutive components. The fit results, along with the Spearman correlation coefficient r and the chance probabilityp are reported in Table  6. This table also gives the number of CIs per emission component and the distribution range of ∆CI.
We find that 263 of 306 (85.9%) pair components satisfy the criterion that the CI does not change. The large fraction of CI without significant time variation for both the Golden and Silver samples suggests that the emission mechanism of the various components in the optical transient may stay the same. This is consistent with the prediction of the external shock afterglow model, for which the emission is dominated by synchrotron radiation, provided that the observation bands do not cover the spectral breaks.
Only 43 of the 306 (14.1%) pair components have a significant difference. For these pairs, we found that 18 include SN components, 5 include revese-shock components, and 5 include flares (see Table 5). The Golden and the Silver samples show similar results. In those pair components that include an SN, 36% show variation. Likewise, 20% of the pair components including flares show variation. This makes the SN and flare components the ones that exhibit most common variations. This is consistent with the observational evidence that the optical transient spectrum typically changes during the associated SN and flares (e.g., Šimon et al. 2001, 2004Zhang et al. 2006;Mu et al. 2016;Geng et al. 2016Geng et al. , 2017, and references therein). We also found a few cases of CI variation for the jet break components. This is note-worthy since a jet-break is believed to be purely dynamical. Therefore, CI evolution is not expected. Here, one should bare in mind the possibility that we could have identified the wrong component. For example, we could have identified a flare as a re-brightening bump if the flare is weak (GRB 021004, GRB 000301C, GRB 060707A). This is because they have similar light curves. Another possibility is that the interpolation procedure for the Silver sample introduces uncertainties affecting the determination of the CI.

Constant Color Indices in the External Shock Model
There is a natural connection between the CIs and the spectral indices. Assuming that the SED of the afterglow in the considered bands is well described by an intrinsically single power law, the spectral flux density (erg cm −2 s −1 Hz −1 ) for the afterglow is described by Equation (2). The relation between the flux and the observed magnitude is given by Equation (2.2), so we can write the CI, which is expressed from the flux density as where C=ZP 2 -ZP 1 , and ZP i = 2.5log 10 (F νi,0 ) is the zeropoint of various bands i. Therefore, the spectral index β o can be connected to the CI as CI = 2.5β o log 10 ν 2 ν 1 According to Equation (9), there exists a linear relation between CIs and spectral indices. Therefore, CIs provide a probe to investigate the spectral properties of the optical afterglow.
Using Equation (9), we derive the index β CI o using constant CIs of Data Set I of the Golden sample (see Appendix Table 7) and the Data Set II of both the Golden and Silver samples during the normal-decay phase. We show the distribution of β CI o in Figure 15. The distribution is well fitted by a Gaussian function, and β CI o =0.68±0.60 for constant CIs of Data Set I and β CI o =0.68±0.68 for the Data Set II during the normal decay phase, which are consistent with each other.
The distributions are consistent with the predictions of the external shock model, with typical β o around 0.75 (Mészáros & Rees 1997;Sari et al. 1998;Zhang et al. 2006). The spectral indices inconsistent with the external shock model are obtained from the high-energy bands, which could point toward a stronger reddening. We note that according to the external shock model, β CI o should cluster between 0.5 and 1.0. However, we found that a small fraction of data deviate from this range, which have large error bars, and thus they cannot accurately constrain the calculation of CI. Then this uncertainty propagates to the calculation of spectral index.
According to the standard synchrotron afterglow model, radiation is produced by electrons distributed with a power-law index p. Therefore, the CIs can be derived from the following equations for both a constant-density interstellar medium (ISM) and a windlike medium (wind): (i) Fast cooling: (ii) Slow cooling: Here ν c is the cooling frequency, ν m is the minimum frequency, and ν 2 and ν 1 are the frequencies of adjacent bands. In the fast cooling regime, electrons at the injection Lorentz factor have time to lose their energy by emitting synchrotron radiation, while they do not in the slow cooling regime. Next, we use the observed CIs (Data Set I of the Golden sample) to derive the electron power-law indices p CI using the Equations (10-11). In Figure 16, we present the distribution of p CI , fit by a Gaussian.
The mean values with their standard deviations are p CI =2.43±1.06 for slow cooling (ν c > ν 2 > ν 1 > ν m ) and p CI =1.43±1.06 for fast/slow cooling (ν 2 > ν 1 > max(ν c , ν m )). We find that the slow cooling scenario is consistent with the theoretical predictions for relativistic shocks, p ∼ 2.5. A wide range of the electron index is obtained, which is consistent with previous studies (Shen et al. 2006;Liang et al. 2008).
Next, we use the closure relations to compute the electron index p from the temporal index α of the normal decay, which was obtained by fitting the light curve, or the spectral index β o . We then further estimate the expected value of the color indices CI th (Data Set III of the Golden sample) using Eq.(10-11). Figure 17 displays the correlation between the CI and CI th for these spectral regimes. Their distributions are also shown with the best Gaussian fits CI obs =0.16±0.14 for observed CI and CI th =0.13±0.07 for theoretical CI.
We find that the data smoothly cluster around the equal line, and both observed CI and theoretical CI th have a similar distributions, consistent with the afterglow synchrotron radiation model with an intrinsic power-law decay.
We also find that a small part of CIs are far away from the equal line, which implies that these data are inconsistent with the model. These CIs are mainly derived from GRB 081029 and 100621A (Fig.17). It is interesting that both GRB 081029 and GRB 100621A present a light curve with re-brightening bumps (see Figure 5). The late re-brightening bumps could come from different emission components. Overlapping bumps could alter the CIs estimation. However, we also find that CIs preceding the bump are still inconsistent with the model. To investigate this further, we searched all the bursts in the Golden sample to find more cases presenting a clear late re-brightening bumps in the light curves. We find that, GRB 071025, GRB 091029, and GRB 100814A are consistent with the theoretical values and present a late re-brightening. The results show that the afterglow emission for both GRB 081029 and GRB 100621A is likely different from other bursts.
In Figure 18, we show color-color diagrams using the Golden sample (Data Set I) to investigate the evolution of CI shift. We calculate the distance of data points to the equal line and the difference value between the two CIs. We find that most of the CIs cluster around the equal line and the typical difference value between two CIs is close to 0. For the distances, one has (gr)-(r-i)=-0.03±0.07, (r-i)-(i-z)=0.00±0.06 and (J-H)-(H-Ks)=0.00±0.14, respectively. For their differences, one has (g-r)-(r-i)=-0.04±0.10, (r-i)-(i-z)=0.01±0.09, and (J-H)-(H-Ks)=0.00±0.20, respectively. This result indicates that the CI shift is not significant, and the intrinsic reddening inside their host galaxies must be quite similar and relatively small for these events. This is also expected from the narrow distributions in wavelength of the ugriz filters. This result also implies that the afterglow spectral shapes are similar from one to another and can be described by a smooth power-law spectral decay, with no bumps or strong lines for these bands. These results are consistent with the finding ofŠimon et al. (2004).

Possible Explanations For The Color Indices Variation
Below we outline three possible explanations for the CI variation, and in §5.2.4 we compare the prediction to the data.

Bands
In the external shock synchrotron model, the cooling frequency, ν c , is likely to cross the optical bands at a time later than approximately 1 hr (Uhm & Zhang 2014). This will cause a change in the spectral indices and thereby in the observed CIs.
According to Eqs. (10)-(11), if the cooling frequency crosses all the optical bands and the spectral regime from a stable phase transition to another stable phase, then the change in the CIs, ∆CI, is given by (i) Fast cooling: (ii) Slow cooling: In the constant environment model (the ISM model), the ν c decreases with time (Sari et al. 1998). The optical band is initially below ν c (bluer spectrum). As ν c decreases in time, the optical band becomes above ν c , and the spectrum becomes redder. This gives rise to a negative ∆CI. On the other hand, in the wind model, ν c increases with time (Chevalier & Li 1999). The optical band is initially above ν c (redder spectrum). As ν c increases in time, the optical band becomes below ν c (the spectrum becomes bluer). This gives a positive ∆CI.
We note that the cooling break in the afterglow spectra may not be sharp (Uhm & Zhang 2014). The transition from one regime to another is rather smooth and requires a very long time, especially in the optical bands. For most of the cases analyzed in this paper the transition has not finished during the period investigated, so what we have tested is an upper limit of ∆CI.

Effects of Host Extinction
The host extinction may also cause the CI change (e.g., Waxman & Draine 2000;Morgan et al. 2014). After being eventually destroyed by the initial intense radiation, dust could reform at an early time near the burst region, changing the host extinction (Perna & Lazzati 2002). This leads to two temporal changes in the CI. During the early afterglow phase (Phase I), dust might be destroyed, and the CI can be described by a red-to-blue change, implying an increase in the CI with ∆CI > 0. On the contrary, at a late time (Phase II), dust could reform and reddening will reappear, and the CI presents a blue-to-red change, which provides a decrease in the CI with ∆CI < 0. We note that Phase II might not occur, since the photodestruction could be irreversible (Perna & Lazzati 2002;Perna et al. 2003;Lazzati & Perna 2003).
The true magnitudes in the optical bands are estimated through considerations of extinctions (in both in our Galaxy and the host galaxy) and the spectral kcorrection. We display the distributions of all correction factors in Figure 19. The factor that affects the achromatic energy band is given by (-A G -K-A H ). This quantity is distributed in [-4.29, 0.22], and, as seen in Fig.19b, it is negligible compared to the brightness (Data Set III).
The factor affecting the CIs can be characterised by the color excess, CE, such that where A H is the magnitude of the host dust extinction and A G is the magnitude of our Galaxy dust extinction.
Here the CE includes contributions from our Galaxy and the host galaxy, while the k-correction factor is the same for both bands. Therefore CE={A G (λ 1 ) − A G (λ 2 ) + A H (λ 1 ) − A H (λ 2 )} is the total correction, ranging from 0.00 to 0.99. Figure 19(c) presents the CE as a function of CI. Since both the CI and CE have similar magnitudes, uncertainties on the corrections have large effects on the value of the CIs.

Effects of the Supernova Emission Component
GRB-associated SNe in some long bursts (Woosley & Bloom 2006;Hjorth et al. 2003;Cano et al. 2017) could produce an additional emission component embedded in the late afterglow synchrotron emission. This results in a CI variation (e.g., Šimon et al. 2004;Gal-Yam et al. 2004), described as a red afterglow bump. From the literature, 13 cases 18 of GRB-SN associations in our sam- 18 GRB 980425/SN 1998bw (Galama et al. 1998), GRB 030329/SN 2003dh (Resmi et al. 2005, GRB 050525A/SN 2005nc ple were reported. We compare the CIs between the GRB-SN emission and their previous components, regardless of the nature of this previous component. We are only concerned here with the variation of the CI (Data Set II). The results displayed by Figure 14 show that a very high fraction (see Table 5) of CIs exhibit a significant variation, indicating a spectral change as expected.

Statistical Analysis
We perform a statistical analysis considering only variable CIs to explore possible physical origins. Figure 20(a) shows the distribution of variable CIs. If the observed variation is due to the crossing of the cooling spectral break, then either positive or negative values of ∆CIs are expected, depending on the density profile of the circumburst medium. As discussed above, a positive value of ∆CIs is expected in the early afterglow phase (earlier than 1 hr) owing to dust destruction. Variable CIs for the SN are defined as the identification of a GRB-SN bump at a later time, which is confirmed in published papers (see §5.2.3). Figure 20(b) presents the temporal evolution of the ratio of varying CIs to the total number of varying CIs in a given time interval. We group the ratios depending, in turn, on (i) the identification of the SN component and whether the ∆CIs are part of it and (ii) the sign of the ∆CIs.
Among the positive ∆CIs (64 % of all cases, labeled II in Figure 20(a)) we find two peaks, at around 500s and ∼ 10 6 s. The early peak is consistent with being due to dust extinction, corresponding to 12% of the ∆CIs. These cases are marked by the red line in the figure. The peak at 10 6 s is consistent with effects of ν c crossing the observed band in a wind environment, causing a red-toblue change, i.e., a shallowing of the spectrum (yellow line in the figure), including 43% of the variable CIs.
Likewise, among the negative ∆CIs (36 % of all cases, labelled I in Figure 20(a)), there are two peaks, one at around 100 s, and one at ∼ 10 5 − 10 6 s. The latter peak contains ∆CIs that are consistent with effects of ν c crossing the observed band in an ISM environment, causing a blue-to-red change, i.e., a steepening of the spectrum (green line in the figure), correspond- ing to 30% of the cases. For the bursts that exhibit a color change at around 100s (blue curve in the figure), the most natural explanation is the transition from reverse-shock (RS) emission to forward-shock (FS) emission. Early on, the emission is likely dominated by the RS component and the optical band is below ν c,r (bluer spectrum). At around 100 s, there is likely a transition from the RS to the FS emission (Type II light curve reported in Zhang et al. 2003). At this time, the FS may be already above ν c,f . This gives a redder color. For the wind model (Wu et al. 2003;Kobayashi & Zhang 2003), a similar situation is expected. This case includes ∼5% of ∆CIs.
Finally, the SN cases dominate the late-time variations (purple curve in the figure), the last one corresponding to 10% of the variable CIs.

CONCLUSION
We studied the multiband optical afterglows of 70 bursts with their CIs. The optical/NIR photometric data are not only based on Swift/UVOT but are also obtained from ground-based instruments, especially the GROND. This catalog of multiwavelength GRB afterglow data provides an opportunity for a statistical study. We corrected the data with the Galactic and host galaxy extinctions. After performing a spectral k -correction, we divided the bursts into two samples depending on whether the different optical bands have a simultaneous observation. Our main results are summarized as follows: • A Golden sample includes 25 out of 70 GRBs with five CIs. The distributions of these CIs are approximately the same. There is no color shift between CIs, and in general they satisfy the afterglow model prediction of a single power-law spectral slope.
• A Silver sample includes 45 out of 70 GRBs with 10 CIs with wider energy ranges. The distributions of these CIs are consistent with the Golden sample in the corresponding energy bands. For most of these CIs, the typical value is around 0.2, and they also, in general, satisfy the afterglow model prediction of a single power-law spectral slope. There are two significantly inconsistent peaks, mainly distributed in the UV bands, with UVM2-UVW1=-0.19 and UVW1-U=0.93. This is consistent with the theoretical prediction that the intrinsic reddening is more significant in the highenergy bands.
• After performing all the corrections for the data, the SED of the optical transient, which is derived from the average CI, presents a deviated single power-law spectral slope. This implies that the traditional extinction models could be not accurate enough to describe the extinction characteristics for some GRBs.
• Most CIs (95.9%) are constant in the Golden sample (Data Set I). They are tightly distributed with standard deviations of ∼ 0.2. The variable CIs (4.1%) are widely distributed, presenting evidence for several different emission mechanisms. We determined that around 30% are due to the crossing of the cooling spectral break in an ISM medium, while 43% are due to a wind-like medium, 10% due to the SN emission, 12% due to early dust extinction, and 5% due to the transition from RS to FS emissions.
• Component-wise, we found that most cases of the varying CIs (Data Set II) are in the transition during the SNe, the Reversed-Shock, or the Flare components.
• We derived the correlations between CIs and spectral indices based on the standard afterglow synchrotron model and obtained β CI o (using the constant CI of data Set I and Data Set II during the normal-decay phase) and p CI (using the constant CI of Data Set I). The typical values are β CI o =0.68±0.60 (Data Set I) and β CI o =0.68±0.68 (Data Set II). We found that these two distributions are consistent with each other. Moreover, p CI =2.43±1.06 for Slow cooling (ν c > ν 2 > ν 1 > ν m ) and p CI =1.43±1.06 for Fast/Slow cooling (ν 2 > ν 1 > max(ν c , ν m )), which are consistent with the theoretical predictions for relativistic shocks. A wide range of p values is obtained, which is consistent with previous findings. In turn, we derived the theoretical CIs using the p values derive from the closure relation and compared them to the observed CIs. We found that they still clustered around the equal line.
• We also investigated the overall behavior of the CIs: (i) The CIs are independent of the redshift.
(ii) More that 90% of CIs are distributed between -1.0 and 1.0. They have negligible variations compared to the decline of the brightness with time.
(iii) The dust extinction correction and the spectral k-correction can significantly affect the values of CIs and spectral indices.
We compiled a sample of multiband observations of optical transients to investigate the temporal evolution of the CIs. The variable CIs are clues to unveiling the origin and details of the optical transient and of GRBs in general, prompting for more prolonged and simultaneous observations.  References-References for z: GRB 970508: Bloom et al. (1998) Table 2 continued  Table 2 continued    Table 3 continued  Table 3 continued  Table 3 continued   b Statistics regardless GRBs and colors.

GRB 050801
Magnitudes Color indices Time (s) Colors

GRB 060218
Magnitudes Color indices Time (s) Colors

GRB 060614
Magnitudes Color indices Time (s) Colors

GRB 061126
Magnitudes Color indices Time (s) Colors

Color indices
Time ( Colors

GRB 080330
Magnitudes   Figure 10. SED of the afterglow, which derives from average values of CIs (Data Set III) and assumes a typical observation magnitudes for r band (Golden sample)/R band (Silver sample). Different colors represents different energy bands. The stars represent the Golden sample (SDSS griz system), while the circles represent the Silver sample (common UBVRI photometry system). Colors Indices     Figure 16. Distributions of observed electron spectral indices p CI , derived from the constant CI of the Golden sample (Data Set I). The dashed lines are the best Gaussian fits giving p CI =2.43±1.06 for the slow cooling case (νc > ν2 > ν1 > νm) and p CI =1.43±1.06 for the fast/slow cooling case (ν2 > ν1 > max(νc, νm)). . Correlation with their distributions between theoretical and observed CI for all GRBs in the three different spectral regimes (Data Set III). The red dashed line is the equal line, and the red solid lines are the best Gaussian fits giving CI obs =0.16±0.14 for observed CI and CI th =0.13±0.07 and for theoretical CI.       Table 7 continued  Table 7 continued     Table 7 continued  Table 7 continued    Table 7 continued      Table 7 continued   Table 7 continued  Table 7 continued     Table 7 continued  Table 7 continued       Table 7 continued