Timing of the accreting millisecond pulsar IGR J17591$-$2342: evidence of spin-down during accretion

NICER (Gendreau et al., 2012) started monitoring the X-ray transient IGR J17591$-$2342 on 2018 August 15 (MJD 58345) up to 2018 October 17 (MJD 58408.2) for a total of 37 dedicated observations. We processed the available NICER observation using the nicerdas pipeline (version V004a) retaining events in the 0.2–12keV energy range, for which the pointing offset was $<54^{\prime\prime}$, the dark Earth limb angle was $>30^{\circ}$, the bright Earth limb angle was $>40^{\circ}$, and the ISS location was outside of the South Atlantic Anomaly (SAA). Moreover, we removed short intervals that showed background flaring in high geomagnetic latitude regions. As a result of the filtering process we retained a total exposure time of $\sim 101$ ks. The spectral background was obtained from a blank-sky observations of the RXTE-6 region. We did not detect any Type-I thermonuclear bursts during the observations analysed in this work. For each observation we modelled the energy spectrum with an absorbed black body component plus a thermally comptonised continuum (Tbabs(bbody+nthcomp) in Xspec) and we extracted the source unabsorbed flux in the energy range 0.5–10 keV (see Fig. 1). Finally, we barycentred the NICER photon arrival times with the BARYCORR tool using the JPL DE-405 Solar System ephemeris and adopting the source coordinates obtained from the ATCA radio detection of the source (Russell et al., 2018).

To perform phase-coherent timing analysis we started by correcting the NICER photon times of arrival for the binary Doppler delay estimated assuming an almost circular orbits (see e.g. Deeter et al., 1981; Burderi et al., 2007; Sanna et al., 2016) with parameters equal to that reported by Sanna et al. (2018d) obtained from the analysis of a small portion of the source outburst (almost 10 out 60 days). We created pulse profiles epoch-folding the available data at the previously determined spin frequency $\nu_{0}=527.42570042$ Hz. To optimise the statistics and increase the number of profiles, we split the datasets in segments with length ranging between 200 and 1500 seconds. We modelled the pulse profile as the combination of two harmonically related sinusoidal functions: the fundamental and the first harmonic components of the spin frequency. To proceed further with the analysis, we selected only statistically significant pulse profiles, i.e. profiles for which the ratio between the sinusoidal amplitude and its 1$\sigma$ uncertainty was greater than three. As shown in Fig. 2, coherent oscillations were clearly detected in the time interval MJD 58345 – MJD 58405, covering almost entirely the outburst duration. The fractional amplitude of the fundamental and the first harmonic components stays overall constant around a mean value of $\sim$11.8% and $\sim$5%, respectively (Fig. 3)

Following standard timing techniques (see e.g. Burderi et al., 2007; Sanna et al., 2018d), we modelled the temporal evolution of the pulse phase delays obtained from the fundamental and first harmonic components as follows:

where $\Delta\nu_{0}=(\nu_{0}-\bar{\nu})$ and $\dot{\nu}$ represent the spin frequency correction and the spin frequency derivative, respectively, estimated with respect to the reference epoch $T_{0}$. $R_{orb}$ models the differential corrections to the ephemeris used to generate the pulse phase delays (see e.g. Deeter et al., 1981). To improve the accuracy on the orbital parameters we modelled simultaneously pulse phase delays from the fundamental and first harmonic components using Eq. 1. More specifically, we fitted the pulse phases linking the $R_{orb}(t)$ component, meaning that the orbital parameters will be the same during the fit of the two datasets. We applied iteratively the method, until no significant differential corrections were found for the parameters of the model. In the right column of Tab. 1, we report the the best-fit parameters, while in Fig. 2 we show the pulse phase delays with the best-fitting model (top panel; black and red points for the fundamental and first harmonic components, respectively), and the corresponding residuals (in phase units) with respect to the best-fitting model to fit the time evolution of the fundamental (middle panel) and first overtone (bottom panel) phase delays, respectively. The value of $\tilde{\chi}^{2}\sim 2.5$ (with 355 degrees of freedom), as well as the distribution of the residuals suggest the presence of timing noise, largely observed in several AMXPs (see e.g. Burderi et al., 2006; Papitto et al., 2007; Riggio et al., 2008; Riggio et al., 2011a). We rescaled the uncertainties on the parameters reported in Tab. 1 by a factor $\sqrt{\chi^{2}}$ to give a realistic estimation of the uncertainties given the relatively poor fit of the pulse phases (see e.g. Finger et al., 1999).

To explore the effect of the positional uncertainties on the spin frequency and spin frequency derivative reported in Tab. 1, we considered the residuals induced by the motion of the Earth for small variations of the source position $\delta_{\lambda}$ and $\delta_{\gamma}$ (ecliptic coordinates) expressed by the relation:

where $y=r_{E}/c$ represents the Earth semi-major axis in light-seconds, $M_{0}=2\pi(T_{0}-T_{v})/P_{\oplus}-\lambda$, $T_{v}$ and $P_{\oplus}$ are the vernal point and the Earth orbital period, respectively, and $\epsilon=2\pi(t-T_{0})/P_{\oplus}$ (see, e.g. Lyne & Graham-Smith, 1990). Degeneracy between phase delay residuals induced by positional uncertainties and those related to spin frequency uncertainties and/or spin frequency derivative is unavoidable on timescales much shorter (i.e. 60 days outburst duration) than Earth’s orbital period. Upper limits on these quantities can be obtained expanding Eq.2 in series of $\epsilon\ll 1$ (see e.g. Burderi et al., 2006, and references therein). More specifically, the systematic error on the spin frequency correction and the spin frequency derivative are $\sigma_{\nu_{pos}}\leq\nu_{0}y\sigma_{v}(1+\sin^{2}\gamma)^{1/2}2\pi/P_{\oplus}$ and $\sigma_{\dot{\nu}_{pos}}\leq\nu_{0}y\sigma_{v}(1+\sin^{2}\gamma)^{1/2}(2\pi/P_{\oplus})^{2}$, respectively, where $\sigma_{v}$ is the positional error circle. Considering the uncertainty on the source position reported in Tab. 1 (Russell et al., 2018), we estimated $\sigma_{\nu_{pos}}\leq 7\times 10^{-9}$ Hz and $\sigma_{\dot{\nu}_{pos}}\leq 2\times 10^{-15}$ Hz/s, respectively. These systematic uncertainties, compatible with the statistical uncertainties reported in this work, have been added in quadrature to the statistical errors of $\nu_{0}$ and $\dot{\nu}$ estimated from the timing analysis.

We carried out studies on the properties of the pulse profile as a function of energy by dividing the NICER energy range between 0.2 keV to 12 keV into 20 intervals. Each background-subtracted pulse profile has been modelled as the superposition of two harmonically related sinusoidal functions (fundamental and first harmonic) for which we calculated the fractional amplitudes and the time lags with respect to the lowest energy selection. In the top panel of Fig. 5 we report the fractional amplitude as a function of energy for the fundamental (black points) and first harmonic (red points) pulse phase component of the pulse profile. The fractional amplitude of the fundamental component increases from $\sim 4\%$ up to $\sim 16\%$ at around 5 keV, decreasing then to $\sim 12\%$ at 12 keV. A similar behaviour, but limited between $\sim 2.6\%$ and $\sim 6.6\%$ is shown by the first harmonic. The bottom panel of Fig. 5 suggests the presence of hard lags (pulsed high energy photons are observed earlier with respect to soft ones) almost linearly increasing with energy for the two components. Interestingly, while the first harmonic component only shows soft lags in the energy range explored, the fundamental component shows no lags or marginally hard lags between 0.5–3.5 keV followed by soft lags above 3.5 keV.

After correcting the photon time of arrivals for the timing solution reported in Tab. 1, we epoch-folded the whole NICER dataset of the source into a 32 bins pulse profile. As shown in Fig. 4, the best-fitting model for the profile requires a combination of four sinusoids, where the fundamental, second, third, and fourth harmonics have fractional amplitudes of $~{}11\%$, $~{}4.8\%$, $~{}0.7\%$, $~{}0.16\%$, respectively. Harmonically rich pulse profiles have been reported in several other AMXPs, such as SWIFT J1749.4$-$2807 (Altamirano et al., 2011), IGR J17511$-$3057 (Riggio et al., 2011a), IGR J17379$-$3747 (Sanna et al., 2018c) and IGR J16597$-$3704 (Sanna et al., 2018b).

As reported in the previous section, during the outburst the pulse phase clearly shows significant evolution relative to a constant spin frequency model. More specifically, when modelled with a second order polynomial function, the pulse phase delays revealed a significant spin frequency derivative describing a deceleration of the NS (spin-down) of the order of $\sim-7\times 10^{-14}$ Hz/s. This result is confirmed independently from the timing analysis of both the fundamental and first harmonic components, for which the values of spin-down derivative are compatible, within errors, with each other. The detection of the second (or higher) harmonic component during the whole outburst has been reported also for other AMXPs such as SWIFT J1749.4$-$2807 (Altamirano et al., 2011), XTE J1814-338 (Papitto et al., 2007), IGR J17511$-$3057 (Riggio et al., 2011a), SAX J1808.4$-$3658 (see e.g. Burderi et al., 2006; Hartman et al., 2008), SWIFT J1756.9-2508 (Bult et al., 2018) and XTE J1807$-$294 (see e.g. Riggio et al., 2008; Patruno et al., 2010).

Four other AMXPs show spin-down during their outbursts: XTE J0929$-$314 with a spin frequency derivative of $\dot{\nu}=-9.2(4)\times 10^{-14}$ Hz s^-1 (Galloway et al., 2002); XTE J1814$-$338 with a spin frequency derivative of $\dot{\nu}=-6.7(7)\times 10^{-14}$ Hz s^-1 (Papitto et al., 2007); IGR J17498$-$2921 with a spin frequency derivative of $\dot{\nu}=-6.3(1.9)\times 10^{-14}$ Hz s^-1 (Papitto et al., 2011) and SAX J1808.4$-$3658 for which evidences of a spin-down derivative have been reported during the decay phase of its 2002 outburst ($\dot{\nu}=-7.6(1.5)\times 10^{-14}$ Hz s^-1; Burderi et al., 2006) and during its 2019 outburst ($\dot{\nu}=-3.0(1)\times 10^{-13}$ Hz s^-1; Bult et al., 2019, however, see the discussion for possible interpretations of the result). It is noteworthy that the observed values of spin-down derivative are similar, despite the different outburst properties in which they have been detected as well as the wide range of NS spin frequencies of the AMXPs showing this phenomenon.

From a very general point of view, the torques acting onto a NS with a significant dipolar magnetic field could be computed by adopting few general prescriptions to describe the interaction of the accreting matter with the NS. In particular: matter is accreting through a Keplerian accretion disc, truncated at the radius $R_{\rm m}$ at which matter is forced, by the action of the magnetic field, to hook onto the magnetic field lines, sharing their rigid rotation motion at the NS angular frequency. This condition implies that accretion is centrifugally inhibited if $R_{\rm m}\gg R_{\rm CO}=(GM/\Omega^{2})^{1/3}$ where $G$ and $M$ are the gravitational constant and the NS mass, respectively, and $\Omega=2\pi\nu$ is the angular frequency of the NS.

Matter accreting onto the NS surface generate a material torque determined by the accretion of its angular momentum onto the NS, $\tau_{\rm+\,\dot{M}}$, where $\dot{M}$ is the mass accretion rate. In particular, the NS is surrounded by a Keplerian accretion disc, whose innermost rim is defined by $R_{\rm m}$ and it is accreting the specific Keplerian angular momentum $\ell=\sqrt{GMR_{\rm m}}$ at $R_{\rm m}$ through a material torque:

In the most general case of a Keplerian disc that is threaded by the magnetic field lines (see e.g. Ghosh & Lamb, 1979; Wang, 1987, 1995, 1996; Rappaport et al., 2004; Kluźniak & Rappaport, 2007a), the interaction of the orbiting matter with the magnetic filed lines (rigidly connected with the NS) determine the rise of two torques of opposite signs. A positive torque $\tau_{\rm+\,B}$ generates (we choose to use the plus sign to represent torques acting concordantly with the material torque) from interaction of the orbiting matter with the magnetic filed lines at radii $R_{\rm m}\leq r\leq R_{\rm CO}$, while a negative torque $\tau_{\rm-\,B}$, arises from the interaction of the orbiting matter with the magnetic filed lines at radii $R_{\rm CO}\leq r\leq R_{\rm LC}$, where $R_{\rm LC}=c/\Omega$ is the light-cylinder radius, beyond which the magneto-static structure of the field is truncated and a radiative solution is present (with $c$ being the speed of light).

To compute the NS torque associated with this radiative solution, we adopted, as a reasonable first order guess, that Larmor’s formula holds for the NS regarded as a magneto-dipole rotator. In this case, the radiative solution beyond $R_{\rm LC}$ implies a negative torque acting onto the NS:

where $\mu=B_{\rm eq}R^{3}$ is the magnetic moment of the NS computed from the value of the surface magnetic field at the NS equator, $B_{\rm eq}$, $R$ is the NS radius and $i$ is the angle between the magnetic dipole moment and the spin axis of the NS.

To compute the torque associated to the interaction of the orbiting matter with the magnetic filed lines, it is reasonable to assume that $\tau_{\rm\pm\,B}$ results from the integration of a torque radial density, $d\tau_{\rm\pm\,B}/dr$, from $R_{\rm CO}$ to $R_{\rm LC}$. The torque radial density is determined by specific prescriptions for the structure of the magnetic field in the mid-plane of the accretion disc. In the following, in line with Kluźniak & Rappaport (2007b), we adopted a quite simple prescription for the ratio of the toroidal to poloidal components of the magnetic field threading the accretion disc at a generic radius $r$:

In line with current wisdom, the poloidal magnetic field is a screened version of the NS dipolar field at the disc mid-plane, $B_{z}(r)=\eta\mu/r^{3}$, with $\eta\sim 1$ ¹¹1On the other hand, several authors (see e.g. Wang, 1987; Rappaport et al., 2004; Kluźniak & Rappaport, 2007b, and references therein) discussed different prescription for the magnetic field structure like the following (computed in detail in Kluźniak & Rappaport, 2007b): $$\frac{B_{\phi}(r)}{B_{z}(r)}=\left\{\begin{array}[]{ll}1-\Omega_{\rm K}(r)/\Omega&{\;\;\;\rm for\,}r>R_{\rm CO}\\ &\\ \Omega/\Omega_{\rm K}(r)-1&{\;\;\;\rm for\,}r\leq R_{\rm CO}\\ \end{array}\right.$$ (6) and other more sophisticated and physically motivated versions of the relations above (see the exhaustive discussion in Kluźniak & Rappaport, 2007b). However, within factors of the order unity, the final results on the torques are quite independent of the detailed prescription adopted for the magnetic field, particularly for $R_{\rm m}\rightarrow R_{\rm CO}$ (Kluźniak & Rappaport, 2007b) which will be the case for IGR J17591$-$2342, as we verified a posteriori (see below)..

Let us consider $\tau_{\rm-\,B}$, resulting from the integration of the torque radial density from $R_{\rm CO}$ to $R_{\rm LC}$. Since the dipolar field decreases as $r^{-3}$, it is reasonable to expect (see e.g. Rappaport et al., 2004) that the integral is dominated by the contribution close to the corotation radius and therefore independent of $R_{\rm LC}$. According to Rappaport et al. (2004) the computation of this term gives the expression

Finally, after some algebraic manipulation, $\tau_{\rm-\,RAD}$ can be expressed in terms of $\tau_{\rm-\,B}$ as

Since $R_{\rm CO}$ and $R_{\rm LC}$ are fixed for a given NS spin frequency, the term $\tau_{\rm-\,RAD}$ can be absorbed in the $\tau_{\rm-\,B}$ term as:

where $\gamma=1+6(\sin i)^{2}(R_{\rm CO}/R_{\rm LC})^{3}\sim 1$, since, typically, $R_{\rm CO}/R_{\rm LC}\ll 1$.

Using an explicit prescription for the torque radial density, as derived from a specific prescription for the magnetic field structure, the positive torque $\tau_{\rm+\,\dot{M}}+\tau_{\rm+\,B}$ can be computed as

In summary, we adopted the general prescription for the total torque exerted on the NS as derived in Kluźniak & Rappaport (2007b, Eq.  36 of the paper), modified by the introduction of the factor $\gamma=1+6(\sin i)^{2}(R_{\rm CO}/R_{\rm LC})^{3}\sim 1$ to take into account the torque from the radiative solution. Defining $\tau_{\rm TOT}=\tau_{+\,\dot{M}}+\tau_{+\,B}+\tau_{-\,B}+\tau_{\rm-\,RAD}$, we find:

which depends on the exact location of the truncation radius $R_{\rm m}$ ²²2Adopting different prescriptions for the torque radial density, as derived from a specific model for the magnetic field structure, the functional prescription for the total torque varies slightly. In particular, (Rappaport et al., 2004) and (Kluźniak & Rappaport, 2007b) computed the total torque acting on the NS for the prescription defined in equation (6). The final torque is reported in Eq. B13 (APPENDIX B; Kluźniak & Rappaport, 2007b) and in Eq. 24 (Rappaport et al., 2004). We note that small typos are probably present in the latter equation, although the limit for $R_{\rm m}\rightarrow R_{\rm CO}$, namely $\tau_{\rm TOT}=\dot{M}\sqrt{GMR_{\rm CO}}-\gamma\mu^{2}/(9R_{\rm CO}^{3})$ is the same for the two expressions, as well as for Eq. 11 adopted in this work. Indeed since we will find a posteriori that the condition $R_{\rm m}\rightarrow R_{\rm CO}$ was applicable during the entire 2018 outburst of IGR J17591$-$2342, we are confident that our results are quite independent of the detailed prescription adopted for the magnetic field structure.. Detailed computations of the truncation radius has been performed by several authors by carefully balancing the torques acting in the disc and requiring that, at the truncation radius, viscous torques vanishes (see e.g. Wang, 1995; Kluźniak & Rappaport, 2007b, and references therein). Defining the Alfvén radius as:

the location of $R_{\rm m}$ depends on the ratio $R_{\rm A}/R_{\rm CO}=\xi$ (see e.g. Wang, 1995; Kluźniak & Rappaport, 2007b) being:

For a reasonable range of Alfvén radii, namely $0.25R_{\rm CO}\leq R_{\rm A}\leq 10R_{\rm CO}$, the results of Eq. 13 can be approximated (from a visual inspection of Fig. 1 of Kluźniak & Rappaport, 2007b) with

In general, in the range $0.25R_{\rm CO}\leq R_{\rm A}\leq 10R_{\rm CO}$, Eq. 14 can be assumed as a reasonable estimate of the truncation radius for several prescriptions of the magnetic field structure. Inserting Eq. 14 in Eq. 11, we find the correct prescription for the torque acting on the NS. In particular, if, throughout the outburst of IGR J17591$-$2342, the condition $R_{\rm A}\gtrsim R_{\rm CO}$ was met, the torque is:

Conservation of angular momentum therefore implies:

where $I_{\rm NS}$ is the NS moment if inertia (for simplicity, here, we assumed that the moment of inertia of the NS does not vary in response to the accretion process, although this small effect can be taken into account if a reliable Equation of State for the NS is assumed³³3More in detail let us assume that the NS moment of inertia has the form $I_{\rm NS}=\zeta(n)(2/5)MR^{2}$, where $\zeta(n)\sim 1$ is a constant that depends, mainly, on the polytropic index $n$ adopted to describe the NS. Logarithmic differentiation and some algebraic manipulation gives for the NS angular momentum: $dL_{\rm NS}/dt=2\pi I_{\rm NS}\dot{\nu}+\dot{M}\Omega R^{2}(2\zeta(n)/5)[1+2(\dot{R}/\dot{M})(M/R)]$. Adopting $2\zeta(n)/5\sim 1$ and $(\dot{R}/\dot{M})(M/R)\sim 0$ (since in most of the Equation of State of the ultra-dense matter the radius is almost independent of the NS mass) we find $dL_{\rm NS}/dt\simeq 2\pi I_{\rm NS}\dot{\nu}+\dot{M}\Omega R^{2}$. This last term has been included in Kluźniak & Rappaport (2007b, Eq. 36 and Eq. B13 of that paper). It is worth noting that Eq. 16 is correct only for negligible variation of $I_{\rm NS}$ during the accretion process, condition that can be easily verified for this work using Eq. 7 in Sanna et al. (2017b), and $\dot{\nu}$ is the NS spin frequency derivative.

If the mass accretion rate is constant, the net torque in Eq. 16 is constant as well, positive or negative depending on the relative strengths of the acting torques. On the other hand, when $\dot{M}(t)$ varies, which is the case during an outburst, variations of the accretion rate must be taken into account when computing the spin frequency derivative. Under the assumption the mass accretion rate is well described by the bolometric accretion luminosity $L_{\rm BOL}$, we can relate the latter quantity to the observed source flux within the NICER energy range. We can then define $L_{\rm BOL}=4\pi\,d^{2}\kappa_{\rm BOL}\Phi_{\rm 0.5-10\,keV}$, where $d$ is the source distance and $\kappa_{\rm BOL}\geq 1$ is a constant, depending on the spectral shape. Finally, the bolometric accretion luminosity is related to the mass accretion rate through the efficiency $\eta_{\rm ACC}$, and therefore:

In summary, we can express the frequency derivative as:

where we used the identity $\sqrt{GMR_{\rm CO}}=\Omega R_{\rm CO}^{2}$.

In recent years, variations of the spin frequency derivative caused by the temporal evolution of the source flux $\phi(t)$ during an outburst were modelled assuming simple analytic prescriptions for the flux evolution (see e.g. IGR J00291+5934 and IGR J17511-3057 Burderi et al., 2007; Riggio et al., 2011a). More recently, a numerical approach in which the torque is rescaled at every time using the instantaneous value of the X-ray flux has been developed and applied successfully to the timing of the X-ray pulsar GRO J1744-28, a slowly spinning ($2.14\,{\rm Hz}$) NS accreting from a low-mass companion (Sanna et al., 2018d). Thanks to the excellent NICER monitoring campaign that furnished a reliable and well sampled (1 point every one or two days) $\rm 0.5-10\,keV$ flux curve of the IGR J17591$-$2342, we were able to apply, for the first time, this technique to an AMXP.

Starting from Eq. 1, we replaced the factor $\frac{1}{2}\dot{\nu}\,(t-T_{0})^{2}$ with the torque-induced phase delay term:

where $\Delta\phi_{\tau M}(t)$ and $\Delta\phi_{\tau B}(t)$ represent the contributions from the material and the threading torques, respectively. To be able to integrate the source flux in the aforementioned relation, we applied cubic spline data interpolation techniques to generate the outburst flux profile starting from that reported in Fig. 1. Using this new prescription, we then modelled simultaneously the pulse phase delays from the fundamental and first harmonic components created following the procedure described in the previous section. More specifically, we linked the model components $R_{orb}(t)$ and $\Delta\phi_{\tau M}(t)$ between the two sets of data. We iterated the method, until no significant differential corrections were found for the parameters of the model. To reduce the possible correlation between the model parameters, we limited the exploration range for the material torque parameter by defining a meaningful interval for the constants included in Eq. 17. More specifically, we limited the accretion efficiency in the range $0.75<\eta_{\rm ACC}<1$ in agreement with some recently proposed estimates (Burderi et al., 2020, submitted), we considered the source located at $d=7.6\pm 0.7$ kpc (obtained from the observation of a photo-spheric radius expansion type-I X-ray burst observed by INTEGRAL during the outburst Kuiper et al., 2020) and we considered the bolometric conversion factor with respect to the unabsorbed flux in the 0.5-10 keV energy range bounded between $2<\kappa_{\rm BOL}<3$. To determine the latter interval, we used the broad-band (0.5-80 keV) energy band spectra obtained combining the simultaneous Swift/XRT, NuSTAR and INTEGRAL IBIS/ISGRI observations of IGR J17591$-$2342 at the beginning of the outburst (Sanna et al., 2018d) to determine the ratio between the low (0.5-10 keV) and high energy (0.5-80 keV) bands obtaining a factor $\sim 3$. Moreover, we estimated the same constant, at different stages of the outburst, by using the 0.5-80 keV unabsorbed flux extrapolated from the best-fitting model of the available NICER energy spectra from which we obtained values ranging between 2 and 3. It is worth to noting that this discrepancy is likely reflecting the limited energy coverage of the NICER data, therefore, we decided to consider the largest interval (2-3) for the parameter.

The best-fit model parameters are reported in Tab. 2, while the top panel of Fig. 6 shows the fundamental (black points) and first harmonic (red points) pulse phase delays, the best-fitting model (solid lines), the spin-down (dashed lines) and spin-up (dotted line) torque components. The medium and bottom panels show the delay residuals (in phase units) with respect to the best-fitting model for the frequency fundamental and first harmonic components, respectively.

As expected, the orbital parameters obtained from this fit are fully consistent with those reported in Tab. 1, since the long baseline of the dataset guarantees to clearly disentangle orbital phase modulation from and phase variation due to the NS spin frequency changes. From the best-fitting parameter of each frequency phase component we infer the dipolar magnetic moment $\mu_{F}\simeq 4.1(4)\times 10^{26}$ G cm³ and $\mu_{1h}\simeq 4.1(5)\times 10^{26}$ G cm³, for the fundamental and the first harmonic, respectively. Adopting the Friedman-Pandharipande-Skyrme (FPS) equation of state (see e.g. Friedman & Pandharipande, 1981; Pandharipande & Ravenhall, 1989) for a 1.4 M_⊙ NS, we estimate a NS radius of $R=1.14\times 10^{6}$ cm, and a moment of inertia $I\simeq 1.45\times 10^{45}$ g cm². This corresponds to an equatorial magnetic field of $B_{eq}=\mu^{2}/R^{3}=2.8(3)\times 10^{8}\simeq 2.8(3)\times 10^{8}$ G. This value is consistent both with the estimation reported in Sanna et al. (2018d) from spin frequency equilibrium considerations for the X-ray pulsar, and the average magnetic field of known AMXPs Mukherjee et al. (2015). It is noteworthy that the value of $B_{eq}$ obtained here is likely overestimated due to a non accurate estimation of mass accretion rate in case of spin-down regime. In fact, as discussed by Kluźniak & Rappaport 2007b, in this regime the angular momentum deposited in the disc by the pulsar torques leads to a substantial additional dissipation of energy in the disc, hence, the mass accretion rate estimated from Eq. 17 results on an upper limit of the quantity. However, since it is not quite clear what fraction of this additional dissipation of energy is released in non-thermal process and what is kept inside the disc, it is complicated to quantify the discrepancy. We started approaching the problem considering the worst case scenario at which the whole rotational energy of the pulsar $I\omega\dot{\omega}$ is transferred and released by the accretion disc. With this hypothesis, we subtract the rotation luminosity from observed luminosity in Eq. 17 and we then fitted the frequency phase delays. We obtained that, taking into account the correction on the mass accretion rate, the NS magnetic field decreases roughly by 15%, reaching the value $B_{eq}=2.4(4)\times 10^{8}$ G. Both the properties of the source during the outburst and the magnetic field estimation obtained from the studies of the spin-down like behaviour of the phase delays do not seem to suggest IGR J17591$-$2342 as an exception among the AMXPs.

We investigated the properties of the pulse profile of IGR J17591$-$2342 as a function of energy (0.2-12 keV) by considering the whole available NICER dataset from its 2018 outburst. Although the overall pulse profile is well described by the superposition of four sinusoidal functions harmonically related (see Fig. 4), while applying energy selections to the events the sensitivity to harmonic components larger than two drops significantly, we therefore, independently studied the fractional amplitude and the lags of the fundamental and first harmonic only.

As shown in the top panel of Fig. 5, the fundamental pulse fractional amplitude clearly increases with energy, varying from 4% to 16% up to 5 keV, followed by a decreasing trend down to $\sim 12$% at 12 keV. The first harmonic fractional amplitude resembles the fundamental component, with a smaller variation between 2.6% and 6.6% up to 5 keV, though. Similar behaviour for the fractional amplitude has been already reported for other AMXPs such as SWIFT J1756.9-2508 (Sanna et al., 2018a), SAX J1808.4-3658 (Patruno et al., 2009a; Sanna et al., 2016), Aql X$-$1 (Casella et al., 2008), XTE J1807$-$294 (Kirsch et al., 2004) and IGR J00291+5934 (Falanga et al., 2005; Sanna et al., 2017a). Mechanisms such as a strong Comptonisation of the beamed radiation have been applied, quite successfully, to describe the properties of a few sources (Falanga & Titarchuk, 2007). However, it is noteworthy that quite different energy distributions of fractional pulse amplitude have been observed in other AMXPs such as XTE J1751$-$305 (Falanga & Titarchuk, 2007), SAX J1808.4$-$3658 (Cui et al., 1998; Falanga & Titarchuk, 2007; Sanna et al., 2017c) and IGR J17511$-$3057 (Papitto et al., 2010; Falanga et al., 2011; Riggio et al., 2011a).

Finally, time lags (Fig. 5, bottom panel) of the fundamental component suggest almost no lags between pulsed photons with energy in the range 0.5-3.5 keV, while hard lags are clearly observed for higher energies. On the other hand, the first harmonic component clearly shows hard lags with increasing energies. Although within a limited energy range, these results are consistent with those reported by Falanga & Titarchuk (2007); Falanga et al. (2011) for the AMXPs XTE J1751$-$305, XTE J17511$-$3057 and SAX J1808.4$-$3658, where the hard lags have been interpreted as the result of soft X-ray photons coming from the disc or the NS surface upscatter off hot electrons in the accretion column.