The detection of polarized x-ray emission from the magnetar 1E 2259+586

We report on IXPE, NICER and XMM-Newton observations of the magnetar 1E 2259+586. We find that the source is significantly polarized at about or above 20% for all phases except for the secondary peak where it is more weakly polarized. The polarization degree is strongest during the primary minimum which is also the phase where an absorption feature has been identified previously (Pizzocaro et al. 2019). The polarization angle of the photons are consistent with a rotating vector model with a mode switch between the primary minimum and the rest of the rotation of the neutron star. We propose a scenario in which the emission at the source is weakly polarized (as in a condensed surface) and, as the radiation passes through a plasma arch, resonant cyclotron scattering off of protons produces the observed polarized radiation. This confirms the magnetar nature of the source with a surface field greater than about 1015G


INTRODUCTION
discovered the X-ray source 1E 2259+586 in observations of the supernova remnant G109.l−1.0 using the Einstein telescope on 17 December 1979.Later identified as one of the Anomalous X-ray Pulsars (AXPs; Mereghetti & Stella 1995), 1E 2259+586 was the first member of this class to be discovered (and the third magnetar after SGR 1806−20 and SGR 0525−66 earlier in 1979).Fahlman & Gregory (1981) identified a periodicity in the X-ray emission from the source of 3.4890 s.Subsequent observations revealed that the pulse profile exhibits two similar peaks over each rotation period.Similar to the other anomalous Xray pulsars the spin period of 1E 2259+586 is gradually increasing (  ≈ 5 × 10 −13 s s −1 , Dib & Kaspi 2014), and this is associated with the presence of a strong magnetic field (≈ 6 × 10 13 G) and the braking of the pulsar through magnetic dipole radiation.Although the dipole magnetic field inferred from the spin down lies at the low end of the magnetar range (in fact several radio pulsars have larger spin-down fields, Manchester et al. 2005), 1E 2259+586 exhibits bursts, glitches and even an anti-glitch similar to other magnetars (Kaspi & Beloborodov 2017), suggesting that a much stronger field of 10 14 -10 15 G, confined in small scale structures close to the surface, is present in 1E 2259+586, similarly to the other "low-field" magnetar SGR 0418+5729 (Tiengo et al. 2013).Pizzocaro et al. (2019) found evidence of phase-dependent spectral features in XMM-Newton observations of 1E 2259+586 that have been consistent in phase and energy over more than a decade while the source has been in quiescence.On the other hand, during an active epoch in 2002, although the spectral feature was still present, its position in energy as a function of phase had changed.Pizzocaro et al. (2019) proposed that this feature could be a signature of resonant cyclotron scattering similar to what has been proposed for a similar feature found in the phase-resolved spectrum of SGR 0418+5729 (Tiengo et al. 2013).If this association is correct in the case of 1E 2259+586, the magnetic field for a proton cyclotron resonance is (3 − 16) × 10 14 G (and below 10 12 G for an electron cyclotron resonance).Here we present polarimetric observations that provide evidence that this spectral feature results from photons scattering off of non-relativistic particles (most likely protons) at the cyclotron resonance.

OBSERVATIONS
IXPE observed 1E 2259+5586 in June-July 2023 and a simultaneous XMM-Newton pointing was carried out on June 30 2023.NICER archival data, partially overlapping the IXPE time window, were also available.Tab. 1 summarizes the observations used for this analysis.

IXPE
IXPE, a NASA mission in partnership with the Italian space agency (ASI; Weisskopf et al. 2022, and references therein) observed 1E 2259+5861 on 2023 June 2-19 and June 30-July 6 for 1.2 Ms in total.The gas-pixel detectors on IXPE register the arrival time, sky position, and energy for each X-ray photon and use the photoelectric effect to provide an estimate of the position angle of each photon (Soffitta et al. 2021).During each observation, photon arrival events registered between energies of 2 and 8 keV, within 48 arcseconds of the position of source (R.A. = 345 • .3,DEC = 58 • .9),were extracted for analysis.The background was estimated from an annular region centered on the source of inner and outer radii of 78 arcseconds and 240 arcseconds, respectively.Background subtraction was applied to the extracted Stokes parameters.Anyway, we observe that the background level for each IXPE detector unit (DU, i.e. telescope) in the 2-8 keV band was less than 2% of the source one.We could, hence, conclude that the energy-integrated analysis we performed is not much affected by background effects.The times of the photon arrivals were finally corrected for the motion of IXPE around the barycenter of the Solar System.

NICER
We used the Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2012;LaMarr et al. 2016;Prigozhin et al. 2016) data collected over 2023 March 19 to 2023 July 24.The data were processed with HEASoft version 6.31 and NICER Data Analysis Software (NICER-DAS) version 10 (2022-12-16_V010a) using the nicerl2 tool with standard filtering criteria, resulting in 6.1 ks of filtered exposure.We performed barycenter corrections in the ICRS reference frame using the JPL DE421 Solar System ephemeris with the barycorr tool in FTOOLS with coordinates R.A. = 345.• 2845, DEC = 58.• 879.

XMM-Newton
A DDT pointing of 1E 2259+586 with XMM-Newton was activated on 2023, June 30, starting at 23:47:46 UTC for an exposure time of ≈ 20 ks.The EPIC-pn (Strüder et al. 2001) as well as the two MOS cameras (Turner et al. 2001) were set in the Small Window mode, with a time resolution of 0.3 s.Row data were processed by means of the SAS version 20.0 and the most updated calibration files.After the subtraction of the intervals in which the background events were dominant, the data were extracted and processed applying standard procedures, for a net exposure of ≈ 19.1 ks for the MOSs and ≈ 18.8 ks for the EPIC-pn.We extracted the source counts from a circular region of radius about 65".Those of the background were extracted from a similar region, within the same CCD where the source lies and ∼2.5' away for the pn, while for the MOSs, due to the use of the small window mode, from another CCD (at a distance of about 9' from the source).The times of the extracted photons were corrected for the barycenter of the Solar System in the ICRS reference frame with barycen tool in the SAS.The background subtracted source count rates were 10.71(3) ct s −1 in the pn and 3.31(1) in the MOSs (1 confidence levels are reported).

Timing analysis
We first started analysing the datasets from each missions.In order to estimate the most reliable timing solution for the 1.2 Ms exposure IXPE dataset, we used a phase-fitting approach (see, for example, Dall'Osso et al. 2003) which gave a period of 6.979281(1) s, or  = 0.14328124(3) Hz, at the reference epoch of 60097.0MJD (a further  component did not significantly improve the fit).The peakto-peak semi-amplitude of the background subtracted light curves folded to the above period resulted to be (33 ± 3)%.For NICER, the same phase-fitting algorithm revealed that a quadratic component was significantly present in the spin period phases as a function of time, resulting in the following timing solution:  = 6.9792783(1)s and  = 5.0(2)×10 −13 s s −1 , reference epoch 60022.0MJD (1 confidence levels are reported), corresponding to  = 0.143281290(3) Hz and  = −1.03(4)× 10 −14 Hz s −1 .Similarly, for XMM-Newton the best timing solution was inferred to be  = 6.97931(2) s or  = 0.1432806(4) Hz at reference epoch 60126.0MJD.The inclusion of a first period derivative  did not significantly improve the fit.The peak-to-peak semi-amplitude of the background subtracted light curves folded to the above period resulted to be (35±3)%.Note that the three timing solutions are in agreement with each other within their uncertainties.Finally, the whole sample of IXPE, NICER and XMM-Newton datasets was used simultaneously to provide the best possible timing solution.A phase-fitting analysis (see Figure 1) gave the following result,  = 6.9792785(1)s and  = 4.7(1) × 10 −13 s s −1 , reference epoch 60022.0MJD, corresponding to  = 0.143281286(2) Hz and  = −9.7(3)× 10 −15 Hz s −1 (1 c.l.; reduced  2 ∼2 for 14 d.o.f. and rms of 0.007 cycles).In Figure 2 we show the light curves of each mission folded to the best solution discussed above.The pulse shape is double peaked and does not change, within uncertainties, considering different energy bands.Figure 3 shows the IXPE pulse profile with the different phase ranges used in the subsequent analysis: Big Dip, Rise, Big Peak, Little Dip and Little Peak (see Section 4 for further details); phase zero was chosen to coincide with that of Pizzocaro et al. (2019).

Spectral analysis
We first performed a phase-integrated spectral analysis of the XMM-Newton observation by fitting simultaneously the EPIC-pn and MOS data in the 0.5-8 keV energy range using XSPEC (Arnaud 1996).Fits with (absorbed 1 ) single-component models (either a blackbody, BB, or a power-law, PL) turned out to be rather unsatisfactory, with a 1 XSPEC model phabs.reduced  2 exceeding 6 for 250 degrees of freedom (dof).A substantial improvement was obtained adding a second component.A purely thermal model (BB+BB), however, still resulted in a poor fit, with  2 = 462.2for 248 dof and a temperature for the hotter blackbody of ∼ 300 keV, difficult to reconcile with the known properties of magnetars (see e.g.Turolla et al. 2015;Kaspi & Beloborodov 2017, for reviews).By adopting a BB+PL decomposition the fit improved, although it was still far from being statistically acceptable, and the best-fitting parameters are compatible with those presented in previous works (Zhu et al. 2008;Pizzocaro et al. 2019).The addition of a Gaussian absorption line, (gabs in XSPEC, as in Pizzocaro et al. 2019), resulted in a further improvement in the quality of the fit (2 = 341.0for 245 dof).By performing an f-test, the probabil-ity that the additional absorption line is unnecessary turns out to be 4.3 × 10 −7 (i.e. the feature is significant at ∼ 5 confidence level).The corresponding BB temperature, PL photon index, line energy and width are in agreement with those reported in Pizzocaro et al. (2019) within 1 confidence level (see Table 2 for the fit parameters).The large  2 for the pn+MOS fit indicates that fitting the merged dataset to the phase-average spectrum may be inadequate.Indeed, by restricting to the EPIC-pn data only resulted in a much better fit ( 2 = 94.1 for 93 dof for the same spectral model, see Figure 4), although the line centroid shifted to slightly higher energies.Results of the phase-resolved spectroscopic analysis are presented in Section 5.
We finally attempted a fit of the IXPE 2-8 keV data using the same spectral decomposition and freezing the column density and line parameters to those obtained from the fit of the EPIC-pn data.The fit is statistically acceptable ( 2 = 138.0for 147 dof) with values of the free parameters in agreement (within statistical uncertainties) with those already obtained in the previous analyses (see again Table 2 and Figure 4).

Polarization analysis
A phase-and energy-integrated study (in the 2-8 keV band) of the polarization properties of the source was performed using the PCUBE algorithm of the ixpeobssim suite (Baldini et al. 2022) 2 .The results for the normalized Stokes parameters / and / are reported in Figure 5 for the single IXPE DUs and for the sum of them.In the figure the loci of constant polarization degree (PD = √︁ (/) 2 + (/) 2 ) and polarization angle (PA = arctan(/)/2) are also shown, together with the value of the minimum detectable polarization at 99% confidence level (MDP 99 ; Weisskopf et al. 2010).We obtained a highly probable detection (significance > 99.9%), with PD = 5.6 ± 1.4% (above the MDP 99 = 4.5%) and PA = −75 • .2± 7 • .4measured East of North (errors at 1).
We performed a phase-integrated, energy-dependent polarimetric analysis as well, by dividing the 2-8 keV band into 3 bins.The only bin with a polarization degree in excess of the MDP 99 is that at low energies (2.0-3.2 keV), with PD = 6.1 ± 1.5% (MDP 99 = 4.6%) and PA = −66 • .4± 7 • .1 (errors at 1).In the rest of the energy range we can derive only upper limits: PD < 14.6% in the 3.2-5.0keV range and PD < 70.0% in the 5.0-8.0 keV one, at 3 confidence level.However, the null hypothesis probability that no polarization is detected is below 3 × 10 −4 .
In order to study the evolution of the polarization with the rotational phase, we divided the counts into 14 equally spaced phase bins taking as phase zero the one reported in Pizzocaro et al. (2019).
In each bin we calculate the normalized Stokes parameters / and / using the likelihood method outlined in González-Caniulef et al. (2023); similar results can be achieved using the standard tools (e.g.ixpeobssim).The results are given in Table 3.There (as well as in the central panel of Figure 7) we listed the values of polarization degree and polarization angle even for the bins where the PD lies below the MDP 99 .This is justified because we measured an overall polarization with high significance (see above), so that the null hypothesis relevant here is not that of unpolarized radiation, but rather of constant polarization.This allows us to probe the agreement between the binned results and the polarization models which are fit directly to the data for individual photons (González-Caniulef et al. 2023).
The Stokes parameters for the bins reported in Table 3 are also

POLARIZATION MODELING
In the highly magnetized environment surrounding a magnetar radiation propagates into two (normal) modes, the ordinary (O) and extraordinary (X) ones.In the former case, the electric field of the wave oscillates in the plane of the (local) magnetic field and of the photon momentum, while in the latter the oscillations are perpendicular to this plane (Gnedin et al. 1978;Pavlov & Shibanov 1979).
Even in the absence of matter, vacuum birefringence will force the polarization vector of photons to follow the direction of the (local) magnetic field until the so-called polarization-limiting radius (Heyl & Shaviv 2000, 2002).For typical magnetars, this radius is estimated to be about 200-300 stellar radii for keV photons (see, e.g., Figure 1 of Taverna et al. 2015, andalso Heyl &Caiazzo 2018), where the field is dominated by the dipole component.The polarization measured at the telescope is, then, expected to be either parallel or perpendicular to the instantaneous projection of the magnetic dipole axis of the star onto the plane of the sky.For this reason, the modulation of the polarization angle with phase is decoupled from the evolution of the polarization degree and intensity (that carry the imprint of the conditions at emission) and most likely should follow the rotating vector model (RVM; Radhakrishnan & Cooke 1969; Poutanen 2020, see also Taverna et al. 2022;González-Caniulef et al. 2023), where  p is the inclination of the spin axis with respect to the lineof-sight,  p the position angle of the spin axis in the plane of the sky (measured East of North),  the inclination of the magnetic axis to the spin axis and  is the spin phase ( 0 is the initial phase).The angle between the dipole axis and the line-of-sight varies between  p −  at  =  0 and  p +  at  =  0 + .Without loss of generality we restrict the parameters as follows: The two best-fitting RVMs for the polarization angle are depicted in Figure 7.The solid curve traces a model where the polarization mode is constant with phase, and the dashed curve shows a model where the polarization mode switches at a phase where the polarization degree is low.This is accomplished by replacing the Stokes parameters ( and ) of the model by their additive inverse over a range of phases  1 <  <  2 where  1 and  2 are two additional parameters.The log-likelihood of the first model (with PD = 11.6%) is 49.6 and that of the second one (PD = 12.7%) is 57.8.The loglikelihoods for 222,043 events, drawn from two models with the observed values of PDs, are distributed approximately normally, with means of 48 and 56 and standard deviations of 10; this indicates that both models are good fits to the data; about 60% of the time random events drawn from these models will yield likelihoods smaller than those measured for the data.Figure 8 depicts posterior distributions of the parameters for these two models.
As the model without mode switching is nested within the modeswitching model, we can calculate the probability to achieve the measured likelihood ratio even if there is no mode switch (the null hypothesis).Twice the difference in likelihoods is distributed as a      2 distribution with two degrees of freedom (for the two additional parameters), yielding a probabilty that the null hypothesis is true of less than 3 × 10 −4 .Furthermore, one of the mode switches can account for the low polarization at phase 0.11 that lies between two high polarization regimes.Figure 9 examines the single-mode RVM in more detail by removing the effect of the motion of the magnetic axis on the plane of the sky from the measured photon polarization angles.To this aim, the  and  Stokes parameters for each photon are rotated into the frame of the best-fitting RVM before measuring the polarization angle and degree.If the low polarization at phase 0.11 results from the smearing of a large intrinsic polarization as the star rotates, the polarization measured after this procedure would be large.However, as Figure 9 shows, the polarization at this phase remains low, so it is indeed a natural time for a mode switch as found in the mode-switching model indicated by dashed lines in Figure 7.Both the single-mode and the mode-switching models allow for solutions with  p > 90 • such that the dipole axis points closest toward the line-of-sight (and so the emission is expected to be brighter, Heyl & Hernquist 1998) at phase 0.13 and 0.02, respectively, landing just before the Big Peak of the light curve.For  p < 90 • , a secondary peak lies at about phase zero.However, it is obvious from Figure 7 that phase zero is the deepest minimum in the light curve.Beyond the mode switching itself, a key difference between the models is that the model without mode switching requires the angle between the magnetic axis and the spin axis () to be larger ( 43 Big Peak.In the mode switching model, this is accomplished with a smaller magnetic obliquity plus a mode switch at phases 0.1 and 0.85, so during the Big Dip the emission is dominated by a different polarization mode than the rest of the time. In the mode-switching model the emission should peak right in the middle of the Big Dip, if it is associated with the orientation of the dipole field.If one ignores the phase region of the Big Dip, the rest of the pulse profile can result from a single hot spot about 10 • in radius located about 20 • from the spin axis, and with the spin axis pointing about 100 • from the line-of-sight (as in Figure 8).These considerations, along with the unique polarization signature of the Big Dip, point toward the hypothesis that the basic emission geometry of the pulsar is straightforward with some sort of obscuration that operates during the Big Dip.1E 2259+586 in 2014, indicating the presence of the spectral feature coincident with the large dip in the pulse profile of the source.This is further highlighted in Figure 10, where the 2-8 keV, phase-resolved spectra for the XMM-Newton 2014 and 2023 observations are shown together with the IXPE one.Although the signal-to-noise ratio of the shorter 2023 observations is lower than that of the 2014 ones, the absorption feature is discernible both in the IXPE and the latest XMM-Newton observations as well.Pizzocaro et al. (2019) argue that this feature may be due to resonant cyclotron scattering of X-rays off of protons in the magnetosphere.In analogy to the results of Tiengo et al. (2013), the scattering plasma should be confined along a magnetic loop above the surface, with proton cyclotron energy ranging from about 2 to 10 keV, corresponding to a magnetic field strength along the loop from 3 × 10 14 G to 2 × 10 15 G, neglecting gravitational redshift.
To examine the implications of this scenario on the polarized flux from the surface of the neutron star, we adapt the expressions for electron resonant cyclotron scattering cross sections in the nonrelativistic limit (Nobili et al. 2008) to the case of scattering off protons, where  ( ′ ) is the photon angle with respect to the magnetic field before (after) scattering,  is the photon frequency,   is the proton cyclotron frequency and  0 is the classical proton radius (a factor of 1,836 smaller than that of the electron).On the left-hand sides, the first (second) subscript indicates the polarization mode of the incoming (scattered) photon.A key feature of this scattering process is that the ordinary mode photons are less likely to scatter (by a factor of three for an isotropic radiation field), and the photon after scattering is three times more likely to be in the extraordinary mode, regardless of the polarization of the incoming photon.Consequently, if the absorption feature evident in the XMM-Newton observations is due to resonant cyclotron scattering, the radiation that passes through the plasma without scattering will be dominated by O-mode photons.
To illustrate this, let us assume that the count rate at phase zero in the absence of scattering is 0.275 Hz, and the radiation before scattering is unpolarized.The rate of scattered photons is 0.150 Hz (from the decrease in flux during the Big Dip as observed by IXPE); so, if the rate of scattering of extraordinary and ordinary photons is 0.090 Hz and 0.060 Hz, respectively (to total 0.150 Hz), the unscattered radiation that we observe at phase zero would have  = 0.030 Hz (dominated by the ordinary mode) and a polarization degree of about 25%, as observed.A modest difference in the scattering probability of 50% in this example is sufficient to account for the observed polarization in the minimum.Two thirds of the scattered photons emerge in the extraordinary mode and one third in the ordinary mode, so over the three phase bins, where the scattering occurs, a net of 0.05 photons per second are scattered into the extraordinary mode (i.e O → X).If we assume that these are visible over the five phase bins from 0.18 to 0.46, the average rate is 0.03 Hz, as observed in these phases at the corresponding polarization angle, if the dominant mode does indeed switch from ordinary to extraordinary at about phase 0.1.The presence of scattered photons in the Big Peak can account for the observed polarization at this phase of the star's rotation.On the other hand, the Little Peak is not appreciably polarized and is somewhat smaller in amplitude; both these features are expected if the scattered photons do not contribute at this phase.In principle the plasma loop is simply hidden behind the star during the Little Peak until its footprints appear over the horizon at the beginning of the Big Dip and the loop begins to block our line-of-sight to the emission region.Clearly, this scenario requires a complicated geometry for the magnetic field near the surface that does not correspond to the large-scale dipole field of the neutron star.This should in principle rule out an explanation of the observed polarization angle in terms of a simple RVM.The fact that the RVM does indeed provide a reliable explanation to the data is further evidence that the polarization direction is determined at a distance far away from the surface (i.e. the polarization-limited radius), where the field and therefore the polarization vectors are aligned with the global dipole direction.

SPECTROPOLARIMETRIC MODELLING
We next examine the extent of polarization averaged over the rotational phase as a function of energy.In order to get a better sense of the polarization properties at the emission, the polarization angles of the photons at each phase are measured relative to the best-fitting RVM (with mode switching, the dashed curve in the third panel of Figure 7).The results are reported in Figure 11.In particular the values of / are consistent with zero which reflects the fact that the polarization states of the photons are conserved as the photons travel out to the polarization-limiting radius (Heyl & Shaviv 2000).
If there were a substantial source of polarized emission outside the polarization-limiting radius, the value of / would not necessarily vanish.The component of polarization along and across the dipole axis (/) ranges from near 20% at 2 keV and drops to be essentially consistent with zero above 4 keV.This dovetails with the hypothe- The polarization angles are generally consistent with zero within the uncertainties, while the polarization degree is similar to that shown in Figure 7, except for phase 0.1, during the rapid swing in the RVM, where the polarization is somewhat washed out in Figure 7.However, even when correcting for the rotation, the polarization degree at this phase lies below MDP 99 .
sis that the observed polarization is generated by resonant cyclotron scattering off of protons, as the proton cyclotron line lies at lower energies in the middle of the Big Dip where the polarization is strongest.
We further examine the energy dependence of the polarization as a function of phase.We focus on just the two regions with large polarized fractions, the Big Dip and the Big Peak and on a relatively narrow range in energy, because the number of photons turns out to be insufficient to reliably determine the polarization above 2.6 keV in these phase intervals.The results are shown in Figure 12.In both the regions considered, we find that the values of / are consistent with zero, as expected.The polarization degree in the Big Dip (upper panel) is essentially constant across the band, indicating that an equal fraction of photons are scattered as a function of energy.Since the spectral feature is broad (see the bottom panel of Figure 7 and Section 5.1), this is not surprising.On the other hand, the lower panel of Figure 12 shows that the polarization degree decreases with  energy and is consistent with zero above 2.3 keV.This feature can also be explained in the scattering picture after we better understand the spectral components as a function of phase (see Section 3.2).
To gain further understanding of the spectral behaviour of the source, we performed the same spectral fitting within each of the phase bins defined in Figure 3 and Table 3.We also fit the 2023 XMM-Newton EPIC-pn spectra within these bins with and without an absorption feature (the results are reported in Table 5).The fits with and without the feature yield acceptable  2 values for the Big Dip and Little Dip, with fits with the spectral feature being preferred.However, for the Big Peak and Little Peak, the fits without the spectral feature are unacceptable.Interestingly, during the Rise, where, according to Figure 7, the line is quickly changing in energy, neither fit is acceptable.In all of the fits, the energy of the spectral line is about 1 keV, with a width about 0.2 keV and a depth ≈ 0.15 keV.The emission during the Big Peak is typically harder than in the Big Dip.We found that the fraction of scattered photons during the Big Dip was approximately constant with energy (upper panel of Figure 12).If these photons are preferentially scattered into the extraordinary mode, the contribution of the scattered photons at another phase will be largest where the spectrum of incoming photons is largest.Because the Big Dip is relatively softer than the Big Peak, the relative contribution of the scattered photons to the observed flux will be larger at lower energies, resulting in a decrease of the polarization degree with energy (lower panel of Figure 12).

The spectral feature
We define the region in energy and phase containing the spectral feature as centred on where and {} denotes the fractional part and  is the rotational phase in radians.We take the width of the feature to be 2 keV.The particular numerical values in equations ( 3) and ( 4) were determined by finding the region of width 2 keV with the smallest mean value of the phase-resolved spectrum normalized by the phase-averaged energy spectrum and then by the energy-integrated pulse profile (the upper middle panel of Figure 13).By design, the mean at each phase is unity.For the XMM-Newton observation, the standard deviation of the mean over a region with the area delineated in the upper panel is 0.01, whereas the value obtained in the XMM-Newton energyphase region is 0.86 (fourteen standard deviations below the expected value).For IXPE in the similarly sized region the standard deviation is 0.005 and the value obtained is 0.88 (twenty-four standard deviations below the expected value), indicating with high confidence that the spectral feature is also present in the recent IXPE observations.The upper most panel of Figure 13 depicts the response of a filter with the shape of the feature against the three datasets.In all three the feature is detected with the matched filter.We can exploit the time and energy-resolution of IXPE to calculate the polarization degree for only the events within the phase and energy domain delineated by equation (3) within the Big Dip phase interval, (0.86, 0.07).We obtain / = 0.05 ± 0.06 and / = 0.28 ± 0.06 (measured relative to North, with MDP 99 = 0.18).If we examine the same phase range and exclude the energy range of the feature, we obtain / = −0.03± 0.04 and / = 0.20 ± 0.04 (with MDP 99 = 0.11), demonstrating that the spectral region defined by the feature in the XMM-Newton observation may be significantly more polarized than the emission at other energies during this phase of the star's rotation.

CONCLUSIONS
Observations of 1E 2259+5586 with IXPE and XMM-Newton support a consistent scenario for the emission from this source.According to our interpretation, the emission is initially only weakly polarized as expected for a condensed surface layer as in 4U 0142+61 (Taverna et al. 2022) and1RXS J170849.0-400910 (Zane et al. 2023), but during particular phases of the star's rotation the radiation that reaches us passes through a loop of plasma, and protons scatter the radiation in the cyclotron resonance (as in SGR 0418+5729, Tiengo et al. 2013).As the scattering is preferentially from the ordinary mode into the extraordinary mode, this results in a net polarization in the ordinary mode during the phase where the radiation passes through the loop towards us (the Big Dip), and in the extraordinary mode during the phases where the loop is visible but does not intersect the line-of-sight to the emission region (the Big Peak and Little Dip), leaving us to observe scattered photons.During the Little Peak, the polarization is weak, presumably because an emission region with its inherently weak polarization is visible, but the plasma loop is hidden from view, so scattered photons do not contribute during this phase.Further observations of 1E 2259+586 with IXPE could confirm the hint that the polarized flux correlates both in energy and phase with the spectral absorption feature found by Pizzocaro et al. (2019).

Figure 1 .
Figure 1.1E 2259+586 phase evolution (in radian units) as a function of time fitted with a linear plus a quadratic component for the whole data sample used in the analysis, that includes IXPE (black filled circles), NICER (red filled squares) and XMM-Newton (blue filled triangles).The dash-dotted line marks the best fit.

Figure 3 .
Figure 3. IXPE Counts as a function of rotational phase and elapsed time for the inferred values of spin frequency and frequency derivative (see text for details).The zero phase and reference time are at MJD 60097.966874079.The spectral analysis phase regions are delineated and named in the upper panel (the Big Dip region spans across phase zero).
1 confidence level. Derived by adopting a 3.2 kpc distance(Kothes & Foster 2012;Pizzocaro et al. 2019). For fitting IXPE data the spectral decomposition was convolved with a constant factor to take into account the different calibration of the 3 DUs (the relative calibration factors obtained from the fit are compatible with those found in previous magnetar analyses, seeTaverna et al. 2022;Zane et al. 2023;Turolla et al. 2023). Frozen to the value obtained from the PN fit.

Figure 4 .
Figure 4. Left: spectral fit of the EPIC-pn XMM-Newton data with the phabs×(bbodyrad+powerlaw)×gabs model in the 0.5-8 keV range.Right: same for the IXPE DU1 (black), DU2 (red) and DU3 (green) data in the 2-8 keV range.The single spectral components are marked by dotted lines (see also Table2) and the background counts by crosses with error bars.

Figure 5 .
Figure 5. Normalized Stokes parameters / and / (filled circles with 1 error bars) averaged over the rotational phase and integrated in the 2-8 keV energy band for the single IXPE DUs (cyan, orange and green, respectively) for the sum (black).The gray-dotted circles represents the levels of constant PD, while the gray-solid lines those of constant PA (measured East of North).The red shaded region represents the MDP at 99% confidence level for the combined measurement.

Figure 6 .
Figure 6.The Stokes Parameters  and  as a function of phase in units of the mean number of counts per second observed with IXPE in the 2-8 keV range.The uncertainties correspond to Δ log  = 1/2 contours of the likelihood.Each cross is labeled by the central value of the corresponding phase bin and coloured according to phase starting with red at phase zero and proceeding along the colour wheel through green then blue.

Figure 7 .
Figure 7. IXPE and XMM-Newton energy-integrated (2-8 keV range) counts (upper), PD (upper-center) and PA (lower and lower-center, blue points with error bars).The uncertainties in the lower panels correspond to Δ log  = 1/2 contours of the likelihood.The orange curve in the upper panel is the XMM-Newton (0.3-12 keV) pulse profile from Pizzocaro et al. (2019) scaled to the mean count rate of IXPE.The red curve in the second panel depicts the value of MDP 99 for each bin.The orange curves in the third panel show the best-fitting RVM model without (solid) or with (dashed) mode switching The lowermost panel shows the phase-resolved spectrum observed in 2014 with XMM-Newton EPIC-pn (as in Fig. 2 of Pizzocaro et al. 2019), normalized to the phase-averaged energy spectrum and energy-integrated pulse profile.The vertical green lines mark the five phase intervals used in our analysis.

Figure 8 .
Figure 8. Posteriors of the RVM for the IXPE observations of 1E 2259+5586The two-dimensional contours correspond to 68%, 95% and 99% confidence levels.The histograms show the normalized one-dimensional distributions for a given parameter derived from the posterior samples.The upper grid depicts constraints on the RVM with a mode switch between phases  1 and  2 .The lower grid depicts the constraints on the RVM without a mode switch.

Figure 9 .
Figure 9. Same as in the upper, upper-center and lower-center panels of Figure7, but with Stokes parameters referred to the frame of the best-fitting RVM without mode switching.The polarization angles are generally consistent with zero within the uncertainties, while the polarization degree is similar to that shown in Figure7, except for phase 0.1, during the rapid swing in the RVM, where the polarization is somewhat washed out in Figure7.However, even when correcting for the rotation, the polarization degree at this phase lies below MDP 99 .

Figure 10 .Figure 11 .
Figure 10.Phase-resolved spectra of 1E 2259+586 in the 2-8 keV range as observed by XMM-Newton EPIC-pn in 2014 (upper panel, as in Fig. 3 of Pizzocaro et al. 2019), by XMM-Newton EPIC-pn and MOS in 2023 (middle panel), and by IXPE in 2023 (lower panel).These are normalized to the phase-averaged energy spectrum but not normalized to the energy-integrated pulse profile.

Figure 12 .
Figure 12.Same as in Figure 11 but for the Big Dip (upper panel) and the Big Peak (lower panel) phase intervals, restricted to the 2-2.6 keV range.Both  and  are measured relative to the best-fitting RVM.

Figure 13 .
Figure 13.Dynamic Phase-and Energy-Normalized Spectra.Upper: Response using a matched filter defined by the line detected in the XMM-Newton 2014 EPIC-pn data when applied to the XMM-Newton EPIC-pn data from 2014, EPIC-pn and MOS data from 2023, and the IXPE data in 2023.The response finds the line feature in all three datasets centered at about the phase 0.97.Upper-Middle: the dynamic spectrum from XMM-Newton EPiC-pn in 2014 with the defined feature region from Eq. 3 superimposed.Lower-Middle: the dynamic spectrum from XMM-Newton EPIC-pn and MOS in 2023.Lower: the dynamic spectrum from IXPE in 2023.Both are normalized by the phase-averaged energy spectrum and then the energy-integrated pulse profile (as Fig. 2 of Pizzocaro et al. 2019, and Fig. 7 above).

Table 4 .
Best-fitting RVM parameters.When the data are consistent with the model, the log-likelihood (log ) is normally distributed.The fit quality in the last column is given in terms of the best-fitting log-likelihood compared with the expected value, where positive values indicate better-than-expected values.