Photometric and Spectroscopic Observations of GRB 190106A: Emission from Reverse and Forward Shocks with Late-time Energy Injection

Early optical observations of gamma-ray bursts can significantly contribute to the study of the central engine and physical processes therein. However, of the thousands observed so far, only a few have data at optical wavelengths in the first minutes after the onset of the prompt emission. Here we report on GRB 190106A, whose afterglow was observed in optical bands just 36 s after the Swift/BAT trigger, i.e., during the prompt emission phase. The early optical afterglow exhibits a bimodal structure followed by a normal decay, with a faster decay after ∼T 0 + 1 day. We present optical photometric and spectroscopic observations of GRB 190106A. We derive the redshift via metal absorption lines from Xinglong 2.16 m/BFOSC spectroscopic observations. From the BFOSC spectrum, we measure z = 1.861 ± 0.002. The double-peak optical light curve is a significant feature predicted by the reverse-forward external-shock model. The shallow decay followed by a normal decay in both the X-ray and optical light curves is well explained with the standard forward-shock model with late-time energy injection. Therefore, GRB 190106A offers a case study for GRB emission from both reverse and forward shocks.


INTRODUCTION
Since their discovery in the 1960s, the understanding of gamma-ray bursts (GRBs) has been greatly advanced by various space satellites and ground-based telescopes (see Zhang 2018, for a review).Based on statistics of their prompt emission, GRBs are divided into two categories defined by their T 90 duration and spectral hardness, i.e., short bursts with hard spectra (also called type I bursts) with T 90 < 2 s and long bursts with soft spectra (also called type II bursts) with T 90 > 2 s (Kouveliotou et al. 1993;Zhang et al. 2009).The association of the short GRB 170817A with the gravitational wave event GW170817 confirmed the longstanding prediction (Eichler et al. 1989) that neutron star mergers can indeed produce short bursts (Abbott et al. 2017a,b;Zhang et al. 2018).Long GRBs are generally believed to originate from massive star collapse (Woosley & Bloom 2006).GRB 980425/SN 1998bw demonstrated that some stellar collapses produce long bursts and broad-lined Type Ic supernovae (Galama et al. 1999).The comparison between longduration GRBs and broad-lined SN rates suggests the ratio GRBs/CC-SNe is a few ∼ 10 −3 (Guetta & Della Valle 2007).
Recent observations suggest that collapsing stars can also produce GRBs that last shorter than 2 s (Ahumada et al. 2021;Zhang et al. 2021;Rossi et al. 2022), compact star mergers seem to be also able to produce long GRBs (Rastinejad et al. 2022), indicating that classification of GRBs is more complex than the simple temporal division.Since its launch in 2004, the Neil Gehrels Swift Observatory (Swift hereafter) has detected more than 1700 GRBs with accurate position, enabling easy follow-up observations with groundbased telescopes (Gehrels et al. 2004).The typical isotropic equivalent energy (since it is the energy equivalent to the one that the GRB should have if the emission were isotropic) of GRBs ranges between 10 50 − 10 54 ergs within a few seconds to minutes (Atteia et al. 2017).The jet break predicted by the collimated jet model has been observed in some GRBs, leading to the estimates of the true energy being 2 − 3 orders of magnitude lower than the isotropic energy release (Frail et al. 2001).
Previous studies have shown that the X-ray light curves of GRB afterglows typically have several stages of power-law decay (F ∝ t −α ): a steep decay phase (typically α ∼ 3), a shallow decay phase also called a plateau (typically α ∼ 0.5), a normal decay phase (typically α ∼ 1.2), a post-jet-break phase (typically α 2), with some bursts also accompanied by flares (Zhang et al. 2006).The typical optical light curve is an early rise followed by a normal decay, then sometimes followed by a re-brightening or a post-jet-break decay.Some long bursts, when at low enough redshift to allow a search, are found to be accompanied by a supernova bump (Zeh et al. 2004;Li et al. 2012).
In the standard model of GRBs, the afterglow is believed to involve a relativistically expanding fireball.Considering a relativistic thin shell with energy E K , initial Lorentz factor Γ 0 , opening angle θ j , and a width in laboratory frame ∆ 0 expanding into the interstellar medium (ISM) with density n.The interaction between the shell and ISM is described by two shocks: 1) a forward shock (FS) propagating into the ISM.The shell first undergoes a coasting phase, where it moves at a nearly constant speed Γ(t) Γ 0 .It starts to decelerate and evolves into the second stage when the mass of the ISM swept by the forward shock is about 1/Γ 0 of the rest mass of the shell.The third phase is the post-jet-break phase when the 1/Γ cone becomes larger than the geometric cone defined by θ j .Finally, the blastwave enters the Newtonian phase (its velocity is much smaller than the speed of light) when it has swept up the ISM with the total rest mass energy comparable to the energy of the shell.During all these three phases, electrons are believed to be accelerated at the forward shock front to a power-law distribution N (γ e ) ∝ γ −p f e .A fraction e,f of the shock energy is distributed into electrons, and a fraction B,f is in the magnetic field generated behind the shock.Accounting for the radiative cooling and the continuous injection of new accelerated electrons at the shock front, one expects a broken power-law energy spectrum of them, which leads to a multi-segment broken power-law radiation spectrum at any epoch (Gao et al. 2013, for a review).The typical afterglow observations are successfully described by the synchrotron emission of electrons from the forward shocks.However, the information of the properties of the fireball is lost in the forward shock, especially when it has been decelerated by the ISM; 2) a reverse shock (RS) propagating into the shell, which occurs in the very early afterglow epoch.In the thin shell case, the reverse shock is too weak to decelerate the shell effectively.The Lorentz factor of the shocked shell material is almost constant while the shock propagates in the shell.The reverse shock crosses the shell at the fireball deceleration time.At this time, the reverse and forward shocked regions have the same Lorentz factor and internal energy density.However, the typical temperature of RS is lower since the mass density of the shell is higher.Consequently, the typical frequency of RS is lower, and RS may play a noticeable role in the low-frequency bands, e.g., optical or radio (Mészáros & Rees 1997).After the RS crosses the shell, the Lorentz factor of the shocked shell may be assumed to have a general power-law decay behavior γ 3 ∝ r −g .The shell expands adiabatically in the shell's comoving frame.The emission from the RS is sensitive to the properties of the fireball.The observations of RS can thus provide important clues on the nature of the GRB source.
Catching very early (a few minutes after the trigger) optical emission from GRBs is not an easy task.There are only a few GRBs with optical observations during the prompt phase (Gao & Mészáros 2015, for a review).Oganesyan et al. (2019) provided a sample of 21 GRBs with at least one observation during the prompt emission phase, and several GRBs in the sample are well observed (eg., GRB 070616, GRB 081008, GRB 121217A).Recently, Oganesyan et al. (2021) reported that they began to observe GRB 210619B ∼ 28s after the Swift trigger with D50 clear band during the prompt emission phase.However, our understanding of the earlier phases of GRBs continues to be quite incomplete compared to the later afterglow phases (occurring tens of minutes after the trigger).
In this paper, we report a special burst, GRB 190106A, for which optical observations began during the prompt emission phase, showing a double-peaked early optical light curve.We present optical photometric and spectroscopic observations of GRB 190106A with the NEXT (Ningbo Bureau Of Education And Xinjiang Observatory Telescope) and the Xinglong 2.16-m Telescopes.With the Xinglong 2.16-m/BFOSC spectroscopic observation, we present the redshift and calculate the equivalent widths of absorption lines in the afterglow spectrum.Combining our data with the Swift/XRT data and other observations from GRB Coordinates Network (GCN) reports, we constrain the parameters in an external shock model of the afterglow.
Our observations and data reduction are described in Section 2. The spectroscopic data analysis and redshift measurements are presented in Section 3. The optical and X-ray afterglow data analysis and the modeling of the afterglow light curve are reported in Section 4. We discuss possible models and present the conclusions in Section 5 and Section 6, respectively.A standard cosmology model is adopted with (Planck Collaboration et al. 2014).

OBSERVATIONS
GRB 190106A triggered Neil Gehrels Swift Observatory (short as Swift) Burst Alert Telescope (BAT, Barthelmy et al. 2005) at 13:34:44 UT on 2019 Jan 06 (Sonbas et al. 2019).It is a long-duration burst with T 90 = 76.8±2.4 s.Konus-Wind also detected the burst with T 90 = 77.1 +1.3  −4.2 s (Tsvetkova et al. 2019).The MASTER-Amur robotic telescope quickly slewed to the burst position, starting observations 36 s after the BAT trigger.An optical transient (with a magnitude of 15.23 in the P-band) was discovered within the Swift/BAT error circle in the second 10 s exposure image (Yurkov et al. 2019;Lipunov et al. 2019) located at equatorial coordinates (J2000.0):RA = 1 h 59 m 31.19 s , Dec. = +23 • 50 44.79, consistent with the Swift/UVOT position RA = 01 h 59 m 31.18s , Dec. = +23 • 50 44.0(Kuin & Sonbas 2019).XRT and UVOT started observations at 81.8 s and 90 s after the BAT trigger, respectively.The BAT data files were downloaded from the Swift data archive1 .We extracted the BAT light curve and spectrum by the HEASoft software package (v6.27.2).The 15-150 keV spectra are analyzed in XSPEC (version 12.12) with power-law approximation over the time interval from 71 s to 81 s after the BAT trigger.We adopted the analysis results of the XRT repository products (Evans et al. 2007(Evans et al. , 2009) ) and downloaded the 0.3-10 keV unabsorbed light curve from the UK Swift Science Data Centre2 .
In our subsequent ground-based optical follow-up observations, the Xinglong 2.16-m obtained spectroscopic observations from which the redshift was determined (Xu et al. 2019), also confirmed by GMG (Mao et al. 2019) and Xshooter (Schady et al. 2019).We also collected observations from the GCNs reported by MASTER (Lipunov et al. 2019) and Sayan Observatory, Mondy (Belkin et al. 2019a,b).The MASTER data are converted to R-band by adding ∆R = 0.51 mag following Xin et al. (2018).The corrected MAS-TER magnitude and Mondy optical observations are shown in Fig. 2 in gray color.

NEXT optical observations
NEXT is a fully automated telescope with a diameter of 60 centimeters, located at Nanshan, Xingming Observatory.NEXT is equipped with a E2V 230-42 back-illuminated CCD, and the CCD controller was made by the Finger Lakes Instrumentation (FLI).The size of the CCD is 2048 × 2048 pixels, with a pixel size of 15µm.The pixel scale is 0.64 /pixel, and the field of view (FOV) is 22 × 22 according to the size of the CCD.The typical gain is 1.85e − /ADU, and typical readout noise is 13 e − with 500 kHz readout speed.The equipped filters were in the standard Johnson-Cousins system.
NEXT is connected to GCN/TAN alert system, but did not observe GRB 190106A immediately as the dome was still closed at the very beginning of the event.NEXT began to obtain the first image at 13:57:21 UT, 22.6 minute after the BAT trigger, and the final one ended at 20:09:26 UT when the airmass was ∼ 9.The NEXT image of the burst location is also shown in Fig. 1.We obtained 88 images, all in the R C filter.Five of them were discarded for poor quality.The observing center time and exposure time of each image used are presented in Table 4.
The data reduction was carried out by standard processes in the IRAF packages (Tody 1986), including bias, flat and dark current corrections.The cosmic rays were also removed by the filtering described in van Dokkum (2001).The final five useful images were combined using the imcombine task in IRAF.The astrometry was calibrated using Astrometry.net(Dey et al. 2019).A weak but obvious source with r ∼ 24.0 mag is located at the GRB position, which could be the host galaxy of the burst.Middle panel: the first NEXT image, obtained 0.39 hours after the BAT trigger, shows a bright source consistent with the UVOT coordinates with R = 16.9 mag.Right panel: the NEXT image obtained at 3.46 hours after the burst.The source in the image center is obviously fainter than the one in the middle panel, with a brightness of R = 18.4 mag.(Evans et al. 2009).The gray circles and gray squares are modified MASTER and Mondy photometric results, respectively.The red diamond points are our observations obtained by NEXT and Xinglong 2.16-m.Most of the error bars are smaller than the marker, and our photometric data are listed in the Table.4. (Lang et al. 2010).The magnitudes were measured using SourceExtractor (Bertin & Arnouts 1996) using a circular aperture with a nine-pixel diameter.The photometric calibration was derived using the Sloan Digital Sky Survey 14 th data release (Abolfathi et al. 2018), with flux/mag conversion of the SDSS system into the Johnson-Cousins system3 .The light curve is shown in Fig. 2

Xinglong 2.16-m Observations
The spectroscopic observations were carried out with the National Astronomical Observatories, Chinese Academy of Sciences (NAOC) 2.16 m telescope (Fan et al. 2016) in Xinglong Observatory on 06 Jan 2019 at 14:02:39, 27.9 minutes after the BAT trigger.The optical spectrum was obtained with the Beijing Faint Object Spectrograph and Camera (BFOSC) equipped with a back-illuminated E2V55-30 AIMO CCD.The optical afterglow magnitude was ∼ 16 mag in R band and the airmass was approximately 1.3 at the beginning of the spectroscopic observations.Due to a seeing of about 2 a slit-width of 1. 8 was used.The slit was oriented in a south-north direction.The grating used was G4 and the order-sorter filter 385LP was also used, leading to a spectral coverage of 3800 to 9000 Å at a resolution of about 10 Å.A total of two spectra were obtained with exposure time of 3600 s and 2400 s, respectively.
The BFOSC spectroscopic data were reduced by the IRAF standard processes to obtain clean images, including bias subtraction, flat-field correction and cosmic-ray removal.With the using of "NOAO" package provided in IRAF, we extracted the one-dimensional spectrum of the burst, and calibrated the wavelength with arc-spectra from an iron-argon lamp.The flux-calibration was performed by the standard star HD+19445 (Oke & Gunn 1983) obtained with the same instrumental setup on the same night.After the processing above, the extracted spectrum is shown in Fig. 3.
The photometric observations were carried out on 7, 8 and 10 of Jan in the R filter again using the BFOSC (Fan et al. 2016).The exposure times were 5 × 360 s, 7 × 360 s, and 10 × 360 s, respectively.The data reduction and processing was the same as NEXT, except dark current, which was negligible.The subsequent magnitude measurements were also performed as described above.The photometric data are pre-sented in Table 4 and shown in Fig. 2. A comparison with a sample of historical GRB optical light curves is shown in Fig. 4.

SPECTROSCOPIC ANALYSIS
As the blue end of the BFOSC spectrum has poor signalto-noise ratio, our spectral analysis starts at 4300 Å.After excluding the A and B band absorption and sky emission lines, we identified the metal absorption lines detected at more than 3σ significance.We detected lines from C IV, Al II, Al III, Fe II, and Mg II.Line centroids were determined by fitting a Gaussian function to the line profiles.We obtain a redshift of z = 1.861 ± 0.002.The identified lines with the restframe equivalent widths (EW obs ) are listed in Table 5 and also shown in Fig. 3.

Line strength analysis
We analyzed the strong features of the GRB host galaxy interstellar medium detected by BFOSC, following de Ugarte Postigo et al. (2012).The Line Strength Parameter (LSP) of the spectrum is LSP = 0.15 ± 0.18, which is very similar to the average value of the sample presented in the paper above.The line strength of the burst compared with the sample is shown in Fig. 5.

Temporal Analysis
The R-band, γ-ray and X-ray at 10 keV light curves are shown in Fig. 2. The optical light curve started at T 0 + 36 s and ended at T 0 + 4.02 d after the burst.According to the decay of the light curve, the optical and X-ray light curves can be divided into five and four stages, respectively.We fit the temporal indices α of the light curves in Fig. 2  Corrected Rc Magnitude in the z=1 system The optical afterglow of the R-band light curve can be divided into five stages with different decay indices, which are listed in Table 1, where the start and end times of different stages are also listed.The first stage is called fast rising, from the beginning of observation to ∼ 67 s after the burst.This is followed by a fast decay until ∼ 1.3 × 10 2 s.In the third stage (slow rising) it is still rising, the index has decreased by about 1 compared with the first stage.The fourth stage, normal decay, it starts from ∼ 2.8 × 10 2 s and continues to ∼ 7.7 × 104 s with a decay index of 0.63.In the final stage, late decay, the decay index increases to 1.29.
The light curve in X-rays seems simpler than in optical, and can be divided into the stages directly.The first stage is the tail of the prompt emission (steep decay) with a decay index of ∼ 4.2, observations started at T 0 + 72.8s.It is followed by a shallow decay from ∼ 2.6 × 10 1 s to ∼ 1.6 × 10 4 s, with an index of 0.36.This stage can also be called a plateau phase, resulting from the small index.In the final stage, which called late decay, the decay index has increased to ∼ 1.3.
The BAT (15 − 350 keV) count rates are also shown in the Fig. 2. The light curve can be simply divided into six episodes.The start, end and peak times of each episode are shown in Table 2.We note that between ∼ 71 s and ∼ 81 s of the prompt emission, the optical light curve has a similar evolutionary behavior: they all peak around 70 s and then fade.However, the peak time is slightly different, i.e., ∼ 76 s for γ − ray, ∼ 76 s for X-ray, and ∼ 68 s for optical.The two high energy bands are almost simultaneous, but the optics peak about 8 s earlier.In addition, the BAT and XRT spectrum of this period can be described using a single power-law F ν ∝ ν −β , and β BAT = 0.77 ± 0.04, β XRT = 0.72 ± 0.10, which is shown in Fig. 6.Although the UVOT optical data are not included in the spectral fitting, they are shown in the best-fit spectra in Fig. 6.As we see, the optical data are much lower than the extrapolation of the best spectral fits for high energies.The SED containing optical and high energy (Xray, γ-ray) can't be fitted by a single power-law.

Afterglow SED analysis
To study the spectral energy properties of the afterglow, we select two epochs to construct its Spectral Energy Distribution (SED) with the best multi-band data: 28.8 ks (epoch 1) and 147.6 ks (epoch 2).The X-ray time sliced spectrum are obtained from the online repository 4 , and the optical data are from Izzo et al. (2019) and Dichiara et al. (2019).We utilized the "grpph" tool of the Xspec package to re-bin the X-ray data to guarantee that each bin contains at least 20 counts to improve the SNR.We fit the spectrum by the power law function from "zdust*zphabs*phabs*powerlaw" model in the Xspec package, where "zdust" represents extinction by dust grains from the host galaxy of the burst,  .The prompt phase SED combined with Optical (black), XRT (red) and BAT (green).The BAT time interval is from 71 s to 81 s, while the XRT is from 85 s to 95 s.The satellite was still slewing before 85 s, so the XRT data before 85 s are ignored.The optical point is converted from the MASTER observation and not included in the fitting.The data points of BAT and XRT were fitted with a single power-law function, and the spectral indices are 0.77± 0.04 and 0.72 ± 0.10, respectively.
"zphabs" and "phabs" are neutral hydrogen photoelectric absorption of the host galaxy and Galaxy, respectively.The redshift and the Galactic hydrogen column is fixed to 1.86 and 7.1 × 10 20 cm −2 (HI4PI Collaboration et al. 2016), respectively.Small Magellanic Clouds extinction law is used as the host galaxy extinction model.The best fitted extinction E(B-V) of the host galaxy is 0.1 for epoch 1, which is also fixed in the fit of epoch 2. The best fitting result gives the photon index Γ = −1.78± 0.03 and Γ = −1.79± 0.01 for epoch 1 and 2, respectively.The results of these two epochs are shown in Fig. 7.

Theoretical Interpretation
The afterglow light curve of GRB 190106A is unusual.There is an optical flash followed by a second bump in the early optical light curve, and then a shallow decay followed by a normal decay.The X-ray light curve shows a steeper decay at the end of the prompt emission, followed by a shallow decay, and finally a normal decay.
To understand the multi-band observational data of GRB 190106A, we consider a relativistic thin shell with energy E K,iso , initial Lorentz factor Γ 0 , opening angle θ j , and a initial width in laboratory frame ∆ 0 expanding into the ISM with density n.The interaction between the shell and ISM is described by two shocks: an RS propagating into the shell and an FS propagating into the ISM (Rees & Meszaros 1992;Mészáros & Rees 1997;Sari et al. 1998;Sari & Piran 1999;Zou et al. 2005).There are four regions separated by the two shocks: Region 1, the ISM with density n 1 ; Region 2, the shocked ISM; Region 3, the shocked shell material; Region 4, the unshocked shell material with density n 4 .First, the pair of shocks (FS and RS) propagating into the ISM and the shell, respectively.After the RS crosses the shell, the blastwave enters the deceleration phase.We now briefly describe these shocks separately.
The FS model we adopt was described in Lei et al. (2016).The dynamical evolution of the shell are calculated numerically using a set of hydrodynamical equations (Huang et al. 2000) where R and t are the radius and time of the event in the burster frame, m is the swept-up mass, M ej = E K,iso (1 − cos θ j )/2(Γ 0 −1)c 2 is the ejecta mass, m p is the proton mass, and β = √ Γ 2 − 1/Γ.The last equation might be replaced with (Geng et al. 2013) when there is energy injection from GRB central engine.For t start < t < t end , the injected power is L inj = L 0 inj (t/t start ) −q , where L 0 inj is the initial injection power, q is the decay power law index, t start and t end are the start and end time for energy injection.By solving these equations with the initial conditions, one can find the evolution of Γ(t) and R(t).
During the dynamical evolution of FS, electrons are believed to be accelerated at the shock front to a power-law distribution N (γ e ) ∝ γ −p f e .Assuming a fraction e,f of the shock energy e 2 = 4Γ 2 n 1 m p c 2 is distributed into electrons, this defines the minimum injected electron Lorentz factor, where m e is electron mass.We also assume that a fraction B,f of the shock energy is in the magnetic field generated behind the shock.This gives the comoving magnetic field The synchrotron power and characteristic frequency from electron with Lorentz factor γ e are ν(γ e ) Γγ 2 e q e B 2πm e c , where σ T is the Thomson cross-section, q e is electron charge.The peak spectra power occurs at ν(γ e ) By equating the lifetime of electron to the time t, one can define a critical electron Lorentz factor γ c the electron distribution shape should be modified for γ e > γ c when cooling due to synchrotron radiation becomes significant.Accounting for the radiative cooling and the continuous injection of new accelerated electrons at the shock front, one expects a broken power-law energy spectrum of them, which leads to a multi-segment broken power-law radiation spectrum separated by three characteristic frequencies at any epoch (Gao et al. 2013, see Equations 15-17 therein).The first two characteristic frequencies ν e and ν c in the synchrotron spectrum are defined by the two electron Lorentz factors γ e and γ c .The third characteristic frequency is the self-absorption frequency ν a , below which the synchrotron photons are self-absorbed (self-absorption optical depth larger than unity).The maximum flux density is F ν,max,f = N e,2 P ν,max /4πD 2 , where N e,2 = 4πR 3 n 1 /3 is the total number of electrons in Region 2 (shocked ISM) and D is the distance of the source.
Besides the numerical calculation of FS as described above, we can also provide an analytical description of the main properties of the evolution and emission of FS.Generally, the evolution includes four phases.The first is a coasting phase, in which we have Γ(t) Γ 0 .In the second phase, the shell starts to decelerate at the deceleration time when the mass m of the ISM swept by FS is about 1/Γ 0 of the rest mass in the ejecta M ej .After t dec , the shell approaches the Blandford & McKee (1976) self-similar evolution Γ(t) (17E K,iso /1024πn 1 m p c 5 t 3 ) 1/8 and R(t) (17E K,iso t/4πn 1 m p c) 1/4 .Later, as the ejecta is decelerated to the post-jet-break phase at the time when the 1/Γ cone becomes larger than θ j .Finally, the blastwave enters the Newtonian phase when it has swept up the ISM with the total rest mass energy comparable to the energy of the ejecta.The dynamics is described by the well-known Sedov-Taylor solution.The synchrotron flux can be described by a series of power-law segments F ν ∝ t −α ν −β (Sari et al. 1998;Gao et al. 2013).For example, in the second phase (deceleration phase), Γ ∝ t −3/8 and R ∝ t 1/4 , one has the scalings for FS spectra parameters, i.e., ν m ∝ t −3/2 , ν c ∝ t −1/2 and F ν,max,f ∝ t 0 (Gao et al. 2013, see Equation 49 therein).In the regime ν a < ν m < ν < ν c , one has (Gao et al. 2013, see Table 13 therein).
The RS model are described in Kobayashi (2000).In the thin shell case, the RS is Newtonian γ34 ∼ 1.The Lorentz factor of the shocked shell γ 3 is nearly constant during the shock crossing the shell.The shell begins to spread around R s = Γ 2 0 ∆ 0 .The scalings before RS crossing time t dec are (Kobayashi 2000).
where N 0 = M ej /m p is the total number of electrons in the ejecta.
We also assume that electrons are accelerated at the RS front to a power-law distribution N (γ e ) ∝ γ −pr e , a fraction e,r of the shock energy e 3 is distributed into electrons and a fraction B,f to the magnetic field in Region 3. Follow the similar procedure as FS, we can obtain the scalings for parameters in the RS synchrotron spectrum (Gao et al. 2013, see Equation 26 therein), After the RS crosses the shell, i.e.,t > t dec , the Lorentz factor of the shocked shell may be assumed to have a general power-law decay behavior γ 3 ∝ r −g .The shell expands adiabatically in the shell's comoving frame.The dynamical behavior in Region 3 are expressed with the scalings (with g ∼ 2), In the same way, we can obtain the scalings for parameters in the RS synchrotron spectrum for t > t dec (Gao et al. 2013, see Equation 31 therein), where ν cut is the cut-off frequency.After RS crossing, no new electrons are accelerated, nu c will be replaced with . In this section, we will show that the double-peak optical light curve as well the X-ray observations of GRB 190106A can be well-explained by synchrotron emission from reverse and forward shocks with late-time energy injection.We adopt a numerical code described above for both FS and RS to model the multi-band observational data.The modeled X-ray and optical light curves corresponding to the above-mentioned parameters (i.e., kinetic energy E K,iso , initial Lorentz factor Γ 0 and opening angle of the relativistic shell, microphysics parameters for FS ( e,f , B,f ,p f ) and RS ( e,r , B,r ,p r ), and energy injection parameters (L inj ,q,t start , t end ) ) are displayed in Fig. 8.However, these model parameters still suffer degeneracy when fitting the optical and X-ray data.In this work, we do not attempt to fit the data across a large parameter space.We present a set of parameter values that interpret the data well, as shown in Table 3.The subscripts "f " and "r" denote the FS and RS emission, respectively.

4.3.1.
Early Double-peak Optical Light curve Zhang et al. (2003) pointed out that, depending on parameters, there are two types of early optical light curve for a fireball interacting with a constant-density ISM, that is, "rebrightening", in which a distinct reverse-shock peak (with α rise = −3p r + 3/2 and α decay = (27p r + 7)/35, see the description of RS model above and also Tables 4 and 5 in Gao et al. (2013)) and forward-shock peak (ν m,f crossing peak) are detectable, and "flattening", in which the FS peak is buried beneath the RS peak.The early optical light curve during the time span ∼ 40 − 1000 s post-trigger shows a double-peak feature, which is clearly a "rebrightening" case.
The first bump (∼ 40 − 133 s) with rising temporal index α O = −1.88±0.57and decay index α O = 1.35±0.32can be interpreted as the combinations of a reverse-shock peak (see the blue dashed line in Fig. 8) and forward-shock emission (see the blue dotted line in Fig. 8).We find the results with p r = 2.1 can roughly explain such bump.
The first peak is also the RS crossing time (or deceleration time) t dec .We can estimate the initial Lorentz factor Γ 0 ∼ 300, if we put E K,iso = 9 × 10 52 erg and n 1 = 0.1cm −1 into Equation ( 11).
The second peak (∼ 133.8 s) of the optical light curve occurs when the typical synchrotron frequency ν m,f crosses the observed frequency (Sari et al. 1998).The slow-rising index α O = −0.58± 0.17 of the second optical bump (∼ 133 − 1000 s) is consistent with emission below ν m in a slow-cooling scenario, i.e., F ν,f ∝ t 1/2 ν 1/3 (Gao et al. 2013, see Table 13 therein).For such FS peak (ν m,f crossing peak), one has (Gao et al. 2013, see Equation 49 therein).The values of jet isotropic kinetic energy E K,iso , ISM density n 1 and the microphysics parameter B,f are chosen to fit this peak.

Shallow Decay Phase
GRB 190106A shows quite a long shallow decay in both X-rays (following the sharp decay) and the optical (following the second bump).In a slow-cooling regime of the forward shock, we expect F ν ∝ ν −(p f −1)/2 t −3(p f −1)/4 for ν m < ν < ν c , which corresponds to α = 3(p f − 1)/4 0.8 using p f = 2.07 (Gao et al. 2013, see Table 13 therein) .Our observed values α O = 0.63±0.01 and α X = 0.36±0.03are shallower than those expected.This phase can be understood within the standard external-shock model with nearly constant energy injection L inj ∝ t −q (see Fig. 8, and Xin et al. 2012 for a similar case).In Table 3, we present the parameter values that interpret the data well, i.e., the index q 0.05, initial injection power L 0 inj 4.5 × 10 49 erg s −1 , and starting time t start > 300/(1 + z) s.The energy injection shuts off at t end ∼ 1500/(1 + z) s.
There are two popular models for the energy injection, i.e., the spin-down of a magnetar and the fall-back accretion onto a stellar mass black hole (BH).First, we consider a fastspinning magnetar as the central engine.The characteristic spin-down luminosity L 0 and the characteristic spin-down timescale τ are related to the magnetar initial parameters as (Zhang & Mészáros 2001): For the NS equation of state (EoS), we adopt EoS GM1 (the radius of the magnetar R = 12.05 km and the rotational inertia I = 3.33 × 10 45 g cm −2 ) as suggested by the recent studies with GRB data (Lü et al. 2015;Gao et al. 2016).Inserting L 0 = L inj = 4.5 × 10 49 erg s −1 and τ 1500/(1 + z) into the above two equations, one can thus infer a magnetar initial period of P 0 ∼ 1.7 ms and a magnetic field B p ∼ 3.5 × 10 15 G.
As we can see, both central engine models can give rise to the energy injection required.However, a fall-back accreting BH tends to produce a giant bump in GRB afterglow light curves (Wu et al. 2013).The plateau phase is a natural expectation from a magnetar central engine (Zhang & Mészáros 2001;Lü et al. 2015).The small injection index q 0.05 (as shown in Table 3) favors the magnetar model.
Assuming that this break is due to the jet break, the predicted temporal index will become α O = 1.5 and α X = 1.7 by adopting p f = 2.07.As shown in Fig. 8, our model can roughly fit the late time optical and X-ray light curves.Using this jet break time, we can estimate the opening angle θ j ∼ 3.3 • if we inset E K,iso = 9 × 10 52 erg, n 1 = 0.1 cm −1 and t 7 × 10 4 /(1 + z) s into the analytical expression Equation ( 12).As shown in Table 3, our numerical modeling gives the opening angle θ j ∼ 4.2 • , which is slightly larger than this analytical estimation.

DISCUSSION
With the redshift measurement in this work, we tested the GRB by using two well-known relations, namely the Amati (Amati et al. 2002;Amati 2006) and Yonetoku (Yonetoku Optical Flux Density (Jy) Figure 8.The modeling of X-ray (thick lines) and optical R band (thin lines) light curves.Observational data are presented with circles for optical (R band) and pluses for X-ray.The blue dashed lines represent emission from FS, while green dot-dashed lines for RS.The combined emission from FS and RS are shown with red solid lines.The parameters adopted are listed in Table 3.The early double-peak optical light curve was expected by Zhang et al. (2003).We show that GRB 190106A is a good case with reverse-forward shocks, just like GRB 041219A (Fan et al. 2005) and GRB 110205A (Zheng et al. 2012).Such double-peak ("rebrightening") light curve is typical for RS+FS emission combinations in a homogeneous ISM case (Zhang et al. 2003;Kobayashi & Zhang 2003).The rising temporal index α O = −1.88< −1/2 favors the RS in thinshell regime.Therefore, we adopted a thin-shell reverseforward shock model in a constant ISM case to investigate GRB 190106A.In this scenario, the first optical peak (RS peak) gives the fireball deceleration time (Zhang et al. 2003), which can be used to estimate the initial Lorentz factor of the GRB, Γ 0 300.The second optical peak appears when ν m crosses the optical band.This peak flux can be used to constrain the values of jet isotropic kinetic energy E K,iso , ISM density n 1 and the microphysics parameter B,f .Based on the fit with this scenario, we infer that the magnetic field strength ratio in RS and FS is B r /B f ∼ 9.5, suggesting a magnetic flux carried by the ejecta from the central engine.a magnetized central engine model, i.e., a BH central engine powered by the BZ process or a magnetar central engine, is preferred in GRB 190106A.
It is worth to note that the first optical peak, like GRB 041219A, can be interpreted with a reverse shock (Fan et al. 2005), internal shocks (Vestrand et al. 2005;Wei 2007), or tail emission (Panaitescu 2020).The last two models assume that the prompt optical and γ-ray emissions have a common origin.However, the radiation mechanism for the prompt γ-ray emission is still poorly understood (Kumar & McMahon 2008).
The shallow decays in both X-ray and optical bands demand a nearly constant late-time energy injection with L inj 4.5 × 10 49 erg s −1 and q ∼ 0.05 lasting from ∼ 300 s to ∼ 1500 s.We investigated several energy injection sources, such as the spin-down of a magnetar, hyperaccreting BH powered by neutrino annihilation, and BH powered by BZ mechanism.The neutrino annihilation mechanism is less efficient for powering the jet at shallow decay phase.Both magnetar and BH central engine with BZ mechanism can give rise to the energy injection required.In general, the flat light curve feature favors the magnetar central engine model.A BH-accreting central engine can also makes a flat light curve, but requires a very small value of viscosity parameter (and thus a very slow accretion) (Lei et al. 2017).On the other hand, the disk will be dominated by advection and contains a strong wind at his afterglow phase.The light curve will then become steeper due to the suppression by the wind (Lei et al. 2017).
In Section 4.1, we divided the X-ray light curve into three stages (steep decay, shallow decay, late decay).As a result, the break time between shallow decay and late decay in X-ray is significantly different with that in optical band (1.6 +0.7 −0.7 × 10 4 s vs 7.7 +1.1 −1.1 × 10 4 s), as show in Table 1.However, the X-ray break time will be roughly consistent with optical if we fit the X-ray light curve with four segments (steep decay, plateau, normal decay, late decay).Considering this uncertainty, we attribute such break to the jet break, revealing an opening angle ∼ 4.2 • .From the observations and modeling, we found that the total jet energy is E total = E γ,iso + E K,iso 1.9 × 10 53 erg.The opening angle-corrected jet energy will be E j ∼ 5 × 10 50 erg, which is well below the maximum rotational energy of 3 × 10 52 erg (Lattimer & Prakash 2016) − 7 × 10 52 erg (Haensel et al. 2009) for a standard neutron star with mass M ∼ 1.4 M .Therefore, our data do not require a BH as the central engine of this GRB.

SUMMARY
We present our optical observation of the afterglow of GRB 190106A with the NEXT telescope, and the Xinglong 2.16-m telescope equipped with BFOSC.With the spectrum obtained from BFOSC, we measured the redshift of the burst and calculated the Line Strength Parameter.The unusual optical light curve obtained by NEXT and Xinglong 2.16-m, combined with the MASTER observations, show a normal decay light curve with early double peak and late break.Therefore, the optical light curve can be divided into five stages: rapid rise in the early stage, followed by decay, then a slow rise reaching a second peak at around 200 s, followed by a shallow decay (plateau), and faster decay from about 1 d onward.The evolution of the X-ray light curve can be divided more simply into three stages: rapid decay, plateau and normal decay.
The multi-band observations, especially the early optical data of GRB 190106A, provide rich information to study the nature of GRB.It is found that GRB 190101A can be successfully explained by the reverse-forward shock model.Our conclusions are summarized as follows: 1.The redshift of the burst is z = 1.861 ± 0.002.The Line Strength Parameter, which is very similar to the average value of the sample.
2. The studies with the Amati and Yonetoku diagrams indicate that GRB 190106A is a typical type II burst.
3. The early double-peak optical light curve can be well interpreted with combination of the RS and FS models.
4. The modeling with RS and FS models indicate that the magnetic field strength ratio in RS and FS is B r /B f ∼ 9.5, suggesting a magnetic flux carried by the ejecta from the central engine.This favors a strongly magnetized central engine model, such as a magnetar central engine model or a BH central engine powered by the BZ process.
5. The initial Lorentz factor can be estimated with the first optical peak time (or deceleration time) as Γ 0 ∼ 300.
6.The multi-band afterglow data can be interpreted with external shock model.From the observations and modeling, we found that the total jet energy is E total = E γ,iso + E K,iso 1.9 × 10 53 erg.
7. The shallow decays in both X-ray and optical bands are explained by nearly constant late-time energy injection with L inj 4.5 × 10 49 erg s −1 and q ∼ 0.05 lasting for ∼ 1000 s.Both magnetar and BH central engine models can give rise to such amount of injected energy, but the flat light curve feature favors a magnetar central engine.
8. The breaks at a few ×10 4 s in both X-ray and optical bands are roughly consistent with the jet break, revealing an opening angle ∼ 4.2 • .The opening angle-corrected jet energy will be E j ∼ 5 × 10 50 erg.

Figure 1 .
Figure 1.The R-band position of GRB 190106A within the field of view 30 × 30 .The location of the burst is circled in red.North is up and East is left.Left panel: the DESI Legacy Imaging Surveys r band image(Dey et al. 2019).A weak but obvious source with r ∼ 24.0 mag is located at the GRB position, which could be the host galaxy of the burst.Middle panel: the first NEXT image, obtained 0.39 hours after the BAT trigger, shows a bright source consistent with the UVOT coordinates with R = 16.9 mag.Right panel: the NEXT image obtained at 3.46 hours after the burst.The source in the image center is obviously fainter than the one in the middle panel, with a brightness of R = 18.4 mag.

Figure 2 .
Figure 2. The light curves of GRB 190106A in γ-ray, X-ray and optical R-band.The X-ray light curve colored in cyan is 0.3-10 keV unabsorbed flux, obtained from UK Swift Science Data Centre (Evans et al. 2009).The gray circles and gray squares are modified MASTER and Mondy photometric results, respectively.The red diamond points are our observations obtained by NEXT and Xinglong 2.16-m.Most of the error bars are smaller than the marker, and our photometric data are listed in the Table.4.

Figure 4 .
Figure 4. Comparison of a large sample of GRB optical afterglow light curves shifted in time and flux to a common redshift of z = 1 following Kann et al. (2010).The gray background are historical data of other GRB light curves.GRB 190106A is shown as a black solid line, which lies at the middle of the distribution in terms of luminosity.

Figure 5 .
Figure 5.The equivalent widths of GRB 190106A obtained from BFOSC are colored in red.These can be compared to the average and standard deviation of strengths of a sample of GRB afterglows in black color, as described by de Ugarte Postigo et al. (2012).
Figure6.The prompt phase SED combined with Optical (black), XRT (red) and BAT (green).The BAT time interval is from 71 s to 81 s, while the XRT is from 85 s to 95 s.The satellite was still slewing before 85 s, so the XRT data before 85 s are ignored.The optical point is converted from the MASTER observation and not included in the fitting.The data points of BAT and XRT were fitted with a single power-law function, and the spectral indices are 0.77± 0.04 and 0.72 ± 0.10, respectively.

Figure 9 .
Figure 9. Left panel: The Amati diagram of GRB 190106A (black star).Dashed and dotted lines represent 1σ and 2σ confidence levels, respectively.The ellipses colored in cyan and green represent type I and type II burst confidence, respectively.The data points for type I and type II GRBs are taken from Minaev & Pozanenko (2020).The burst lies within the 1σ confidence region of type II bursts, consistent with the Ep − Eiso correlation.Right panel: The Yonetoku diagram of GRB 190106A (black star).The data points for type I and type II GRBs are taken from Nava et al. (2012) and Zitouni et al. (2022).
and the photometric results are presented in Table4.The spectrum obtained with the Xinglong 2.16-m/BFOSC.Grey is the original spectrum and the blue spectrum has smoothed for display purposes.The identified metal absorption lines are indicated with vertical lines in the figure.Telluric features are marked with the telluric symbol.

Table 1 .
List of X-ray and optical light curve decay indices.

Table 2 .
The start, end and peak times of each episode.

Table 3 .
Values of parameters adopted for interpreting the broadband data of GRB 190106A.

Table 5 .
List of features in the spectra and their observatory frame equivalent widths.