# The rotational phase dependence of magnetar bursts

The rotational phase dependence of magnetar bursts Abstract The trigger for the short bursts observed in γ-rays from many magnetar sources remains unknown. One particular open question in this context is the localization of burst emission to a singular active region or a larger area across the neutron star. While several observational studies have attempted to investigate this question by looking at the phase dependence of burst properties, results have been mixed. At the same time, it is not obvious a priori that bursts from a localized active region would actually give rise to a detectable phase dependence, taking into account issues such as geometry, relativistic effects, and intrinsic burst properties such brightness and duration. In this paper, we build a simple theoretical model to investigate the circumstances under which the latter effects could affect detectability of dependence of burst emission on rotational phase. We find that even for strongly phase-dependent emission, inferred burst properties may not show a rotational phase dependence, depending on the geometry of the system and the observer. Furthermore, the observed properties of bursts with durations short as 10–20 per cent of the spin period can vary strongly depending on the rotational phase at which the burst was emitted. We also show that detectability of a rotational phase dependence depends strongly on the minimum number of bursts observed, and find that existing burst samples may simply be too small to rule out a phase dependence. stars: magnetars, magnetic fields, gamma-rays: stars, X-rays: bursts 1 INTRODUCTION Magnetars, the most highly magnetized neutron stars (with dipole fields ≳1013 G), are isolated stars powered primarily by magnetic field decay (Thompson & Duncan 1993, 1995; Kouveliotou et al. 1998, 1999). One of their key characteristics is the sporadic emission of soft γ-ray bursts (for reviews of this specific aspect, see Woods & Thompson 2006; Turolla, Zane & Watts 2015; Kaspi & Beloborodov 2017). Durations and fluences vary, but most of these bursts are short, lasting ∼0.01–1 s, less than a typical magnetar spin period P ∼ 6 s. The slow decay of the strong magnetic field is assumed to build up stresses in the system: stress release must involve rapid reconfiguration of the external magnetic field, particle acceleration, and γ-ray emission. However, what triggers the occurrence of individual bursts, the way in which a burst progresses, and the associated emission processes, remain very poorly understood (Turolla et al. 2015). The failure point could be internal, within the crust of the star, or in the external magnetosphere itself. One question is whether the bursts are triggered within a specific active region, fixed in the rotating frame of the star. Some crust zones, for example, are expected to be particularly prone to magnetically induced faulting and yielding (see e.g. Gourgouliatos et al. 2015; Lander et al. 2015; Thompson, Yang & Ortiz 2017). Certain regions of the magnetosphere could also be more active than others, in which case there may also be a preferred height above the neutron star surface. Being able to identify whether this is the case would certainly help in efforts to determine the burst mechanism. One way to determine this, suggested by Lyutikov (2002), is to look at whether there is any rotational phase dependence to the bursts. A number of observational studies have attempted to investigate this, using various different measures such as the phase dependence of the time at which the burst peak is recorded, or the phase distribution of all of the burst photons. The evidence for phase dependence using these measures is mixed (see Section 2 for more details). What has never been done is to determine from a theoretical perspective the circumstances under which bursts from a localized active region would actually give rise to a detectable phase dependence. This will depend on geometry, gravitational light-bending, any beaming factor associated with the burst emission, the size of the burst sample for a given source, and the intrinsic burst properties (e.g. brightness, duration). Under some circumstances, bursts may well be visible throughout most of the rotational phase cycle even if they do originate from a specific active region. In this paper we address this fundamental question of the circumstances under which emission from a localized bursting region would be detectable as a rotational phase dependence according to the measures used in the literature. We then revisit the observational studies carried out to date to see what constraints they actually place on the degree to which bursting might be localized. We also consider the degree to which the rotational phase at which a burst is emitted might affect the properties measured by an observer. First, we provide an overview of magnetar burst phase-dependence studies in the literature in Section 2. Different methods have been applied to various sources in an effort to assess the (non-)phase dependence of magnetar bursts. Next in Section 3, we briefly outline the method through which we aim to answer the aforementioned questions. We choose to simulate sequences of elementary bursts of which we can control the input parameters and study any phase-dependent effects we may observe. The light-curve model is treated in Section 4 and the simulations are described in Section 5. We discuss the results of the simulations and assess the claims made in the literature in Sections 6 and 7. We find that under certain conditions the properties of the observed bursts may become significantly phase-dependent. However, we also find that for a large range of input burst parameters and configurations, a guaranteed detection of phase dependence requires many more bursts than have commonly been observed. 2 OVERVIEW OF PUBLISHED BURST PHASE-DEPENDENCE ANALYSIS Here we provide a review of previous work where, given the acquired data, the phase dependence of magnetar bursts has been evaluated. We focus on how the data were obtained and processed, and what method was used to determine the absence or presence of a phase dependence in the burst occurrences or properties. In practice, three methods have been applied: (i) searching for any significant deviations from uniformity of burst occurrence and photon arrival-time distributions against phase, (ii) searching for any correlation between the phase at which bursts occur and the pulse maxima of the (underlying) pulsed emission, and (iii) Fourier analysis on the burst occurrence times in an effort to search for significant periodicities. The latter has been applied only once; the first two are far more common. It is worth mentioning that the first two methods depend on the accuracy of the ascertained timing ephemeris. The longer the time baseline spanned by the bursts, the greater the risk of undetected time anomalies, such as glitches or spin-down deviations, that may undermine the inference of the phase. Table 1 provides a summary of the references that have carried out phase-dependence analysis of magnetar bursts. Table 1. Summary of burst phase-dependence results in the literature. Reference  Source  Dates of bursts  Satellite/  nbursts  Methodb      [dd/mm/‘yy]  Instrumenta      Gavriil et al. (2004)  1E 2259+586  18/06/‘02  RXTE/PCA  80  (i)  Savchenko et al. (2010)  SGR J1550–5428  22/01/‘09  INTEGRAL/ACS  84  (i)  Scholz & Kaspi (2011)    (22/01–30/09)/‘09  Swift/XRT  303  (i)  Lin et al. (2012)    22/01/‘09  Swift/XRT  31  (i)      30/01/‘09        Collazzi et al. (2015)    03/10/‘08–17/04/‘09  Fermi/GBM  354  (i)  Mus et al. (2015)    22/01/‘09  RXTE/PCA  4  (ii)      06/02/‘09            30/03/‘09            11/01/‘10        Gavriil et al. (2002)  1E 1048.1–5937  29/10/‘01  RXTE/PCA  2  (ii)      14/11/‘01        Gavriil et al. (2006)    29/06/‘04  RXTE/PCA  1  (ii)  Dib et al. (2009)    29/10/‘01  RXTE/PCA  4  (ii)      14/11/‘01            29/06/‘04            28/04/‘08        An et al. (2014)    (17–27)/07/‘13  NuSTAR  8  (ii)  Woods et al. (2005)  XTE J1810–197  22/07/‘03  RXTE/PCA  6  (ii)      16/02/‘04            19/04/‘04            19/05/‘04        Gavriil et al. (2011)  4U 0142+61  (06/04–26/06)/‘06  RXTE/PCA  6  (ii)  Palmer (1999)  SGR 1806–20  (10–15)/11/‘83  ICE  33  (iii)  Palmer (2002)  SGR 1900+14  -  -  -  -  Reference  Source  Dates of bursts  Satellite/  nbursts  Methodb      [dd/mm/‘yy]  Instrumenta      Gavriil et al. (2004)  1E 2259+586  18/06/‘02  RXTE/PCA  80  (i)  Savchenko et al. (2010)  SGR J1550–5428  22/01/‘09  INTEGRAL/ACS  84  (i)  Scholz & Kaspi (2011)    (22/01–30/09)/‘09  Swift/XRT  303  (i)  Lin et al. (2012)    22/01/‘09  Swift/XRT  31  (i)      30/01/‘09        Collazzi et al. (2015)    03/10/‘08–17/04/‘09  Fermi/GBM  354  (i)  Mus et al. (2015)    22/01/‘09  RXTE/PCA  4  (ii)      06/02/‘09            30/03/‘09            11/01/‘10        Gavriil et al. (2002)  1E 1048.1–5937  29/10/‘01  RXTE/PCA  2  (ii)      14/11/‘01        Gavriil et al. (2006)    29/06/‘04  RXTE/PCA  1  (ii)  Dib et al. (2009)    29/10/‘01  RXTE/PCA  4  (ii)      14/11/‘01            29/06/‘04            28/04/‘08        An et al. (2014)    (17–27)/07/‘13  NuSTAR  8  (ii)  Woods et al. (2005)  XTE J1810–197  22/07/‘03  RXTE/PCA  6  (ii)      16/02/‘04            19/04/‘04            19/05/‘04        Gavriil et al. (2011)  4U 0142+61  (06/04–26/06)/‘06  RXTE/PCA  6  (ii)  Palmer (1999)  SGR 1806–20  (10–15)/11/‘83  ICE  33  (iii)  Palmer (2002)  SGR 1900+14  -  -  -  -  aSpacecraft/instrument acronyms: Rossi X-ray Timing Explorer (RXTE), Proportional Counter Array (PCA), Nuclear Spectroscopic Telescope Array (NuSTAR), International Cometary Explorer (ICE), Anti-Coincidence Shield (ACS), X-ray Telescope (XRT), Nuclear Spectroscopic Telescope Array (NuSTAR), and Gamma-ray Burst Monitor (GBM). bThe methods are specified in Section 2. View Large The active phase of 1E 2259+586 on 2002 June 18 consisted of 80 bursts and was studied using method (i) (Gavriil, Kaspi & Woods 2004); it was claimed that the burst peak phase occurrences tended to correlate with the intensity of the pulsed emission, yet no phase dependencies were observed for the burst durations, fluences, peak fluxes, and rise/fall times. In excess of 300 bursts were observed from SGR J1550–5428 between 2008 March and 2010 January and many of those bursts were detected by multiple space-based telescopes simultaneously. Savchenko et al. (2010) found that the burst start times (the moment the burst exceeds 5σ above background) of 84 bursts, observed with the Anti-Coincidence Shield (ACS) aboard the INTEGRAL spacecraft, appear to be distributed randomly across phase, i.e. no significant departure from the mean bursts per phase bin was identified. Scholz & Kaspi (2011) and Collazzi et al. (2015) studied the burst peak times of, respectively, 303 and 354 bursts, and both found no significant (>3σ) deviations from the mean number of burst peaks per phase bin. Scholz & Kaspi (2011) however do show that the phase-folded photon times of arrival of the bursts exhibit an apparent pulse which has an offset with respect to maximum of the associated quiescent pulse profile. Lin et al. (2012) study a sample of 31 bursts and similarly find that the burst count distribution is not uniform across phase. Moreover, they find that the phase probability density anti-correlates with the phase profile of the persistent emission (with a correlation factor of −0.5 and chance probability of 3.4 × 10−2), which may suggest that the burst emission region is distinct to that of the persistent emission. Contrary to these results however, Collazzi et al. (2015) do not find a significant (>3σ) pulse shape in the epoch folded burst emission light curves. Note that the data set used by Lin et al. (2012) constitutes a subset of the data used by Collazzi et al. (2015). A total of 12 bursts from 1E 1048.1–5937 were analysed using method (ii); 4 of which were observed with RXTE/PCA between 2001 October 29 and 2008 April 28 (Gavriil, Kaspi & Woods 2002; Gavriil, Kaspi & Woods 2006; Dib, Kaspi & Gavriil 2009) and 8 of which were observed with NuSTAR in 2013 July 17–27 (An et al. 2014). It was determined that the majority of the bursts1 observed with RXTE/PCA had a probable chance alignment with pulse maxima of less than 0.01. For the latter eight bursts from the same source however there is no evidence for a preferred phase occurrence. Six bursts from XTE J1810–197 observed with RXTE/PCA were also studied with method (ii) (Woods et al. 2005). These bursts consisted of individual burst spikes, which in turn occurred near the corresponding pulse maxima of the source, either leading or trailing. A chance alignment of these spikes with the pulse maxima was estimated at roughly 0.004. Gavriil, Dib & Kaspi (2011) studied six bursts of 4U 0142+61 observed with RXTE/PCA in 2006 from April to June. They found that several bursts appear to occur near the maxima of contemporaneous folded pulse profiles (no significance criteria are specified in the reference). They argue that this may indicate that the bursts comprise extreme episodes of local transient emission sites. Palmer (1999, 2002) studied the burst properties of SGR 1806–20 and SGR 1900+14, where for the former source it was found that active bursting episodes emerge from local active regions characterized as ‘relaxation systems’. From a larger burst sample, a group of 33 bursts was identified as belonging to a single relaxation system. Method (iii), i.e. Fourier analysis on the burst occurrences of this group, revealed no apparent modulation at the rotation frequency of the NS, indicating the lack of a phase dependence.2 Here we will focus mainly on the (non-)uniform phase occurrence of burst peaks as a proxy for the phase dependency of magnetar bursts. We briefly discuss the use of alternative methods in Section 7. 3 METHODOLOGY To understand how the observed emission may depend on the rotational phase, we intentionally introduce a phase dependency by fixing the burst location to a certain region or burst patch on the magnetar surface and then set out to describe and simulate the process from emission, where we control the input parameters, to detection and characterization. Subsequently, we can study the effects of a certain configuration on the burst parameters by investigating the phase distributions of the observed burst properties. Moreover, we can establish detectability criteria for the phase dependency for certain input values/distributions and system configurations. In order to do so, we require a light-curve model that describes how the burst emission is modified depending on the location of the bursts and additional system parameters, e.g. the inclination angle to the observer and compactness of the source. The latter parameter will reshape the trajectory of emitted photons through gravitational light bending. In Section 4 we ascertain an expression that describes the fraction of rays, i.e. paths along which the emitted photons propagate, that extend out from the burst location and intersect with an observer at infinity. Subsequently, in Section 5, we simulate sequences of bursts and investigate how the burst properties are modified through the correction of the burst intensity by the aforementioned expression. 4 LIGHT CURVE MODEL In the following we adopt natural units, i.e. G = c = 1, and the spatial spherical coordinates (r, φ, θ), where φ is the polar angle to the y-axis and θ is azimuthal angle to the z-axis. Since magnetars rotate slowly (typically $$|\boldsymbol{\Omega }|\sim 10^{-1}$$ rad s−1), they can be considered to be almost spherically symmetric. Accordingly, we may assume that the metric external to the star is approximately given by the Schwarzschild space–time solution,   $$\text{d}s^2=-\mathcal {A}(r)\text{d}t^2+\mathcal {A}^{-1}(r)\text{d}r^2+r^2 \left(\text{d}\varphi ^2+\sin ^2\varphi \,\text{d}\theta ^2\right)$$ (1)where $$\mathcal {A}(r)=(1-R_S/r)$$, and RS = 2M the Schwarzschild radius, with M corresponding to the gravitational mass of the compact object. To model the effect of gravitational light bending on the burst emission, we consider the configuration illustrated in Fig. 1, which is based on the work done by Pechenick, Ftaclas & Cohen (1983). The stellar surface is located at a distance R from the origin; for neutron stars, R lies roughly in the range 2.5 − 4RS. For now we assume that the burst emission originates at3r = R over a circular patch of angular radius ψ centred at point $$\boldsymbol{p}(R, \varphi , \theta )$$ with total intensity I. Depending on the burst emission mechanism, I may be anisotropic and depends on δ ∈ [0, π), i.e. the angle between the normal vector to the stellar surface $$\hat{\boldsymbol{n}}$$ and the outgoing emission vector $$\hat{\boldsymbol{k}}_R$$. The latter lies at r = R along the associated null geodesic $$\mathcal {G}$$ of the outgoing emission, which in turn intersects with the observer at r = r0, where $$\hat{\boldsymbol{k}}_R\rightarrow \hat{\boldsymbol{k}}$$. We presume that the region r < R is opaque and r > R is entirely transparent. Moreover, due to the comparatively long rotation period of magnetars, we may neglect certain corrections, such as oblateness of the stellar surface, light travel-time delays, and Doppler effects, which become significant for NSs with $$|\boldsymbol{\Omega }|\gg 1$$ rad s−1 (e.g. Morsink et al. 2007). Figure 1. View largeDownload slide Schematic representation of null geodesic $$\mathcal {G}$$ connecting the burst patch, centred at point $$\boldsymbol{p}(\theta _0)$$ of angular width ψ, with the observer. The observer is set in the +z-direction, where the z-axis makes an angle χ with the rotational axis of the neutron star (NS) $$\boldsymbol{\Omega }$$. The normal vector to the NS surface $$\hat{\boldsymbol{n}}$$ at $$\boldsymbol{p}$$ is at an angle α with $$\boldsymbol{\Omega }$$ and at an angle θ0 with the observer's line of sight, where the latter depends on the rotational phase of the NS, i.e. θ0 = θ0(ϕ) as prescribed by equation (16). In a stationary frame at r = R, photons radiate from the burst location (along $$\hat{\boldsymbol{k}}_R$$) at and angle δ to the normal. The intensity of the source may be isotropic or beamed I = I(δ), depending on the underlying emission mechanism and properties of the emitting region. Figure 1. View largeDownload slide Schematic representation of null geodesic $$\mathcal {G}$$ connecting the burst patch, centred at point $$\boldsymbol{p}(\theta _0)$$ of angular width ψ, with the observer. The observer is set in the +z-direction, where the z-axis makes an angle χ with the rotational axis of the neutron star (NS) $$\boldsymbol{\Omega }$$. The normal vector to the NS surface $$\hat{\boldsymbol{n}}$$ at $$\boldsymbol{p}$$ is at an angle α with $$\boldsymbol{\Omega }$$ and at an angle θ0 with the observer's line of sight, where the latter depends on the rotational phase of the NS, i.e. θ0 = θ0(ϕ) as prescribed by equation (16). In a stationary frame at r = R, photons radiate from the burst location (along $$\hat{\boldsymbol{k}}_R$$) at and angle δ to the normal. The intensity of the source may be isotropic or beamed I = I(δ), depending on the underlying emission mechanism and properties of the emitting region. The Schwarzschild solution admits four Killing vectors associated with conserved quantities that emerge from symmetries inherent in the solution for which,   $$K_\mu \dot{x}^\mu =\text{constant},$$ (2)where $$\dot{x}^{\mu }=\text{d}x^{\mu }/\text{d}\lambda$$, with λ is some affine parameter. Considering the orbital motion of a photon with tangent 4-vector $$V^\mu =\dot{x}^{\mu }$$ in the equatorial plane, i.e. φ = π/2, the Schwarzschild solution admits two Killing vectors   \begin{eqnarray} \epsilon _\mu =(\mathrm{\partial} _t)_\mu = (-\mathcal {A}(r),0,0,0), \end{eqnarray} (3)  \begin{eqnarray} J_\mu =(\mathrm{\partial} _\varphi )_\mu = (0,0,0,r^2), \end{eqnarray} (4)associated, respectively, with conservation of energy and the magnitude of angular momentum. Accordingly, we may define   \begin{eqnarray} \epsilon _\mu \dot{x}^\mu =-\mathcal {A}(r)V^t\equiv -1, \end{eqnarray} (5)  \begin{eqnarray} J_\mu \dot{x}^\mu =r^2V^\theta \equiv b, \end{eqnarray} (6)where b ≥ 0 denotes the impact parameter of the photon trajectory, i.e. the null geodesic $$\mathcal {G}$$. Since $$\dot{x}^\theta =0$$ and $$g_{\mu \nu } \dot{x}^\mu \dot{x}^\nu =0$$ for massless particles, we find that the tangent 4-vector of an outgoing photon is given by,   $$V^{\mu }=\left(\mathcal {A}^{-1}(r), \sqrt{1-\frac{b^2}{r^2}\mathcal {A}(r)},0, \frac{b}{r^2}\right).$$ (7)Setting b = 0, we obtain the tangent 4-vector of a radially outgoing photon,   $$W^{\mu }=(\mathcal {A}^{-1}(r),1,0,0).$$ (8)A stationary observer with 4-velocity   $$U^{\mu }=(\mathcal {A}^{-1/2}(r),0,0,0),$$ (9)will observe an angle   $$\cos \xi =\sqrt{1-\frac{b^2}{r^2}\mathcal {A}(r)},$$ (10)between the photons prescribed by Vμ and Wμ (Pechenick et al. 1983). Note that $$\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{k}}_R=\cos \delta = \cos \,(\xi |_{r=R})$$, such that we may write   $$\delta (b)=\arcsin \left(\frac{b}{b_{\rm max}}\right), \ \ \ \ \text{with}\ \ b_{\rm max}=\frac{R}{\mathcal {A}^{1/2}(r=R)}.$$ (11)The total angular deflection of $$\mathcal {G}$$, which determines the ‘bending’ of the photon trajectory from the surface patch to the observer, is given by   $$\theta _*(b) =\int _R^{r_0}\frac{b}{r^2}\left[1-\frac{b^2}{r^2}\mathcal {A}(r)\right]^{-1/2}\text{d}r.$$ (12)Incidentally, the total coordinate light travel time along $$\mathcal {G}$$ is given by   $$T_*(b) =\int _R^{r_0}\mathcal {A}^{-1}(r)\left[1-\frac{b^2}{r^2}\mathcal {A}(r)\right]^{-1/2}\text{d}r.$$ (13) The difference in travel time between radially emitted photons (with b = 0) and those with an arbitrary impact parameter can be estimated accordingly,   $$\Delta t_*(b) =\int _R^{r_0}\mathcal {A}^{-1}(r)\left\lbrace \left[1-\frac{b^2}{r^2}\mathcal {A}(r)\right]^{-1/2}-1\right\rbrace \text{d}r.$$ (14)The maximum travel time delay then for a typical NS with R = 2.5 RS (R = 106 cm, M = 1.5 M⊙), is Δt*(bmax) ≃ 6.7 × 10−2 ms ≪ P ∼ 6 s. To an observer in the +z-direction, the system is axisymmetric around the z-axis, such that the location of the burst patch $$\boldsymbol{p}$$ can be uniquely described by θ0. The angle between the observer's line of sight and the rotation axis of the neutron star $$\boldsymbol{\Omega }$$ is denoted by χ. Furthermore, the angle between the location of the burst patch and $$\boldsymbol{\Omega }$$ is given by α. Accordingly, depending on the rotational phase of the neutron star,   $$\phi (t)=\frac{2\pi t}{P},$$ (15)the angle between the observer and burst patch θ0 is given by the relation   $$\cos [\theta _0(\phi )] = \cos \chi \cos \alpha + \sin \chi \sin \alpha \cos \phi .$$ (16) The observed brightness is the integral of the intensity at the observer,   $$\text{d}I_{\rm obs}=I(r_0, \Omega ^{\prime })\text{d}\Omega ^{\prime }=\mathcal {A}^2(r=R)I(R, \Omega ^{\prime })\text{d}\Omega ^{\prime },$$ (17)over the solid angle subtended in the observer's sky,   $$\text{d}\Omega ^{\prime }=\sin \theta ^{\prime } \text{d}\theta ^{\prime } \text{d}\varphi ^{\prime }\simeq \theta ^{\prime } \text{d}\theta ^{\prime } \text{d}\varphi ^{\prime }.$$ (18)Due to the axisymmetry dφ΄ = dφ ≡ Φ(θ*, θ0), where we define the polar differential distance as the function Φ(θ*, θ0), which under the conditions θ0 + ψ ≤ θ* ≤ π and θ0 − ψ ≥ 0 is given by the following expression:   \begin{eqnarray} \Phi (\theta _*, \theta _0) = \left\lbrace \begin{array}{ll}2\arccos \left(\frac{\cos \psi -\cos \theta _0\cos \theta _*}{\sin \theta _0\sin \theta _*}\right) & \text{if } \theta _0-\psi \le \theta _* \\ &\le \theta _0 +\,\psi ,\\ 0 & \text{otherwise.} \end{array} \right. \end{eqnarray} (19)Consequently, together with θ΄ = ξ and equation (10) evaluated at r0 → ∞, we obtain   $$d\Omega ^{\prime }\simeq \Phi [\theta _*(b), \theta _0]\, \xi \text{d}\xi \simeq \frac{\Phi (b, \theta _0)}{r_0^2}b\,\text{d}b.$$ (20)We write the brightness of the source as   $$I(R, \Omega ^{\prime })=I_0f[\delta (b)],$$ (21)with the beaming functions given by f(b). Currently, we do not have a physical model for the shape of the beaming function, which will most likely depend on the radiative transfer properties of the local magnetic field. Accordingly, the influence of the magnetic field on the light trajectories is ignored for now. An example of a more realistic model was considered by van Putten et al. (2016) in the case of fireball beaming. Here, for descriptive purposes, we consider a Gaussian shape for the beaming function,   $$f[\delta (b)]= \left\lbrace \begin{array}{ll}1 & \text{isotropic}, \\ \sqrt{\frac{\pi }{2\sigma _{\rm b}^2}}\ {\rm erf}\left(\frac{\pi }{2\sqrt{2}\sigma _{\rm b}}\right)^{-1}\exp \left(-\frac{\delta (b)^2}{2\sigma _{\rm b}^2}\right)& \text{beamed}, \end{array} \right.$$ (22)where σb parameterizes the beam width. We neglect rotational aberration of light effects, since the star rotates slowly. Finally, we find the expression for the observed intensity of the source   $$I_{\rm obs}(\theta _0)=I_0\left(\frac{R}{r_0}\right)^2\kappa (\theta _0),$$ (23)where we define   $$\kappa (\theta _0)\equiv \left(\frac{R^{1/2}}{b_{\rm max}}\right)^4\int ^{b_{\rm max}}_0f[\delta (b)]\, \Phi (b, \theta _0)\,b\,\text{d}b.$$ (24) Fig. 2 shows the observed burst emission Iobs as a function of burst patch location θ0 for a compact object with R = 2.5 RS. In this case the size of the patch is ψ = 1°. We consider the observed intensity for the case of isotropic and beamed emission. Note that the location of the terminator lies ‘behind’ the star, i.e. beyond θ0 = π/2, at θ0 ≃ 0.72π. At this angle, only photons with an impact parameter of bmax ≃ 3.23RS reach the observer. Note that the beamed emission is more prominent at θ0 = 0 and drops off faster than the isotropic emission with increasing θ0. For comparison, we plotted the emission profiles of sources with R = 1.6 RS and R ≫ RS, where the former is close to the most extreme case, i.e. R > 1.5 RS and the latter approximates flat spacetime. Figure 2. View largeDownload slide Burst emission as a function of the angle between the observer and the burst location, parametrized by θ0 (see Fig. 1) for a burst spot area of ψ = 1°. The solid (dashed) curves denote emission from an isotropic (a beamed) source. A beam width of σb = π/6 was used. The black curves represent emission from a compact object with R = 1.6 RS (close to the most extreme case, i.e. R > 1.5 RS). The magenta curves denote the emission from a typical neutron star R = 2.5 RS (R = 106 cm, M = 1.5 M⊙). The blue curves denote the emission from a non-relativistic source R ≫ RS, i.e. if we neglect the GR light bending effects. The dotted lines indicate the location of the terminator, respectively, at θ0 = 0.5π and θ0 ∼ 0.72π for R ≫ RS and R = 2.5 RS, beyond which the burst patch is invisible. Note that emission is always visible in the R = 1.6 RS case, regardless the location of the burst, i.e. ∀θ0. Moreover, in this extreme case the intensity of the burst patch is greatly amplified at θ0 → π due to gravitational lensing effects. This plot is similar to fig. 4 in Pechenick et al. (1983). Figure 2. View largeDownload slide Burst emission as a function of the angle between the observer and the burst location, parametrized by θ0 (see Fig. 1) for a burst spot area of ψ = 1°. The solid (dashed) curves denote emission from an isotropic (a beamed) source. A beam width of σb = π/6 was used. The black curves represent emission from a compact object with R = 1.6 RS (close to the most extreme case, i.e. R > 1.5 RS). The magenta curves denote the emission from a typical neutron star R = 2.5 RS (R = 106 cm, M = 1.5 M⊙). The blue curves denote the emission from a non-relativistic source R ≫ RS, i.e. if we neglect the GR light bending effects. The dotted lines indicate the location of the terminator, respectively, at θ0 = 0.5π and θ0 ∼ 0.72π for R ≫ RS and R = 2.5 RS, beyond which the burst patch is invisible. Note that emission is always visible in the R = 1.6 RS case, regardless the location of the burst, i.e. ∀θ0. Moreover, in this extreme case the intensity of the burst patch is greatly amplified at θ0 → π due to gravitational lensing effects. This plot is similar to fig. 4 in Pechenick et al. (1983). In the simulations we concentrate on the relative changes in intensity between the input and observed burst. Accordingly, from equations (16) and (24), we define   $$\kappa _*(\phi )\equiv \frac{\kappa [\theta (\phi )]}{\kappa _{\rm max}},$$ (25)which depends on the angles χ, α, and the phase of the neutron star, and describes the fraction of rays that intersect with the observer at infinity, from the entire ray-bundle that extends outwards from the burst patch. In the following section, we use this expression as our measure for how the burst intensity is modulated. Varying χ separately from α, or vice versa, acts as a multiplicative factor to the absolute intensity. Since, we only consider the fractional intensity, we may explore the parameter space of these angles by setting χ = α. Fig. 3 illustrates the shape of κ*(ϕ) for R = 2.5 RS in four different angle configurations, in the case of both isotropic and beamed emissions. Figure 3. View largeDownload slide Fraction of rays that intersect with the observer at infinity, from the entire ray-bundle that extends outwards from the burst patch, i.e. κ*(ϕ). The solid (dashed) curves represent isotropic (beamed, with σb = π/6) emission from a source with R = 2.5 RS. The colours represent four different geometries, given by χ and α. Figure 3. View largeDownload slide Fraction of rays that intersect with the observer at infinity, from the entire ray-bundle that extends outwards from the burst patch, i.e. κ*(ϕ). The solid (dashed) curves represent isotropic (beamed, with σb = π/6) emission from a source with R = 2.5 RS. The colours represent four different geometries, given by χ and α. 5 SIMULATIONS Per simulation run we produce a sequence of n bursts where we control the input parameters of the bursts and the system. Nonetheless, we treat the detection of individual photons, which are entirely described by their times-of-arrival (TOA), in a probabilistic fashion. We assume for simplicity that a single magnetar burst can be modelled with a exponential rise/exponential decay profile that allows for asymmetry, i.e. the rise time and fall time may be different. When simulating photons emitted at the source, this corresponds to drawing N random photon times-of-emission (TOE) from a skewed Laplace distribution,   $$p_{\rm TOE}(t)=\frac{1}{(1+s)\tau }\left\lbrace \begin{array}{ll}\exp \left[(t-t_0)/\tau \right] & \text{if }\ t < t_0, \\ \exp \left[-(t-t_0)/(s\tau )\right] & \text{if }\ t \ge t_0, \end{array} \right.$$ (26)where τ denotes the exponential rise time of the burst, (sτ) the decay time, with s the skewness factor, and t0 the peak time of the burst. The profile of a simulated burst is shown in Fig. 4. We deliberately adopt an oversimplified burst profile to better understand the differences between the input and output data. Huppenkothen et al. (2015) decompose complex magnetar bursts, observed from SGR J1550-5418, into several spike-like components, which in turn are modelled with a similar profile as in equation (26). In this paper, we assume that a burst can be represented as a single spike, as a simple model that lets us explore the relevant effects. Note that we define the duration of the burst, T90, as the time it takes for the fluence to increase from 0.05 to 0.95 of the total burst fluence. We fix the compactness of the source to R = 2.5 RS, corresponding to a typical neutron star with R = 106 cm and M = 1.5M⊙, and set the rotation period to P = 6 s. We choose a light-curve bin width, δt, and background count level, b, such that the background count rate approximates that from Fermi/GBM data (Huppenkothen et al. 2015), i.e. ζGBM ∼ 318 counts s−1. A list of the simulation parameters is given in Table A1. Figure 4. View largeDownload slide Profile of an input burst, i.e. the photon emission rate at the burst patch, where t0 represents the time at which the burst peaks, A denotes the burst amplitude, τ denotes the rise time, τs parameterizes the decay time, determined by s the skewness parameter for a given τ. For s > 1 (s < 1) the burst rises faster (slower) than it decays. The burst duration is given by T90, which is defined as the time the fluence increases from 0.05 to 0.95 of the total burst fluence; the interval T90 contains 0.90N photons. Figure 4. View largeDownload slide Profile of an input burst, i.e. the photon emission rate at the burst patch, where t0 represents the time at which the burst peaks, A denotes the burst amplitude, τ denotes the rise time, τs parameterizes the decay time, determined by s the skewness parameter for a given τ. For s > 1 (s < 1) the burst rises faster (slower) than it decays. The burst duration is given by T90, which is defined as the time the fluence increases from 0.05 to 0.95 of the total burst fluence; the interval T90 contains 0.90N photons. 5.1 General simulation procedure Here we proceed to describe in more detail the general form of a simulation run step by step: We decide on the number of bursts n we wish to produce per simulation run and set the inclination angle of the source χ. Next, we assign values to the burst parameters ψ, α, N, t0, τ, s, and T90, where the latter three parameters cannot be defined independently of each other. In the following simulation runs, described in Sections 5.2, 5.3, and 5.4, we only consider symmetric input bursts, s = 1, with a size of ψ = 1°, and draw their peak time from a uniform distribution, i.e. t0 ∼ Uniform(0, P), where P is the rotation period of the magnetar which we set to P = 6. We generate a single burst by drawing N photon emission times (TOEs) from pTOE(t): We draw random numbers from a uniform distribution Uniform(0, 1) and transform these values to follow the required skewed Laplace distribution by using the inverse cumulative distribution function of the latter, i.e. the percent point function,   \begin{eqnarray} {\rm TOE}(x)=\left\lbrace \begin{array}{ll}t_0 + \tau \ln \left[\left(1+s\right)x\right] & \text{if }\ x < (1+s)^{-1}, \\ t_0 - s\tau \ln \left[\left(1+s^{-1}\right)\left(1-x\right)\right] & \text{if }\ x \ge (1+s)^{-1}. \end{array} \right. \nonumber\\ \end{eqnarray} (27) Using equation (15), we determine the phase of each TOEi, i.e. ϕi. Whether an emitted photon reaches the observer depends on whether the ray, along which the photon propagates, intersects with the detector, given by κ*(ϕi), which denotes the fraction of photons directed into our line of sight. In order to decide whether a given photon intersects with the detector, we use rejection sampling: for each TOEi we generate a latent variable z drawn from p(z) = Uniform(0, 1). We only keep the TOEi if z < κ*(ϕi).4 The TOEs that we save are detected by the observer. A detected photon is recorded as a count with a corresponding TOAi, where the TOAi = TOEi of the respective photon, since we consider a perfect detector and may neglect the distance to the source and gravitational time-delay effects (see Section 4). Furthermore, we add background counts or TOAs uniformly to the detected TOA data from the burst (with length Ndet ≤ N), such that the mean background count rate becomes approximately ζ. We bin the total TOA data in $$\mathcal {N}_{\rm bins}$$ time bins of length δt, whereby the counts in each bin follow Poisson statistics. We proceed by applying a similar burst identification algorithm as used by Gavriil et al. (2004), assuming that we can infer the background count rate to be ζ, yet have no prior knowledge of the burst and system input parameters. The probability of the number of counts ki in the ith bin occurring is given by the Poisson distribution,   $$P_i=\frac{\mu ^{k_i}\text{e}^{-\mu }}{k_i!},$$ (28)where μ represents the mean count level, which in our case will be b = ζ δt. Bins for which   $$P_i \le 3\times 10^{-3}\mathcal {N}_{\rm bins}^{-1}$$ (29)are recorded as significant departures from the mean, where we have corrected for the number of trails by dividing by $$\mathcal {N}_{\rm bins}$$, i.e. the total number of bins searched over. From these, the time bin containing the maximum departure ysig is labeled as $$t_0^{\rm sig}$$. The burst edges, labeled as tin and tout, are ascertained by making use of a running mean, i.e. when the mean count level of an interval μ* of ΔTinterval = 0.25 s, moving outwards in steps of δt on both sides of $$t_0^{\rm sig}$$, falls below b* = 1.1b, the burst edges are then given by the centre time of the respective intervals, tin before and tout after the burst. The duration of this interval is denoted ΔT = |tin − tout|. We fit the light curve of any identified bursts with the following burst model,   $$m(t) = N p_{\rm TOE}(t\,|\,t_0, \tau ,s) + \zeta ,$$ (30)using an L-BFGS-B constrained optimizer (Zhu et al. 1997) to determine the maximum (Poisson) likelihood, whereby we fix the background parameter ζ and provide initial guesses for the remaining parameters:   \begin{eqnarray} t_0^{\rm init} = t_0^{\rm sig}, \end{eqnarray} (31)  \begin{eqnarray} \tau ^{\rm init} = \left(t_0^{\rm sig} - t_{\rm in}\right)\left[\ln \left(\frac{y^{\rm sig}}{b^*}\right)\right]^{-1}, \end{eqnarray} (32)  \begin{eqnarray} s^{\rm init} = 1, \end{eqnarray} (33)and Ninit is defined as the number of counts in the interval ΔT minus the background counts,5 i.e. ζΔT. After the fit, we delete the bins in ΔT from the light curve, and repeat steps (V) and (VI) until no significant departures from the mean, i.e. μ = b, are recorded. We return to step (II) until we have generated the pre-defined number of bursts n. Note that the number of observed bursts might be different, since bursts might go undetected, or be interpreted as multiple separate bursts. Moreover, some identified ‘bursts’ may simply be significant statistical deviations from the background level. However, according to the condition stated in equation (29), we only expect this to be the case in ∼0.3 per cent of the input bursts. 5.2 Run 1: initial simulation run We start with the simplest scenario, where we consider simulations of sequences of identical bursts, referred to as Run 1. Per simulation we fix the values for χ, α, ψ, s, and T90. The latter two parameters determine the value of τ. Subsequently, we define N using the condition that the input burst amplitude A is 104 photons s−1. We run simulations for three separate burst durations (see Fig. 5) and vary the angles χ and α, to study their effects on the observed quantities. Figure 5. View largeDownload slide Folded profiles of input bursts of the initial simulation run. We consider symmetric bursts of three separate durations. Here the rise-time τ is determined by the values for s and T90. We choose N such that the burst amplitude is A = 104 photons s−1 at t0. Figure 5. View largeDownload slide Folded profiles of input bursts of the initial simulation run. We consider symmetric bursts of three separate durations. Here the rise-time τ is determined by the values for s and T90. We choose N such that the burst amplitude is A = 104 photons s−1 at t0. We concentrate on the difference in input and observed best-fitting value for the time of the burst peak (respectively, t0 and $$t_0^{\rm bf}$$), where the difference is parametrized as $$\Delta t_0\equiv t_0^{\rm bf}-t_0$$, the rise-time τ, skewness factor s, and burst duration T90. The input values of the latter three parameters are denoted as τ0, s0, and T90,0. All input parameters of Run 1 are listed in Table 2. Table 2. Input parameters for Run 1, consisting of 12 separate simulations. We consider a constant input burst profile with peak times distributed uniformly across phase. Per simulation we vary the burst duration T90, and angles χ, α, where we set χ = α (see Section 4). Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  0.15, 1.5, 3.0  δt (s)  200−1  Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  0.15, 1.5, 3.0  δt (s)  200−1  View Large 5.3 Run 2: T90-distribution Next we perform a simulation run, Run 2, where in step (I) of the general simulation procedure (Section 5.1), we draw the burst duration T90 for each individual burst from a lognormal distribution centred at $$\overline{T}_{90}=0.1$$ s, with a width of $$\sigma _{T_{90}}=1$$ (e.g. Göğüş et al. 2001), and lower and upper burst duration cut-off at, respectively, $$T_{90}^{\rm min}=300^{-1}$$ s and $$T_{90}^{\rm max}=3$$ s < P. Fixing the input burst amplitude at A = 104 photons s−1, as done in Run 1, we find that the rise time of the shortest admissible burst duration, i.e. 300−1 s, is τmin ∼ 7.2 × 10−4 s. Accordingly, we set δt = 1400−1 s for this simulation run. The input parameters are summarized in Table 3 and the results are presented in Section 6.2.2. Table 3. Input parameters for Run 2, consisting of four separate simulations. The input burst durations T90 are drawn from a lognormal distribution with lower and upper cut-off, respectively, at $$T^{\rm min}_{90}=300^{-1}$$ s and $$T^{\rm max}_{90}=3$$ s. Per simulation run we vary the angles χ, α, where we set χ = α. Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  $${\sim }{\rm LogNormal(\overline{{\it T}}_{90}, \sigma ^2_{{\it T}_{90}})}$$  $$\overline{T}_{90}$$ (s)  0.1  $$\sigma _{T_{90}}$$ (s)  1  δt (s)  1400−1  Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  $${\sim }{\rm LogNormal(\overline{{\it T}}_{90}, \sigma ^2_{{\it T}_{90}})}$$  $$\overline{T}_{90}$$ (s)  0.1  $$\sigma _{T_{90}}$$ (s)  1  δt (s)  1400−1  View Large 5.4 Run 3: burst amplitude distribution In this simulation run we fix the burst duration to T90 = 1 s and draw an amplitude A for each individual burst from a power-law distribution,   $$\frac{\text{d}n}{\text{d}A}\propto A^{-\Gamma }, \ \ \ \text{with}\ \ \ A^{\rm min} < A < A^{\rm max},$$ (34)where Amin and Amax represent the limits of the distribution, and Γ denotes the power-law index. Note that the number of emitted photons during the burst N are linearly proportional to A, such that these are distributed in a similar fashion. In accordance with the observation of the energy distributions of magnetar bursts, we choose Γ = 5/3 (e.g. Cheng et al. 1996). The input parameters are summarized in Table 4 and the results are presented in Section 6.2.3. Table 4. Input parameters for Run 3, consisting of four separate simulations. The input burst amplitudes A are drawn from a power-law distribution [equation (34)], with Amin = 5 × 102 photons s−1, and Amax = 106 photons s−1. Per simulation run we vary the angles χ, α, where we set χ = α. Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  ∼Powerlaw(Γ)  Γ  5/3  T90 (s)  1  δt (s)  200−1  Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  ∼Powerlaw(Γ)  Γ  5/3  T90 (s)  1  δt (s)  200−1  View Large 6 RESULTS 6.1 Predictions for run 1 To better understand the results of the simulations, we first examine how κ*(ϕ) will affect the burst parameters. In Fig. 6, we plot the predicted phase distributions of the burst parameters (rows) for the three separate burst durations (columns) of Run 1. These curves were obtained by fitting the burst model to the theoretical light curve that results when the input model (see Fig. 5) is modulated, for a given phase, by the appropriate κ*(ϕ) but without taking into account any photon noise or detectability effects (which are treated properly in the full simulations). It gives an idea of the general trends expected, but no idea of the scatter. Furthermore, we only fit modulated burst profiles with a peak rate of ≳600 counts s−1, since ones with lower peak rates will likely go undetected in the simulations. Figure 6. View largeDownload slide Predicted phase distributions of burst parameters (top to bottom) for the three bursts (left to right) studied in the initial simulation run, Run1. From top to bottom: best-fitting burst counts, time difference between burst input t0 and the best-fitting value (Δt0), best-fitting burst rise time τ, best-fitting skewness factor s, and burst duration T90, inferred from the latter τ and s (the input parameters are given by N0, τ0, s0, and T90,0). The distinct colours represent different values for the angles, χ and α. These curves were obtained by fitting the burst model to theoretical burst profiles (Fig. 5) that are modulated by κ* (equation 25; Fig. 3), depending on their phase occurrence. We only proceed to fit the modulated profile if the peak rate is ≳600 counts s−1. Note that for longer burst durations and larger angles, the best-fitting parameters deviate more from their input values. Figure 6. View largeDownload slide Predicted phase distributions of burst parameters (top to bottom) for the three bursts (left to right) studied in the initial simulation run, Run1. From top to bottom: best-fitting burst counts, time difference between burst input t0 and the best-fitting value (Δt0), best-fitting burst rise time τ, best-fitting skewness factor s, and burst duration T90, inferred from the latter τ and s (the input parameters are given by N0, τ0, s0, and T90,0). The distinct colours represent different values for the angles, χ and α. These curves were obtained by fitting the burst model to theoretical burst profiles (Fig. 5) that are modulated by κ* (equation 25; Fig. 3), depending on their phase occurrence. We only proceed to fit the modulated profile if the peak rate is ≳600 counts s−1. Note that for longer burst durations and larger angles, the best-fitting parameters deviate more from their input values. Based on the predicted curves, we expect that the parameter distributions that we obtain from Run 1 will deviate from their input parameters more strongly for longer burst durations T90 and larger angles χ and α. Approaching ϕ = π from below (above), we find that the bursts will appear to occur earlier (later), rise slower (faster), to become more skewed, and last longer, than their input counterparts. Note furthermore that, in contrast to the predicted phase distributions of N, Δt0, s, and T90, the phase distribution of τ is neither symmetric nor perfectly anti-symmetric about ϕ = π. The results of Run 1 are presented in Section 6.2.1. 6.2 Burst properties from simulations 6.2.1 Run 1: initial simulation run In Figs 7, 8, and 9, we plot the phase distributions (left) and parameter densities (right) of the obtained bursts parameters for three separate input burst durations T90,0, respectively, 0.15 s, 1.5 s, and 3.0 s. Table 5 lists the amount of bursts that were identified per configuration. Figure 7. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 0.15 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 7. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 0.15 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 8. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 1.5 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 8. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 1.5 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 9. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 3.0 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 9. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 3.0 s. The theoretical predictions for these distributions are shown in Fig. 6. Table 5. Number of identified bursts nid for Run 1, per configuration. The input number of bursts for each simulation run was n = 104. We expect that ∼30 of the identified ‘bursts’ simply constitute statistical deviations that exceed the burst identification threshold (given by equation 29). T90 (s)  χ = α (°)  nid  0.15  30  10 034    45  10 022    60  8875    90  6592  1.50  30  10 012    45  10 017    60  9940    30  7438  3.00  30  10 017    45  10 012    60  10 007    90  9553  T90 (s)  χ = α (°)  nid  0.15  30  10 034    45  10 022    60  8875    90  6592  1.50  30  10 012    45  10 017    60  9940    30  7438  3.00  30  10 017    45  10 012    60  10 007    90  9553  View Large As predicted, we find that especially for longer duration bursts and larger angles, the phase dependence of the burst parameters becomes more pronounced. Evidently this is much less the case for bursts with T90 ≪ P – the parameter densities remain strongly peaked around their input values (e.g. Fig. 7). Nevertheless, in those cases around ϕ ∼ π the bursts still go undetected for large values of χ and α, because either no rays extending from the burst patch intersect with the detector during the burst (i.e. the bursts are invisible) or they do not significantly stand out from the background level. The results confirm that when approaching ϕ = π from below (above), the bursts will appear to occur earlier (later), rise slower (faster), and last longer than their input counterparts. Moreover, the predicted asymmetric profile of the rise-time phase distribution (most notably in the parameter densities of Figs 8 and 9) is clearly observed. Looking at the parameter densities, it appears that the rise times of bursts going out of view are more spread out, yet those of bursts coming into view are more clustered. This is in accordance with the predictions; the initial slope of the rise-time phase distribution (from ∼0 − 4π/5) is steeper compared to the final slope (from ∼6π/5 − 2π). The predicted values between ∼4π/5 − 6π/5 are produced less well in the simulations, since the amount of detected photons is minimal around ϕ = π, complicating burst identification and characterization. We find that both in the predictions and, even more so, in the simulations that the majority of observed bursts have τ/τ0 < 1. Since, T90 ∝ τ we also find for most observed bursts that T90/T90,0 < 1, i.e. the bursts seem to last shorter than their input counterparts. In general, the observed scatter is likely due to photon noise effects, which become most significant near ϕ = π. These effects influence the efficacy of the burst-identification algorithm and the observed burst morphology. 6.2.2 Run 2: T90 distribution The results of Run 2, where we draw the burst duration for each individual burst from a lognormal distribution, are shown in Fig. 10. Table 6 lists the number of identified bursts per configuration. The results closely resemble those of Run 1, with T90 = 0.15 s (see Fig. 7). We only find weak phase dependencies that only become noticeable for large values of the angles, χ and α. In Fig. 11, the input (dotted histogram) and best-fitting (solid histogram) T90 distributions are shown for the four separate configurations. Notice the small dearth of short duration bursts in each histogram; short duration bursts contain fewer counts and may therefore be missed by the burst-identification algorithm (step (V) in Section 5.1). We furthermore find that there is a slight excess at T90 ∼ 0.6 s (although not apparent when χ = α = 90°), which is due to the fact that for most observed bursts τ/τ0 < 1 and T90 ∝ τ. Figure 10. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 2, where the input burst durations are drawn from a lognormal distribution (see Fig. 11). Figure 10. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 2, where the input burst durations are drawn from a lognormal distribution (see Fig. 11). Figure 11. View largeDownload slide Burst duration T90 distributions of four separate simulations (with different values for χ and α) of Run 2. The dotted (solid) histograms represent the input (observed) burst duration distributions. The slight dearth of observed short duration bursts, present in each simulation, is simply due to the fact that they consist of fewer counts and are therefore less likely to be identified by the burst-identification algorithm. Figure 11. View largeDownload slide Burst duration T90 distributions of four separate simulations (with different values for χ and α) of Run 2. The dotted (solid) histograms represent the input (observed) burst duration distributions. The slight dearth of observed short duration bursts, present in each simulation, is simply due to the fact that they consist of fewer counts and are therefore less likely to be identified by the burst-identification algorithm. Table 6. Number of identified bursts nid for Run 2, per configuration. The input number of bursts for each configuration was n = 104. χ = α (°)  nid  30  10 028  45  9871  60  7382  90  5971  χ = α (°)  nid  30  10 028  45  9871  60  7382  90  5971  View Large 6.2.3 Run 3: burst amplitude distribution The results of Run 3, where we draw the burst amplitude/number of emitted burst photons for each individual burst from a power-law distribution, are presented in Fig. 12. Table 7 lists the number of identified bursts per configuration. We find much less spread in the phase distributions of the parameters, compared to e.g. the results from the 1.5 s burst in Run 1 (Fig. 8), however we do observe a considerable amount of scatter. The latter is likely due to the fact that the majority of input bursts (∼0.87) are low-amplitude bursts, i.e. A ≲ 104 photons s−1, which are more difficult to characterize, i.e. their morphology is relatively heavily affected by Poisson noise. Fig. 13 displays the input and observed burst amplitude distributions. Despite a slight offset at larger angles, we find that the slope of the distributions is reproduced by the observed bursts. Input bursts with an amplitude ≲103 photons s−1 may go unidentified as they will likely fall below the significance threshold of the burst identification algorithm. Figure 12. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 3, where the input burst amplitudes are drawn from a power-law distribution (see Fig. 13). Figure 12. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 3, where the input burst amplitudes are drawn from a power-law distribution (see Fig. 13). Figure 13. View largeDownload slide Burst amplitude distributions of four separate simulations (with different values for χ and α) of Run 3. The dotted (solid) histograms represent the input (observed) burst amplitude distributions. The cut-off of observed low-amplitude bursts (at ≲103 photons s−1) is due to the fact that the amplitude of these bursts likely occurs below the threshold of the burst-identification algorithm. Figure 13. View largeDownload slide Burst amplitude distributions of four separate simulations (with different values for χ and α) of Run 3. The dotted (solid) histograms represent the input (observed) burst amplitude distributions. The cut-off of observed low-amplitude bursts (at ≲103 photons s−1) is due to the fact that the amplitude of these bursts likely occurs below the threshold of the burst-identification algorithm. Table 7. Number of identified bursts nid for Run 3, per configuration. The input number of bursts for each configuration was n = 104. χ = α (°)  nid  30  6645  45  5831  60  4734  90  3717  χ = α (°)  nid  30  6645  45  5831  60  4734  90  3717  View Large 6.3 Detectability of burst phase dependence Here we set out to test the main method used in studies of burst phase dependence to date, see Section 2. During the simulations we determine and record the phase occurrence of the burst peak $$\phi _0^{\rm bf}$$, i.e. the phase occurrence of the best-fitting burst peak time $$t_0^{\rm bf}$$. After each burst we compile a distribution of the values for $$\phi _0^{\rm bf}$$ of all previous bursts up to the most recent one, and compare this burst phase occurrence distribution to a uniform distribution using a Kolmogorov–Smirnov (K–S) test, from which we obtain a p-value. We set the significance threshold at a p-value of 0.003,corresponding to a 0.3 per cent probability that the observed burst peak phase occurrences are distributed uniformly across phase. To be clear, we are simulating emission from a fixed point on the NS surface, from which bursts are being emitted at random rotational phase. The naive expectation that this would result in an observable phase dependence most from the expected modulation of the intensity (see Fig. 3) that results, for example, in missing some bursts emitted on the dark side of the star. In Fig. 14 we plot the evolution of the p-value against the number of observed bursts up to that point in the simulation for the three separate burst durations (from top to bottom) of Run 1. The different curves per subplot correspond to different values of χ and α. The horizontal dashed lines denote the threshold level and the vertical dotted lines indicate the number of bursts at which the p-value of a given simulation drops below the threshold level. Note that of these 12 simulations, the p-value does not fall below the threshold before 103 bursts for χ = α ≤ 45°. Nevertheless, we do find a decreasing trend for T90 = 3.0 s after ∼500 bursts, reaching the threshold at ∼104 bursts (the length of the simulation), for those angles. Remarkably, the p-value associated with the simulation where T90 = 1.50 s and χ = α = 60° does not show a decreasing trend before 104 bursts. The remaining configurations do drop below the threshold fairly soon, i.e. after ∼20 − 150 bursts. Figure 14. View largeDownload slide Evolution of the p-value against the number of bursts for three separate burst durations (top to bottom) of Run 1. The distinct curves per subplot represent different values for the angles χ and α. The horizontal dashed line denotes the threshold level. The vertical dotted lines denote the number of bursts at which the p-value drops below the threshold (see the text for more details). Figure 14. View largeDownload slide Evolution of the p-value against the number of bursts for three separate burst durations (top to bottom) of Run 1. The distinct curves per subplot represent different values for the angles χ and α. The horizontal dashed line denotes the threshold level. The vertical dotted lines denote the number of bursts at which the p-value drops below the threshold (see the text for more details). To determine the minimum number of bursts required to guarantee that the p-value drops below the threshold, we ran the simulation per configuration $$\mathcal {N}_{\rm s}$$ times, each time until the threshold was reached, and recorded the number of bursts. Fig. 15 shows the evolution of the p-value for $$\mathcal {N}_{\rm s}$$ = 10, 100, and 400 simulations, with T90 = 0.15 s and χ = α = 90°. We found that the maximum obtained p-value (out of all $$\mathcal {N}_{\rm s}$$ simulations for a given configuration) after a specific number of bursts decreases at a certain rate (denoted by the black markers). In Fig. 16, we plot the maximum p-value attained, over $$\mathcal {N}_{\rm s}=400$$ simulations per configuration, against the number of bursts. Subsequently, we fit a straight line to the decreasing trends of the log of the p-value and record the number of bursts at which these lines intersect with the threshold level. Accordingly, we find an estimate for the minimum number of bursts nmin at which, assuming a certain configuration, the observed $$\phi _0^{\rm bf}$$ distribution should deviate significantly from a uniform distribution. If the $$\phi _0^{\rm bf}$$ distribution does not significantly deviate from a uniform distribution after nmin bursts, then the configuration will likely be such that the modulation in intensity is less strong than assumed. The latter is dependent on assumptions on the parameters that determine the shape of κ(θ0) (equation 24). Figure 15. View largeDownload slide Evolution of the p-value for $$\mathcal {N}_{\rm s}$$ simulations for T90 = 0.15 s and χ = α = 90°. From left to right we increase the value of $$\mathcal {N}_{\rm s}$$ from 10, to 100, to 400. The distinct coloured curves denote the evolution of the p-value for separate simulations. The black markers denote the maximum p-value that was attained (out of all simulations) for a given number of bursts. The horizontal dashed line denotes the threshold value, the cyan dash-dotted curve represents the fit to the decreasing trend of the log of the maximum p-values, and the vertical dotted curve denotes the intersection of the fit with the threshold level. The latter occurs, from left to right, at nmin = 86, 95, and 104, respectively. The maximum p-values in the panel on the right are also plotted in the top panel of Fig. 16. Figure 15. View largeDownload slide Evolution of the p-value for $$\mathcal {N}_{\rm s}$$ simulations for T90 = 0.15 s and χ = α = 90°. From left to right we increase the value of $$\mathcal {N}_{\rm s}$$ from 10, to 100, to 400. The distinct coloured curves denote the evolution of the p-value for separate simulations. The black markers denote the maximum p-value that was attained (out of all simulations) for a given number of bursts. The horizontal dashed line denotes the threshold value, the cyan dash-dotted curve represents the fit to the decreasing trend of the log of the maximum p-values, and the vertical dotted curve denotes the intersection of the fit with the threshold level. The latter occurs, from left to right, at nmin = 86, 95, and 104, respectively. The maximum p-values in the panel on the right are also plotted in the top panel of Fig. 16. Figure 16. View largeDownload slide Maximum p-value attained, out of $$\mathcal {N}_{\rm s}$$ simulations per configuration, against the number of bursts for the three separate burst durations (top to bottom) of Run 1. The horizontal dashed line denotes the threshold level, the blue dash-dotted curves represent the fits to the decreasing trends, and the vertical dotted lines denote the number of bursts where the fits intersect with the threshold level. Figure 16. View largeDownload slide Maximum p-value attained, out of $$\mathcal {N}_{\rm s}$$ simulations per configuration, against the number of bursts for the three separate burst durations (top to bottom) of Run 1. The horizontal dashed line denotes the threshold level, the blue dash-dotted curves represent the fits to the decreasing trends, and the vertical dotted lines denote the number of bursts where the fits intersect with the threshold level. For the burst durations and configurations that we study in Run 1, we find for χ = α = 90° that nmin ∼ 100 bursts. For χ = α = 60°, we find nmin = 1446 bursts and nmin = 296 bursts, for T90 = 0.15 s and T90 = 3.0 s, respectively. Yet, we do not find an nmin for T90 = 1.50 s, since the attained maximum p-value does not exhibit a decreasing trend before 103 bursts; consistent with the simulation run displayed in the middle panel of Fig. 14. This is because the burst spot remains (partially) visible throughout the NS's rotation, such that enough counts can be detected for the duration of the burst, and the fact that the Δt0 remains comparatively small, i.e. the corresponding parameter density comprises a narrow peak (Fig. 8), in contrast to e.g. the parameter density of the 3.0 s burst, which is much more spread out (Fig. 9). 7 DISCUSSION AND CONCLUSION We have studied, from a theoretical perspective, the conditions under which magnetar bursts from a predefined localized active region or burst patch on the NS surface would give rise to a detectable phase dependence. By adopting a straightforward input burst model, we were able to examine the changes in the observed bursts after they were modulated by the phase-dependent function κ*(ϕ), which takes into account the effects of gravitational light bending and depends on the configuration of the system. We found that the degree to which the inferred burst properties become phase dependent is strongly contingent on the duration of the bursts and geometry of the system; we find a stronger phase dependency of the burst properties for longer duration bursts and larger values of the angles χ and α. The former is because longer bursts sample a wider range of photon trajectories and the latter is due to the fact that for larger values of χ and α the GR effects become more significant. Furthermore, the majority of observed bursts turn out to have τ/τ0 < 1 and T90/T90,0 < 1, i.e. they rise faster and appear shorter than their input counterparts. Attempts to infer the properties of individual bursts with durations greater than ∼10−20  per cent of the spin period should certainly take into account potential distortion due to phase-dependent effects. Adopting a lognormal burst duration distribution that peaks at $$\overline{T}_{90}=0.1$$ s (as observed for well-sampled sources), from which we draw the input duration for each individual burst, we found that phase distributions of the parameters closely resembled those of Run 1, for which T90 = 0.15 s. When considering a power-law distribution for the input burst amplitudes and burst duration of T90 = 1 s, we observed a weak phase dependency of the burst parameters and a considerable amount of scatter, which in turn is caused by the large fraction of low-amplitude input bursts, which are more affected by Poisson noise. We conclude that the observed distributions of burst properties from well-sampled sources are likely not strongly distorted due to phase-dependent effects, by virtue of being dominated by short bursts. We studied the detectability of phase dependence, using the most commonly-used measure (see Section 2) whereby one concentrates on the phase occurrence of the burst peaks. In our set-up all bursts originate at a specific small active region – in some respects the most extreme phase-dependent scenario. However, rotational phase dependence of the peak occurrences was not always apparent. We found that one would require a minimum number of bursts for certain input burst properties and a given system configuration, to guarantee observing a phase dependence. Only in the case of the most extreme geometries, i.e. χ = α > 60°, does this approach the burst sample sizes that were examined in the literature, which range from tens to several hundred bursts. Studies that have not found a phase dependence in the distribution of the burst peak occurrences as yet (e.g. Savchenko et al. 2010; Scholz & Kaspi 2011; Collazzi et al. 2015), might simply require a larger burst sample in order to rule it out for certain geometries. For other geometries, however, it will never be possible to rule out the presence of a burst phase dependence. In our study we have considered only a restricted range of scenarios, where the emission region is tied to the stellar surface. One factor that we have not simulated in detail is that of any potential beaming of the burst emission. To offer brief insights for such influences, we show in Fig. 17 the theoretical phase distributions of the observed burst properties in the case of beamed emission (equation 22 with σb = π/6) from a burst with T90 = 1.5 s; the corresponding shape of κ*(ϕ) is shown by the dashed curves in Fig. 3. We compare it to its isotropic counterpart and find that the phase dependency of the burst properties will be enhanced in the presence of beaming. Figure 17. View largeDownload slide Predicted phase distributions of burst parameters in the presence of beamed emission for a burst of T90 = 1.5 s (solid curves). Beaming is described by equation (22), where we have set the beaming width σb to π/6. For comparison, the dotted curves represent the burst parameter distributions for an isotropic burst with the same duration. Figure 17. View largeDownload slide Predicted phase distributions of burst parameters in the presence of beamed emission for a burst of T90 = 1.5 s (solid curves). Beaming is described by equation (22), where we have set the beaming width σb to π/6. For comparison, the dotted curves represent the burst parameter distributions for an isotropic burst with the same duration. The detectability of a burst phase dependence depends strongly on the shape of κ*(ϕ): the stronger the variation with ϕ the greater the modification to the input burst profiles. Introducing additional bursts patches or allowing for active regions to occur at a certain height above the surface will cause the phase dependence of κ* to decrease. A burst phase dependence in those cases may then only become detectable if the emission is also strongly beamed. In this paper we have not studied the method whereby a phase dependence is searched for in the epoch-folded photon times of arrival. This method can, and should, be subjected to the same level of scrutiny. An additional challenge with this method, however, is to determine a proper false alarm rate. Straightforwardly looking for deviations from uniformity of the times of arrival does not work, since a single burst already consists a significant departure. One must instead quantify the conditions under which one would detect a burst photon phase dependence, even if the bursts originated at random locations on or above the NS surface. We defer this topic to future studies. ACKNOWLEDGEMENTS C.E. and A.L.W. acknowledge support from NOVA, the Dutch Top Research School for Astronomy. D.H. was partially supported by the Moore-Sloan Data Science Environment at New York University, and the James Arthur Postdoctoral Fellowship at New York University. D.H. acknowledges support from the DIRAC Institute in the Department of Astronomy at the University of Washington. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. We thank Matthew Baring for useful comments. Footnotes 1 Dib et al. (2009) discuss the four bursts from AXP 1E 1048.1–5937 and note that only three of them occur near pulse maximum, whereas the fourth burst does not. 2 No further details of the analysis procedure on the SGR 1806–20, such as the applied nominal threshold to determine significance, are given in the article. The phase-dependence analysis procedure of the SGR 1900+14 data is also not described. 3 Note that here we only consider the case where burst emission escapes from the system at the stellar surface (as it would from a trapped fireball, e.g., due to the reduced scattering opacity close to the surface; Thompson & Duncan 1995). It is for surface emission that the effects of GR will be most significant. We argue that bursts that occur high-up in the magnetosphere will be much less affected by the effects of GR or occultation of the star itself, and thus may exhibit weak to no phase-dependent properties. In effect, we are considering the most optimistic case for the detection of phase-dependent effects. 4 If no burst photons are detected, i.e. Ndet = 0, we move on to the next burst [step (II)]. 5 If $$N^{\rm init}<y_0^{\rm sig}$$, we set $$N^{\rm init}=(1+s^{\rm init})\tau ^{\rm init}y_0^{\rm sig}/\delta t$$. REFERENCES An H. et al.  , 2014, ApJ , 790, 60 CrossRef Search ADS   Cheng B., Epstein R. I., Guyer R. A., Young A. C., 1996, Nature , 382, 518 CrossRef Search ADS   Collazzi A. C. et al.  , 2015, ApJS , 218, 1 CrossRef Search ADS   Dib R., Kaspi V. M., Gavriil F. P., 2009, ApJ , 702, 614 CrossRef Search ADS   Gavriil F. P., Kaspi V. M., Woods P. M., 2002, Nature , 419, 142 CrossRef Search ADS PubMed  Gavriil F. P., Kaspi V. M., Woods P. M., 2004, ApJ , 607, 959 CrossRef Search ADS   Gavriil F. P., Kaspi V. M., Woods P. M., 2006, ApJ , 641, 418 CrossRef Search ADS   Gavriil F. P., Dib R., Kaspi V. M., 2011, ApJ , 736, 138 CrossRef Search ADS   Gourgouliatos K. N., Kondić T., Lyutikov M., Hollerbach R., 2015, MNRAS , 453, L93 CrossRef Search ADS   Göğüş E., Kouveliotou C., Woods P. M., Thompson C., Duncan R. C., Briggs M. S., 2001, ApJ , 558, 228 CrossRef Search ADS   Huppenkothen D. et al.  , 2015, ApJ , 810, 66 CrossRef Search ADS   Kaspi V. M., Beloborodov A. M., 2017, ARA&A , 55, 261 CrossRef Search ADS   Kouveliotou C. et al.  , 1998, Nature , 393, 235 CrossRef Search ADS   Kouveliotou C. et al.  , 1999, ApJ , 510, L115 CrossRef Search ADS   Lander S. K., Andersson N., Antonopoulou D., Watts A. L., 2015, MNRAS , 449, 2047 CrossRef Search ADS   Lin L. et al.  , 2012, ApJ , 756, 54 CrossRef Search ADS   Lyutikov M., 2002, ApJ , 580, L65 CrossRef Search ADS   Morsink S. M., Leahy D. A., Cadeau C., Braga J., 2007, ApJ , 663, 1244 CrossRef Search ADS   Mus S. S, Göǧüş E., Kaneko Y., Chakraborty M., Aydin B., 2015, arXiv.org, p. 42 Palmer D. M., 1999, ApJ , 512, L113 CrossRef Search ADS   Palmer D. M., 2002, Mem. Soc. Astron. Italiana , 73, 578 Pechenick K. R., Ftaclas C., Cohen J. M., 1983, ApJ , 274, 846 CrossRef Search ADS   Savchenko V., Neronov A., Beckmann V., Produit N., Walter R., 2010, A&A , 510, A77 CrossRef Search ADS   Scholz P., Kaspi V. M., 2011, ApJ , 739, 94 CrossRef Search ADS   Thompson C., Duncan R. C., 1993, The ApJ , 408, 194 CrossRef Search ADS   Thompson C., Duncan R. C., 1995, MNRAS , 275, 255 CrossRef Search ADS   Thompson C., Yang H., Ortiz N., 2017, ApJ , 841, 54 CrossRef Search ADS   Turolla R., Zane S., Watts A. L., 2015, Rep. Progr. Phys. , 78, 116901 CrossRef Search ADS   van Putten T., Watts A. L., Baring M. G., Wijers R. A. M. J., 2016, MNRAS , 461, 877 CrossRef Search ADS   Woods P. M., Thompson C., 2006, in Lewin W., van der Klis M., eds, Cambridge Astrophysics Series, Compact Stellar X-ray Sources. Cambridge Univ. Press, p. 547 Woods P. M. et al.  , 2005, ApJ , 629, 985 CrossRef Search ADS   Zhu C., Byrd R. H., Lu P., Nocedal J., 1997, ACM Trans. Math. Softw. , 23, 550 CrossRef Search ADS   APPENDIX: PARAMETER TABLE Here we list brief descriptions of the simulation parameters and their associated symbols (see Table A1). Table A1. A table of the simulation parameters appearing in Sections 5 and 6. Symbol  Description  n  Number of input bursts  nid  Number of identified bursts  nmin  Number of bursts at which the p-value is    Guaranteed to drop below the threshold  R  Neutron star radius  P  Neutron star rotation period  ψ  Size of the burst patch  χ  Angle between NS axis of rotation and    the line-of-sight  α  Colatitude of the burst patch  ζ  Background rate  δt  Light-curve bin width  b  Background level  N  Number of emitted photons  Ndet  Number of detected photons  t0  Burst peak time  ϕ0  Burst peak phase occurrence  τ  Rise-time  s  Skewness factor  A  Burst amplitude  $$\mathcal {N}_{\rm bins}$$  Number of time bins  ysig  Maximum departure amplitude  $$t_0^{\rm sig}$$  Time bin with ysig  ΔT  Time interval of an observed burst  tin, tout  Limits of ΔT  μ  Mean count level ∼b  μ*  Running mean  ΔTinterval  Time interval over which μ* is estimated  Δt0  Difference input and best-fitting burst peak    time: Δt0 ≡ tbf − t0  T90  Burst duration  $$\overline{T}_{90}$$  Mode of the duration distribution  $$\sigma _{T_{90}}$$  Width of the duration distribution  $$T_{90}^{\rm min}$$  Minimum burst duration  $$T_{90}^{\rm max}$$  Maximum burst duration  Γ  Power-law index of the amplitude    distribution  Amin  Minimum burst amplitude  Amax  Maximum burst amplitude  $$\mathcal {N}_{\rm s}$$  Number of simulations  subscript ‘0’  Input parameter  superscript ‘init’  Initial guess  superscript ‘bf’  Best-fitting parameter  Symbol  Description  n  Number of input bursts  nid  Number of identified bursts  nmin  Number of bursts at which the p-value is    Guaranteed to drop below the threshold  R  Neutron star radius  P  Neutron star rotation period  ψ  Size of the burst patch  χ  Angle between NS axis of rotation and    the line-of-sight  α  Colatitude of the burst patch  ζ  Background rate  δt  Light-curve bin width  b  Background level  N  Number of emitted photons  Ndet  Number of detected photons  t0  Burst peak time  ϕ0  Burst peak phase occurrence  τ  Rise-time  s  Skewness factor  A  Burst amplitude  $$\mathcal {N}_{\rm bins}$$  Number of time bins  ysig  Maximum departure amplitude  $$t_0^{\rm sig}$$  Time bin with ysig  ΔT  Time interval of an observed burst  tin, tout  Limits of ΔT  μ  Mean count level ∼b  μ*  Running mean  ΔTinterval  Time interval over which μ* is estimated  Δt0  Difference input and best-fitting burst peak    time: Δt0 ≡ tbf − t0  T90  Burst duration  $$\overline{T}_{90}$$  Mode of the duration distribution  $$\sigma _{T_{90}}$$  Width of the duration distribution  $$T_{90}^{\rm min}$$  Minimum burst duration  $$T_{90}^{\rm max}$$  Maximum burst duration  Γ  Power-law index of the amplitude    distribution  Amin  Minimum burst amplitude  Amax  Maximum burst amplitude  $$\mathcal {N}_{\rm s}$$  Number of simulations  subscript ‘0’  Input parameter  superscript ‘init’  Initial guess  superscript ‘bf’  Best-fitting parameter  View Large © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# The rotational phase dependence of magnetar bursts

Monthly Notices of the Royal Astronomical Society, Volume 476 (1) – May 1, 2018
15 pages

/lp/ou_press/the-rotational-phase-dependence-of-magnetar-bursts-ancX3DDvTA
Publisher
Oxford University Press
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty321
Publisher site
See Article on Publisher Site

### Abstract

Abstract The trigger for the short bursts observed in γ-rays from many magnetar sources remains unknown. One particular open question in this context is the localization of burst emission to a singular active region or a larger area across the neutron star. While several observational studies have attempted to investigate this question by looking at the phase dependence of burst properties, results have been mixed. At the same time, it is not obvious a priori that bursts from a localized active region would actually give rise to a detectable phase dependence, taking into account issues such as geometry, relativistic effects, and intrinsic burst properties such brightness and duration. In this paper, we build a simple theoretical model to investigate the circumstances under which the latter effects could affect detectability of dependence of burst emission on rotational phase. We find that even for strongly phase-dependent emission, inferred burst properties may not show a rotational phase dependence, depending on the geometry of the system and the observer. Furthermore, the observed properties of bursts with durations short as 10–20 per cent of the spin period can vary strongly depending on the rotational phase at which the burst was emitted. We also show that detectability of a rotational phase dependence depends strongly on the minimum number of bursts observed, and find that existing burst samples may simply be too small to rule out a phase dependence. stars: magnetars, magnetic fields, gamma-rays: stars, X-rays: bursts 1 INTRODUCTION Magnetars, the most highly magnetized neutron stars (with dipole fields ≳1013 G), are isolated stars powered primarily by magnetic field decay (Thompson & Duncan 1993, 1995; Kouveliotou et al. 1998, 1999). One of their key characteristics is the sporadic emission of soft γ-ray bursts (for reviews of this specific aspect, see Woods & Thompson 2006; Turolla, Zane & Watts 2015; Kaspi & Beloborodov 2017). Durations and fluences vary, but most of these bursts are short, lasting ∼0.01–1 s, less than a typical magnetar spin period P ∼ 6 s. The slow decay of the strong magnetic field is assumed to build up stresses in the system: stress release must involve rapid reconfiguration of the external magnetic field, particle acceleration, and γ-ray emission. However, what triggers the occurrence of individual bursts, the way in which a burst progresses, and the associated emission processes, remain very poorly understood (Turolla et al. 2015). The failure point could be internal, within the crust of the star, or in the external magnetosphere itself. One question is whether the bursts are triggered within a specific active region, fixed in the rotating frame of the star. Some crust zones, for example, are expected to be particularly prone to magnetically induced faulting and yielding (see e.g. Gourgouliatos et al. 2015; Lander et al. 2015; Thompson, Yang & Ortiz 2017). Certain regions of the magnetosphere could also be more active than others, in which case there may also be a preferred height above the neutron star surface. Being able to identify whether this is the case would certainly help in efforts to determine the burst mechanism. One way to determine this, suggested by Lyutikov (2002), is to look at whether there is any rotational phase dependence to the bursts. A number of observational studies have attempted to investigate this, using various different measures such as the phase dependence of the time at which the burst peak is recorded, or the phase distribution of all of the burst photons. The evidence for phase dependence using these measures is mixed (see Section 2 for more details). What has never been done is to determine from a theoretical perspective the circumstances under which bursts from a localized active region would actually give rise to a detectable phase dependence. This will depend on geometry, gravitational light-bending, any beaming factor associated with the burst emission, the size of the burst sample for a given source, and the intrinsic burst properties (e.g. brightness, duration). Under some circumstances, bursts may well be visible throughout most of the rotational phase cycle even if they do originate from a specific active region. In this paper we address this fundamental question of the circumstances under which emission from a localized bursting region would be detectable as a rotational phase dependence according to the measures used in the literature. We then revisit the observational studies carried out to date to see what constraints they actually place on the degree to which bursting might be localized. We also consider the degree to which the rotational phase at which a burst is emitted might affect the properties measured by an observer. First, we provide an overview of magnetar burst phase-dependence studies in the literature in Section 2. Different methods have been applied to various sources in an effort to assess the (non-)phase dependence of magnetar bursts. Next in Section 3, we briefly outline the method through which we aim to answer the aforementioned questions. We choose to simulate sequences of elementary bursts of which we can control the input parameters and study any phase-dependent effects we may observe. The light-curve model is treated in Section 4 and the simulations are described in Section 5. We discuss the results of the simulations and assess the claims made in the literature in Sections 6 and 7. We find that under certain conditions the properties of the observed bursts may become significantly phase-dependent. However, we also find that for a large range of input burst parameters and configurations, a guaranteed detection of phase dependence requires many more bursts than have commonly been observed. 2 OVERVIEW OF PUBLISHED BURST PHASE-DEPENDENCE ANALYSIS Here we provide a review of previous work where, given the acquired data, the phase dependence of magnetar bursts has been evaluated. We focus on how the data were obtained and processed, and what method was used to determine the absence or presence of a phase dependence in the burst occurrences or properties. In practice, three methods have been applied: (i) searching for any significant deviations from uniformity of burst occurrence and photon arrival-time distributions against phase, (ii) searching for any correlation between the phase at which bursts occur and the pulse maxima of the (underlying) pulsed emission, and (iii) Fourier analysis on the burst occurrence times in an effort to search for significant periodicities. The latter has been applied only once; the first two are far more common. It is worth mentioning that the first two methods depend on the accuracy of the ascertained timing ephemeris. The longer the time baseline spanned by the bursts, the greater the risk of undetected time anomalies, such as glitches or spin-down deviations, that may undermine the inference of the phase. Table 1 provides a summary of the references that have carried out phase-dependence analysis of magnetar bursts. Table 1. Summary of burst phase-dependence results in the literature. Reference  Source  Dates of bursts  Satellite/  nbursts  Methodb      [dd/mm/‘yy]  Instrumenta      Gavriil et al. (2004)  1E 2259+586  18/06/‘02  RXTE/PCA  80  (i)  Savchenko et al. (2010)  SGR J1550–5428  22/01/‘09  INTEGRAL/ACS  84  (i)  Scholz & Kaspi (2011)    (22/01–30/09)/‘09  Swift/XRT  303  (i)  Lin et al. (2012)    22/01/‘09  Swift/XRT  31  (i)      30/01/‘09        Collazzi et al. (2015)    03/10/‘08–17/04/‘09  Fermi/GBM  354  (i)  Mus et al. (2015)    22/01/‘09  RXTE/PCA  4  (ii)      06/02/‘09            30/03/‘09            11/01/‘10        Gavriil et al. (2002)  1E 1048.1–5937  29/10/‘01  RXTE/PCA  2  (ii)      14/11/‘01        Gavriil et al. (2006)    29/06/‘04  RXTE/PCA  1  (ii)  Dib et al. (2009)    29/10/‘01  RXTE/PCA  4  (ii)      14/11/‘01            29/06/‘04            28/04/‘08        An et al. (2014)    (17–27)/07/‘13  NuSTAR  8  (ii)  Woods et al. (2005)  XTE J1810–197  22/07/‘03  RXTE/PCA  6  (ii)      16/02/‘04            19/04/‘04            19/05/‘04        Gavriil et al. (2011)  4U 0142+61  (06/04–26/06)/‘06  RXTE/PCA  6  (ii)  Palmer (1999)  SGR 1806–20  (10–15)/11/‘83  ICE  33  (iii)  Palmer (2002)  SGR 1900+14  -  -  -  -  Reference  Source  Dates of bursts  Satellite/  nbursts  Methodb      [dd/mm/‘yy]  Instrumenta      Gavriil et al. (2004)  1E 2259+586  18/06/‘02  RXTE/PCA  80  (i)  Savchenko et al. (2010)  SGR J1550–5428  22/01/‘09  INTEGRAL/ACS  84  (i)  Scholz & Kaspi (2011)    (22/01–30/09)/‘09  Swift/XRT  303  (i)  Lin et al. (2012)    22/01/‘09  Swift/XRT  31  (i)      30/01/‘09        Collazzi et al. (2015)    03/10/‘08–17/04/‘09  Fermi/GBM  354  (i)  Mus et al. (2015)    22/01/‘09  RXTE/PCA  4  (ii)      06/02/‘09            30/03/‘09            11/01/‘10        Gavriil et al. (2002)  1E 1048.1–5937  29/10/‘01  RXTE/PCA  2  (ii)      14/11/‘01        Gavriil et al. (2006)    29/06/‘04  RXTE/PCA  1  (ii)  Dib et al. (2009)    29/10/‘01  RXTE/PCA  4  (ii)      14/11/‘01            29/06/‘04            28/04/‘08        An et al. (2014)    (17–27)/07/‘13  NuSTAR  8  (ii)  Woods et al. (2005)  XTE J1810–197  22/07/‘03  RXTE/PCA  6  (ii)      16/02/‘04            19/04/‘04            19/05/‘04        Gavriil et al. (2011)  4U 0142+61  (06/04–26/06)/‘06  RXTE/PCA  6  (ii)  Palmer (1999)  SGR 1806–20  (10–15)/11/‘83  ICE  33  (iii)  Palmer (2002)  SGR 1900+14  -  -  -  -  aSpacecraft/instrument acronyms: Rossi X-ray Timing Explorer (RXTE), Proportional Counter Array (PCA), Nuclear Spectroscopic Telescope Array (NuSTAR), International Cometary Explorer (ICE), Anti-Coincidence Shield (ACS), X-ray Telescope (XRT), Nuclear Spectroscopic Telescope Array (NuSTAR), and Gamma-ray Burst Monitor (GBM). bThe methods are specified in Section 2. View Large The active phase of 1E 2259+586 on 2002 June 18 consisted of 80 bursts and was studied using method (i) (Gavriil, Kaspi & Woods 2004); it was claimed that the burst peak phase occurrences tended to correlate with the intensity of the pulsed emission, yet no phase dependencies were observed for the burst durations, fluences, peak fluxes, and rise/fall times. In excess of 300 bursts were observed from SGR J1550–5428 between 2008 March and 2010 January and many of those bursts were detected by multiple space-based telescopes simultaneously. Savchenko et al. (2010) found that the burst start times (the moment the burst exceeds 5σ above background) of 84 bursts, observed with the Anti-Coincidence Shield (ACS) aboard the INTEGRAL spacecraft, appear to be distributed randomly across phase, i.e. no significant departure from the mean bursts per phase bin was identified. Scholz & Kaspi (2011) and Collazzi et al. (2015) studied the burst peak times of, respectively, 303 and 354 bursts, and both found no significant (>3σ) deviations from the mean number of burst peaks per phase bin. Scholz & Kaspi (2011) however do show that the phase-folded photon times of arrival of the bursts exhibit an apparent pulse which has an offset with respect to maximum of the associated quiescent pulse profile. Lin et al. (2012) study a sample of 31 bursts and similarly find that the burst count distribution is not uniform across phase. Moreover, they find that the phase probability density anti-correlates with the phase profile of the persistent emission (with a correlation factor of −0.5 and chance probability of 3.4 × 10−2), which may suggest that the burst emission region is distinct to that of the persistent emission. Contrary to these results however, Collazzi et al. (2015) do not find a significant (>3σ) pulse shape in the epoch folded burst emission light curves. Note that the data set used by Lin et al. (2012) constitutes a subset of the data used by Collazzi et al. (2015). A total of 12 bursts from 1E 1048.1–5937 were analysed using method (ii); 4 of which were observed with RXTE/PCA between 2001 October 29 and 2008 April 28 (Gavriil, Kaspi & Woods 2002; Gavriil, Kaspi & Woods 2006; Dib, Kaspi & Gavriil 2009) and 8 of which were observed with NuSTAR in 2013 July 17–27 (An et al. 2014). It was determined that the majority of the bursts1 observed with RXTE/PCA had a probable chance alignment with pulse maxima of less than 0.01. For the latter eight bursts from the same source however there is no evidence for a preferred phase occurrence. Six bursts from XTE J1810–197 observed with RXTE/PCA were also studied with method (ii) (Woods et al. 2005). These bursts consisted of individual burst spikes, which in turn occurred near the corresponding pulse maxima of the source, either leading or trailing. A chance alignment of these spikes with the pulse maxima was estimated at roughly 0.004. Gavriil, Dib & Kaspi (2011) studied six bursts of 4U 0142+61 observed with RXTE/PCA in 2006 from April to June. They found that several bursts appear to occur near the maxima of contemporaneous folded pulse profiles (no significance criteria are specified in the reference). They argue that this may indicate that the bursts comprise extreme episodes of local transient emission sites. Palmer (1999, 2002) studied the burst properties of SGR 1806–20 and SGR 1900+14, where for the former source it was found that active bursting episodes emerge from local active regions characterized as ‘relaxation systems’. From a larger burst sample, a group of 33 bursts was identified as belonging to a single relaxation system. Method (iii), i.e. Fourier analysis on the burst occurrences of this group, revealed no apparent modulation at the rotation frequency of the NS, indicating the lack of a phase dependence.2 Here we will focus mainly on the (non-)uniform phase occurrence of burst peaks as a proxy for the phase dependency of magnetar bursts. We briefly discuss the use of alternative methods in Section 7. 3 METHODOLOGY To understand how the observed emission may depend on the rotational phase, we intentionally introduce a phase dependency by fixing the burst location to a certain region or burst patch on the magnetar surface and then set out to describe and simulate the process from emission, where we control the input parameters, to detection and characterization. Subsequently, we can study the effects of a certain configuration on the burst parameters by investigating the phase distributions of the observed burst properties. Moreover, we can establish detectability criteria for the phase dependency for certain input values/distributions and system configurations. In order to do so, we require a light-curve model that describes how the burst emission is modified depending on the location of the bursts and additional system parameters, e.g. the inclination angle to the observer and compactness of the source. The latter parameter will reshape the trajectory of emitted photons through gravitational light bending. In Section 4 we ascertain an expression that describes the fraction of rays, i.e. paths along which the emitted photons propagate, that extend out from the burst location and intersect with an observer at infinity. Subsequently, in Section 5, we simulate sequences of bursts and investigate how the burst properties are modified through the correction of the burst intensity by the aforementioned expression. 4 LIGHT CURVE MODEL In the following we adopt natural units, i.e. G = c = 1, and the spatial spherical coordinates (r, φ, θ), where φ is the polar angle to the y-axis and θ is azimuthal angle to the z-axis. Since magnetars rotate slowly (typically $$|\boldsymbol{\Omega }|\sim 10^{-1}$$ rad s−1), they can be considered to be almost spherically symmetric. Accordingly, we may assume that the metric external to the star is approximately given by the Schwarzschild space–time solution,   $$\text{d}s^2=-\mathcal {A}(r)\text{d}t^2+\mathcal {A}^{-1}(r)\text{d}r^2+r^2 \left(\text{d}\varphi ^2+\sin ^2\varphi \,\text{d}\theta ^2\right)$$ (1)where $$\mathcal {A}(r)=(1-R_S/r)$$, and RS = 2M the Schwarzschild radius, with M corresponding to the gravitational mass of the compact object. To model the effect of gravitational light bending on the burst emission, we consider the configuration illustrated in Fig. 1, which is based on the work done by Pechenick, Ftaclas & Cohen (1983). The stellar surface is located at a distance R from the origin; for neutron stars, R lies roughly in the range 2.5 − 4RS. For now we assume that the burst emission originates at3r = R over a circular patch of angular radius ψ centred at point $$\boldsymbol{p}(R, \varphi , \theta )$$ with total intensity I. Depending on the burst emission mechanism, I may be anisotropic and depends on δ ∈ [0, π), i.e. the angle between the normal vector to the stellar surface $$\hat{\boldsymbol{n}}$$ and the outgoing emission vector $$\hat{\boldsymbol{k}}_R$$. The latter lies at r = R along the associated null geodesic $$\mathcal {G}$$ of the outgoing emission, which in turn intersects with the observer at r = r0, where $$\hat{\boldsymbol{k}}_R\rightarrow \hat{\boldsymbol{k}}$$. We presume that the region r < R is opaque and r > R is entirely transparent. Moreover, due to the comparatively long rotation period of magnetars, we may neglect certain corrections, such as oblateness of the stellar surface, light travel-time delays, and Doppler effects, which become significant for NSs with $$|\boldsymbol{\Omega }|\gg 1$$ rad s−1 (e.g. Morsink et al. 2007). Figure 1. View largeDownload slide Schematic representation of null geodesic $$\mathcal {G}$$ connecting the burst patch, centred at point $$\boldsymbol{p}(\theta _0)$$ of angular width ψ, with the observer. The observer is set in the +z-direction, where the z-axis makes an angle χ with the rotational axis of the neutron star (NS) $$\boldsymbol{\Omega }$$. The normal vector to the NS surface $$\hat{\boldsymbol{n}}$$ at $$\boldsymbol{p}$$ is at an angle α with $$\boldsymbol{\Omega }$$ and at an angle θ0 with the observer's line of sight, where the latter depends on the rotational phase of the NS, i.e. θ0 = θ0(ϕ) as prescribed by equation (16). In a stationary frame at r = R, photons radiate from the burst location (along $$\hat{\boldsymbol{k}}_R$$) at and angle δ to the normal. The intensity of the source may be isotropic or beamed I = I(δ), depending on the underlying emission mechanism and properties of the emitting region. Figure 1. View largeDownload slide Schematic representation of null geodesic $$\mathcal {G}$$ connecting the burst patch, centred at point $$\boldsymbol{p}(\theta _0)$$ of angular width ψ, with the observer. The observer is set in the +z-direction, where the z-axis makes an angle χ with the rotational axis of the neutron star (NS) $$\boldsymbol{\Omega }$$. The normal vector to the NS surface $$\hat{\boldsymbol{n}}$$ at $$\boldsymbol{p}$$ is at an angle α with $$\boldsymbol{\Omega }$$ and at an angle θ0 with the observer's line of sight, where the latter depends on the rotational phase of the NS, i.e. θ0 = θ0(ϕ) as prescribed by equation (16). In a stationary frame at r = R, photons radiate from the burst location (along $$\hat{\boldsymbol{k}}_R$$) at and angle δ to the normal. The intensity of the source may be isotropic or beamed I = I(δ), depending on the underlying emission mechanism and properties of the emitting region. The Schwarzschild solution admits four Killing vectors associated with conserved quantities that emerge from symmetries inherent in the solution for which,   $$K_\mu \dot{x}^\mu =\text{constant},$$ (2)where $$\dot{x}^{\mu }=\text{d}x^{\mu }/\text{d}\lambda$$, with λ is some affine parameter. Considering the orbital motion of a photon with tangent 4-vector $$V^\mu =\dot{x}^{\mu }$$ in the equatorial plane, i.e. φ = π/2, the Schwarzschild solution admits two Killing vectors   \begin{eqnarray} \epsilon _\mu =(\mathrm{\partial} _t)_\mu = (-\mathcal {A}(r),0,0,0), \end{eqnarray} (3)  \begin{eqnarray} J_\mu =(\mathrm{\partial} _\varphi )_\mu = (0,0,0,r^2), \end{eqnarray} (4)associated, respectively, with conservation of energy and the magnitude of angular momentum. Accordingly, we may define   \begin{eqnarray} \epsilon _\mu \dot{x}^\mu =-\mathcal {A}(r)V^t\equiv -1, \end{eqnarray} (5)  \begin{eqnarray} J_\mu \dot{x}^\mu =r^2V^\theta \equiv b, \end{eqnarray} (6)where b ≥ 0 denotes the impact parameter of the photon trajectory, i.e. the null geodesic $$\mathcal {G}$$. Since $$\dot{x}^\theta =0$$ and $$g_{\mu \nu } \dot{x}^\mu \dot{x}^\nu =0$$ for massless particles, we find that the tangent 4-vector of an outgoing photon is given by,   $$V^{\mu }=\left(\mathcal {A}^{-1}(r), \sqrt{1-\frac{b^2}{r^2}\mathcal {A}(r)},0, \frac{b}{r^2}\right).$$ (7)Setting b = 0, we obtain the tangent 4-vector of a radially outgoing photon,   $$W^{\mu }=(\mathcal {A}^{-1}(r),1,0,0).$$ (8)A stationary observer with 4-velocity   $$U^{\mu }=(\mathcal {A}^{-1/2}(r),0,0,0),$$ (9)will observe an angle   $$\cos \xi =\sqrt{1-\frac{b^2}{r^2}\mathcal {A}(r)},$$ (10)between the photons prescribed by Vμ and Wμ (Pechenick et al. 1983). Note that $$\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{k}}_R=\cos \delta = \cos \,(\xi |_{r=R})$$, such that we may write   $$\delta (b)=\arcsin \left(\frac{b}{b_{\rm max}}\right), \ \ \ \ \text{with}\ \ b_{\rm max}=\frac{R}{\mathcal {A}^{1/2}(r=R)}.$$ (11)The total angular deflection of $$\mathcal {G}$$, which determines the ‘bending’ of the photon trajectory from the surface patch to the observer, is given by   $$\theta _*(b) =\int _R^{r_0}\frac{b}{r^2}\left[1-\frac{b^2}{r^2}\mathcal {A}(r)\right]^{-1/2}\text{d}r.$$ (12)Incidentally, the total coordinate light travel time along $$\mathcal {G}$$ is given by   $$T_*(b) =\int _R^{r_0}\mathcal {A}^{-1}(r)\left[1-\frac{b^2}{r^2}\mathcal {A}(r)\right]^{-1/2}\text{d}r.$$ (13) The difference in travel time between radially emitted photons (with b = 0) and those with an arbitrary impact parameter can be estimated accordingly,   $$\Delta t_*(b) =\int _R^{r_0}\mathcal {A}^{-1}(r)\left\lbrace \left[1-\frac{b^2}{r^2}\mathcal {A}(r)\right]^{-1/2}-1\right\rbrace \text{d}r.$$ (14)The maximum travel time delay then for a typical NS with R = 2.5 RS (R = 106 cm, M = 1.5 M⊙), is Δt*(bmax) ≃ 6.7 × 10−2 ms ≪ P ∼ 6 s. To an observer in the +z-direction, the system is axisymmetric around the z-axis, such that the location of the burst patch $$\boldsymbol{p}$$ can be uniquely described by θ0. The angle between the observer's line of sight and the rotation axis of the neutron star $$\boldsymbol{\Omega }$$ is denoted by χ. Furthermore, the angle between the location of the burst patch and $$\boldsymbol{\Omega }$$ is given by α. Accordingly, depending on the rotational phase of the neutron star,   $$\phi (t)=\frac{2\pi t}{P},$$ (15)the angle between the observer and burst patch θ0 is given by the relation   $$\cos [\theta _0(\phi )] = \cos \chi \cos \alpha + \sin \chi \sin \alpha \cos \phi .$$ (16) The observed brightness is the integral of the intensity at the observer,   $$\text{d}I_{\rm obs}=I(r_0, \Omega ^{\prime })\text{d}\Omega ^{\prime }=\mathcal {A}^2(r=R)I(R, \Omega ^{\prime })\text{d}\Omega ^{\prime },$$ (17)over the solid angle subtended in the observer's sky,   $$\text{d}\Omega ^{\prime }=\sin \theta ^{\prime } \text{d}\theta ^{\prime } \text{d}\varphi ^{\prime }\simeq \theta ^{\prime } \text{d}\theta ^{\prime } \text{d}\varphi ^{\prime }.$$ (18)Due to the axisymmetry dφ΄ = dφ ≡ Φ(θ*, θ0), where we define the polar differential distance as the function Φ(θ*, θ0), which under the conditions θ0 + ψ ≤ θ* ≤ π and θ0 − ψ ≥ 0 is given by the following expression:   \begin{eqnarray} \Phi (\theta _*, \theta _0) = \left\lbrace \begin{array}{ll}2\arccos \left(\frac{\cos \psi -\cos \theta _0\cos \theta _*}{\sin \theta _0\sin \theta _*}\right) & \text{if } \theta _0-\psi \le \theta _* \\ &\le \theta _0 +\,\psi ,\\ 0 & \text{otherwise.} \end{array} \right. \end{eqnarray} (19)Consequently, together with θ΄ = ξ and equation (10) evaluated at r0 → ∞, we obtain   $$d\Omega ^{\prime }\simeq \Phi [\theta _*(b), \theta _0]\, \xi \text{d}\xi \simeq \frac{\Phi (b, \theta _0)}{r_0^2}b\,\text{d}b.$$ (20)We write the brightness of the source as   $$I(R, \Omega ^{\prime })=I_0f[\delta (b)],$$ (21)with the beaming functions given by f(b). Currently, we do not have a physical model for the shape of the beaming function, which will most likely depend on the radiative transfer properties of the local magnetic field. Accordingly, the influence of the magnetic field on the light trajectories is ignored for now. An example of a more realistic model was considered by van Putten et al. (2016) in the case of fireball beaming. Here, for descriptive purposes, we consider a Gaussian shape for the beaming function,   $$f[\delta (b)]= \left\lbrace \begin{array}{ll}1 & \text{isotropic}, \\ \sqrt{\frac{\pi }{2\sigma _{\rm b}^2}}\ {\rm erf}\left(\frac{\pi }{2\sqrt{2}\sigma _{\rm b}}\right)^{-1}\exp \left(-\frac{\delta (b)^2}{2\sigma _{\rm b}^2}\right)& \text{beamed}, \end{array} \right.$$ (22)where σb parameterizes the beam width. We neglect rotational aberration of light effects, since the star rotates slowly. Finally, we find the expression for the observed intensity of the source   $$I_{\rm obs}(\theta _0)=I_0\left(\frac{R}{r_0}\right)^2\kappa (\theta _0),$$ (23)where we define   $$\kappa (\theta _0)\equiv \left(\frac{R^{1/2}}{b_{\rm max}}\right)^4\int ^{b_{\rm max}}_0f[\delta (b)]\, \Phi (b, \theta _0)\,b\,\text{d}b.$$ (24) Fig. 2 shows the observed burst emission Iobs as a function of burst patch location θ0 for a compact object with R = 2.5 RS. In this case the size of the patch is ψ = 1°. We consider the observed intensity for the case of isotropic and beamed emission. Note that the location of the terminator lies ‘behind’ the star, i.e. beyond θ0 = π/2, at θ0 ≃ 0.72π. At this angle, only photons with an impact parameter of bmax ≃ 3.23RS reach the observer. Note that the beamed emission is more prominent at θ0 = 0 and drops off faster than the isotropic emission with increasing θ0. For comparison, we plotted the emission profiles of sources with R = 1.6 RS and R ≫ RS, where the former is close to the most extreme case, i.e. R > 1.5 RS and the latter approximates flat spacetime. Figure 2. View largeDownload slide Burst emission as a function of the angle between the observer and the burst location, parametrized by θ0 (see Fig. 1) for a burst spot area of ψ = 1°. The solid (dashed) curves denote emission from an isotropic (a beamed) source. A beam width of σb = π/6 was used. The black curves represent emission from a compact object with R = 1.6 RS (close to the most extreme case, i.e. R > 1.5 RS). The magenta curves denote the emission from a typical neutron star R = 2.5 RS (R = 106 cm, M = 1.5 M⊙). The blue curves denote the emission from a non-relativistic source R ≫ RS, i.e. if we neglect the GR light bending effects. The dotted lines indicate the location of the terminator, respectively, at θ0 = 0.5π and θ0 ∼ 0.72π for R ≫ RS and R = 2.5 RS, beyond which the burst patch is invisible. Note that emission is always visible in the R = 1.6 RS case, regardless the location of the burst, i.e. ∀θ0. Moreover, in this extreme case the intensity of the burst patch is greatly amplified at θ0 → π due to gravitational lensing effects. This plot is similar to fig. 4 in Pechenick et al. (1983). Figure 2. View largeDownload slide Burst emission as a function of the angle between the observer and the burst location, parametrized by θ0 (see Fig. 1) for a burst spot area of ψ = 1°. The solid (dashed) curves denote emission from an isotropic (a beamed) source. A beam width of σb = π/6 was used. The black curves represent emission from a compact object with R = 1.6 RS (close to the most extreme case, i.e. R > 1.5 RS). The magenta curves denote the emission from a typical neutron star R = 2.5 RS (R = 106 cm, M = 1.5 M⊙). The blue curves denote the emission from a non-relativistic source R ≫ RS, i.e. if we neglect the GR light bending effects. The dotted lines indicate the location of the terminator, respectively, at θ0 = 0.5π and θ0 ∼ 0.72π for R ≫ RS and R = 2.5 RS, beyond which the burst patch is invisible. Note that emission is always visible in the R = 1.6 RS case, regardless the location of the burst, i.e. ∀θ0. Moreover, in this extreme case the intensity of the burst patch is greatly amplified at θ0 → π due to gravitational lensing effects. This plot is similar to fig. 4 in Pechenick et al. (1983). In the simulations we concentrate on the relative changes in intensity between the input and observed burst. Accordingly, from equations (16) and (24), we define   $$\kappa _*(\phi )\equiv \frac{\kappa [\theta (\phi )]}{\kappa _{\rm max}},$$ (25)which depends on the angles χ, α, and the phase of the neutron star, and describes the fraction of rays that intersect with the observer at infinity, from the entire ray-bundle that extends outwards from the burst patch. In the following section, we use this expression as our measure for how the burst intensity is modulated. Varying χ separately from α, or vice versa, acts as a multiplicative factor to the absolute intensity. Since, we only consider the fractional intensity, we may explore the parameter space of these angles by setting χ = α. Fig. 3 illustrates the shape of κ*(ϕ) for R = 2.5 RS in four different angle configurations, in the case of both isotropic and beamed emissions. Figure 3. View largeDownload slide Fraction of rays that intersect with the observer at infinity, from the entire ray-bundle that extends outwards from the burst patch, i.e. κ*(ϕ). The solid (dashed) curves represent isotropic (beamed, with σb = π/6) emission from a source with R = 2.5 RS. The colours represent four different geometries, given by χ and α. Figure 3. View largeDownload slide Fraction of rays that intersect with the observer at infinity, from the entire ray-bundle that extends outwards from the burst patch, i.e. κ*(ϕ). The solid (dashed) curves represent isotropic (beamed, with σb = π/6) emission from a source with R = 2.5 RS. The colours represent four different geometries, given by χ and α. 5 SIMULATIONS Per simulation run we produce a sequence of n bursts where we control the input parameters of the bursts and the system. Nonetheless, we treat the detection of individual photons, which are entirely described by their times-of-arrival (TOA), in a probabilistic fashion. We assume for simplicity that a single magnetar burst can be modelled with a exponential rise/exponential decay profile that allows for asymmetry, i.e. the rise time and fall time may be different. When simulating photons emitted at the source, this corresponds to drawing N random photon times-of-emission (TOE) from a skewed Laplace distribution,   $$p_{\rm TOE}(t)=\frac{1}{(1+s)\tau }\left\lbrace \begin{array}{ll}\exp \left[(t-t_0)/\tau \right] & \text{if }\ t < t_0, \\ \exp \left[-(t-t_0)/(s\tau )\right] & \text{if }\ t \ge t_0, \end{array} \right.$$ (26)where τ denotes the exponential rise time of the burst, (sτ) the decay time, with s the skewness factor, and t0 the peak time of the burst. The profile of a simulated burst is shown in Fig. 4. We deliberately adopt an oversimplified burst profile to better understand the differences between the input and output data. Huppenkothen et al. (2015) decompose complex magnetar bursts, observed from SGR J1550-5418, into several spike-like components, which in turn are modelled with a similar profile as in equation (26). In this paper, we assume that a burst can be represented as a single spike, as a simple model that lets us explore the relevant effects. Note that we define the duration of the burst, T90, as the time it takes for the fluence to increase from 0.05 to 0.95 of the total burst fluence. We fix the compactness of the source to R = 2.5 RS, corresponding to a typical neutron star with R = 106 cm and M = 1.5M⊙, and set the rotation period to P = 6 s. We choose a light-curve bin width, δt, and background count level, b, such that the background count rate approximates that from Fermi/GBM data (Huppenkothen et al. 2015), i.e. ζGBM ∼ 318 counts s−1. A list of the simulation parameters is given in Table A1. Figure 4. View largeDownload slide Profile of an input burst, i.e. the photon emission rate at the burst patch, where t0 represents the time at which the burst peaks, A denotes the burst amplitude, τ denotes the rise time, τs parameterizes the decay time, determined by s the skewness parameter for a given τ. For s > 1 (s < 1) the burst rises faster (slower) than it decays. The burst duration is given by T90, which is defined as the time the fluence increases from 0.05 to 0.95 of the total burst fluence; the interval T90 contains 0.90N photons. Figure 4. View largeDownload slide Profile of an input burst, i.e. the photon emission rate at the burst patch, where t0 represents the time at which the burst peaks, A denotes the burst amplitude, τ denotes the rise time, τs parameterizes the decay time, determined by s the skewness parameter for a given τ. For s > 1 (s < 1) the burst rises faster (slower) than it decays. The burst duration is given by T90, which is defined as the time the fluence increases from 0.05 to 0.95 of the total burst fluence; the interval T90 contains 0.90N photons. 5.1 General simulation procedure Here we proceed to describe in more detail the general form of a simulation run step by step: We decide on the number of bursts n we wish to produce per simulation run and set the inclination angle of the source χ. Next, we assign values to the burst parameters ψ, α, N, t0, τ, s, and T90, where the latter three parameters cannot be defined independently of each other. In the following simulation runs, described in Sections 5.2, 5.3, and 5.4, we only consider symmetric input bursts, s = 1, with a size of ψ = 1°, and draw their peak time from a uniform distribution, i.e. t0 ∼ Uniform(0, P), where P is the rotation period of the magnetar which we set to P = 6. We generate a single burst by drawing N photon emission times (TOEs) from pTOE(t): We draw random numbers from a uniform distribution Uniform(0, 1) and transform these values to follow the required skewed Laplace distribution by using the inverse cumulative distribution function of the latter, i.e. the percent point function,   \begin{eqnarray} {\rm TOE}(x)=\left\lbrace \begin{array}{ll}t_0 + \tau \ln \left[\left(1+s\right)x\right] & \text{if }\ x < (1+s)^{-1}, \\ t_0 - s\tau \ln \left[\left(1+s^{-1}\right)\left(1-x\right)\right] & \text{if }\ x \ge (1+s)^{-1}. \end{array} \right. \nonumber\\ \end{eqnarray} (27) Using equation (15), we determine the phase of each TOEi, i.e. ϕi. Whether an emitted photon reaches the observer depends on whether the ray, along which the photon propagates, intersects with the detector, given by κ*(ϕi), which denotes the fraction of photons directed into our line of sight. In order to decide whether a given photon intersects with the detector, we use rejection sampling: for each TOEi we generate a latent variable z drawn from p(z) = Uniform(0, 1). We only keep the TOEi if z < κ*(ϕi).4 The TOEs that we save are detected by the observer. A detected photon is recorded as a count with a corresponding TOAi, where the TOAi = TOEi of the respective photon, since we consider a perfect detector and may neglect the distance to the source and gravitational time-delay effects (see Section 4). Furthermore, we add background counts or TOAs uniformly to the detected TOA data from the burst (with length Ndet ≤ N), such that the mean background count rate becomes approximately ζ. We bin the total TOA data in $$\mathcal {N}_{\rm bins}$$ time bins of length δt, whereby the counts in each bin follow Poisson statistics. We proceed by applying a similar burst identification algorithm as used by Gavriil et al. (2004), assuming that we can infer the background count rate to be ζ, yet have no prior knowledge of the burst and system input parameters. The probability of the number of counts ki in the ith bin occurring is given by the Poisson distribution,   $$P_i=\frac{\mu ^{k_i}\text{e}^{-\mu }}{k_i!},$$ (28)where μ represents the mean count level, which in our case will be b = ζ δt. Bins for which   $$P_i \le 3\times 10^{-3}\mathcal {N}_{\rm bins}^{-1}$$ (29)are recorded as significant departures from the mean, where we have corrected for the number of trails by dividing by $$\mathcal {N}_{\rm bins}$$, i.e. the total number of bins searched over. From these, the time bin containing the maximum departure ysig is labeled as $$t_0^{\rm sig}$$. The burst edges, labeled as tin and tout, are ascertained by making use of a running mean, i.e. when the mean count level of an interval μ* of ΔTinterval = 0.25 s, moving outwards in steps of δt on both sides of $$t_0^{\rm sig}$$, falls below b* = 1.1b, the burst edges are then given by the centre time of the respective intervals, tin before and tout after the burst. The duration of this interval is denoted ΔT = |tin − tout|. We fit the light curve of any identified bursts with the following burst model,   $$m(t) = N p_{\rm TOE}(t\,|\,t_0, \tau ,s) + \zeta ,$$ (30)using an L-BFGS-B constrained optimizer (Zhu et al. 1997) to determine the maximum (Poisson) likelihood, whereby we fix the background parameter ζ and provide initial guesses for the remaining parameters:   \begin{eqnarray} t_0^{\rm init} = t_0^{\rm sig}, \end{eqnarray} (31)  \begin{eqnarray} \tau ^{\rm init} = \left(t_0^{\rm sig} - t_{\rm in}\right)\left[\ln \left(\frac{y^{\rm sig}}{b^*}\right)\right]^{-1}, \end{eqnarray} (32)  \begin{eqnarray} s^{\rm init} = 1, \end{eqnarray} (33)and Ninit is defined as the number of counts in the interval ΔT minus the background counts,5 i.e. ζΔT. After the fit, we delete the bins in ΔT from the light curve, and repeat steps (V) and (VI) until no significant departures from the mean, i.e. μ = b, are recorded. We return to step (II) until we have generated the pre-defined number of bursts n. Note that the number of observed bursts might be different, since bursts might go undetected, or be interpreted as multiple separate bursts. Moreover, some identified ‘bursts’ may simply be significant statistical deviations from the background level. However, according to the condition stated in equation (29), we only expect this to be the case in ∼0.3 per cent of the input bursts. 5.2 Run 1: initial simulation run We start with the simplest scenario, where we consider simulations of sequences of identical bursts, referred to as Run 1. Per simulation we fix the values for χ, α, ψ, s, and T90. The latter two parameters determine the value of τ. Subsequently, we define N using the condition that the input burst amplitude A is 104 photons s−1. We run simulations for three separate burst durations (see Fig. 5) and vary the angles χ and α, to study their effects on the observed quantities. Figure 5. View largeDownload slide Folded profiles of input bursts of the initial simulation run. We consider symmetric bursts of three separate durations. Here the rise-time τ is determined by the values for s and T90. We choose N such that the burst amplitude is A = 104 photons s−1 at t0. Figure 5. View largeDownload slide Folded profiles of input bursts of the initial simulation run. We consider symmetric bursts of three separate durations. Here the rise-time τ is determined by the values for s and T90. We choose N such that the burst amplitude is A = 104 photons s−1 at t0. We concentrate on the difference in input and observed best-fitting value for the time of the burst peak (respectively, t0 and $$t_0^{\rm bf}$$), where the difference is parametrized as $$\Delta t_0\equiv t_0^{\rm bf}-t_0$$, the rise-time τ, skewness factor s, and burst duration T90. The input values of the latter three parameters are denoted as τ0, s0, and T90,0. All input parameters of Run 1 are listed in Table 2. Table 2. Input parameters for Run 1, consisting of 12 separate simulations. We consider a constant input burst profile with peak times distributed uniformly across phase. Per simulation we vary the burst duration T90, and angles χ, α, where we set χ = α (see Section 4). Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  0.15, 1.5, 3.0  δt (s)  200−1  Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  0.15, 1.5, 3.0  δt (s)  200−1  View Large 5.3 Run 2: T90-distribution Next we perform a simulation run, Run 2, where in step (I) of the general simulation procedure (Section 5.1), we draw the burst duration T90 for each individual burst from a lognormal distribution centred at $$\overline{T}_{90}=0.1$$ s, with a width of $$\sigma _{T_{90}}=1$$ (e.g. Göğüş et al. 2001), and lower and upper burst duration cut-off at, respectively, $$T_{90}^{\rm min}=300^{-1}$$ s and $$T_{90}^{\rm max}=3$$ s < P. Fixing the input burst amplitude at A = 104 photons s−1, as done in Run 1, we find that the rise time of the shortest admissible burst duration, i.e. 300−1 s, is τmin ∼ 7.2 × 10−4 s. Accordingly, we set δt = 1400−1 s for this simulation run. The input parameters are summarized in Table 3 and the results are presented in Section 6.2.2. Table 3. Input parameters for Run 2, consisting of four separate simulations. The input burst durations T90 are drawn from a lognormal distribution with lower and upper cut-off, respectively, at $$T^{\rm min}_{90}=300^{-1}$$ s and $$T^{\rm max}_{90}=3$$ s. Per simulation run we vary the angles χ, α, where we set χ = α. Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  $${\sim }{\rm LogNormal(\overline{{\it T}}_{90}, \sigma ^2_{{\it T}_{90}})}$$  $$\overline{T}_{90}$$ (s)  0.1  $$\sigma _{T_{90}}$$ (s)  1  δt (s)  1400−1  Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  104  T90 (s)  $${\sim }{\rm LogNormal(\overline{{\it T}}_{90}, \sigma ^2_{{\it T}_{90}})}$$  $$\overline{T}_{90}$$ (s)  0.1  $$\sigma _{T_{90}}$$ (s)  1  δt (s)  1400−1  View Large 5.4 Run 3: burst amplitude distribution In this simulation run we fix the burst duration to T90 = 1 s and draw an amplitude A for each individual burst from a power-law distribution,   $$\frac{\text{d}n}{\text{d}A}\propto A^{-\Gamma }, \ \ \ \text{with}\ \ \ A^{\rm min} < A < A^{\rm max},$$ (34)where Amin and Amax represent the limits of the distribution, and Γ denotes the power-law index. Note that the number of emitted photons during the burst N are linearly proportional to A, such that these are distributed in a similar fashion. In accordance with the observation of the energy distributions of magnetar bursts, we choose Γ = 5/3 (e.g. Cheng et al. 1996). The input parameters are summarized in Table 4 and the results are presented in Section 6.2.3. Table 4. Input parameters for Run 3, consisting of four separate simulations. The input burst amplitudes A are drawn from a power-law distribution [equation (34)], with Amin = 5 × 102 photons s−1, and Amax = 106 photons s−1. Per simulation run we vary the angles χ, α, where we set χ = α. Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  ∼Powerlaw(Γ)  Γ  5/3  T90 (s)  1  δt (s)  200−1  Parameter  Value  n  104  χ, α (°)  30, 45, 60, 90  A (photons s−1)  ∼Powerlaw(Γ)  Γ  5/3  T90 (s)  1  δt (s)  200−1  View Large 6 RESULTS 6.1 Predictions for run 1 To better understand the results of the simulations, we first examine how κ*(ϕ) will affect the burst parameters. In Fig. 6, we plot the predicted phase distributions of the burst parameters (rows) for the three separate burst durations (columns) of Run 1. These curves were obtained by fitting the burst model to the theoretical light curve that results when the input model (see Fig. 5) is modulated, for a given phase, by the appropriate κ*(ϕ) but without taking into account any photon noise or detectability effects (which are treated properly in the full simulations). It gives an idea of the general trends expected, but no idea of the scatter. Furthermore, we only fit modulated burst profiles with a peak rate of ≳600 counts s−1, since ones with lower peak rates will likely go undetected in the simulations. Figure 6. View largeDownload slide Predicted phase distributions of burst parameters (top to bottom) for the three bursts (left to right) studied in the initial simulation run, Run1. From top to bottom: best-fitting burst counts, time difference between burst input t0 and the best-fitting value (Δt0), best-fitting burst rise time τ, best-fitting skewness factor s, and burst duration T90, inferred from the latter τ and s (the input parameters are given by N0, τ0, s0, and T90,0). The distinct colours represent different values for the angles, χ and α. These curves were obtained by fitting the burst model to theoretical burst profiles (Fig. 5) that are modulated by κ* (equation 25; Fig. 3), depending on their phase occurrence. We only proceed to fit the modulated profile if the peak rate is ≳600 counts s−1. Note that for longer burst durations and larger angles, the best-fitting parameters deviate more from their input values. Figure 6. View largeDownload slide Predicted phase distributions of burst parameters (top to bottom) for the three bursts (left to right) studied in the initial simulation run, Run1. From top to bottom: best-fitting burst counts, time difference between burst input t0 and the best-fitting value (Δt0), best-fitting burst rise time τ, best-fitting skewness factor s, and burst duration T90, inferred from the latter τ and s (the input parameters are given by N0, τ0, s0, and T90,0). The distinct colours represent different values for the angles, χ and α. These curves were obtained by fitting the burst model to theoretical burst profiles (Fig. 5) that are modulated by κ* (equation 25; Fig. 3), depending on their phase occurrence. We only proceed to fit the modulated profile if the peak rate is ≳600 counts s−1. Note that for longer burst durations and larger angles, the best-fitting parameters deviate more from their input values. Based on the predicted curves, we expect that the parameter distributions that we obtain from Run 1 will deviate from their input parameters more strongly for longer burst durations T90 and larger angles χ and α. Approaching ϕ = π from below (above), we find that the bursts will appear to occur earlier (later), rise slower (faster), to become more skewed, and last longer, than their input counterparts. Note furthermore that, in contrast to the predicted phase distributions of N, Δt0, s, and T90, the phase distribution of τ is neither symmetric nor perfectly anti-symmetric about ϕ = π. The results of Run 1 are presented in Section 6.2.1. 6.2 Burst properties from simulations 6.2.1 Run 1: initial simulation run In Figs 7, 8, and 9, we plot the phase distributions (left) and parameter densities (right) of the obtained bursts parameters for three separate input burst durations T90,0, respectively, 0.15 s, 1.5 s, and 3.0 s. Table 5 lists the amount of bursts that were identified per configuration. Figure 7. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 0.15 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 7. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 0.15 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 8. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 1.5 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 8. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 1.5 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 9. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 3.0 s. The theoretical predictions for these distributions are shown in Fig. 6. Figure 9. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 1 for an input burst of T90,0 = 3.0 s. The theoretical predictions for these distributions are shown in Fig. 6. Table 5. Number of identified bursts nid for Run 1, per configuration. The input number of bursts for each simulation run was n = 104. We expect that ∼30 of the identified ‘bursts’ simply constitute statistical deviations that exceed the burst identification threshold (given by equation 29). T90 (s)  χ = α (°)  nid  0.15  30  10 034    45  10 022    60  8875    90  6592  1.50  30  10 012    45  10 017    60  9940    30  7438  3.00  30  10 017    45  10 012    60  10 007    90  9553  T90 (s)  χ = α (°)  nid  0.15  30  10 034    45  10 022    60  8875    90  6592  1.50  30  10 012    45  10 017    60  9940    30  7438  3.00  30  10 017    45  10 012    60  10 007    90  9553  View Large As predicted, we find that especially for longer duration bursts and larger angles, the phase dependence of the burst parameters becomes more pronounced. Evidently this is much less the case for bursts with T90 ≪ P – the parameter densities remain strongly peaked around their input values (e.g. Fig. 7). Nevertheless, in those cases around ϕ ∼ π the bursts still go undetected for large values of χ and α, because either no rays extending from the burst patch intersect with the detector during the burst (i.e. the bursts are invisible) or they do not significantly stand out from the background level. The results confirm that when approaching ϕ = π from below (above), the bursts will appear to occur earlier (later), rise slower (faster), and last longer than their input counterparts. Moreover, the predicted asymmetric profile of the rise-time phase distribution (most notably in the parameter densities of Figs 8 and 9) is clearly observed. Looking at the parameter densities, it appears that the rise times of bursts going out of view are more spread out, yet those of bursts coming into view are more clustered. This is in accordance with the predictions; the initial slope of the rise-time phase distribution (from ∼0 − 4π/5) is steeper compared to the final slope (from ∼6π/5 − 2π). The predicted values between ∼4π/5 − 6π/5 are produced less well in the simulations, since the amount of detected photons is minimal around ϕ = π, complicating burst identification and characterization. We find that both in the predictions and, even more so, in the simulations that the majority of observed bursts have τ/τ0 < 1. Since, T90 ∝ τ we also find for most observed bursts that T90/T90,0 < 1, i.e. the bursts seem to last shorter than their input counterparts. In general, the observed scatter is likely due to photon noise effects, which become most significant near ϕ = π. These effects influence the efficacy of the burst-identification algorithm and the observed burst morphology. 6.2.2 Run 2: T90 distribution The results of Run 2, where we draw the burst duration for each individual burst from a lognormal distribution, are shown in Fig. 10. Table 6 lists the number of identified bursts per configuration. The results closely resemble those of Run 1, with T90 = 0.15 s (see Fig. 7). We only find weak phase dependencies that only become noticeable for large values of the angles, χ and α. In Fig. 11, the input (dotted histogram) and best-fitting (solid histogram) T90 distributions are shown for the four separate configurations. Notice the small dearth of short duration bursts in each histogram; short duration bursts contain fewer counts and may therefore be missed by the burst-identification algorithm (step (V) in Section 5.1). We furthermore find that there is a slight excess at T90 ∼ 0.6 s (although not apparent when χ = α = 90°), which is due to the fact that for most observed bursts τ/τ0 < 1 and T90 ∝ τ. Figure 10. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 2, where the input burst durations are drawn from a lognormal distribution (see Fig. 11). Figure 10. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 2, where the input burst durations are drawn from a lognormal distribution (see Fig. 11). Figure 11. View largeDownload slide Burst duration T90 distributions of four separate simulations (with different values for χ and α) of Run 2. The dotted (solid) histograms represent the input (observed) burst duration distributions. The slight dearth of observed short duration bursts, present in each simulation, is simply due to the fact that they consist of fewer counts and are therefore less likely to be identified by the burst-identification algorithm. Figure 11. View largeDownload slide Burst duration T90 distributions of four separate simulations (with different values for χ and α) of Run 2. The dotted (solid) histograms represent the input (observed) burst duration distributions. The slight dearth of observed short duration bursts, present in each simulation, is simply due to the fact that they consist of fewer counts and are therefore less likely to be identified by the burst-identification algorithm. Table 6. Number of identified bursts nid for Run 2, per configuration. The input number of bursts for each configuration was n = 104. χ = α (°)  nid  30  10 028  45  9871  60  7382  90  5971  χ = α (°)  nid  30  10 028  45  9871  60  7382  90  5971  View Large 6.2.3 Run 3: burst amplitude distribution The results of Run 3, where we draw the burst amplitude/number of emitted burst photons for each individual burst from a power-law distribution, are presented in Fig. 12. Table 7 lists the number of identified bursts per configuration. We find much less spread in the phase distributions of the parameters, compared to e.g. the results from the 1.5 s burst in Run 1 (Fig. 8), however we do observe a considerable amount of scatter. The latter is likely due to the fact that the majority of input bursts (∼0.87) are low-amplitude bursts, i.e. A ≲ 104 photons s−1, which are more difficult to characterize, i.e. their morphology is relatively heavily affected by Poisson noise. Fig. 13 displays the input and observed burst amplitude distributions. Despite a slight offset at larger angles, we find that the slope of the distributions is reproduced by the observed bursts. Input bursts with an amplitude ≲103 photons s−1 may go unidentified as they will likely fall below the significance threshold of the burst identification algorithm. Figure 12. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 3, where the input burst amplitudes are drawn from a power-law distribution (see Fig. 13). Figure 12. View largeDownload slide Phase distributions (left) and parameter densities (right) of burst parameters of Run 3, where the input burst amplitudes are drawn from a power-law distribution (see Fig. 13). Figure 13. View largeDownload slide Burst amplitude distributions of four separate simulations (with different values for χ and α) of Run 3. The dotted (solid) histograms represent the input (observed) burst amplitude distributions. The cut-off of observed low-amplitude bursts (at ≲103 photons s−1) is due to the fact that the amplitude of these bursts likely occurs below the threshold of the burst-identification algorithm. Figure 13. View largeDownload slide Burst amplitude distributions of four separate simulations (with different values for χ and α) of Run 3. The dotted (solid) histograms represent the input (observed) burst amplitude distributions. The cut-off of observed low-amplitude bursts (at ≲103 photons s−1) is due to the fact that the amplitude of these bursts likely occurs below the threshold of the burst-identification algorithm. Table 7. Number of identified bursts nid for Run 3, per configuration. The input number of bursts for each configuration was n = 104. χ = α (°)  nid  30  6645  45  5831  60  4734  90  3717  χ = α (°)  nid  30  6645  45  5831  60  4734  90  3717  View Large 6.3 Detectability of burst phase dependence Here we set out to test the main method used in studies of burst phase dependence to date, see Section 2. During the simulations we determine and record the phase occurrence of the burst peak $$\phi _0^{\rm bf}$$, i.e. the phase occurrence of the best-fitting burst peak time $$t_0^{\rm bf}$$. After each burst we compile a distribution of the values for $$\phi _0^{\rm bf}$$ of all previous bursts up to the most recent one, and compare this burst phase occurrence distribution to a uniform distribution using a Kolmogorov–Smirnov (K–S) test, from which we obtain a p-value. We set the significance threshold at a p-value of 0.003,corresponding to a 0.3 per cent probability that the observed burst peak phase occurrences are distributed uniformly across phase. To be clear, we are simulating emission from a fixed point on the NS surface, from which bursts are being emitted at random rotational phase. The naive expectation that this would result in an observable phase dependence most from the expected modulation of the intensity (see Fig. 3) that results, for example, in missing some bursts emitted on the dark side of the star. In Fig. 14 we plot the evolution of the p-value against the number of observed bursts up to that point in the simulation for the three separate burst durations (from top to bottom) of Run 1. The different curves per subplot correspond to different values of χ and α. The horizontal dashed lines denote the threshold level and the vertical dotted lines indicate the number of bursts at which the p-value of a given simulation drops below the threshold level. Note that of these 12 simulations, the p-value does not fall below the threshold before 103 bursts for χ = α ≤ 45°. Nevertheless, we do find a decreasing trend for T90 = 3.0 s after ∼500 bursts, reaching the threshold at ∼104 bursts (the length of the simulation), for those angles. Remarkably, the p-value associated with the simulation where T90 = 1.50 s and χ = α = 60° does not show a decreasing trend before 104 bursts. The remaining configurations do drop below the threshold fairly soon, i.e. after ∼20 − 150 bursts. Figure 14. View largeDownload slide Evolution of the p-value against the number of bursts for three separate burst durations (top to bottom) of Run 1. The distinct curves per subplot represent different values for the angles χ and α. The horizontal dashed line denotes the threshold level. The vertical dotted lines denote the number of bursts at which the p-value drops below the threshold (see the text for more details). Figure 14. View largeDownload slide Evolution of the p-value against the number of bursts for three separate burst durations (top to bottom) of Run 1. The distinct curves per subplot represent different values for the angles χ and α. The horizontal dashed line denotes the threshold level. The vertical dotted lines denote the number of bursts at which the p-value drops below the threshold (see the text for more details). To determine the minimum number of bursts required to guarantee that the p-value drops below the threshold, we ran the simulation per configuration $$\mathcal {N}_{\rm s}$$ times, each time until the threshold was reached, and recorded the number of bursts. Fig. 15 shows the evolution of the p-value for $$\mathcal {N}_{\rm s}$$ = 10, 100, and 400 simulations, with T90 = 0.15 s and χ = α = 90°. We found that the maximum obtained p-value (out of all $$\mathcal {N}_{\rm s}$$ simulations for a given configuration) after a specific number of bursts decreases at a certain rate (denoted by the black markers). In Fig. 16, we plot the maximum p-value attained, over $$\mathcal {N}_{\rm s}=400$$ simulations per configuration, against the number of bursts. Subsequently, we fit a straight line to the decreasing trends of the log of the p-value and record the number of bursts at which these lines intersect with the threshold level. Accordingly, we find an estimate for the minimum number of bursts nmin at which, assuming a certain configuration, the observed $$\phi _0^{\rm bf}$$ distribution should deviate significantly from a uniform distribution. If the $$\phi _0^{\rm bf}$$ distribution does not significantly deviate from a uniform distribution after nmin bursts, then the configuration will likely be such that the modulation in intensity is less strong than assumed. The latter is dependent on assumptions on the parameters that determine the shape of κ(θ0) (equation 24). Figure 15. View largeDownload slide Evolution of the p-value for $$\mathcal {N}_{\rm s}$$ simulations for T90 = 0.15 s and χ = α = 90°. From left to right we increase the value of $$\mathcal {N}_{\rm s}$$ from 10, to 100, to 400. The distinct coloured curves denote the evolution of the p-value for separate simulations. The black markers denote the maximum p-value that was attained (out of all simulations) for a given number of bursts. The horizontal dashed line denotes the threshold value, the cyan dash-dotted curve represents the fit to the decreasing trend of the log of the maximum p-values, and the vertical dotted curve denotes the intersection of the fit with the threshold level. The latter occurs, from left to right, at nmin = 86, 95, and 104, respectively. The maximum p-values in the panel on the right are also plotted in the top panel of Fig. 16. Figure 15. View largeDownload slide Evolution of the p-value for $$\mathcal {N}_{\rm s}$$ simulations for T90 = 0.15 s and χ = α = 90°. From left to right we increase the value of $$\mathcal {N}_{\rm s}$$ from 10, to 100, to 400. The distinct coloured curves denote the evolution of the p-value for separate simulations. The black markers denote the maximum p-value that was attained (out of all simulations) for a given number of bursts. The horizontal dashed line denotes the threshold value, the cyan dash-dotted curve represents the fit to the decreasing trend of the log of the maximum p-values, and the vertical dotted curve denotes the intersection of the fit with the threshold level. The latter occurs, from left to right, at nmin = 86, 95, and 104, respectively. The maximum p-values in the panel on the right are also plotted in the top panel of Fig. 16. Figure 16. View largeDownload slide Maximum p-value attained, out of $$\mathcal {N}_{\rm s}$$ simulations per configuration, against the number of bursts for the three separate burst durations (top to bottom) of Run 1. The horizontal dashed line denotes the threshold level, the blue dash-dotted curves represent the fits to the decreasing trends, and the vertical dotted lines denote the number of bursts where the fits intersect with the threshold level. Figure 16. View largeDownload slide Maximum p-value attained, out of $$\mathcal {N}_{\rm s}$$ simulations per configuration, against the number of bursts for the three separate burst durations (top to bottom) of Run 1. The horizontal dashed line denotes the threshold level, the blue dash-dotted curves represent the fits to the decreasing trends, and the vertical dotted lines denote the number of bursts where the fits intersect with the threshold level. For the burst durations and configurations that we study in Run 1, we find for χ = α = 90° that nmin ∼ 100 bursts. For χ = α = 60°, we find nmin = 1446 bursts and nmin = 296 bursts, for T90 = 0.15 s and T90 = 3.0 s, respectively. Yet, we do not find an nmin for T90 = 1.50 s, since the attained maximum p-value does not exhibit a decreasing trend before 103 bursts; consistent with the simulation run displayed in the middle panel of Fig. 14. This is because the burst spot remains (partially) visible throughout the NS's rotation, such that enough counts can be detected for the duration of the burst, and the fact that the Δt0 remains comparatively small, i.e. the corresponding parameter density comprises a narrow peak (Fig. 8), in contrast to e.g. the parameter density of the 3.0 s burst, which is much more spread out (Fig. 9). 7 DISCUSSION AND CONCLUSION We have studied, from a theoretical perspective, the conditions under which magnetar bursts from a predefined localized active region or burst patch on the NS surface would give rise to a detectable phase dependence. By adopting a straightforward input burst model, we were able to examine the changes in the observed bursts after they were modulated by the phase-dependent function κ*(ϕ), which takes into account the effects of gravitational light bending and depends on the configuration of the system. We found that the degree to which the inferred burst properties become phase dependent is strongly contingent on the duration of the bursts and geometry of the system; we find a stronger phase dependency of the burst properties for longer duration bursts and larger values of the angles χ and α. The former is because longer bursts sample a wider range of photon trajectories and the latter is due to the fact that for larger values of χ and α the GR effects become more significant. Furthermore, the majority of observed bursts turn out to have τ/τ0 < 1 and T90/T90,0 < 1, i.e. they rise faster and appear shorter than their input counterparts. Attempts to infer the properties of individual bursts with durations greater than ∼10−20  per cent of the spin period should certainly take into account potential distortion due to phase-dependent effects. Adopting a lognormal burst duration distribution that peaks at $$\overline{T}_{90}=0.1$$ s (as observed for well-sampled sources), from which we draw the input duration for each individual burst, we found that phase distributions of the parameters closely resembled those of Run 1, for which T90 = 0.15 s. When considering a power-law distribution for the input burst amplitudes and burst duration of T90 = 1 s, we observed a weak phase dependency of the burst parameters and a considerable amount of scatter, which in turn is caused by the large fraction of low-amplitude input bursts, which are more affected by Poisson noise. We conclude that the observed distributions of burst properties from well-sampled sources are likely not strongly distorted due to phase-dependent effects, by virtue of being dominated by short bursts. We studied the detectability of phase dependence, using the most commonly-used measure (see Section 2) whereby one concentrates on the phase occurrence of the burst peaks. In our set-up all bursts originate at a specific small active region – in some respects the most extreme phase-dependent scenario. However, rotational phase dependence of the peak occurrences was not always apparent. We found that one would require a minimum number of bursts for certain input burst properties and a given system configuration, to guarantee observing a phase dependence. Only in the case of the most extreme geometries, i.e. χ = α > 60°, does this approach the burst sample sizes that were examined in the literature, which range from tens to several hundred bursts. Studies that have not found a phase dependence in the distribution of the burst peak occurrences as yet (e.g. Savchenko et al. 2010; Scholz & Kaspi 2011; Collazzi et al. 2015), might simply require a larger burst sample in order to rule it out for certain geometries. For other geometries, however, it will never be possible to rule out the presence of a burst phase dependence. In our study we have considered only a restricted range of scenarios, where the emission region is tied to the stellar surface. One factor that we have not simulated in detail is that of any potential beaming of the burst emission. To offer brief insights for such influences, we show in Fig. 17 the theoretical phase distributions of the observed burst properties in the case of beamed emission (equation 22 with σb = π/6) from a burst with T90 = 1.5 s; the corresponding shape of κ*(ϕ) is shown by the dashed curves in Fig. 3. We compare it to its isotropic counterpart and find that the phase dependency of the burst properties will be enhanced in the presence of beaming. Figure 17. View largeDownload slide Predicted phase distributions of burst parameters in the presence of beamed emission for a burst of T90 = 1.5 s (solid curves). Beaming is described by equation (22), where we have set the beaming width σb to π/6. For comparison, the dotted curves represent the burst parameter distributions for an isotropic burst with the same duration. Figure 17. View largeDownload slide Predicted phase distributions of burst parameters in the presence of beamed emission for a burst of T90 = 1.5 s (solid curves). Beaming is described by equation (22), where we have set the beaming width σb to π/6. For comparison, the dotted curves represent the burst parameter distributions for an isotropic burst with the same duration. The detectability of a burst phase dependence depends strongly on the shape of κ*(ϕ): the stronger the variation with ϕ the greater the modification to the input burst profiles. Introducing additional bursts patches or allowing for active regions to occur at a certain height above the surface will cause the phase dependence of κ* to decrease. A burst phase dependence in those cases may then only become detectable if the emission is also strongly beamed. In this paper we have not studied the method whereby a phase dependence is searched for in the epoch-folded photon times of arrival. This method can, and should, be subjected to the same level of scrutiny. An additional challenge with this method, however, is to determine a proper false alarm rate. Straightforwardly looking for deviations from uniformity of the times of arrival does not work, since a single burst already consists a significant departure. One must instead quantify the conditions under which one would detect a burst photon phase dependence, even if the bursts originated at random locations on or above the NS surface. We defer this topic to future studies. ACKNOWLEDGEMENTS C.E. and A.L.W. acknowledge support from NOVA, the Dutch Top Research School for Astronomy. D.H. was partially supported by the Moore-Sloan Data Science Environment at New York University, and the James Arthur Postdoctoral Fellowship at New York University. D.H. acknowledges support from the DIRAC Institute in the Department of Astronomy at the University of Washington. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. We thank Matthew Baring for useful comments. Footnotes 1 Dib et al. (2009) discuss the four bursts from AXP 1E 1048.1–5937 and note that only three of them occur near pulse maximum, whereas the fourth burst does not. 2 No further details of the analysis procedure on the SGR 1806–20, such as the applied nominal threshold to determine significance, are given in the article. The phase-dependence analysis procedure of the SGR 1900+14 data is also not described. 3 Note that here we only consider the case where burst emission escapes from the system at the stellar surface (as it would from a trapped fireball, e.g., due to the reduced scattering opacity close to the surface; Thompson & Duncan 1995). It is for surface emission that the effects of GR will be most significant. We argue that bursts that occur high-up in the magnetosphere will be much less affected by the effects of GR or occultation of the star itself, and thus may exhibit weak to no phase-dependent properties. In effect, we are considering the most optimistic case for the detection of phase-dependent effects. 4 If no burst photons are detected, i.e. Ndet = 0, we move on to the next burst [step (II)]. 5 If $$N^{\rm init}<y_0^{\rm sig}$$, we set $$N^{\rm init}=(1+s^{\rm init})\tau ^{\rm init}y_0^{\rm sig}/\delta t$$. REFERENCES An H. et al.  , 2014, ApJ , 790, 60 CrossRef Search ADS   Cheng B., Epstein R. I., Guyer R. A., Young A. C., 1996, Nature , 382, 518 CrossRef Search ADS   Collazzi A. C. et al.  , 2015, ApJS , 218, 1 CrossRef Search ADS   Dib R., Kaspi V. M., Gavriil F. P., 2009, ApJ , 702, 614 CrossRef Search ADS   Gavriil F. P., Kaspi V. M., Woods P. M., 2002, Nature , 419, 142 CrossRef Search ADS PubMed  Gavriil F. P., Kaspi V. M., Woods P. M., 2004, ApJ , 607, 959 CrossRef Search ADS   Gavriil F. P., Kaspi V. M., Woods P. M., 2006, ApJ , 641, 418 CrossRef Search ADS   Gavriil F. P., Dib R., Kaspi V. M., 2011, ApJ , 736, 138 CrossRef Search ADS   Gourgouliatos K. N., Kondić T., Lyutikov M., Hollerbach R., 2015, MNRAS , 453, L93 CrossRef Search ADS   Göğüş E., Kouveliotou C., Woods P. M., Thompson C., Duncan R. C., Briggs M. S., 2001, ApJ , 558, 228 CrossRef Search ADS   Huppenkothen D. et al.  , 2015, ApJ , 810, 66 CrossRef Search ADS   Kaspi V. M., Beloborodov A. M., 2017, ARA&A , 55, 261 CrossRef Search ADS   Kouveliotou C. et al.  , 1998, Nature , 393, 235 CrossRef Search ADS   Kouveliotou C. et al.  , 1999, ApJ , 510, L115 CrossRef Search ADS   Lander S. K., Andersson N., Antonopoulou D., Watts A. L., 2015, MNRAS , 449, 2047 CrossRef Search ADS   Lin L. et al.  , 2012, ApJ , 756, 54 CrossRef Search ADS   Lyutikov M., 2002, ApJ , 580, L65 CrossRef Search ADS   Morsink S. M., Leahy D. A., Cadeau C., Braga J., 2007, ApJ , 663, 1244 CrossRef Search ADS   Mus S. S, Göǧüş E., Kaneko Y., Chakraborty M., Aydin B., 2015, arXiv.org, p. 42 Palmer D. M., 1999, ApJ , 512, L113 CrossRef Search ADS   Palmer D. M., 2002, Mem. Soc. Astron. Italiana , 73, 578 Pechenick K. R., Ftaclas C., Cohen J. M., 1983, ApJ , 274, 846 CrossRef Search ADS   Savchenko V., Neronov A., Beckmann V., Produit N., Walter R., 2010, A&A , 510, A77 CrossRef Search ADS   Scholz P., Kaspi V. M., 2011, ApJ , 739, 94 CrossRef Search ADS   Thompson C., Duncan R. C., 1993, The ApJ , 408, 194 CrossRef Search ADS   Thompson C., Duncan R. C., 1995, MNRAS , 275, 255 CrossRef Search ADS   Thompson C., Yang H., Ortiz N., 2017, ApJ , 841, 54 CrossRef Search ADS   Turolla R., Zane S., Watts A. L., 2015, Rep. Progr. Phys. , 78, 116901 CrossRef Search ADS   van Putten T., Watts A. L., Baring M. G., Wijers R. A. M. J., 2016, MNRAS , 461, 877 CrossRef Search ADS   Woods P. M., Thompson C., 2006, in Lewin W., van der Klis M., eds, Cambridge Astrophysics Series, Compact Stellar X-ray Sources. Cambridge Univ. Press, p. 547 Woods P. M. et al.  , 2005, ApJ , 629, 985 CrossRef Search ADS   Zhu C., Byrd R. H., Lu P., Nocedal J., 1997, ACM Trans. Math. Softw. , 23, 550 CrossRef Search ADS   APPENDIX: PARAMETER TABLE Here we list brief descriptions of the simulation parameters and their associated symbols (see Table A1). Table A1. A table of the simulation parameters appearing in Sections 5 and 6. Symbol  Description  n  Number of input bursts  nid  Number of identified bursts  nmin  Number of bursts at which the p-value is    Guaranteed to drop below the threshold  R  Neutron star radius  P  Neutron star rotation period  ψ  Size of the burst patch  χ  Angle between NS axis of rotation and    the line-of-sight  α  Colatitude of the burst patch  ζ  Background rate  δt  Light-curve bin width  b  Background level  N  Number of emitted photons  Ndet  Number of detected photons  t0  Burst peak time  ϕ0  Burst peak phase occurrence  τ  Rise-time  s  Skewness factor  A  Burst amplitude  $$\mathcal {N}_{\rm bins}$$  Number of time bins  ysig  Maximum departure amplitude  $$t_0^{\rm sig}$$  Time bin with ysig  ΔT  Time interval of an observed burst  tin, tout  Limits of ΔT  μ  Mean count level ∼b  μ*  Running mean  ΔTinterval  Time interval over which μ* is estimated  Δt0  Difference input and best-fitting burst peak    time: Δt0 ≡ tbf − t0  T90  Burst duration  $$\overline{T}_{90}$$  Mode of the duration distribution  $$\sigma _{T_{90}}$$  Width of the duration distribution  $$T_{90}^{\rm min}$$  Minimum burst duration  $$T_{90}^{\rm max}$$  Maximum burst duration  Γ  Power-law index of the amplitude    distribution  Amin  Minimum burst amplitude  Amax  Maximum burst amplitude  $$\mathcal {N}_{\rm s}$$  Number of simulations  subscript ‘0’  Input parameter  superscript ‘init’  Initial guess  superscript ‘bf’  Best-fitting parameter  Symbol  Description  n  Number of input bursts  nid  Number of identified bursts  nmin  Number of bursts at which the p-value is    Guaranteed to drop below the threshold  R  Neutron star radius  P  Neutron star rotation period  ψ  Size of the burst patch  χ  Angle between NS axis of rotation and    the line-of-sight  α  Colatitude of the burst patch  ζ  Background rate  δt  Light-curve bin width  b  Background level  N  Number of emitted photons  Ndet  Number of detected photons  t0  Burst peak time  ϕ0  Burst peak phase occurrence  τ  Rise-time  s  Skewness factor  A  Burst amplitude  $$\mathcal {N}_{\rm bins}$$  Number of time bins  ysig  Maximum departure amplitude  $$t_0^{\rm sig}$$  Time bin with ysig  ΔT  Time interval of an observed burst  tin, tout  Limits of ΔT  μ  Mean count level ∼b  μ*  Running mean  ΔTinterval  Time interval over which μ* is estimated  Δt0  Difference input and best-fitting burst peak    time: Δt0 ≡ tbf − t0  T90  Burst duration  $$\overline{T}_{90}$$  Mode of the duration distribution  $$\sigma _{T_{90}}$$  Width of the duration distribution  $$T_{90}^{\rm min}$$  Minimum burst duration  $$T_{90}^{\rm max}$$  Maximum burst duration  Γ  Power-law index of the amplitude    distribution  Amin  Minimum burst amplitude  Amax  Maximum burst amplitude  $$\mathcal {N}_{\rm s}$$  Number of simulations  subscript ‘0’  Input parameter  superscript ‘init’  Initial guess  superscript ‘bf’  Best-fitting parameter  View Large © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations