Broadband X-Ray Timing and Spectral Characteristics of the Accretion-powered Millisecond X-Ray Pulsar MAXIJ1816-195

We studied the broadband X-ray timing and spectral behaviors of the newly confirmed accreting millisecond X-ray pulsar MAXI J1816−195 during its 2022 outburst. We used data from the Insight-HXMT Medium Energy (ME) and High Energy (HE) telescopes, NICER, and NuSTAR that cover the energy range between 0.8 and 210 keV. A coherent timing analysis of solely the Insight-HXMT HE data across the full outburst revealed a complex behavior of the timing residuals, also prominently visible in the independent Insight-HXMT ME and NICER data, particularly at the rising part of the outburst and at the very end in the NICER data. Therefore, we broke down the full outburst into a (noisy) rising part, covering only about five days, from MJD 59737.0 to 59741.9, and a decaying part, lasting for 19 days across MJD 59741.9–59760.6. Fitting for the decaying part, a timing model including a frequency ν and frequency time derivative ν̇ component yielded a value of (+9.0 ± 2.1) × 10−14 Hz s−1 for ν̇ , which could be interpreted as a spinup under our model assumptions. We detected X-ray pulsations up to ∼95 keV in a combination of Insight-HXMT HE observations. The pulse profiles were quite stable over the whole outburst and could be well described by a truncated Fourier series using two harmonics, the fundamental and the first overtone. Both components kept alignment in the range 0.8–64 keV. The joint and time-averaged NICER and Insight-HXMT spectra in the energy range 1–150 keV were well fitted by the absorbed Comptonization model compps plus disk blackbody with two additional Gaussian components. Using the bolometric flux and spinup values both evaluated during the decay phase, we determined a magnetic field strength of (0.2–2) × 108 G for MAXI J1816−195.

The new X-ray transient MAXI J1816−195 was discovered by MAXI/Gas Slit Camera on June 7 during its 2022 outburst (Negoro et al. 2022).Using the Swift localization, NICER detected pulsations and thermonuclear type I X-ray bursts from MAXI J1816−195, confirming that the source is an AMXP with a spin frequency of 528 Hz, an orbital period of 4.83 hr and a projected semi-major axis of 0.26 lt-s (Kennea et al. 2022;Bult et al. 2022b,c,a).The mass of the companion star was determined to be in the range 0.10-0.55M ⊙ (Bult et al. 2022a).Meanwhile, Insight-HXMT found X-ray pulsations from MAXI J1816−195 in the hard Xray/soft γ-ray band (Li et al. 2022).Bult et al. (2022a) proposed the flux-bias model, which accounts for the effects of the accretion torque and/or the wandering of the hot spot on the NS surface, to explain the structures left in the timing residuals across the outburst.Chen et al. (2022) reported the detection of 73 type I X-ray bursts by Insight-HXMT from MAXI J1816−195, and obtained an upper limit of the distance as 6.3 kpc.Wang et al. (2023, in preparation) found the type I Xray bursts rate decreasing for the bolometric flux higher than ∼ 1.04 × 10 −8 erg cm −2 s −1 and estimated the distance to the source at ∼ 3.4 kpc (Cavecchi et al. 2020).
Since the onset of the 2022 outburst a search for radio, optical and near-infrared emission from MAXI J1816−195 had been carried out.The radio emission of MAXI J1816−195 showed significant evolution, at 1.28 GHz with a flux density of 4.2 ± 0.4 mJy on 2022 June 9 and a 3σ upper limit of 70 µJy on 2022 June 11 by MeerKAT, at 5.5 GHz with a flux density of 135 ± 35 µJy on 2022 June 12 by the Australia Telescope Compact Array (ATCA), and at 6 GHz with a flux density of 457 ± 20 µJy on 2022 June 18 by VLA (Beauchamp et al. 2022;Bright et al. 2022).However, radio pulsations were not detected.The likely optical and near-infrared counterpart of MAXI J1816−195 had been identified (de Martino et al. 2022;Kong 2022).When MAXI J1816−195 evolved towards the quiescent state about one month after the start of the outburst, X-ray emission from a dust scattered halo around the source position was detected by Swift (Beardmore et al. 2022).
In this work, we present the broad band timing and spectral analysis of MAXI J1816−195.We introduce the instruments, Insight-HXMT, NICER and NuSTAR, and the performed observations in Sect. 2. The timing and the broad-band spectral results are reported in Sect. 3 and Sect. 4, respectively, and our findings are discussed in Sect. 5.

OBSERVATIONS AND DATA REDUCTION
We analyzed the data from MAXI J1816−195 collected by Insight-HXMT, NICER, and NuSTAR.

Insight-HXMT
Insight-HXMT (Insight Hard X-ray Modulation Telescope, Zhang et al. 2020) is the first Chinese X-ray telescope, and is equipped with three slat-collimated instruments: the Low Energy X-ray telescope (LE, 1-12 keV; Chen et al. 2020), the Medium Energy X-ray telescope (ME, 5-35 keV;Cao et al. 2020) and the High Energy X-ray telescope (HE, 20-350 keV;Liu et al. 2020), providing capabilities for broadband X-ray timing and spectroscopy.Insight-HXMT carried out highcadence observations of MAXI J1816−195 starting on MJD 59738.17,∼0.5 days after the detection by MAXI.The set of 79 ME and 77 HE observations includes runs P0404275001 -P0404275024.The HE and ME data were studied to investigate the timing and spectral properties of the source because of their broadband Xray coverage and good time resolution ( ∼2 µs for HE, ∼20 µs for ME).We analyzed the data using the Insight-HXMT Data Analysis Software (HXMTDAS) version 2.05.The ME and HE data were calibrated by using the scripts mepical and hepical, respectively.The good time intervals were individually selected from the scripts megtigen and hegtigen with the standard criteria, that is, ELV > 0, the satellite located outside the South Atlantic Anomaly region and the offset angle from the pointing direction is smaller than 0. ′′ 04.Totally, 73 type I X-ray bursts have been observed by ME (Chen et al. 2022), which were removed in the timing and spectral analysis.No bursts were found in HE because the burst spectra are usually soft.Background subtracted light curves for the ME and HE were generated (see Figure 2.3).The outburst light curve of Insight-HXMT ME showed a fast rise during the early days of the outburst and a slow decay in the next 25 d.However, the Insight-HXMT HE light curve for the 20-50 keV band dropped continuously from ∼ 100 cts s −1 since the start of the observations to the quiescent level during the outburst.The spectra and their response matrix files are produced by the tools hespecgen and herspgen for HE, and mespecgen and merspgen for ME, respectively.Finally, we obtained the cleaned events using mescreen and hescreen and barycentered with the tool hxbary.

NICER
started regular observations of MAXI J1816−195 on 2022 June 7 (MJD 59737.6;obs.ID 5202820101), several hours after the detection by MAXI, and ended these on 2022 August 3 (MJD 59794.01;obs.ID 5533013504) when the source was already several weeks in its quiescent state.We performed the standard data processing using the NICER Data Analysis Software (NICERDAS).The default filtering criteria were applied to extract the cleaned event data.Next, we extracted 1-s light curves by xselect for the 0.5-10 and 12-15 keV bands to search for type I X-ray bursts and flaring background, respectively.Time intervals containing X-ray bursts or background flares were removed in further analysis.The 0.5-10 keV light curve showed a profile similar to Insight-HXMT ME, i.e., a fast rise in first few days, and a decay to quiescent level during the next ∼ 25 days (see the top panel of Figure 2.3).
The background spectra were produced from the nibackgen3C50 tool (Remillard et al. 2022).The redistribution matrix file and ancillary response file were generated from the tools nicerarf and nicerrmf, respectively.
In the pulse-profile analysis of the NICER data (see Sect. 3.2) we used obs.IDs 5533013001-5533013502 ; cleaned exposure time 14.388 ks) performed when the source had already reached its quiescent state with undetectable (pulsed) flux levels, for the estimation of the local background to obtain background corrected fractional amplitudes.

NuSTAR
NuSTAR observed MAXI J1816−195 on 2022 June 23 (obs.ID 90801315001; MJD ∼ 59753.456− 59754.413)for about 35.7 ks.We cleaned the event file using the NuSTAR pipeline tool nupipeline for both FPMA and FPMB.The light curves and spectra were extracted from a circular region with a radius of 100 ′′ centered on the source location by using nuproducts, and the response and ancillary response files were produced simultaneously.To extract the background spectra, we chose a source-free circular background region centered at (α 2000 , δ 2000 ) = (18 h 16 m 28.s 4048, −19 • 39 ′ 56.′′ 308) with a radial aperture of 100 ′′ .From the light curves, four type I X-ray bursts were identified (see also Mandal et al. 2023).For the timing and spectral analysis the time intervals during these bursts are discarded resulting in a total exposure time of 29.6 ks.
In the timing analysis, we barycentered event data from the source region, adopting a 90 ′′ extraction radius, using HEASOFT multi-mission tool barycorr v2.16,The light curves of MAXI J1816−195 during its 2022 outburst.From top to bottom: the 100 s binned light curves from NICER, Insight-HXMT ME, and Insight-HXMT HE are displayed, respectively.The energy range of each light curve is indicated in each panel.The grey area marks the time of NuSTAR observation.Time intervals covering X-ray bursts are removed.The blue dashed lines represent the time interval that we used for the derivation of our timing models, see Table 1 and Sect.3.
with NuSTAR fine clock-correction file #142, yielding time tags accurate at 60-100 µs level (Bachetti et al. 2021) in absolute time.To obtain background corrected timing characteristics, such as fractional amplitudes, we used the above mentioned background region, however, now with a 90 ′′ extraction radius.

TIMING ANALYSIS
Irrespective the instrument in timing analyses we converted the Terrestial Time (TT) arrival times of the selected events to arrival times at the solar system barycenter (in a Barycentric Dynamical Time (TDB) time scale).For this process, throughout in this work we used: (1) the JPL DE405 solar system ephemeris; (2) the instantaneous spacecraft position with respect to the Earth's center; (3) the (most accurate) VLA position of MAXI J1816−195, α 2000 = 18 h 16 m 52.s 41168 and δ 2000 = −19 • 37 ′ 57.′′ 40138 (Beauchamp et al. 2022), which differ by 0. ′′ 4 in declination from the Swift location (Kennea et al. 2022).Moreover, in the selection procedures we ignored time intervals during which type I X-ray bursts or high-background flaring occurred.

Insight-HXMT and NICER Timing analysis
Initially, we constrained the timing analysis to Insight-HXMT HE data using only events with measured energies between 20 and 35 keV.Before generating phasecoherent timing models (ephemerides) describing accurately the spin frequency evolution across the outburst we first need to correct the barycentered event arrival times for the periodic orbital motion effects.Early mission data covering only a part (MJD 59738-59752) of the full dataset available for Insight-HXMT HE analysis had been used to improve the orbital parameters of MAXI J1816−195 applying an optimisation scheme based on a SIMPLEX algorithm (see De Falco et al. 2017a;Li et al. 2021, for more details) by iteratively improving the Z 2 2 (ϕ)-test statistics (Buccheri et al. 1983) with respect to the pulse frequency ν, the orbital period P orb , and time of the ascending node T asc .The optimised values for the (most critical) orbital parameters P orb and T asc , while fixing the eccentricity to zero and the projected semi-major axis a x sin i to the value quoted in Bult et al. (2022b) are listed in Table 1.The obtained values are fully consistent within 1σ with the best orbital parameters obtained later by Bult et al. (2022a) using the full NICER database.Therefore, we kept the values for the orbital parameters as derived from our three parameter optimisation scheme using 20-35 keV Insight-HXMT HE data.
After correcting for the orbital motion effects we performed a time-of-arrival (ToA; see Kuiper & Hermsen 2009, for more details) analysis to obtain accurate estimates for the pulse frequency ν, and (if required) its first time derivative ν over a certain time interval.We obtained 69 ToA's for Insight-HXMT HE (20-35 keV) covering MJD 59738-59761.Similar analyses yielded 73 ToA's for Insight-HXMT ME restricted to the 10-20 keV band (MJD 59738-59761), and 26 ToA's for NICER using 1.8-10 keV events across MJD 59737-59761.
Fitting the complete Insight-HXMT HE ToA dataset with a constant frequency model resulted in a very poor fit with an unacceptable χ 2 .Ignoring 22 ToA's from the first ∼ 5 days of the outburst yielded reasonable fits assuming a constant frequency model and a model including a frequency derivative term.Both are listed in Table 1.The constant frequency model is fully consistent at 1σ confidence with the corresponding model, based on solely NICER data, shown in Table 1 of Bult et al. (2022a) for 46 degrees of freedom.
The model with both the ν and ν terms, however, improves upon the constant frequency model by reducing the χ 2 from 86.78 to 68.18 adding one fit parameter.Applying a maximum likelihood ratio test to assess the model improvement yields an improvement of (86.78 − 68.18) 0.5 = 4.3σ, adopting one degree of freedom (=the additional fit parameter). 1 The measured ν of (9.0 ± 2.1) × 10 −14 Hz s −1 can be interpreted as a spin-up across the time interval MJD 59741.9-59760.6.We have also investigated the effects of employing a flux-bias model in the ToA fit procedure by adding a (flux) bias term in the form, Φ bias = α(C/C max ) Γ with C max = 895.37c s −1 (the maximum observed count rate in the 0.5-10 keV NICER bandpass), to the pulse-phase residual relation in analogy with Bult et al. (2022a) (see their Eq.1).Note, that in Table 1 of Bult et al. (2022a) both optimized flux-bias model parameters, their b and Γ, are listed abusively with a plus sign.Taken this sign inconsistency into account we could reproduce their findings and obtained Γ = −1.05± 0.15, consistent within one sigma confidence with their value of −1.20 ± 0.20 employing the full NICER ToA dataset running from MJD 59737.748-59762.610.Constraining the NICER ToA dataset to MJD 59742.370-59762.610, and so ignoring the rising and maximum part of the light curve (see Figure 2.3), we obtained a flux-bias index of Γ = −0.77± 0.09, however, as with the full range NICER model shown before, with a statistically unac- 1.8-10 keV) ToAs with respect to a two-parameter timing model (including a quadratic / spin-up term) valid for the range MJD 59741.9-59760.6.These models are based on solely Insight-HXMT HE ToAs (black data points).The interval bounded by the two blue dashed lines represents this time interval that has been used to derive the timing models given in Table 1.Note, that there are no ToA measurements beyond MJD 59760.6 for Insight-HXMT HE (and ME) contrary to NICER.
ceptable fit quality of χ 2 r ∼ 3.54 compatible with the results shown in Bult et al. (2022a).Next, we fitted a constant frequency-plus flux-bias model with a fixed index Γ = −0.77 to the Insight-HXMT-HE ToA dataset covering MJD 59741.9-59760.6,the decaying part of the outburst.This model improves upon the constant frequency model by reducing the χ 2 from 86.78 to 83.10 adding one free parameter (the flux-bias model scale α), which represents a (86.78 − 83.10) 0.5 = 1.92σ improvement with respect to the constant frequency model, applying a maximum likelihood ratio-test.The best fit flux-bias scale parameter is α = (−7.1 ± 3.7) × 10 −3 , while ν = 528.611105 851(11) Hz.Therefore, the spinup model for the decaying part of the outburst (covering MJD 59741.9-59760.6)provides a better data description than the flux-bias model for the Insight-HXMT-HE ToA dataset.
The pulse-phase residuals of the complete sets of ToA's from all three instruments combined (each determined for a different energy band) after folding upon the (Insight-HXMT HE 20-35 keV based) spin-up model are shown in Figure 1.Irrespective the energy band each instrument showed a quite similar behaviour: during the first ∼ 5 days of the outburst -the rising flux part -there is a structured residual feature which can be described by a (short-duration) spin-down model.This feature is the cause of the very poor χ 2 when the full ToA datasets are considered (see also Table 1 of Bult et al. 2022a, indicating unacceptable χ 2 fit values).Also, the steep positive excursions for the NICER ToAs at the early beginning and very end (not covered by Insight-HXMT observations) of the outburst are striking in Figure 1.It is clear that even a (bolometric) flux bias model as well as the constant frequency and quadratic models considered in this work cannot describe the full set of residuals, consistent with Bult et al. (2022a).

Pulse profile analysis
With the accurate orbital and (spin-up) timing models (see Table 1) we converted the event arrival times from each high-energy instrument to pulse phases.To obtain high-statistics pulse-phase distributions (pulse profiles including background) we combined the phase information for each instrument across the validity interval (MJD 59741.9-59760.9)and sorted these on (measured) energy.We verified for each instrument that the pulse shape for each individual observation constituting the combination remained the same within statistical uncertainties.The measured pulse-phase distributions for the four different high-energy instruments are shown in Figure 2 for various different energy bands covering the 0.8-210 keV range.
Panels (a)-(c) of Figure 2 show the NICER pulse profiles for the energy intervals 0.8-1.8,1.8-3, and 3-10 keV, respectively, while panels (d)-(f) show the Insight-HXMT ME profiles for the 5-10, 10-20 and 20-35 keV bands.The NuSTAR profiles are displayed in panels (g)-(k) covering the 3-79 keV energy range.Noteworthy is that comparing the NICER and NuSTAR 3-10 keV distributions, shown in panels (c) and (g), respectively, a small phase shift of 0.065 ± 0.020 is detectable which is just compatible with the absolute timing accuracy of NuSTAR (Bachetti et al. 2021).Pulsed emission is detected significantly up to ∼60 keV, with a nonuniformity significance of 5.4σ for the 35-60 keV band (panel j) adopting a Z 2 2 -test (Buccheri et al. 1983).This test applied to the 60-79 keV band (panel k) yielded an insignificant result of 1.3σ.Finally, panels (l)-(p) show the Insight-HXMT HE profiles for energies of 20-30, 30-50, 50-65, 65-95 and 95-210 keV, respectively.Interestingly, pulsed emission is detected significantly up to ∼95 keV.The 65-95 keV band (panel o) shows a nonuniform significance of 4.4σ applying a Z 2 2 -test with indications for a second weaker pulse near phase 0.15 that seems to pop-up at energies above ∼50 keV (panels n- o).The non-uniformity significance for the distribution shown in panel (p) for the 95-210 keV band is about 1.1σ, and so insignificant.
To obtain quantitative information about morphology changes of the pulse profiles as a function of energy, we fit the pulse profiles with a truncated Fourier series given by a formula where A 1 and A 2 are the amplitudes, ϕ 1 and ϕ 2 are the phase angles (in radians/2π), of the fundamental and the first overtone, respectively, and A 0 is the constant level of the profile.We verified that higher harmonics are not required to accurately describe the data.The left panel of Figure 3 shows the background corrected fractional amplitudes A k /A 0 of the fundamental (black) and first overtone (red) for all four instruments covering the 0.8-95 keV band, across which we detected significant pulsed emission.
The background correction determination was straightforward for NICER and NuSTAR.We used for the former non-imaging instrument the data from observations when the source was completely Off (see Sect. 2.2) with a combined Good Time Interval (GTI) exposure of 14.388 ks.Along with the On-time GTI exposure of 51.735 ks this yielded a scale factor s bg of 3.5957 to be applied in the background subtraction.For the latter (imaging) instrument the data from a differ- ent source-free location in the field-of-view with equal area as source region (s bg ≡ 1.0) was used.
For the non-imaging Insight-HXMT ME/HE (with GTI exposures of 186.171 ks and 132.093 ks, respectively, for the validity interval during the On source period) the last observations of MAXI J1816−195 performed on MJD 59765 (2022 July 5) totaling exposure times of 13.474 ks and 9.022 ks for ME and HE, respectively, were used to estimate the underlying background, when the source did not show pulsed emission anymore, but was not completely Off (see, e.g., Beardmore et al. 2022, who reported the detection of the source at a 4.1σ level in 5.2 ks Swift XRT data taken a day later on MJD 59766).
Thus, the ME/HE event data of MJD 59765 still contain a tiny contribution from MAXI J1816−195, that would yield an over-subtraction when applied uncorrected in the background subtraction procedure.Demanding that the background subtraction did not yield a negative value for energy bands at the high end of the sensitivity window of both instruments (ME ≳ 30 keV; HE ≳ 95 keV) we could estimate the genuine background contribution in both samples.For the ME/HE, we found that ∼9% and ∼5% of the events from the background sample was originating from the 'almost' Off MAXI J1816−195, respectively.Taking these fractions into account in the background subtraction pro-cess showed that the A 1 /A 0 (and also A 2 /A 0 ) values obtained for ME/HE matched with NuSTAR values obtained in the overlapping energy range.
It is clear from the left panel of Figure 3 that the pulsed fraction of the fundamental (black) gradually increases with energy from ∼0.5% to 10%-12% going from ∼0.8 to ∼10 keV, where a plateau is reached.Above ∼35 keV a decrease sets in to values of 3%-5% above 50 keV.The first overtone (red) saturates at ≳10 keV to a ∼2% value, without a clear decrease beyond.This means that the A 2 /A 1 ratio increases the higher the energy, implying that the first overtone contributes more significantly at higher energies.
The right panel of Figure 3 shows the phases of the maxima of the fundamental (black) and first overtone (red) for the four instruments, specifying the alignment of the components constituting the pulse profile.It is noteworthy that these quantities, ϕ 1 and ϕ 2 , are independent on the used background method.
For NuSTAR, which provided a snapshot during the full outburst, we have subtracted 0.065 (in phase units) -still within its absolute timing accuracy -from the measured ϕ 1 and ϕ 2 values for alignment purposes with NICER.It is clear that irrespective the instrument both quantities are stable within statistical uncertainties till ∼50 keV, beyond which the first overtone moves closer to the fundamental.Combined with the increase of the A 2 /A 1 ratio at those energies this means that a second bump/maximum in the pulse profile in front of the main maximum becomes more prominent at ≳50 keV energies.

SPECTRAL ANALYSIS
The spectral analysis was carried out using xspec version 12.12.1 (Arnaud 1996).We first performed the joint spectral evolution analysis for NICER and Insight-HXMT observations (Sect.4.1).And then, we studied the quasi-simultaneous broadband spectra of NICER, NuSTAR and Insight-HXMT (Sect.4.2).All spectra are grouped by the tool ftgrouppha.For NICER and NuS-TAR spectra, we applied the optimal binning method (Kaastra & Bleeker 2016).The Insight-HXMT spectra were grouped with a minimal signal-to-noise ratio of 3 related to their backgrounds.All uncertainties of the spectral parameters are provided at a 1σ confidence level for a single parameter.

The spectral evolution of joint NICER and Insight-HXMT observations
There were many quasi-simultaneous NICER and Insight-HXMT observations of MAXI J1816−195, which allows us to perform a detailed broadband spectral analysis for the whole outburst.We found that each NICER observation in the MJD 59738-59761 range can be (partially) covered by 1-6 Insight-HXMT observations.We combined the Insight-HXMT spectra, and then carried out joint NICER and Insight-HXMT spectral fitting.We adopted the energy ranges of 1-10 keV for NICER spectra, 10-35 keV for Insight-HXMT ME and 25-150 keV for Insight-HXMT HE.We fit all the spectra by using the combination of the thermal Comptonization model, compps, in the slab geometry (Poutanen & Svensson 1996), the disk blackbody, diskbb, modified by the interstellar absorption described by the model, tbabs.This model has been used to fit the broadband spectra of many AXMPs (see, e.g., Gierliński & Poutanen 2005;Falanga et al. 2005aFalanga et al. ,b, 2011Falanga et al. , 2012;;De Falco et al. 2017b,a;Kuiper et al. 2020;Li et al. 2018Li et al. , 2021)).To account for cross-calibration uncertainties between different instruments a multiplication factor is included in the fits.The factor is fixed at unity for NICER and set free for Insight-HXMT ME/HE.The main parameters are the Thomson optical depth across the slab, τ T , the electron temperature, kT e , the temperature of the soft seed photons, kT bb , the normalization factor for the seed black body photons, K bb = (R km /d 10 ) 2 (with d 10 being distance in units of 10 kpc), and the inclination angle θ (fixed at 60 • ) between the slab normal and the line of sight to the observer for the compps model; the disk temperature, kT diskbb and normalization for the multi-black body model, and the hydrogen column density, N H .Most of NICER spectra showed emission lines around ∼1.6 and ∼6.5 keV, therefore, we added two Gaussian components to account for these.The full model is tbabs (gaussian+gaussian+diskbb+compps) in xspec.The ∼1.6 keV component could have an instrumental origin, while the ∼6.5 keV Gaussian component may come from the iron Kα line emission.
The fitted parameters are shown in Figure 4. Most of the spectra can be well fitted with reduced χ 2 < 1.3.For two observations, the reduced χ 2 are larger than 1.5, which show no features in the NICER residuals and possibly caused by the spectral evolution during the long exposures.We calculated the unabsorbed bolometric flux in the 1-250 keV range by using the tool cflux.During the outburst, the disk blackbody temperature increased from 0.4 to 0.8 keV in the first few days and then dropped back to ∼0.4 keV.The hydrogen column density did not change much, and the mean value is (2.31 ± 0.07) × 10 22 cm −2 .The optical depth was in the range 1.7-3.0,which explains the visibility of hard Xray emissions during the whole outburst.The electron temperature was around 10-25 keV.The bolometric flux reached a peak value of 1.02×10 −8 erg cm −2 s −1 during the first few days, and in the next 20 days, it decayed to 8.5 × 10 −10 erg cm −2 s −1 at MJD 59760.9.
We also replaced the compps component by a cutoff power-law model.This model provided compatible results compared with the compps model, i.e., N H , the temperature and normalization of disk black body component, the χ 2 , and the unabsorbed flux.

Broadband spectra including NuSTAR observation
We noted that the only NuSTAR observation was carried out quasi-simultaneously with some of NICER and Insight-HXMT observations.We fitted the spectra, including the Insight-HXMT ME and HE data between MJD 59753.46-59754.40(obs.ID P0404275015) with a total exposure of 20.99 ks (ME) and 17.51 ks (HE), NuSTAR FPMA (29.55 ks), FPMB (29.70 ks) and NICER data starting on MJD 59753.30(5.21 ks, obs.ID 5533011501).We kept the default energy range of 3-79 keV for NuSTAR FPMA/FPMB, and for the other instruments the same ranges as in Sect.4.1.The multiplication factors were also included to account for the cross-calibration uncertainties between the instruments.The factor was fixed at unity for the NuS-TAR FPMA, and set free for the other instruments.We first fitted the spectra by using an absorbed cutoff power-law model in combination with a disk black From top to bottom, the bolometric flux (in units of 10 −8 erg cm −2 s −1 ), the hydrogen column density (in 10 22 cm −2 ), the electron temperature (in keV), the optical depth, the temperature of disk blackbody (in keV), and the reduced χ 2 are shown.
body model and two additional Gaussians.The best fitted model yielded N H ≈ 2.30 × 10 22 cm −2 , a photon index Γ ≈ 1.54, a cutoff energy E cut ≈ 24.1 keV, the Gaussian at 5.44 keV with a width of 1.32 keV, a disk temperature kT diskbb ≈ 0.61 keV, and a reduced χ 2 ∼ 1.15 for a total d.o.f. of 2300.We then replaced the phenomenological cutoff power-law model by the thermal Comptonization model compps, which is tbabs (gaussian+gaussian+diskbb+compps).We obtained a best-fit reduced χ 2 of 1.17.We noticed that a reflection feature has been reported by Chauhan et al. (2022), therefore the reflection parameter, R, of compps was set free.The best-fit reduced χ 2 decreased to 1.11 for a total d.o.f. of 2298, better than the cut-off powerlaw model.The best-fit parameters are listed in Table 2 and the broad-band spectrum is shown in Figure 5.The NICER spectrum contributed the most significant residuals.Considering such a large data set, this model provides an acceptable fit to the spectra.We obtained a normalization factor for the seed photons, K bb = 65.9±6.8,which implies a size of blackbody emitting region R bb = 0.34 √ K bb ≈ 3 km for a distance of 3.4 kpc (Wang et al. submitted, see also Chen et al. 2022).This region is smaller than the whole NS surface with a typical radius of 12 km.However, the true emitting size should be larger since part of photons were reflected by the accretion disk.It is worthy to note that it is difficult to constrain the refection contribution by only using the joint NICER and Insight-HXMT/ME/HE spectra, i.e., without NuSTAR data, the χ 2 r is 1.0, see Sect.4.1.We also changed the compps model to a relativistic reflection model, relxillCp (Dauser et al. 2016), in which the reflection spectrum is calculated using the Comptonization continuum nthcomp.After many attempts, during the fit, we set the iron abundance, A Fe , the density of the accretion disk, log N , the ionization parameter, log ξ, the power-law index of the incident spectrum, Γ, the reflection fraction, R f , the outer emissivity index, q 2 , free to vary.The spin parameter was fixed at 0.248 for the NS spin frequency of 528 Hz.The outer disk radius was fixed at 1000 R ISCO .The inclination angle, i, was also fixed at 60 degree.Other parameters were fixed at their default values.This model gave a best-fit reduced χ 2 of 1.11 for a total d.o.f. of 2295, slightly larger than the compps model.The results are listed in Table 2.

SUMMARY AND DISCUSSION
In this paper we reported on detailed timing and spectral analyses of Insight-HXMT, NICER and NuS-TAR data collected from the newly confirmed AMXP MAXI J1816−195 during its 2022 outburst.The outburst of MAXI J1816−195 reached its peak flux within the first few days at ∼ 1 × 10 −8 erg cm −2 s −1 and de- F bol (10 −9 erg s −1 cm −2 ) a 3.97±0.01

3.96±0.01
Note-Best parameters determined from the fits to the average broad-band spectrum of MAXI J1816−195 performed with the NICER/NuSTAR/Insight-HXMT data.The multiplication factor for all instruments are provided.
a Unabsorbed flux in the 1-250 keV energy range.
cayed into quiescent state in the next ∼ 25 days.The outburst profile observed from MAXI J1816−195 is quite similar to that of other AMXPs, showing a fast rise (∼1-2 days), an exponential followed by a linear decay, which can be explained by the disk instability model (King & Ritter 1998;Powell et al. 2007).
During the full outburst hard X-ray/soft γ-ray emission was detected up to ∼ 200 keV.The joint NICER and Insight-HXMT spectra in the energy of 1-150 keV were well fitted by a combination of a multi-color disk black body and Comptonization model compps, with two additional Gaussian lines, exposed to interstellar extinction.A joint spectral fit including quasisimultaneous NuSTAR, NICER and Insight-HXMT ME/HE data showed emission reflected from the accretion disk, which can be modelled by the reflected compps model or the relxillCP model.
Usually, the hard X-ray emission from AMXPs is relatively weak and has typically pulsation amplitudes of less than 10% (see e.g., Papitto et al. 2020, and references therein; see also Falanga et al. 2005a for IGR J00291+5934 with the pulsed fraction of ∼ 12 − 20% at ∼ 100 keV).Therefore, the pulsations of most AMXPs are mostly detected at soft X-rays with instruments aboard e.g RXTE, XMM-Newton or NICER.At hard Xrays, Insight-HXMT observed Swift J1756.9-2508 during its 2018 outburst, and detected the X-ray pulsations up to ∼ 100 keV using an ephemeris determined by NICER at soft X-rays (Li et al. 2021).In the case of MAXI J1816−195, at hard X-rays both the ME and HE aboard of Insight-HXMT were capable of deriving the pulsar ephemeris directly from their event data, obtaining parameters well consistent with those derived by NICER at soft X-rays.Particularly, our timing analysis using the DE405 solar system ephemeris and the most accurate VLA (Beauchamp et al. 2022) source location for MAXI J1816−195 in the barycentering process, contrary to Bult et al. (2022a) who used DE430 and the less accurate (∼ 0. ′′ 4 difference) Swift location (Kennea et al. 2022), yielded consistent orbital and spin parameters with Bult et al. (2022a) at 1σ-confidence.
From Insight-HXMT observations hard X-ray pulsations have been detected significantly up to ∼ 95 keV.The pulse profiles in 0.8-95 keV band are well described by a truncated Fourier series.The phase angles of the fundamental and the first overtone are nearly constant for energies below ∼ 60 keV, implying that most photons are likely emitted from the same region at the NS surface with a very similar emission pattern.The observed pulse profiles could be used to measure the NS mass and radius (Poutanen & Gierliński 2003) of MAXI J1816−195 in a future study.
The stability of the pulse profiles across a broad energy range allows us to determine the absolute timing accuracy of Insight-HXMT in detail.Cross-correlating Insight-HXMT ME and NICER pulse profiles for the overlapping 5-10 keV energy band yielded a phase shift δΦ of only −0.008 ± 0.005, consistent with zero, suggesting that both instruments are well aligned within 14 µs, with the Insight-HXMT ME events arriving a little bit later than NICER ones.A similar correlation analysis now between Insight-HXMT ME and HE pulse profiles for the 20-35 keV band resulted in a difference of only 0.9 µs with an uncertainty of 13.9 µs.Bult et al. (2022a) performed a coherent timing analysis adopting a constant frequency-and a flux bias model using NICER ToA's covering the full outburst period.They found unacceptably high χ 2 values for both models assumptions resulting in complex timing residuals for the first 4-5 days during the rising part of the outburst, a flat shape during the next ∼ 20 days in the decay phase, and finally, (linearly) increasing phase deviations at the end of outburst, also visible at the very beginning of the outburst.
In this work we analysed the timing data from Insight-HXMT ME and HE observations as well as from NICER, and found that the ToAs from the different (independent) instruments were well consistent with each other in spite their different bandpass.The timing residuals showed erratic behavior (or V-shaped wedge) during the rising part of the outburst between MJD 59737.0 and 59741.9,similar to the results shown in Bult et al. (2022a).
For the decaying part of the outburst between MJD 59741.9-59760.6we detected a significant positive quadratic term ν of (9.0 ± 2.1) × 10 −14 Hz s −1 fitting a two-parameter timing model, including ν and ν, yielding an acceptable χ 2 , which could interpreted as a spin-up.
We realize that such spin-up or spin-down terms can arise when an inaccurate source location is used in the barycentering process.However, in our work we have used the sub-arcsecond accurate VLA-location for MAXI J1816−195, with an accuracy better than 0. ′′ 005.This uncertainty on the location can account for an uncertainty in the spin-up rate (see e.g., Sanna et al. 2020, for the formulae) of only ∼ 2.0 × 10 −16 Hz s −1 , 2.7 order of magnitude smaller than the measured ν value, and much smaller than our error estimate of 2.1 × 10 −14 Hz s −1 .On the contrary, if we should have used the Swift-XRT location in the barycentering process with a centroid 0. ′′ 4 from the VLA-location and an estimated error radius of 2. ′′ 2 then the uncertainty in the quadratic term would be ∼80 times larger, i.e. of the order of 1.6 × 10 −14 Hz s −1 , comparable to our measurement error.Therefore, we assume that the measured value for ν during the decay phase is accurate and genuine, and could be interpreted as a gain of angular momentum by the NS from the accreted matter.
During the outburst X-ray pulsations have been detected at bolometric flux states with a flux between ∼ (0.1 − 1) × 10 −8 erg cm −2 s −1 .Adopting a source distance of 3.4 kpc, as constrained by Wang et al. (2023, submitted, see also Chen et al. 2022), these bolomet-ric fluxes correspond to accretion rates in the range (0.012 − 0.12) ṀEdd , where ṀEdd = 2 × 10 −8 M ⊙ yr −1 is the Eddington rate, assuming that L = η Ṁ c 2 and a typical value of η = 0.1 for NS (Frank et al. 2002).By adopting Equations 11 and 12 and the same assumptions in Psaltis & Chakrabarty (1999), these limits allow us to constrain the surface dipole magnetic field strength of (0.25 − 74) × 10 8 G (see also Hartman et al. 2008).If 20% uncertainty of distance is considered, the magnetic field is in the range of (0.2 − 85) × 10 8 G.The possible bolometric correction factor of the flux would not change the magnetic field significantly because the measured fluxes are covered a broadband energy range.
We appreciate the referee for all valuable comments and suggestions that led to the improvement of the manuscript.This work was supported the Major Science and Technology Program of Xinjiang Uygur Autonomous Region (No. 2022A03013-3).Z.S.L., Y.Y.P. and L.J. were supported by National Natural Science Foundation of China (12103042, U1938107, 12273030, 12173103).S.Z. is supported by the National Key R&D Program of China (2021YFA0718500).J.P. acknowledges support from the Academy of Finland grant 333112.
This work is also supported by International Partnership Program of Chinese Academy of Sciences (Grant No.113111KYSB20190020).This work made use of data from the Insight-HXMT mission, a

Figure 1 .
Figure 1.Pulse phase residuals of Insight-HXMT HE (open black triangles; 20-35 keV), Insight-HXMT ME (open orange squares; 10-20 keV) and NICER (open aqua diamonds;1.8-10keV) ToAs with respect to a two-parameter timing model (including a quadratic / spin-up term) valid for the range MJD 59741.9-59760.6.These models are based on solely Insight-HXMT HE ToAs (black data points).The interval bounded by the two blue dashed lines represents this time interval that has been used to derive the timing models given in Table1.Note, that there are no ToA measurements beyond MJD 59760.6 for Insight-HXMT HE (and ME) contrary to NICER.

Figure 3 .
Figure3.Left: The fractional amplitude of various harmonics: 1 -fundamental (black crosses) and 2 -first overtone (red crosses) as a function of energy using background-corrected data from NICER (filled squares), NuSTAR (triangles), Insight-HXMT ME (filled circles) and Insight-HXMT HE (open circles).Right: The phase angle (in units radians/2π) as a function of energy for the two Fourier components, the fundamental (in black) and the first overtone (in red), which are independent of the background subtraction method.

Figure 4 .
Figure4.The fitted results of joint NICER and Insight-HXMT spectra by the model constant×tbabs(gaussian + gaussian + diskbb + compps).From top to bottom, the bolometric flux (in units of 10 −8 erg cm −2 s −1 ), the hydrogen column density (in 10 22 cm −2 ), the electron temperature (in keV), the optical depth, the temperature of disk blackbody (in keV), and the reduced χ 2 are shown.

Figure 5 .
Figure 5.The unabsorbed spectra of MAXI J1816−195 from NICER, NuSTAR FPMA, NuSTAR FPMB, Insight-HXMT ME and Insight-HXMT HE in the 1-150 keV energy range.The fitting model is constant×tbabs (gaussian+gaussian+diskbb+compps).The black dashed lines represent the two gaussian lines.The best model is plotted as black solid line.The disk blackbody and compps components are shown in dotted and dash-dot-dot-dot lines, respectively.

Table 1 .
Positional, orbital and spin parameters derived in this work and used as fixed values from literature for MAXI J1816−195.