Study of the X-ray properties of the neutron-star binary 4U 1728$-$34 from the soft to hard state

We studied five XMM-Newton observations of the neutron-star binary 4U 1728$-$34 covering the hard, intermediate and soft spectral states. By jointly fitting the spectra with several reflection models, we obtained an inclination angle of 25$-$53$\deg$ and an iron abundance up to 10 times the solar. From the fits with reflection models, we found that the fluxes of the reflection and the Comptonised components vary inconsistently; since the latter is assumed to be the illuminating source, this result possibly indicates the contribution of the neutron star surface/boundary layer to the disc reflection. As the source evolved from the relatively soft to the intermediate state, the disc inner radius decreased, opposite to the prediction of the standard accretion disc model. We also explore the possible reasons why the supersolar iron abundance is required by the data and found that this high value is probably caused by the absence of the hard photons in the XMM-Newton data.

A reflection spectrum, as the result of the hard coronal radiation illuminating an accretion disc, has been observed in several accreting black hole (BH; e.g. George & Fabian 1991;Magdziarz & Zdziarski 1995;Nowak, Wilms & Dove 2002;Miller et al. 2013) and neutron star (NS; e.g. Bhattacharyya & Strohmayer 2007;Cackett et al. 2008Cackett et al. , 2010Wang et al. 2017) systems. The combination of the high fluorescent yield and large cosmic abundance makes the iron emission line at 6.4-7 keV the most prominent feature in the reflection spectrum of these systems (see the Monte Carlo simulation results in Reynolds 1996). As the energy of some incident X-ray photons is much larger than the binding energy of the atomic electron in the disc, where those photons are scattered, Compton recoil occurs. This leads to a hump at high energies (e.g. George E-mail: yanan@astro.rug.nl & Fabian 1991;Matt, Perola & Piro 1991), known as the Compton hump, peaking at 30 keV in the reflection spectrum.
Unlike in BH and faint NS low-mass X-ray binaries (LMXBs), where the illuminating source of the disc is assumed to be a hot corona of highly energetic elements, in bright NS-LMXBs the NS boundary layer could contribute significantly to the reflection spectrum as well (Cackett et al. 2010;Miller et al. 2013;Ludlam et al. 2017). Regardless of the nature of the illuminating source, when reflection occurs in the vicinity of the compact object, the reflection spectrum can be modified by Doppler effects, light bending, and gravitational redshift; the combination of all these effects produce a broadened and skewed line profile with a red wing extending to low energies (e.g. Fabian et al. 2000;Reynolds & Nowak 2003;Miller, Turner & Reeves 2008). Therefore, by studying the asymmetrically broadened profile of such lines, we can investigate the geometry and the extension of the accretion disc.
Even though modelling the reflection spectrum has so far provided one of the best methods to estimate the spin parameter in BH The X-ray properties of 4U 1728-34 systems, the derived high iron abundance (several times the solar value; e.g. Cyg X-1, Parker et al. 2015;GX 339-4, Fürst et al. 2015;García et al. 2015) of the disc rises concerns about the accuracy of the spin estimates. Currently, there is no plausible physical explanation for these systems to be so iron rich. Fürst et al. (2015) found that the high iron abundance in GX 339-4 is model dependent. Once they allowed the photon indices of the direct power-law component and the power-law component that illuminates the disc to be different, the best-fitting iron abundance decreased and the fit statistically improved. Alternatively, Tomsick et al. (2018) reported that by applying high-density (up to 10 22 cm −3 ) reflection models, the fit no longer required a supersolar iron abundance in Cyg X-1.
The source states in atoll-type NSs are called the 'island' and 'banana' states, based on the shape of the colour-colour diagram (CD) and the timing properties of these sources, which correspond to the 'hard' and 'soft' states in other X-ray binaries, respectively. We used the latter nomenclature hereafter in this paper. The source states in these systems, which is likely related to changes in the mass accretion rate (Hasinger & van der Klis 1989), are usually associated with the evolution of the accretion flow. For instance, as a source evolves from the soft to the hard state, the edge of the disc moves outwards, from the innermost stable circular orbit (ISCO) to a larger radius (e.g. Esin, McClintock & Narayan 1997;Done, Gierliński & Kubota 2007). However, Sanna et al. (2014) found that the inner radius of the accretion disc was uncorrelated with the spectral state for the NS 4U 1636-53.
A broad iron emission line has been detected in the X-ray spectra of 4U 1728-34 with several instruments, e.g. BEPPOSAX (Di Salvo et al. 2000;Piraino, Santangelo & Kaaret 2000), XMM-Newton (Ng et al. 2010;Egron et al. 2011), ASTROSAT/LAXPC (Verdhan Chauhan et al. 2017), NUSTAR and SWIFT (Sleator et al. 2016;Mondal et al. 2017). Both Sleator et al. (2016) and Mondal et al. (2017) fitted the spectra of the SWIFT and NUSTAR data of this source using a reflection model. Sleator et al. (2016) found a disc inclination angle of 37 o , an iron abundance of the accretion disc of <1 times solar and an upper limit for the inner disc radius of ≤2 R g , where R g = GM/c 2 . Mondal et al. (2017) reported that the inclination angle in this system is 22 o -40 o , the disc iron abundance is 2-5 times solar and, as the source evolved from the soft to the hard state, the inner radius changed from 2.3 +2.1 −1.0 to 3.7 +2.2 −0.7 R ISCO , consistent with being constant.
In this paper, we conduct timing and spectral analysis of the NS LMXB 4U 1728-34 with XMM-Newton data and the (quasi) simultaneous RXTE data to study how the accretion flow changed while the source evolved from the soft-to-hard state. The paper is organized as follows: in Section 2 we describe the observations and the data reduction; our results of the spectral analysis are presented in Section 3; we discuss our results in Section 4, and we summarize our conclusions in Section 5.

O B S E RVAT I O N S A N D DATA R E D U C T I O N
The XMM-Newton observatory (Jansen et al. 2001) carries three high-throughput X-ray telescopes, each of them containing an European Photon Imaging Camera (EPIC, 0.1-12 keV). Two of these cameras are equipped with Metal Oxide Semi-conductor (MOS) CCDs (Turner et al. 2001) and one carries PN CCDs (Strüder et al. 2001). Reflection grating spectrometers (RGS, 0.35-2.5 keV;den Herder et al. 2001) are installed behind two of these telescopes. The five XMM-Newton observations of 4U 1728-34 used here were taken between 2011 August and October 7. We show the details of the observations in Table 1 and refer to them as Obs. 1-5 according to the observing time. We used data obtained with the EPIC-PN in Timing mode and with the RGS in Standard spectroscopy. To reduce and analyse the raw data we used version 16.1.0 of the XMM-Newton Scientific Analysis Software (SAS) package. Using the command epatplot, we found that the PN data were affected by pile up and we hence excluded the central region of the point-spread function source to mitigate this effect.
There were 14 type-I X-ray bursts in the PN light curves; we excluded these periods when we produced the PN spectra. We extracted all the PN background spectra from the outer columns of the central CCD (RAWX in 4-10) and found that the extracted background spectra are contaminated with the source (see also Ng et al. 2010;Hiemstra et al. 2011). We hence used the PN observation (ObsID 0085680601) of GX 339-4, which is on similar sky coordinates and column density along the line of sight, when this source was in the quiescent state, as a blank field to extract background spectra for all the five PN observations. We re-binned the PN spectra to have a minimum of 25 counts or to oversample the instrumental energetic resolution by a maximum factor of 3 in each bin. We fitted the PN spectra between 2.5 and 11 keV, avoiding the detector Si K-edge at 1.8 keV and the mirror Au M-edge at 2.3 keV (Egron et al. 2011).
We extracted the RGS data using the SAS tool rgsproc to produce calibrated event files, spectra and response matrices. The RGS data were grouped to provide a minimum of 25 counts per bin. We fitted the RGS spectra between 1 and 2 keV to constrain models in the soft band. We fitted the X-ray spectra using XSPEC (12.9.1a). To account for the interstellar absorption, in all fits we used the component TBABS with solar abundances from Wilms, Allen & McCray (2000) and cross-sections from Verner et al. (1996). Unless explicitly mentioned, we quote all errors at 1σ confidence level and at 95 per cent confidence for upper limits.
There were also 22 Rossi X-ray Timing Explorer (RXTE) observations of 4U 1728-34 (quasi) simultaneous with our XMM-Newton data. To search for the presence of QPOs, we first generated standard good-time interval files (GTIs) to remove instrumental drop-outs and other technical anomalies from the Proportional Counting Array (PCA) observations as suggested by the RXTE documentation. 1 Type I bursts have been detected and removed as well. We then divided each observation into segments of 16 s and extracted power spectra using the full energy band with a Nyquist frequency of 2048 Hz and averaged all the segments to obtain a single power spectrum for each observation.

Timing analysis
According to Zhang et al. (2016), some of the RXTE observations of 4U 1728-34 are contaminated by the nearby active transient 4U 1730-335 (the Rapid Burster). Both of the sources are in the PCA field of view and this transient was in outburst at the same time with the RXTE observations. Because the Rapid Burster only displayed significant power in the low-frequency range (Rutledge et al. 1995;Stella et al. 1988), we ignored the frequency range, <200 Hz, of the power spectra of 4U 1728-34 to avoid the contamination from the Rapid Burster. We linearly rebinned the power spectra by a factor of 200 to a frequency resolution of 12.5 Hz to improve the signal-tonoise ratio and fitted the power spectra with a constant to represent the Poisson noise and one or two Lorentzians to represent the kHz QPO(s).
We found significant kHz QPOs only in two observations: ObsIDs 96322-01-03-00 and 96322-01-03-01, both of them corresponding to Obs. 3 of the XMM-Newton data. The kHz QPO in observation 96322-01-03-00 has a frequency of 604 ± 17 Hz and a fractional rms amplitude of 7.3 ± 1.8 per cent, at a level of significance of 2.8σ , calculated as the ratio of the integral power of the fitted Lorentzian with 1σ negative error. Another kHz QPO in observation 96322-01-03-01 has a frequency of 583 ± 19 Hz and a fractional rms amplitude of 9.8 ± 1.6 per cent, at a level of significance of 4σ .

Colour-colour diagram and long-term light curve
To explore the source state of 4U 1728-34, we took the data from Zhang et al. (2016) and plotted the CD of the RXTE data in the upper panel of Fig. 1. As the definition of the colours in their work, the soft and hard colours are the 3.5-6.0/2.0-3.5 keV and 9.7-16.0/6.0-9.7 keV count rate ratios, respectively. Type I bursts have been removed from the RXTE data and the colours of 4U 1728-34 are normalized to the colours of Crab. Zhang et al. (2016) parametrized the position of the source on the CD through the value of the parameter S a , that gives quantitatively the position of the source along the path traced by the source in the CD . They fixed the values of S a = 1 and S a = 2 at the top right and the bottom left vertex of the CD, respectively.
We assigned an S a value to each XMM-Newton observation as the average S a value of the simultaneous RXTE data and indicated them with the red, green, blue, magenta, and olive squares in the upper panel of Fig. 1. During our observations, the source evolved from  Table 1 from top to bottom correspond to the simultaneous RXTE/SWIFT data in red, green, blue, magenta, and olive squares/lines. the left bottom to the right top on the CD as S a decreased. As some of the RXTE observations are contaminated by the Rapid Burster, this prevented us from using the simultaneous RXTE data to do spectral analysis and the emission from the Rapid Burster may have also affected the colours of these observations. For instance, Obs. 2 and 3 are off the main track in the CD and both of them are entirely contaminated; the other RXTE observations are only partially contaminated. To check whether the evolution of the source in the RXTE CD is reliable, we created a SWIFT/BAT long-term light curve in the energy of 15-50 keV at around the time of the observations with XMM-Newton of 4U 1728-34; we show the SWIFT/BAT light curve in the lower panel of Fig. 1. Corresponding to the five XMM-Newton observations, the count rate of the SWIFT/BAT light curve increased from Obs. 1-3, remained constant within errors during Obs. 3 and 4, and increased again from Obs. 4 to 5.
Both the source evolution on the CD in the upper panel of Fig. 1 and of the light curve in the lower panel of Fig. 1 indicate that the source indeed transited from the relatively soft to the hard state.

Spectral analysis
We initially fitted the PN spectra of the five XMM-Newton observations in the energy range 2.5-11 keV with a thermally Comptonized component (NTHCOMP in XSPEC; Zdziarski, Johnson & Magdziarz 1996;Życki, Done & Smith 1999) plus a single temperature blackbody component (BBODYRAD in XSPEC). The fit was bad, χ 2 = 2199.9 for ν = 652, where ν is the number of degrees of freedom, and the fit showed prominent residuals at 5-8.5 keV.
We then fitted the spectra with the same components, but only in the energy ranges of 2.5-5 and 8.5-11 keV; we show the datato-model ratio of individual observation in Fig. 2. A strong broad asymmetric emission feature appears at around 5-9 keV in each spectrum in this plot. During this fit, we also found that: (1) the seed photon temperature of the NTHCOMP component, kT seed , in all the spectra, except for Obs. 1, is consistent with 0; we therefore linked this parameter across the spectra of Obs. 2-4 and got an upper limit at 0.4 keV. However, in order to use a value that was consistent with the one used in the models that we applied in the following sections, we fixed this parameter at kT seed = 0.05 keV. This improved the constraints on other parameters without extra effect on the fit. The electron temperature of the NTHCOMP component, kT e , pegged at its upper limit, 1000 keV, in the spectra of Obs. 1-4 and we thus fixed kT e at 300 keV in these observations to be consistent with the value that is required by the other models (see details in the following sections). Both the seed photon temperature in Obs. 1 and the electron temperature in Obs. 5 of the NTHCOMP component were free to vary. If we change NTHCOMP to CUTOFFPL to describe the hard part of the spectrum, we obtain a worse fit in this case; the χ 2 increased from χ 2 = 469.2 with ν = 374 to χ 2 = 516 for ν = 375. Using this continuum model, we fitted the broad emission feature with a simple Gaussian component. In this case, we fitted the data over 2.5-11 keV the full energy range. We call this model M1 gau; the fit yields χ 2 = 770.1 for ν = 638 and with a null hypothesis probability of 2.4 × 10 −4 . The seed temperature kT seed , in NTHCOMP is 0.70 ± 0.05 keV in Obs. 1 and the electron temperature kT e , in NTHCOMP, is 3.4 ± 0.1 keV in Obs. 5. As an example, we show the individual components and model residuals in terms of sigmas for Obs. 1 and 5 in the upper panels Fig. 3, since Obs. 1 and 5 represent, respectively, the softest and hardest spectra of the source in our samples.
We report the best-fitting parameters and the individual flux of each component of M1 gau in Table 2 and show the evolution of the parameters and flux of each individual component as a function of S a in Figs 4 and 5, respectively. The photon index , in NTHCOMP, and the blackbody temperature kT bb , increased monotonically with S a , consistent with the softening of the spectrum as the source evolved in the CD. The centroid energy of the Gaussian component decreased first and then increased, while the width of the line changed in the opposite way. The fluxes of the NTHCOMP, F Compt , and the Gaussian, F gau , components change in correlation with each other except in Obs. 2, indicating that the corona was probably the main illuminating source of the reflection component, here represented   Note. In this and the following tables, the symbol l indicates that the parameters are linked across the observations, f means that the parameter is fixed during the fit, p denotes that the parameter pegged at its limit, and u stands for 95% confidence upper limit. All the unabsorbed fluxes are in units of 10 −10 erg cm −2 s −1 in the 2.5-11 keV range. Errors are quoted at the 1σ confidence level.
by the GAUSSIAN line. Even though the flux of the soft component, F bb , in Obs. 1 in the 2.5-11 keV energy band is almost four times higher than that in other observations, the hard NTHCOMP component dominates the emission during the entire evolution; the total flux, F tt , peaks in Obs. 1 and does not appear to change in a simple manner with the source state.
In order to test if adding the RGS data to the fits can help constraining the column density, we fitted the RGS spectra in the energy range between 1 and 2 keV, together with PN spectra in the energy range 2.5-11 keV. For the two RGS and one PN spectra of the same observation, we tied all parameters to each other, with two multiplicative factors, one for each RGS instrument, left free to vary; for the same instrument, among different observations, this multiplicative factors were linked. The best-fitting value of the column density, N H = 4.1 ± 0.03 × 10 22 cm −2 , obtained for the joint fit of the RGS and PN spectra, is similar with the value, N H = 4.5 ± 0.1 × 10 22 cm −2 , that we derived from the fit to the PN spectra only. In the end, we found that adding the RGS data did not improve the value of N H significantly, and therefore we continued using the PN spectra only.
In previous studies, the best-fitting hydrogen column density along the line of sight was 2.4 − 2.6 × 10 22 cm −2 (e.g. D' Aí et al. 2006;Egron et al. 2011). However, since in these papers they used different cross-sections and solar-abundance tables from ours, it is no surprising that the column density in our case is not consistent with theirs. Mondal et al. (2017) and Sleator et al. (2016) analysed NUSTAR and SWIFT data of 4U 1728-34 using the same cross-sections and solar-abundance tables as ours to calculate the column density, and found the column density as N H ∼ 3.9 − 4.5 × 10 22 cm −2 .

Relativistic reflection model
Since the plots in Fig. 2 suggest that the broad profile of the emission line at 7 keV is not symmetric and it has been argued in the past that this may be due to Doppler and relativistic effects, we fitted the spectra with the self-consistent reflection model RELXILLCP 2 v0.5b 2 http://www.sternwarte.uni-erlangen.de/∼dauser/research/relxill/ (Dauser et al. 2014;García et al. 2014). This component includes the thermal Comptonization model NTHCOMP as the illuminating continuum. To limit the number of the free parameters, we set the inner and outer emissivity indices of this component to be the same time, we set both of them to be the same, q in = q out , and let q in free to vary across observations. Following Braje, Romani & Rauch (2000) and assuming a 1.4 M NS, we adopted a dimensionless spin parameter a * = 0.47/P ms , where P ms is the spin period in ms. Since the spin frequency of 4U 1728-34 is 363 Hz (Strohmayer et al. 1996), we fixed a * = 0.17.
Our result of the fit of the model M1 gau in the previous section showed that the best-fitting values of the electron temperature were much larger than the upper bound of the PN energy range in all observations except for Obs. 5. To make a fair comparison between models, as the high-energy rollover is fixed at 300 keV by default in RELXILLD (and also in REFLIONX-based models, which we will apply in the following sections), we fixed the electron temperature in RELXILLCP at 300 keV in Obs. 1-4. Once we got a good fit with the model TBABS * (BBODYRAD+RELXILLCP), we froze the parameter reflection fraction, refl frac, at its negative value to force this component to only account for the reflected part, we added the direct NTHCOMP component to the model and linked the common parameters, the photon index and the electron temperature, of both the components NTHCOMP and RELXILLCP. This procedure allows us to get the fluxes of the individual components separately. We followed the same procedure when we used other RELXILL-based models in this paper.
The overall model, hereafter M2 Cp, became TBABS * (BBODYRAD+RELXILLCP+NTHCOMP), which yields χ 2 /ν = 685.4/631. Compared with the fit with M1 gau, the χ 2 of this fit decreased by χ 2 = 84.7 for 7 degrees of freedom less. The emissivity index was not well constrained and was marginally consistent within errors in all the observations. We therefore linked this parameter across observations to improve the constraint on other parameters, which yields χ 2 /ν = 694.6/635 and with a null hypothesis probability of 0.05 (see the unfolded spectra and models in Fig. 3).
We show the relevant parameters of this model in Table 3 and plot some of the parameters of each component as a function From the top to the bottom panels, the parameters are the blackbody temperature (keV) and the normalization (R 2 km /D 2 10 , where R km is the source radius in km and D 10 is the distance to the source in units of 10 kpc), the line energy (keV)/the disc ionization (erg cm s −1 ), the line width (keV)/the disc inner radius (R ISCO ), the Gaussian normalization (10 −3 )/the RELXILLCP/RELXILLD normalization (10 −3 ), the photon index and the NTHCOMP normalization, respectively. The green arrow indicates the 95 per cent confidence upper limit of the RELXILLCP normalization of Obs. 2. of S a in Fig. 4. In Obs. 1-3, the spectrum is dominated by the reflection component, whereas in Obs. 4 and 5 the fluxes of the reflection and the Comptonized components are comparable. There are, however, two issues with this fit: (1) the best-fitting value of the iron abundance is very high, A Fe = 10 times solar abundance, which is the upper bound of this parameter (see the contour plot for the iron abundance versus the inclination in Fig. 6); (2) some of the best-fitting parameters in this model are not consistent with the same parameters in M1 gau, e.g. both the blackbody temperature and the photon index in M1 gau monotonously increase with S a whereas the same parameters in M2 Cp first increase and then decrease or remain more or less constant. If we forced the iron abundance to be 1, the fit becomes worse and χ 2 increases by χ 2 / ν = 100.9/1. In all the Cp-type versions of the RELXILL-based models, the seed temperature is fixed at 0.05 keV by default, which is more than 10 times smaller than the best-fitting value of kT seed that we obtained from M1 gau in Obs. 1. This discrepancy may partly cause the inconsistent results between M1 gau and M2 Cp.
We also tried to fit the data with other types of the RELXILL-based reflection models: the fit with the model RELXILL was almost as good as the one with M2 Cp, χ 2 /ν = 697.1/635; as for the lamppost model, RELXILLLPCP, the fit yielded a similar χ 2 , χ 2 /ν = 696.2/635. However, as the Compton hump is not covered by XMM-Newton, we cannot constrain the key parameter, the height of the corona, in this model well. It is worth noting that the iron abundance, A Fe = 10 times solar abundance, pegs at its upper limit in the fits with all these reflection models.

Reflection model with high density
A high iron abundance using these reflection models has been reported in previous works and, in some cases, the authors argued that this was the result of the low density of the accretion disc, n e = 10 15 cm −3 used in the calculation of the models (e.g. García et al. 2016;Tomsick et al. 2018). To test the possible effect of the density on the disc iron abundance, we fitted the spectra with the extended reflection model RELXILLD (García et al. 2016) that allows the electron density parameter to vary between n e = 10 15 and 10 19 cm −3 , which we call M2 hd. We replaced RELXILLCP by RELXILLD, NTHCOMP by CUTOFFPL with the cut-off energy, E cut , fixed at 300 keV since this is required by the RELXILLD component, and applied the same fixed parameters as for M2 Cp. The electron density, log (n e ), was linked to be the same in Obs. 1-5.
The best-fitting parameters and individual unabsorbed flux are given in Table 4 and are added as red triangles to Figs 4 and 5. Compared to the fit with M2 Cp, the new fit slightly improved,  Note. All the symbols and units are the same as in Table 2. &x266E x0266E; The parameter, refl frac, has been frozen at its negative value to force the component RELXILLCP to only account for the reflection part. such that χ 2 decreased by χ 2 = 7.1 for the same ν and with a null hypothesis probability of 0.073. The iron abundance still pegged at 10, with a high density of log (n e ) up to 19; the evolution of the photon index in M2 hd is similar to that in M1 gau. The righthand panel in Fig. 6 shows the contour plot for the iron abundance versus the inclination for M2 hd. If we force the iron abundance to be 1, the fit becomes worse, with χ 2 = 131.9 for ν = 1. Similar to M2 Cp, the flux of M2 hd in Obs. 1-3 is dominated by the reflection component, whereas in Obs. 4 and 5 it is dominated by the Comptonized component.

Alternative reflection model
To check the robustness of the values we obtained from the fits with the RELXILL-based models, and especially to explore the issue of the supersolar iron abundance, we also fitted the data with the model REFLIONX (Ross & Fabian 2005) that characterizes the emergent reflection spectrum arising from an illuminating powerlaw spectrum, including a high-energy exponential cutoff with e-folding energy fixed at 300 keV; we convolved this component with the relativistic convolution model KERRCONV (Brenneman & Reynolds 2006). The model that we fitted in XSPEC was TBABS * (BBODYRAD+KERRCONV * REFLIONX+CUTOFFPL), which we call M3 pl. As in the previous fits, in REFLIONX we tied the inner and outer emissivity indices, q in = q out , and set the spin parameter a * = 0.17. The value of q in in M3 pl was consistent with being the same in all observations, so we linked this parameter across all the observations. Additionally, we linked the photon index and the cut-off energy in REFLIONX to those in CUTOFFPL. The final fit is worse than M2 Cp, χ 2 /ν = 727.3/637 and with a null hypothesis probability of 6.8 × 10 −3 ; the parameters are listed in Table A1 in the appendix. The best-fitting iron abundance is 5.8 +0.7 −0.04 and the inclination angle of the accretion disc with respect to the line of sight is 24.6 ± 1.2. The average blackbody temperature for M3 pl is smaller than for M2 Cp, but the trends of the photon index and the inner radius are similar to those for M2 Cp. The reflection flux for M3 pl in Obs. 1-3 is larger than in the rest of the observations; the reflection and Comptonized fluxes in Obs. 4 are almost equal, and the Comptonized flux in Obs. 5 is dominant.
We also applied the high electron density version of this model, REFLIONX HD (M3 hd; Tomsick et al. 2018), in which the density in the reflector can go up to 10 22 cm −3 ; the iron abundance is fixed at the solar abundance. The overall fit is worse than M3 pl, χ 2 /ν = 763.5/637 and with a null hypothesis probability of 3.7 × 10 −4 ; the corresponding best-fitting parameters are shown in Table A2. Compared to M2 hd, the evolution of the inner radius of both models is analogous. The spectrum fitted with M3 hd is dominated by the reflection component all the time.
Since 4U 1728-34 is a NS, we tried to fit the reflection spectrum with another version of REFLIONX, RE-FLIONX BB (Ludlam et al. 2017), in which the illuminating source is the blackbody component. In this model,  Note. All the symbols and units are the same as in Table 2. &x266E x0266E; The parameter, refl frac, has been frozen at its negative value to force the component RELXILLD to only account for the reflection part.
TBABS * (BBODYRAD+KERRCONV * REFLIONX BB+CUTOFFPL), which we call M3 bb, we linked the blackbody temperature, kT in REFLIONX BB, to the blackbody temperature, kT bb , in BBODYRAD in all observations. The fit is worse than M3 pl, χ 2 /ν = 753.3/637 and with a null hypothesis probability of 3.9 × 10 −5 (see Table A3), suggesting that the illuminating source in 4U 1728-34 cannot only be the blackbody component. However, different from the above results, the iron abundance derived from M3 bb is 0.78 +0.01 −0.09 and the inclination angle is 52.9 +1.6 −0.5 . Two other differences are that the column density and the overall blackbody temperature in M3 bb are higher than those in M3 pl. The spectrum for this model is dominated by the Comptonized component in all observations.
To further identify whether the illuminating source is the corona or the NS surface/the boundary, we combined the Comptonized REFLIONX and the blackbody REFLIONX BB versions together in a model, M3 pl bb. We assumed that the iron abundance and the ionization of the disc in both components are the same. In Obs. 5 the normalization of the REFLIONX BB component is negligible, as well as the normalization of the REFLIONX component in Obs. 2. We show the parameters in Table A4 and the fit yields χ 2 /ν = 730.2/634 and with a null hypothesis probability of 3.2 × 10 −3 , similar with M3 pl, although the inner radius in Obs. 1 is very large, R in = 57 R ISCO , and the iron abundance is consistent with the one in M3 pl. The flux for Obs. 1 and 3 is dominated by the REFLIONX component; the flux for Obs. 2 and 5 is dominated by the Comptonized component. In Obs. 4, the fluxes of the REFLIONX and the Comptonized components are almost the same. Except in Obs. 2, the REFLIONX flux is always larger than that of the REFLIONX BB component.

Tests with NUSTAR data
As previously mentioned in Section 1, 4U 1728-34 was also observed with NUSTAR. Mondal et al. (2017) analysed two NUSTAR observations (ObsIDs: 80001012002 and 80001012004) plus two simultaneous SWIFT/XRT observations (ObsIDs: 00080185001 and 00080185002) and found an iron abundance A Fe = 2-5 times solar. To test if the discrepancy in the iron abundance derived from our and their models is due to the lack of coverage of the 11 keVhigh-energy range (above ) in our data, we re-analysed NUSTAR observation 80001012002 and the simultaneous XRT observation 00080185001 in which the source was in the hard state. We used M2 Cp to jointly fit the NUSTAR observation in the energy ranges 3.5-50.0/3.5-11.0 keV and the XRT observation in the energy range 1.0-7.5 keV. Type I bursts were detected and removed from the NUSTAR spectra. Although M2 Cp and M2 hd fit the XMM-Newton data equally well, a cut-off energy is required by the NUSTAR spectra (Mondal et al. 2017) and the cut-off energy is frozen at 300 keV as default in RELXILLD, therefore here we chose M2 Cp to do this test.
In Table 5, we show the best-fitting results when the photons above 11 keV are either included or excluded in the NUSTAR spectra. The results show that most of the parameters are marginally consistent with each other no matter whether the hard photons are included in the spectra or not; as expected, the parameters that are affected the most are the photon index, , and the electron temperature, kT e , from the NTHCOMP component. On the other hand, both the reflection and Comptonized components are less well constrained when we exclude the hard photons. Another significant difference is that the iron abundance, A Fe , increases from two times solar when we include the NUSTAR data above 11 keV, to eight times solar when we fit the NUSTAR spectra only in the 3.5-11 keV range.

D I S C U S S I O N
We analysed five XMM-Newton observations of the NS LMXB source 4U 1728-34 obtained in 2011, when the source evolved from Note. The energy range of the SWIFT/XRT data used here is always between 1.0 and 7.5 keV; only the energy range of the NUSTAR data changes. All the symbols and units are the same as in Table 2. the soft-to-hard state, to explore how the accretion flow changed between those states. A broad emission line at 6.5-6.7 keV in the spectrum of this source indicates the presence of a reflection component in this system. By jointly fitting all the five spectra with several reflection models, we obtained an inclination angle of 25 o -53 o and an iron abundance of up to 10 times the solar abundance. In what follows, we compare the spectral parameters derived from the fits with different models, identify the kHz QPOs that we observed in the power spectra, and discuss the possible reasons why a supersolar iron abundance appears to be required by the data.

Comparisons of all applied the models
In this paper, we fitted the continuum spectrum of 4U 1728-34 with a single temperature blackbody BBODYRAD, plus a Comptonized component, NTHCOMP/CUTOFFPL depending on the requirement of the model, to account for the soft and hard photons in the spectra, respectively. As shown in Fig. 2, a strong emission feature appears to be present in the 5-9 keV energy range of each spectrum. We used several components to fit this emission: a Gaussian component in M1 gau and the reflection components RELXILLCP in M2 Cp, RELXILLD in M2 hd, KERRCONV * REFLIONX in M3 pl, KERRCONV * REFLIONX HD in M3 hd, KERRCONV * REFLIONX BB in M3 bb, KERRCONV * (REFLIONX+REFLIONX BB) in M3 pl bb, respectively. In M1 gau, as shown in Table 2 and Fig. 5, the line flux is the same in Obs. 1 and 5, which represent, respectively, the softest and the hardest state observations in this work; even though the spectrum was dominated by the hard component all the time, as the total flux decreased, the blackbody flux in Obs. 1 dramatically dropped to one fourth of that in Obs. 5.
When the emission feature was fitted with the reflection component RELXILLCP in M2 Cp, the spectrum was dominated by the reflection component in Obs. 1-3 and the fluxes of the reflection and the Comptonized components were equally strong in the last two observations. When we took the high-density effect (> 15 cm −3 ) into account in the reflection component, the fit with M2 hd was slightly better than that of M2 Cp. Similar to M2 Cp in Obs. 1-3, the dominant component in M2 hd are the reflection component but in Obs. 4 and 5, the dominant component in M2 hd turns to be the Comptonized component (see Tables 3 and 4).
Ever though the trends and values of the parameters derived from M2 Cp and M2 hd are consistent within errors in Obs. 2-4, these parameters in Obs. 1 and 5 are different. For instance, the blackbody flux, in units of 10 −10 erg cm −2 s −1 , of Obs. 1 increased from 1.1 ± 0.3 in M2 Cp to 4.1 ± 0.3 in M2 hd and that of Obs. 5 decreased from 2.6 ± 0.4 in M2 Cp to 1.3 ± 0.1 in M2 hd. On the contrary, the Comptonized flux, in units of 10 −10 erg cm −2 s −1 , of Obs. 1 decreased from 7.3 ± 0.9 in M2 Cp to 2.1 +7.5 −0 p in M2 hd and that of Obs. 5 increased from 6.2 ± 0.9 in M2 Cp to 10.0 ± 0.5 in M2 hd.
The iron abundance derived from M2 Cp and M2 hd both pegs at the upper limit, A Fe = 10 in solar units. When we replaced the self-consistent reflection models RELXILLCP and RELXILLD with the REFLIONX-based components convolved with the relativistic blurring kernel KERRCONV in M3 pl, M3 hd, and M3 bb, the fit became worse, with χ 2 increasing 40-76 for 2 degrees of freedom more (see Tables A1-A3). The iron abundance in M3 pl and M3 bb were 5.8 +0.04 −0.8 and 0.78 +0.01 −0.09 times solar, respectively. When A Fe was forced to be 1 and the density of the disc was allowed to be as high as 10 22 cm −3 in M3 hd, there was not improvement on the fit.
The inclination derived from all the models above was around 30 o , except in M3 bb in which the inclination was 52.9 +1.6 −0.5 but the χ 2 of the fit was very large. Although the version of the combination of the REFLIONX and REFLIONX BB components, M3 pl bb, improves the fit compared to the version of the REFLIONX BB component alone, this fit does not yield an iron abundance as low as in M3 bb.

Identification of the kHz QPO
As we mentioned in Section 3.1, two single kHz QPOs with a frequency of, respectively, 604 ± 17 Hz and 583 ± 19 Hz have been detected in the RXTE data, corresponding to our Obs. 3 with XMM-Newton. As reported by , Di Salvo et al. (2001) andvan Straaten et al. (2002), the frequencies of the upper and lower kHz QPOs in 4U 1728-34 fall in the range 500-880 Hz and 380-1160 Hz, respectively. Di Salvo et al. (2001) and van Straaten et al. (2002) studied the fractional rms amplitude of both kHz QPOs as a function of the QPO frequency. In order to tell whether we have detected the lower or the upper kHz QPOs, we compared the rms amplitude and frequency of our kHz QPOs to the ones in their papers, and found that both our detected QPOs are more likely the upper kHz QPO.  found that, as a function of S a , the frequencies of the upper and lower QPOs follow well-defined separate tracks. If both kHz QPOs that we detected here were the lower kHz QPO, according to Fig. 4 in , the corresponding S a would be 1.9, which is much higher than the one of Obs. 3, S a = 1.3 and, different from what we observe, it would put the source in the transitional intermediate state, close to the vertex of the CD. If the QPOs that we detected were the upper kHz QPO, the corresponding S a from the same figure would be between 1.3 and 1.4, consistent with our value of S a for Obs. 3. We therefore conclude that the two detections of kHz QPOs in the RXTE data that are simultaneous with our XMM-Newton Obs. 3 correspond to the upper kHz QPOs in 4U 1728-34.

Inner radius uncorrelated with source states
The evolution of the source on the RXTE CD and in the SWIFT/BAT light curve in Fig. 1 give an idea of the spectral evolution of 4U 1728-34 during the XMM-Newton observations presented here, from a relatively soft to the hard state. The evolution of the spectral parameters of M1 gau support this idea: the blackbody temperature and the photon index of NTHCOMP increase as S a increases, even though the spectra are dominated by the hard component, NTHCOMP, at all times (see Figs 4 and 5). The flux of the Gaussian component followed the same trend as that of the NTHCOMP component, except in Obs. 2.
When we fitted the data with the relativistic reflection models M2 Cp and M2 hd, the model parameters follow a similar trend to that of model M1 gau, except the RELXILLD normalization in model M2 hd. In the standard truncated accretion disc model, as mass accretion rate increases the inner disc radius moves inwards (Esin et al. 1997;Done et al. 2007). However, Fig. 4 shows that the inner radius derived from both models first decreased from Obs. 1 and 2, it then remains constant in Obs. 2-4 and it finally increased from Obs. 4 and 5. The evolution of the inner radius in Obs. 2-5 supports the truncated disc model above, indicating that the inner radius moves outwards with decreasing mass accretion disc. However, going from Obs. 1 and 2, with an apparently decreasing mass accretion rate, the inner radius moves inwards.
In the standard accretion disc model, gas pressure dominates when both the accretion rate and X-ray luminosity (L x < 10 36 ergs s −1 ) are low (Shakura & Sunyaev 1973). On the contrary, when the luminosity is high, radiation pressure should dominate. Popham & Sunyaev (2001) showed that when the luminosity approaches the Eddington limit, the radiation feedback from the NS surface leads to an increase of the inner radius. As the flux of Obs. 1 is the largest one in our samples, this process may result in the inner radius variation that we observe.

Iron abundance deduced from XMM-Newton and NUSTAR data
As we showed in Fig. 6, the best-fitting value of the iron abundance in 4U 1728-34 from the fits to the XMM-Newton data is 10 times solar or higher, which differs from what Mondal et al. (2017) found with NUSTAR and SWIFT data. Mondal et al. (2017) analysed two simultaneous NUSTAR and SWIFT observations carried out in 2013, and inferred that during these two observations the source was in the hard and soft state, respectively. Similar to what they did, we assumed that the spin parameter is 0.17, and applied similar models to fit the reflection spectrum: they used RELXILL and we used RELXILLCP; the inclination angle in their and our work are consistent, around 30 o , but the iron abundance they obtained is 2-4 times the solar, about half to one-fifth of the value that we find.
A high electron density of the accretion disc has been suggested as a potential solution of the supersolar disc iron abundance (e.g. García et al. 2016;Tomsick et al. 2018). Tomsick et al. (2018) explained that a high density produces more soft emission, resulting in a harder power law, which provides a better match to the hard spectrum, as well as an extra soft excess below 1 keV. However, compared to the fit with model M2 Cp, the fit with model M2 hd that allows for higher density than M2 Cp, only improved slightly, with the iron abundance pegging at 10 times solar and the density pegging at 10 19 cm −3 . Allowing for a higher density of the disc only increased the column density of the interstellar medium and the disc temperature in our fits.
As the iron abundance derived from model M2 hd pegged at its upper limit, we tested another model, REFLIONX HD, with an electron density that can go up to 10 22 cm −3 . Unfortunately, this model, M3 hd, did not return a good fit and the density still pegged at the upper limit (see Table A2).
Since the iron abundance reported by Mondal et al. (2017) is very different from ours, we did another test with NUSTAR and SWIFT data, as Mondal et al. (2017) used, to see if the lack of the data at energies above 11 keV plays a role in this result. As we described in Section 3.3, the iron abundance, A Fe , increases from 2 in solar units when we fit the full NUSTAR data up to 50 keV, to 8 in solar units when we fit the NUSTAR spectra only in the 3.5-11 keV range. At the same time, the NTHCOMP becomes negligible if we ignore the NUSTAR data above 11 keV. A possible explanation for this result is that in order to produce a similar significant reflection spectrum, more iron is required. Even though neither the RELXILLCP nor the NTHCOMP components are well constrained when the hard photons are ignored, we cannot exclude this hypothesis.

The possible illuminating source of the reflection component
Most of the fits show that reflection makes a significant contribution to the entire spectrum. For instance, the reflection component in Obs. 1-3 dominated the total emission in all the models, except in M3 bb. The reflection fraction, refl frac, remains constant within errors among the observations, and the reflection flux is independent of the Comptonized flux in the fits with models M2 Cp and M2 hd, in both of which we assume that the corona is responsible for the disc reflection. We identify two possible explanations for our finding that the changes of the reflection flux and the Comptonized flux are uncorrelated. The first possibility is that both the NS surface/boundary and corona irradiated the disc and contributed to the reflection spectrum. Alternatively, light bending may play a role in the reflection process as well. As the reflection happens in the vicinity of a compact object, due to the strong gravitational light-bending effect, more of the Comptonized photons would be bent towards the disc, which results in less of the Comptonized photons being observed directly at infinity. Miniutti & Fabian (2004) identified three different regimes in which the reflection-dominated component (and the iron line) is correlated, anti-correlated or almost independent with respect to the direct continuum, and they concluded that the relation between the reflection and direct component is correlated to the source state and the height of the illuminating source. Cackett et al. (2010) studied broad iron emission lines in 10 NS-LMXBs and concluded that the boundary layer is the illuminating source irradiating the accretion disc in these systems. In Section 3.2.3, we explored the relative contribution of the corona and the NS surface/boundary layer to the reflection spectrum. Comparing the fits with models M3 pl and M3 bb, the former gives a better fit, χ 2 = 25.9 with the same ν, which suggests that the boundary layer might not be the only contributor to the reflection spectrum in all observations. Thanks to model M3 pl bb, we can make a direct comparison of the contribution to the reflection spectrum between the corona and the NS surface/boundary layer. In Table A4, we show that the flux of the REFLIONX component is much larger than that of the REFLIONX BB component except in Obs. 2. The boundary layer contributed 4 per cent-43 per cent of the total flux to the reflection component in Obs. 1-4, not strong but still required by the data; the contribution of the corona to the reflection component is considerable, 25-63 per cent of the total flux in Obs. 1 and 3-5. This suggests that most of the time the disc is mainly illuminated by the corona, and the contribution of the illuminating source is not affected by the source state. It is worthwhile to emphasize that neither the changes of CUTOFFPL and the REFLIONX fluxes nor these of the BBODYRAD and the REFLIONX BB fluxes are correlated.

Some caveats
Note that even though compared to other models, models M2 Cp and M2 hd statistically give the best fits, low χ 2 and null hypothesis probabilities, the iron abundance derived from these two models pegged at the upper allowed limit. If we forced the iron abundance to be 1, the fits with the reflection models, RELXILLCP and RELXILLD, are worse than the fit with Gaussian. This fact may affect the other best-fitting parameters derived from both reflection models. However, the relative evolution of these parameters should be still reliable.

C O N C L U S I O N S
The NS LMXBs 4U 1728-34 has been jointly observed by XMM-Newton and RXTE in 2011. We carried out the spectral and timing analysis with both instruments, and found that the source evolved from the soft-to-hard state during a period of 40 d. We fitted the PN spectra with several reflection models; the fits yield a disc inclination angle of 25 o -53 o and an iron abundance as high as 10 times solar, which is probably the result of the lack of high-energy coverage of the XMM-Newton instruments. Besides that, when the source evolved from the soft to intermediate state, we found that the changes in the inner radius of the accretion disc do not support the standard accretion model. We finally concluded that during the entire evolution, both the corona and NS surface/boundary layer contributed to the reflection component, but the former was dominant most of the time.

AC K N OW L E D G E M E N T S
This work is partly supported by China Scholarship Council (CSC), under the grant number 201404910530. TMB acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0. DA acknowledges support from the Royal Society. EMR acknowledges the support from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-Brazil). GZ acknowledges funding support from the CAS Pioneer Hundred Talent Program Y7CZ181001. This work has made use of data from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA/Goddard Space Flight Center (GSFC).  Note. In this and the following tables, the symbol l indicates that the parameters are linked to vary across the observations, f means that the parameter is fixed during the fit, p denotes that the parameter pegs at its limit and u stands for 95% confidence upper limit. All the fluxes are in units of 10 −10 erg cm −2 s −1 in the 2.5-11 keV range. Errors are quoted at 1σ confidence level. Note. All the symbols and units are the same as in Table A1. Note. All the symbols and units are the same as in Table A1. Note. All the symbols and units are the same as in Table A1.

A P P E N D I X A : A D D I T I O NA L B E S T-F I T T I N G PA R A M E T E R S F O R T H E 4 U 1 7 2 8 -3 4
This paper has been typeset from a T E X/L A T E X file prepared by the author.