The study of multi-peaked type-I X-ray bursts in the neutron-star low mass X-ray binary 4U 1636$-$536 with RXTE

We have found and analysed 16 multi-peaked type-I bursts from the neutron-star low mass X-ray binary 4U 1636$-$53 with the Rossi X-ray Timing Explorer (RXTE). One of the bursts is a rare quadruple-peaked burst which was not previously reported. All 16 bursts show a multi-peaked structure not only in the X-ray light curves but also in the bolometric light curves. Most of the multi-peaked bursts appear in observations during the transition from the hard to the soft state in the colour-colour diagram. We find an anti-correlation between the second peak flux and the separation time between two peaks. We also find that in the double-peaked bursts the peak-flux ratio and the temperature of the thermal component in the pre-burst spectra are correlated. This indicates that the double-peaked structure in the light curve of the bursts may be affected by enhanced accretion rate in the disc, or increased temperature of the neutron star.


INTRODUCTION
Thermonuclear (Type I) X-ray bursts show a sudden increase in X-ray intensity, becoming ∼ 10 − 100 times brighter than the persistent level, triggered by unstable ignition of accreted fuel on the surface of an accreting neutron star (NS) in low mass X-ray binaries (LMXBs) (Galloway et al. 2008;Galloway & Keek 2017). Type I X-ray bursts were first detected in 1975 in the binary 3 1820 − 30 in the globular cluster NGC 6624 (Grindlay 1976); subsequently a growing population of bursters has been observed by different Xray satellites (Galloway & Keek 2017). In a typical X-ray burst, the light curve shows a single-peaked profile with a fast rise (∼ 1 − 5 s) and an exponential decay within 10 − 100 s (Lewin et al. 1993;Strohmayer & Bildsten 2006;Galloway et al. 2008).
Besides the single-peaked normal bursts, multi-peaked bursts have also been reported in previous studies. Double-peaked bursts have been reported in several NS-LMXBs, e.g., 4U 1608−52 (Penninx et al. 1989), GX 17+2 (Kuulkers et al. 2002), 4U 1709−267 (Jonker et al. 2004) and MXB 1730−335 (Bagnoli et al. 2014). With the Rossi X-ray Timing Explorer (RXTE), Watts & Maurer (2007) analyzed 4 double-peaked bursts in 4U 1636−53. Still more rare, ★ E-mail: lichen@ynao.ac.cn † E-mail: zhangguobao@ynao.ac.cn bursts with triple-peaked structure have also been observed in 4U 1636−53 (van Paradĳs et al. 1986;Zhang et al. 2009). While investigating the cooling phase of X-ray bursts in 4U 1636−53, Zhang et al. (2011) reported 12 double-peaked bursts and found that most of them appeared at the vertex of colour-colour diagram. Recently, there were two new observations of double-peaked burst, one in the soft spectral state in 4U 1608−52 (Jaisawal et al. 2019), and another one in SAX J1808.4−3658 (Bult et al. 2019), both using the Neutron Star Interior Composition Explorer (NICER). The double-peaked structures in the burst light curve can be separated into two groups. The first one consists of bursts with a double-peaked profile in X-rays but a single-peaked profile in the bolometric lightcurve, generally accompanied by photospheric radius expansion (PRE), where the flux of the burst reaches the Eddington luminosity. In this case, the temperature of the photosphere temporarily shifts out of the instrument passband, causing an apparent dip in the observed X-ray light curve (Paczynski 1983). The other group consists of bursts that have a double-peaked profile both in X-rays and the bolometric lightcurve. Most of these bursts have low peak flux, although PRE bursts with double-peaked profiles both in X-rays and bolometric luminosity have been recently observed with NICER (Jaisawal et al. 2019;Bult et al. 2019).
Several theoretical models have been proposed to explain the double-peaked bursts. Fujimoto et al. (1988) proposed a model of stepped thermonuclear energy generation due to shear instabilities in the fuel on the NS surface. Melia & Zylstra (1992) suggested that the double-peaked bursts are due to the scattering of the X-ray emission by material evaporated from the disk during the burst. These models, however, can not reproduce the observed doublepeaked profiles in the light curve, black-body temperature and radius (Bhattacharyya & Strohmayer 2006a). Fisker et al. (2004) suggested that a waiting point impedes the nuclear reaction flow and causes a stepped release of thermonuclear energy, but this idea has difficulties in explaining the large dips observed between the two peaks (Bhattacharyya & Strohmayer 2006a). The thermonuclear flame spreading model provided by Bhattacharyya & Strohmayer (2006a) suggested that the double-peaked structure is caused by high latitude ignition and stalling approaching the equator. This model qualitatively explains the essential features of the light curve and reproduces the spectral evolution of two double-peaked bursts in 4U 1636−53. However, one of the problems of the flame spreading model is that it can not explain the triple-peaked bursts (Zhang et al. 2009). Lampe et al. (2016) in their simulations found that low accretion rate and high metallicity could affect the burst morphology and produce twin-peaked structure when a large amount of hydrogen has been depleted. Recently, Bult et al. (2019) suggested that the bright double-peaked bursts are due to the local Eddington limits associated with the hydrogen and helium layers of the NS envelope. Understanding these mechanisms is important, because current models for multi-peaked X-ray bursts have met with only partial success in explaining their light curves and temperature profiles.
Burst properties in an individual system depend mainly on accretion rate (Fujimoto et al. 1981;Bildsten 2000;Zhang et al. 2011;Galloway & Keek 2017). For a specific source, given a certain global accretion rate, the local accretion rate varies with latitude, being higher at the equator and lower at high latitude (Cooper & Narayan 2007). The ignition latitude depends on the column depth related to accretion rate to trigger a burst. As the increasing global accretion rate, the lower latitude region firstly reach the critical local accretion rate and the fuel become stable burning, the ignition should occur at higher latitude even pole. Concerning the research of Cooper & Narayan (2007) with the thermonuclear spreading model of Bhattacharyya & Strohmayer (2006a), Watts & Maurer (2007) expected to find more double-peaked bursts at higher global accretion rates than the single-peaked bursts. However, Watts & Maurer (2007) presented the analysis of the accretion rate limited of 4 double-peaked bursts and posed a challenge to the above expectation. In this paper, we collect a large sample to provide a more complete description of the observational features of multi-peaked bursts and further discuss the relation between accretion rate and multi-peaked structure.
In the time-resolved spectral analysis, the standard approach is to fit the X-ray burst spectra by assuming a constant persistent emission (non-burst component) during the burst (Galloway et al. 2008). However, recent studies provide evidence of enhanced accretion during type-I X-ray bursts (Worpel et al. 2013(Worpel et al. , 2015. The method gives improvements in the quality of spectral fit compared to the standard approach, and the analysis is sensitive to changes in the persistent spectrum in the 2.5 − 20 keV (Worpel et al. 2015). In this paper, we adopt both the standard approach and the method to analyze our sample of multi-peaked bursts.
The LMXB 4U 1636−53 is one of the best-studied sources of X-ray bursts. The NS is in a binary system in a 3.8 hr orbit (van Paradĳs et al. 1990) with an 18th magnitude blue star companion (Galloway et al. 2008), and the spin period of the NS is 581 Hz (Strohmayer et al. 1998a,b). 4U 1636−53 is an Atoll source, and as the source moves in the colour-colour diagram (hereafter CCD) from the top right to bottom right, the accretion rate gradually increases, with a transition from the Island to the Banana state (Hasinger & van der Klis 1989). The single peak bursts show a uniform distribution in the CCD (Zhang et al. 2011). About a dozen double-peaked and two triple-peaked X-ray bursts have been discovered from this source using different satellites (Sztajno et al. 1985;van Paradĳs et al. 1986;Lewin et al. 1987;Bhattacharyya & Strohmayer 2006a,b;Watts & Maurer 2007;Galloway et al. 2008;Zhang et al. 2009Zhang et al. , 2011. The large multi-peaked bursts sample makes 4U 1636−53 an ideal source to study the properties and evolution of this kind of bursts.
The structure of this paper is organised as follows. In Section 2 we describe the data analysis of our sample. In Section 3 we show the results of light curves and spectra. In Section 4 we discuss our findings in the context of previous theoretical work.

DATA ANALYSIS
The Rossi X-ray Timing Explorer (RXTE) was launched in 1995, and operated until 2012 with a circular orbit at an altitude of 580 km, correspoding to an orbital period of about 96 min (Bradt et al. 1993). We analysed all archived data from the proportional Counter Array (PCA) which is the main instrument onboard RXTE. The PCA consists of five collimated proportional counter units (PCUs), which are sensitive in the 2−60 keV energy range with an energy resolution of ∼ 1 keV at 6 keV (Jahoda et al. 2006). For each observation we used the Standard2 data (16-s time resolution and 129 energy channels) to calculate X-ray colours. We used the Standard1 mode (only the PCU2) to produce the burst light curves. For the timeresolved spectral analysis of the bursts, we extracted spectra in 64 channels from the Event data of all available PCUs.
We studied 336 type I X-ray bursts (as in Zhang et al. 2013) in LMXB 4U 1636−53 with RXTE and discover 16 multi-peaked bursts. The 0.125s bin light curves of these 16 bursts in the 2 − 60 keV are shown in Figure 1. Following the procedure in Zhang et al. (2011) and Zhang et al. (2013), we searched the 1s Standard1 light curve for burst visually, and considered that the start time of a burst is when the flux is larger than 3 times the 1 error of the average persistent flux. We find 14 double-peaked bursts (4 bursts has investigated by Watts & Maurer 2007), one triple-peaked burst (Zhang et al. 2009) and one quadruple-peaked burst which was not reported in previous work. We divided the 16 bursts into two classes according to the number of peaks. All the double-peaked bursts are classified as Class 1 bursts and the bursts with more than two peaks are classified as Class 2 bursts.
To study these bursts in detail, we introduce several parameters to characterise the burst light curve (see Table 1 and Table 2). In Table 1, the column is the burst number sorted by the value of peak flux ratio, 1,2 (see details later). Because most of the data around the peak show a symmetric distribution, we used a Gaussian function to fit the data around each single peak in each burst light curve to get the peaking time. The quantity ,1 is the peak time of the first peak, is the separation time between the first and second peak, and ,2 is the peak time of the second peak. In Table 2, we list the characteristics of the burst light curve in Class 2 bursts. The fitting process of triple-peaked burst is similar to that of the double-peaked bursts. About the quadruple-peaked burst, we have a more detailed discussion in Sec 3.5. The quantities ,3 and ,4 are the peak time of the third and fourth peak, respectively, 23 and 34 are the separation time between the second and the third peak, as well as between the third and the fourth peak, respectively. For all of these parameters we give the 1 error.
To trace the spectral state of the source when these multipeaked bursts appear, we made a CCD as shown in Figure 2. We defined the soft colour as the ratio of the count rate in the 3.5 − 6.0 keV to the count rate in the 2.0 − 3.5 keV bands, and the hard colour as the ratio of the count rate in the 9.7 − 16.0 keV to the count rate in the 6.0 − 9.7 keV bands (Zhang et al. 2009). The colours of the source are normalised by the Crab. For the colours of the source before the burst, we used 64-s of the pre-burst spectrum. In Figure  2, the grey points represent all available observations. The black crosses represent all the bursts in this source. The red filled circles stand for Class 1 bursts, and the blue filled triangles indicate the Class 2 bursts. We find that most of the multiple-peaked bursts are located close to the vertex of the CCD.
We extracted the 64-s interval spectrum prior to a burst as persistent emission. We generated the instrument response matrix using the tool pcarsp and the instrumental background using the tool pcabackest in HEAsoft for each spectrum. In this work, we fitted the spectrum in the 3.0 − 20 keV band using XSPEC version 12.10.1 (Arnaud 1996). We added a 0.5% systematic error to the pre-burst spectra because of calibration uncertainties. During the fitting process, we included the effect of interstellar absorption using the cross-sections of Balucinska-Church & McCammon (1992) and solar abundances from Anders & Grevesse (1989), with a fixed hydrogen column density, H , of 0.36 × 10 22 cm −2 (Pandel et al. 2008). Table 1. Properties of the light curves of Class 1 bursts in 4U 1636−53. The column gives burst number sorted by the value of the peak flux ratio ( 1,2 ) which is the ratio of the first peak to the second peak flux. We used a Gaussian function to fit each peak of every burst. The quantify ,1 gives the peak time of the first peak, is the separation time between the first and the second peak, ,2 is the peak time of the second peak.  Table 2. Properties of the light curves of Class 2 bursts in 4U 1636−53. The quantities ,3 and ,4 are the peak time of the third and fourth peak, respectively, 23 and 34 are the separation time between the second and the third peak, as well as between the third and the fourth peak, respectively. For the triple-peaked burst, we use the same method that we used in double-peaked bursts. For the quadruple-peaked burst, we use two GAUSSIAN functions adding two BURS models to fit all data of the light curve (see detailed discussion in Sec 3.5).

Start time (UTC)
Obsid   Grey dashed line shows the ratio 1,2 = 1. All these data are given in Table  1.
We used all available PCUs during the X-ray bursts to produce time-resolved spectrum. We corrected every spectrum for dead time according to the methods supplied by the RXTE team. Since the light curve of bursts decay is quite smooth, to compensate the lower count rates, we extracted spectrum over longer intervals in the tail of the bursts. For each multiple-peaked burst, we generated one instrument response matrix using the pcarsp and the instrumental background using the pcabackest in HEAsoft.
For the burst spectral analysis, we initially used a singletemperature blackbody model, TBabs*bbodyrad, to fit the net burst spectra, which is well established as a standard procedure in X-ray burst analysis (Kuulkers et al. 2002;Galloway et al. 2008). The model provides the blackbody colour temperature, bb , and the normalization, bb , proportional to the square of the blackbody radius of the burst emission surface, and it allows us to estimate the bolometric luminosity as a function of time assuming a distance of 5.95 kpc (Fiocchi et al. 2006). The bursts bolometric flux are calculated as where bb = 2 km / 2 10 , km is the effective radius of the emitting surface in km, and 10 is distance to the source in units of 10 kpc. We defined 1,2 as the ratio of the first peak to the second peak flux. We note that, to reduce the effect of the instrument, we used the bolometric flux (not the net burst count rate ratio in Figure 1) from the standard approach to calculate the ratio 1,2 .
After subtracting the instrumental background for the timeresolved spectra, we used another model that allows us to vary the pre-burst spectrum by a free scaling factor (Worpel et al. 2013). We used two models to describe the persistent emission: TBabs*(bbodyrad + powerlaw) and TBabs*(diskbb + powerlaw) in XSPEC. The fit results for the above two models are shown in Tables 3 and 4, respectively. In Table 3, bb is the blackbody temperature and bb is the normalization of the balckbody, is the power-law index and is the normalization of the power law. The spectral parameters of our best-fitting model are given with 1 error. In Table 4, dbb is the disk blackbody temperature and dbb is the normalization of the disk balckbody. Comparing the fitting results of the two pre-burst spectra models (see Table 3 and Table 4), we selected the TBabs*(bbodyrad + powerlaw) as the best pre-burst model, so we used TBabs*( *(bbodyrad + powerlaw) + bbodyrad) in our analysis (the so-called method) to re-fit the net burst spectra. The parameter was allowed to vary between −100 and 100 during the fits.
In Tables 3 and 4, p represents the persistent unabsorbed flux in the 2.5 − 25 keV band and is the luminosity of the source before the burst in Eddington units. We used the unabsorbed 2.5 − 25 keV flux, p , a bolometric correction factor bol = 1.2 (Galloway et al. 2008), and the source distance of 5.95 kpc to estimate the X-ray persistent luminosity, p . In the calculation of Eddington luminosity we assume that the ratio of the colour temperature to the effective temperature, / , is 1.4 (Madej et al. 2004), the NS mass is 1.4 and the hydrogen mass fraction X = 0.7. We note that, in the model of TBabs*(bbodyrad+powerlaw) corresponds to bb and in the model of TBabs*(diskbb+powerlaw) corresponds to dbb . Figure 1 shows the 14 bursts with double-peaked profiles. The bursts display large or gentle dips in their light curves. The first peak can be either weaker/shorter or stronger/longer than the second peak.

Multi-peaked bursts light curves
In Figure 3 we show the rise time, ,1 , of the first peak and the separation time, , of the two peaks as a function of the peak-flux ratio, 1,2 . The grey dashed line represents 1,2 = 1. The rise time of the first peak is in the range of ∼ 1 − 6s, and the separation in time between the two peaks is in the range of ∼ 3 − 7s. We find that: (1) in the case of 1,2 < 1, as the peak flux ratio increases, both the rise time of the first peak and the separation in time increase; (2) in the case of 1,2 > 1, the above parameters do not depend upon peak flux ratio. Besides the triple-peaked burst (burst #15) reported by Zhang et al. (2009), we also discovered a ∼ 40 s long quadruple-peaked burst (burst #16; date of observation: 2006 August 17), which was not reported in previous studies. Figure 2 shows the position of all multi-peaked bursts in the CCD of 4U 1636−53. The normal, single-peaked, bursts observed with RXTE are distributed more or less uniformly across the CCD. However, except for burst #1 and #8, most of the multiple-peaked bursts are located close to the vertex of the CCD. Table 3. Best-fitting parameters of the persistent emission before the 16 multi-peaked bursts in 4U 1636−53 with the model TBabs*(bbodyrad+powerlaw). The quantity bb is the blackbody temperature, bb is the normalisation of the blackbody, is the power-law index, is the normalisation of the power law, p represents the persistent unabsorbed flux in the 2.5 − 25 keV band and is the luminosity of the source before the burst in Eddington units. The spectral parameters of our best-fitting model are given with 1 error.  Table 4. Best-fitting parameters of the persistent emission before the 16 multi-peaked bursts in 4U 1636−53 with the model TBabs*(diskbb+powerlaw). The quantity dbb is the disc blackbody temperature, dbb is the normalisation of the disk blackbody. is the power-law index, is the normalisation of the power law, p represents the persistent unabsorbed flux in the 2.5 − 25 keV band and is the luminosity of the source before the burst in Eddington units. The spectral parameters of our best-fitting model are given with 1 error. With the lowest peak flux ratio ( 1,2 ∼ 0.3) and highest peak flux (∼ 5.2 × 10 −8 ergs cm −2 s −1 using the standard approach), burst #1 is located at the bottom right, when the source was in the so-called upper banana branch in the CCD. In the light curve of burst #1, the first peak is relative weak ( 1,2 ∼ 0.3) and the whole burst is dominated by the second peak. Burst #1 has also the second shortest rise time ( ,1 ∼ 1.6 s) of the first peak and shortest separation time ( ∼ 4.0 s).

Multi-peaked bursts in the CCD
Burst #8 is located at the position between the Banana and the Island states in the CCD. Different from burst #1, burst #8 has the longest rise time ( ,1 ∼ 5.0 s) of the first peak and the longest separation time ( ∼ 6.6 s), and a high 1,2 (∼ 0.8).

The relation between pre-burst spectrum and burst profile
To investigate the relation between the pre-burst spectra and the multi-peaked burst light curve profiles, we compared the spectral parameters of the persistent emission of Class 1 bursts with two properties of the light curve: the peak flux ratio, 1,2 , and the separation time of two peaks, . Figure 4 shows the spectral parameters of the persistent emission (in two models) against the peak flux ratio. The red dots represent the Class 1 bursts, and the blue triangles represent Class 2 bursts. The left panels (a) and (b) in Figure 4 display, respectively, the blackbody temperature, bb , and power-law index, , against the peak flux ratio, 1,2 , for the model TBabs*(bbodyrad+powerlaw). The blackbody temperature ranges from 1.57 keV to 1.87 keV. There appears to be a positive correlation between the blackbody temperature and the peak flux ratio. In order to check that, we firstly fit these data with a constant model and get 2 = 24.6 for 13 d.o.f. We then fit these data with a line function, and get a slope of 0.07 ± 0.02 with a 2 = 13.8 for 12 d.o.f. The F-test probability for these two fits is 0.009, indicating that a linear function is slightly better than a constant. The right panels (c) and (d) in Figure 4 show, respectively, the disk-blackbody temperature, dbb , and power-law index, , against the peak flux ratio, 1,2 for the model TBabs*(diskbb+powerlaw). For panel (c) of Figure 4, we do the same analysis as in panel (a), obtaining a 2 = 31.7 for 13 d.o.f. for a constant model and a slope of 0.09 ± 0.03 with a reduced chi-square of 1.64 ( 2 / = 19.7/12) for the line function. The F-test probability for these two fits is 0.019. The above analysis indicates that the peak flux ratio, 1,2 , is marginally correlated with the temperature of the thermal component in the pre-burst spectra.

Time-resolved burst spectra
We adopted both the standard approach and the method to fit the time-resolved burst spectra. An example (burst #5) of the bestfitting parameters using standard approach is shown in the left panel of Figure 5. From the top to the third panel, we display the blackbody temperature, blackbody radius and bolometric flux as a function of time. We used a 0.5-s bin light curve in the 2 − 60 keV range in the fourth panel. The double-peaked structures occur simultaneously in the X-ray and bolometric light curves, which indicates that the double-peaked profiles during the X-ray bursts is not due to a passband effect of the instrument. The time-resolved temperature shows two local maxima, as is also the case in the bolometric flux curve, however, the peaking time in temperature is different from the peaking time in bolometric flux curve. After initially growing, the radius continues increasing following a dip, and then remains more or less constant. The double-peaked structures are also apparent in the bolometric flux in the rest of the bursts of Class 1 bursts. We then measured the local peak temperature of the doublepeaked bursts in the standard approach. The peak temperature ratio, represents the ratio of the first peak to the second peak temperature. Figure 6 shows the peak temperature ratios against the peak flux ratios in the 14 double-peaked bursts. The vertical and horizontal dashed lines represent 1,2 =1.0 and =1.2, respectively. When 1,2 > 1, the peak temperature ratios are always larger than 1.2. In the case of 1,2 < 1, the first peak temperature can be higher or lower than the second peak temperature and is always less than 1.2. There appears to be a bimodal distribution in the temperature ratios. Figure 7 shows the first and second peak flux in the standard approach against the duration time between two peaks, respectively. The typical PRE peak flux is in the range of (6 − 8) × 10 −8 ergs cm −2 s −1 in 4U 1636−53 (Lyu et al. 2015). In our present sample, there are no PRE events. There is no significant trend between the first peak flux and the duration time, however, an anti-correlation is present between the second peak flux and the separation time between two peaks.
We show the fits results of burst #5 with the method in the right-hand panel of Figure 5. As in the standard approach, both the bolometric flux and the X-ray light curve show a double-peaked structure. The evolution of black body temperature and radius in the method is similar to that in the standard approach. The increased black body radius in the cooling phase may be explained in two ways. One would be the influence of persistent emission (accretion geometry) in the soft state (Kajava et al. 2014), the other would be different canonical composition in the NS atmosphere (Suleimanov et al. 2011;Zhang et al. 2011). The plotted in the bottom panel shows values larger than one during the whole burst. Due to the new parameter , the bolometric flux is lower than in the standard approach and has larger error bars. The horizontal solid line in the bottom panels stands for =1. In general, the reduced 2 distribution in the method is lower than that in the standard approach.
Finally, we investigated the time evolution of these 16 multipeaked bursts in Figure 8 using the same order as in Figure 1. The horizontal black dashed lines correspond to = 1. In all multipeaked bursts, the values vary in the range of ∼ 1 − 7 during the bursting period, which is similar to the range observed by Worpel et al. (2013) in other bursts ( ∼ 2 − 10). In general, high values of coincide with large flux values during the burst. We study the correlation between and bolometric flux (e.g. blackbody temperature and radius) for the 16 multi-peaked bursts using the method of cross-correlation lags (Peterson et al. 1998;Sun et al. 2018). We do not find any obvious correlation between and burst spectral parameters, and also we do not find any evidence of a relation between the centroid of the cross-correlation function (time delay between and bolometric flux) and the double-peaked structure (e.g. the peak-flux ratio). .

The quadruple-peaked X-ray burst
In Figure 9, we show the light curve of the quadruple-peaked X-ray burst. We use two GAUSSIAN functions adding two BURS models to fit the quadruple-peaked burst light curve data, where BURS represents a model to describe the burst component in QDP 1 .
where is time, is the photons per second, and represent the start and peak time of each peak, respectively, and and stand for the normalisation and decay factor, respectively.
In the top panel of Figure 9, we show the fitting results, getting a 2 = 120.64 for 62 d.o.f. In the bottom, we show the residual of the fit in units of the error, and the black horizontal line represents the (data-fit)/err = 0. To verify the existence of the third peak, we also use three components (one GAUSSIAN function adding two BURS models) to re-fit the light curve, getting a 2 = 386.79 for 65 d.o.f. The F-test probability for these two fits is 1.1 × 10 −15 , which indicates that the third peak is significant. This is the first time that a quadruple-peaked burst is reported in 4U 1636−53.
There are three relatively strong peaks (∼ 3s, ∼ 7s, ∼ 22s for the first, second, fourth peak) in the light curve and one weak peak (∼ 19s for the third peak). There is a very long time separation (∼ 12s) between the second and the third peak. The separation time between the first and last peak in this quadruple-peaked burst is ∼ 18s, which is similar to the separation (∼ 17s) in the triplepeaked burst reported by van Paradĳs et al. (1986), but about two times longer than that in the triple-peaked burst (∼ 8s) in Zhang et al. (2009).
We investigated the time-resolved spectra of the quadruplepeaked burst with two models. In the standard approach (the lefthand panel of Figure 10), both the bolometric light curve and X- (d) Figure 4. Parameter of the persistent spectrum before the multi-peaked burst in 4U 1636−53. Left panels: blackbody temperature ( bb ), and power-law index ( ) for the model TBabs*(bbodyrad+powerlaw) against the peak flux ratio ( 1,2 ); right panels: the disk blackbody temperature ( dbb ) and the power-law index ( ) for the model TBabs*(diskbb+powerlaw) against the peak flux ( 1,2 ). Red dots and blue triangles represent Class1 and Class 2 bursts, respectively. We use a linear (red line) and constant (black line) function to fit only the Class 1 bursts data.
ray light curve show four peaks, three pronounced and small one. The total burst fluence is (26.9 ± 3.5) × 10 −8 ergs cm −2 and the bolometric first peak flux is (2.0±0.6) ×10 −8 ergs cm −2 s −1 . There is a ∼ 5s long waiting period from 13s to 18s in the light curve, and the bolometric flux is always higher than the persistent emission. The effective temperature profile, after an initially rising, decreases with time and shows three noticeable local maxima. The blackbody radius increases steadily during the whole burst, which is similar to what was observed in the triple-peaked burst in Zhang et al. (2009). In the method (the right-hand panel of Figure 10), both the bolometric flux and X-ray light curve display four peaks as in the standard approach. Similar to the double-peaked bursts, the bolometric flux in the method is lower than that of the standard approach due to the variation being larger than one during the first two peaks. The evolution of the blackbody radius and the blackbody temperature in the method is similar to that in the standard approach.

DISCUSSION
We analyzed all available RXTE data of the LMXB 4U 1636−53, and we found 14 double-peaked bursts with complex profiles, one triple-peaked burst, and one quadruple-peaked burst that had not been reported before. The bolometric flux light curve of the multiple-peaked bursts exhibits a similar profile as the X-ray light curve, which indicates that the multi-peaked structures in 4U 1636−53 are not due to passband instrumental effect (Jaisawal et al. 2019). Between the peaks, the flux of these multi-peaked bursts never drops near or below the pre-burst level, which is different from the triple bursts in Boirin et al. (2007) and Beri et al. (2019).
We find a marginal positive correlation between the peakseparation time and the peak-flux ratio in Class 1 bursts (bottom panel of Figure 3). This suggests that double-peaked bursts with high peak-flux ratio value have a longer peak separation time. This phenomenon may be explained if the burst consumes less fuel during the first peak, such that it is easier to release subsequently the fuel that is left to trigger the second peak. Therefore, the bursts with higher peak flux ratio need more time to accumulate fuel to release the second peak energy. As we know, in the bottom panel of Figure  7, we also find an anti-correlation between the second peak flux and the separation time between two consecutive peaks. This indicates that the double-peaked bursts with longer peak separation time have a weaker second peak. If we assume that the mass accretion is negligible during the burst, this anti-correlation indicates that when the separation time is long, there is less fuel for the second peak. These results indicate that a single peak in these double-peaked bursts is not isolated, and the double-peaked profile is affected by the strength and separation time of the two single peaks.
Most of the multi-peaked bursts in our sample appear during the transition from the hard to the soft state in the CCD where PRE bursts are also present (Watts & Maurer 2007;Zhang et al. 2009). That the PRE and multi-peaked bursts appear in the same state of the source may provide an important clue to understanding the origin of the multi-peaked bursts, indicating that the appearance of the multi-peaked bursts could be affected by the transition of the source from the hard to the soft state. If mass accretion rate increases from the upper right to the lower right in the CCD, there is no apparent relation between mass accretion rate and the parameter 1,2 . We have checked the correlation between the persistent flux and the double-peaked burst parameters (peak ratio 1,2 ; peak separation ), but do not find a clear trend between them either. To further investigate if mass accretion affects the burst light curve structure, we studied the pre-burst spectra. We find that there is a positive correlation between the peak-flux ratio and the temperature of thermal component in the pre-burst spectra for the doublepeaked bursts, where we use a blackbody model to fit the burst time-resolved spectra. The soft thermal component of X-ray spectra in accreting NS-LMXBs is generally explained by the emission from the NS surface and the accretion disc. Sanna et al. (2013)  analyzed the X-ray spectra of six observations of 4U 1636−53 in different spectral states taken with XMM-Newton and RXTE simultaneously. They find that the temperature at the inner disc radius is ∼ 0.1 − 0.8 keV, and the temperature at the NS surface is ∼ 1.4 − 2.0 keV. We used a black-body or a disk black-body to depict the enhancement of the thermal component. The emission from the NS surface and inner region of the accretion disc are degenerate in our persistent spectrum. The high blackbody (or disc blackbody) temperature could be due to an enhanced accretion rate in the disc or an increased temperature of the NS surface. Since the above discus- kT Figure 6. The peak temperature ratio ( ) against the peak flux ratio ( 1,2 ) for the double-peaked bursts in 4U 1636−53 obtained from fits using the standard procedure to fit time-resolved spectra of X-ray bursts. The vertical and horizontal dashed lines correspond to the 1,2 =1.0, =1.2 respectively. We note that the local maxima of temperature and flux do not necessarily occur at the same time.   sion indicates that the double-peaked properties are not affected by accretion rate, we suggest that the double-peaked bursts with high peak flux ratio might appear when the NS surface temperature is high.
To date, there are several models to explain the unusual doublepeaked structure in X-ray bursts. One of them is the thermonuclear flame spreading model (Bhattacharyya & Strohmayer 2006a,b), which succeeded in explaining the large dip in the X-ray light curve and reproducing the spectral evolution of double-peaked bursts. Cooper & Narayan (2007) suggested that ignition latitude has a positive correlation with the accretion rate. Combing these ideas with the results in Figure 4 and assuming that the temperature of the NS is correlated to the accretion rate, the high peak flux ratio of double-peaked bursts would correspond to high latitude ignition (Bhattacharyya & Strohmayer 2006a,b).
Alternatively, Lampe et al. (2016) investigated the possibility of the nuclear origin of the double-peaked bursts based on the model of nuclear waiting points in the rp-process explaining the double-peaked profile (Fisker et al. 2004). They find that for certain metallicities and low accretion rate, there is an anti-correlation between the peak-flux ratio ( 1,2 ) and accretion rate. We do not find that the peak-flux ratio of double-peaked burst decreases when the source moves in the CCD (see Figure 2) from the top right to the bottom right, as the accretion rate gradually increases. If the temperature before the burst is propotional to mass accretion rate, the correlation between the peak flux and the accretion rate does not agree (see Figure 4) with the simulation of Lampe et al. (2016) . In their simulation, as accretion rate increases, the double-peaked structure shows a distinct two stage, different from the large dip in the bolometric flux profile of our sample (See Figure 1).
The newly discovered quadruple-peaked burst in 4U 1636−53 poses a problem to the thermonuclear spreading model of Bhattacharyya & Strohmayer (2006a,b). If this model is applied to ex- plain the four peaks, the burning front not only needs to stall three time, but it also needs a longer "waiting" period between the second and third peak than the separation time in the double-peaked burst. Similar to normal single-peaked bursts, we find that the value is large when the flux is high during the multi-peaked bursts. We investigated the relation between and the burst spectral parameters, but do not find any obvious correlation between them. The enhanced can be interpreted as changes in the accretion flow rate by Poynting-Robertson drag (Walker 1992). However, it is difficult to explain the dip between peaks by reducing the accretion.
The double-peaked structure could be due to the influence of the variation of the accretion geometry, but we should observe more double-peaked bursts to test this hypothesis. We note also that it is hard to explain the existence of triple or quadruple-peaked profiles with the variation of the geometry.

SUMMARY
We found 16 bursts with multi-peaked structure by investigating 336 type I X-ray bursts in LMXB 4U 1636−53 with RXTE. Our sample contains 14 double-peaked, one triple-peaked and one quadruplepeaked burst; the latter had never been reported before.
(i) Most of the multi-peaked bursts in our sample appear during the transition from the hard to the soft state in the CCD.
(ii) We find that double-peaked bursts with high peak-flux ratio appear when the pre-burst temperature is high; the high NS temperature may be due to enhanced accretion rate on to the NS surface.
(iii) We find an anti-correlation between the second peak flux and the peak-separation time for double-peaked bursts.
(iv) The quadruple-peaked burst shows a long separation time between the second and the third peak which is difficult to explain with current models.
(v) We use the method to re-analyse these 16 bursts and we find no evidence that the multi-peaked structure is due to enhanced accretion during the bursts.

APPENDIX A: DOUBLE-PEAKED EVENTS IN OUR SAMPLE
Some bursts appear to be flat after the onset followed by a clear peak, such as burst #6 and #8. In order to verify the significance of the double-peaked structure of these bursts, we use either one or two Gaussian functions to fit the bursts around the peaking time data.
For burst #6, at the top left panel of A1, we show the fitting results of the two Gaussian functions, yielding 2 = 113.3 for 17 d.o.f, and at the top right panel of A1, we show the fitting results with one Gaussian function, yielding 2 = 304.8 for 20 d.o.f. The F-test probability for these two fits is 0.0006, which indicates that the double-peaked structure is significant.
For burst #8, at the bottom left panel of A1, we show the fitting results of the two Gaussian functions, we yielding 2 = 59.8 for 21 d.o.f, and at the bottom right panel of A1, we show the fitting results with one Gaussian function, yielding 2 = 241.9 for 24 d.o.f. The F-test probability for these two fits is 1.4 × 10 −6 , which also indicates that the double-peaked structure is significant. This paper has been typeset from a T E X/L A T E X file prepared by the author.  Figure A1. The light curves of burst #6 and #8 of 4U 1636−53. The peaking time data of these two bursts can be well described by two Gaussian function.