Soft excess in the quiescent Be/X-ray pulsar RX J0812.4-3114

We report a 72 ks XMM-Newton observation of the Be/X-ray pulsar (BeXRP) RX J0812.4-3114 in quiescence ($L_X \approx 1.6 \times 10^{33}~\mathrm{erg~s^{-1}}$). Intriguingly, we find a two component spectrum, with a hard power-law ($\Gamma \approx 1.5$) and a soft blackbody-like excess below $\approx 1~\mathrm{keV}$. The blackbody component is consistent in $kT$ with a prior quiescent Chandra observation reported by Tsygankov et al. and has an inferred blackbody radius of $\approx 10~\mathrm{km}$, consistent with emission from the entire neutron star (NS) surface. There is also mild evidence for an absorption line at $\approx 1~\mathrm{keV}$ and/or $\approx 1.4~\mathrm{keV}$. The hard component shows pulsations at $P \approx 31.908~\mathrm{s}$ (pulsed fraction $0.84 \pm 0.10$), agreeing with the pulse period seen previously in outbursts, but no pulsations were found in the soft excess (pulsed fraction $\lesssim 31\%$). We conclude that the pulsed hard component suggests low-level accretion onto the neutron star poles, while the soft excess seems to originate from the entire NS surface. We speculate that, in quiescence, the source switches between a soft thermal-dominated state (when the propeller effect is at work) and a relatively hard state with low-level accretion, and use the propeller cutoff to estimate the magnetic field of the system to be $\lesssim 8.4 \times 10^{11}~\mathrm{G}$. We compare the quiescent thermal $L_X$ predicted by the standard deep crustal heating model to our observations and find that RX J0812.4-3114 has a high thermal $L_X$, at or above the prediction for minimum cooling mechanisms. This suggests that RX J0812.4-3114 either contains a relatively low-mass NS with minimum cooling, or that the system may be young enough that the NS has not fully cooled from the supernova explosion.

is forced to move along magnetic field lines when the magnetic pressure equals the ram pressure of the infalling matter, defining the magnetospheric radius within which a hot disc will be disrupted. High magnetic fields in systems with short spin periods may halt or largely suppress accretion by forming centrifugal barriers when the magnetospheric radius is larger than the corotation radius (where the orbital period in the disc matches the NS spin period). In this situation, material threading on to magnetic field lines must be accelerated to higher velocities, which moves it outwards through the disc, inhibiting accretion; this is known as the 'propeller regime' (Illarionov & Sunyaev 1975). As the mass accretion rate falls during an outburst decline, allowing the magnetosphere to expand, an abrupt drop in X-ray luminosity has often been observed, suggesting the system has entered the propeller regime (Stella, White & Rosner 1986;Campana et al. 2001;Tsygankov et al. 2016). However, signs of continued accretion are still observed in some low-luminosity systems that have seemingly transitioned to the propeller regime. For example, pulsations from 3A 0535+26 were detected when the source was at L X ≈ 2-4 × 10 33 erg s −1 (Negueruela et al. 2000;Mukherjee & Paul 2005). Multiple scenarios have been proposed to explain the mechanisms of matter penetrating the centrifugal barrier (e.g. Romanova et al. 2004;Doroshenko, Santangelo & Suleimanov 2011). Tsygankov et al. (2017b) recently proposed that in systems with sufficiently long spin period, below a certain accretion rate, the disc temperature may fall below the hydrogen ionization temperature (∼6500 K) rendering a recombined neutral 'cold disc', which can penetrate through the centrifugal barrier of the magnetosphere. Several sources have been observed to maintain quasi-stable accretion at an intermediate luminosity, higher than the limiting luminosity for the propeller regime (e.g. Rouco Escorial, van den Eijnden & Wijnands 2018).
RX J0812.4-3114 (hereafter J0812) was first identified by the ROSAT Galactic Plane Survey (Motch et al. 1991) as an X-ray source that positionally coincides with a Be star (LS 992, B0.2IVe; see Motch et al. 1997;Reig et al. 2001). Corbet (1999) reported that the source entered an active state in early 1998, with a series of prominent outbursts. Reig & Roche (1999) reported two observations in February 1998 with the Proportional Counter Array (PCA) on the Rossi X-ray Timing Explorer (RXTE), with which they detected strong X-ray pulsations at a period of 31.885 s. The X-ray pulsar (XRP) nature of J0812 was hence corroborated. Using data from the All Sky Camera (ASM) on board RXTE, the orbital period was then found to be ≈81 d by Corbet & Peele (2000), who used the RXTE/PCA observation on 1999 March 25 to again confirm the strong pulsations at a period of 31.88 s. Fig. 1 shows the full ASM light curve, from which we see that the source stayed in an active state until early July of 2000, and returned to a relatively low-count state ever since.
The Chandra X-ray Observatory observed J0812 in 2013 July, as part of a campaign to systematically study quiescent BeXRPs (PI: Wijnands, ObsIDs: 14635-14650), during which it had a relatively low L X of ∼ 2 × 10 33 erg s −1 and a very soft spectrum. A fit with an absorbed power-law gave a photon index of ≈5.6, which suggests a blackbody-like fit would be more appropriate. Intriguingly, the blackbody fit gave an unusually low temperature of kT ≈ 0.1 keV, suggesting thermal emission from a large but poorly constrained inferred emission radius (up to ∼10 km; see Tsygankov et al. 2017b for more details).
In this work, we report results from our recent XMM-Newton observation of this BeXRP. The paper is organized as follows: in Section 2, we show observational information, the methodologies of our data reduction, and the results of spectral and temporal analyses. In Section 3, we present our discussions on the possible nature of the source and some relevant calculations, and in Section 4, we draw conclusions.

O B S E RVAT I O N A N D A NA LY S I S
We use data from our 72 ks XMM-Newton observation on 2018 October 9 using the European Photon Imaging Camera (EPIC; ObsID: 0822050101). Both PN and MOS detectors used full frame mode, and the source was covered by both detectors (see Fig. 2). To avoid optical contamination, medium filters were applied to both PN and MOS cameras. We also made use of the 4.6 ks Chandra ACIS-S observation (ObsID: 14637; PI: Wijnands) taken on 2013 July 13.
For the XMM-Newton observation, we used the event files from the Processing Pipeline System products, as derived from the Observation Data Files for further reduction and analyses. We cleaned the flaring particle background by first generating high energy (10-12 keV), single event (PATTERN = 0) light curves for both PN and MOS cameras, using the evselect task from the latest XMM SCIENCE ANALYSIS SOFTWARE (SAS; version 17.0). 1 Based on the light curves, we identified a period of low and steady background, with count rates ≤ 0.9 counts s −1 for PN, ≤ 0.2 counts s −1 for MOS1, and ≤ 0.4 counts s −1 for MOS2. These thresholds were then applied to the tabgtigen task to find good time intervals (GTIs), rendering effective exposures of ≈50 ks for PN, and of ≈63 ks for MOS1 and MOS2. The GTI files are then used to create filtered event files that are used to further generate spectra and time series. The Chandra data set was reprocessed with the chandra repro task from CIAO 4.11 (CALDB 4.8.2) 2 to align it with the up-to-date calibration. The reprocessed level-2 event file was used for further analyses.
To study the source's long-term accretion history, we also obtained the RXTE/ASM light curve from the ASM Data Product page. 3 The light-curve spans from 1996 January 5 to 2011 December 29.

Spectral analysis
Spectra extracted from the XMM-Newton and Chandra data sets were analysed with the HEASOFT/XSPEC software. The PN and MOS spectra were rebinned to at least 20 counts per bin for simultaneous fits using χ 2 -statistics. For each fit, we report the reduced χ 2 (χ 2 ν hereafter) together with the corresponding degrees of freedom (d.o.f.) as χ 2 ν (d.o.f.). Uncertainties and upper (or lower) limits of parameters are reported at the 90 per cent confidence level. The distance used in all analyses is the distance to the Be star (LS 992; Motch et al. 1997)  where distances are estimated from Gaia parallaxes by Bayesian analysis with a weak distance prior (Bailer-Jones et al. 2018;Gaia Collaboration 2018). Note that the distance estimate might be different using different types of measurement. For example, d = 8.6 ± 1.8 kpc according to Coleiro & Chaty (2013). For all analyses of XMM-Newton spectra, we use the energy channels between 0.2 and 10 keV while channels between 0.5 and 10 keV are noticed for Chandra fits. We accounted for interstellar absorption by convolving our models with the Tuebingen-Boulder ISM absorption model (tbabs in XSPEC) using wilm abundances (Wilms, Allen & McCray 2000). Since we have a moderately large number of spectral counts in the low energies, we tried fits with a free n H (hydrogen column density), and fits with n H fixed to the expected (based on H I) Galactic value (≈ 0.48 × 10 22 cm −2 ; see Kalberla et al. 2005). We first tried a simple absorbed power-law (powerlaw in XSPEC, pow hereafter) fit to the XMM-Newton spectra, finding a photon index ( ) of 1.25 ± 0.10. However, the fit exhibits strong residuals below ∼1 keV that lead to a χ 2 ν of 1.68 (95), suggestive of an additional soft component. The n H from this fit is a factor of ∼5 below the Galactic value. However, if we force n H to the Galactic value, the model gives a significantly worse fit with χ 2 ν = 2.66 (96). We proceeded by adding a blackbody component (pow+bbodyrad) to account for the soft excess. The fit was significantly improved to a χ 2 ν of 1.003 (93), with a kT = 0.06 ± 0.01 keV, a softer power-law index = 1.74 +0.18 −0.17 , and an enhanced absorption column density (n H = 1.25 +0.24 −0.23 × 10 22 cm −2 ). The inferred radius of the blackbody emission region (R bb ) reached 440.2 +1223.6 −317.3 kmtoo large to be consistent with the scale of an NS. The fit is significantly worse (χ 2 ν = 1.31 (94)) when n H is fixed to the Galactic value; however, this gives a much smaller emission region (R bb = 10.3 +5.6 −3.6 km) and a slightly higher blackbody temperature (kT bb = 0.10 ± 0.01 keV), which is consistent with the results from the 2013 Chandra observation. Similarly, a substitute for the soft component with a thermal bremsstrahlung model (pow+bremss) also suggests an enhanced n H = 1.36 +0.24 −0.23 × 10 22 cm −2 , giving a good fit (χ 2 ν = 1.02 (93)), while fixing n H to the Galactic value rendered a worse fit (χ 2 ν = 1.42 (94)). The soft component can also be modelled by a Gaussian emission line (pow+gauss), with the line energy located at 0.63 +0.06 −0.08 keV and a broad line width of 0.13 +0.04 −0.03 keV when n H is fixed. This model is acceptable either when n H is free (χ 2 ν = 0.99 (92)) or fixed (χ 2 ν = 1.02 (93)); however, no strong and broad emission features are expected in XRPs around this energy.
We also tried fits with a magnetic NS atmosphere model (pow+nsmaxg; see Ho, Potekhin & Chabrier 2008;Potekhin, Chabrier & Ho 2014), assuming a hydrogen atmosphere on a 1.4 M , 12 km NS, for different choices of magnetic field (10 10-13 G). Initially, we assumed emission from the entire NS (normalisation=(R em /R NS ) 2 = 1.0, where R em is the radius of the emission region). These fits were not superior to those using a blackbody. For example, when B = 10 12 G, the fit yielded χ 2 ν = 1.61 (94) if n H is freed (n H = 0.26 +0.05 −0.06 × 10 22 cm −2 ), and a worse χ 2 ν = 2.00 (95) if n H is fixed. The former resulted in a surface temperature (kT s ) of 0.07 ± 0.01 keV, while the latter gives a slightly higher surface temperature at 0.082 ± 0.003 keV. We then tried fits with a free normalisation and found that they improved (χ ν = 1.13 (93) for n H -free versus χ ν = 1.45 (94) for n H -fixed); however, the inferred size of the emission region is too large to be plausible (R em /R NS = 812 +120 −401 for n H -free fit versus R em /R NS = 6.2 +2.9 −2.4 for n H -fixed fit). We obtained better fits with B = 10 11 G. When the normalisation is fixed to unity, χ 2 ν = 1.46 (94) for free n H (= 0.29 +0.05 −0.05 × 10 22 cm −2 ) while χ 2 ν = 1.86 (95) for fixed n H . The n H -free fit gives a kT = 0.072 +0.005 −0.007 keV while the n H -fixed fit results in a similar kT = 0.082 +0.003 −0.003 keV. With free normalisations, the fits are somewhat improved to χ 2 ν = 1.13 (93) when n H is free and to χ 2 ν = 1.22 (94) when n H is fixed. The problem with these fits, however, is still that the normalisations infer larger emission regions than the NS surface. The n H -free fit (n H = 0.75 +0.17 −0.15 × 10 22 cm −2 ) has quite a low kT ≤ 0.04 keV that exceeds the allowed lower limit, so the corresponding inferred emission radius is much greater than the NS radius (R em /R NS = 36 +40 −23 ); the n H -fixed fit gives a slightly higher kT = 0.04 +0.01 −0.01 keV yet still yields an emission radius greater than the NS radius (R em /R NS = 8 +7 −3 ). We also tried to fix the normalisation to values smaller than unity, modelling surface hotspots, but did not find any significant improvement in the fits.
We also substituted other physically motivated models for the power-law component. We first tried a Comptonization model (comptt+bbodyrad), assuming that soft photons from the blackbody component are up-scattered to form the hard component. We obtained a fair fit (χ 2 ν = 1.02 (91)), with an optical depth (τ ) of 0.10 +2.94 −0.08 , but an unconstrained plasma temperature kT e ≥ 36.80 keV. The fit also suggests a higher n H = 1.22 +0.28 −0.22 × 10 22 cm −2 , and gave a poorer fit (χ 2 ν = 1.19 (92)) when n H was frozen at the Galactic value.
We found an equivalently good fit (χ 2 ν = 1.01 (92) when n H is free) using a power-law with an exponential highenergy cut-off (cutoffpl+bbodyrad). When n H was free, we found = 1.73 +0.18 −0.74 but the cut-off energy was unconstrained (E cut ≥ 5.33 keV). The cut-off energy was better constrained to 2.97 +2.93 −1.06 keV with n H fixed to the Galactic value, but the fit itself became worse (χ 2 ν = 1.19 (93)). The pow+bbodyrad, comptt+bbodyrad, cutoffpl+bbodyrad, and some of the pow+nsmaxg fits all statistically suggest an n H above the Galactic value, but, with the high n H , infer a very high intrinsic luminosity of ∼ 10 35 erg s −1 , which seems unlikely for a quiescent system. For example, the pow+bbodyrad fit would give a blackbody component with L X (0.4-1 keV) ∼ 10 35 erg s −1 , while the power-law component has a much smaller contribution (L X (1-10 keV) ∼ 10 33 erg s −1 ). As the soft component must come from either reprocessed accretion energy (through either the NS surface, or an accretion disc), or other stored heat in the NS, it seems quite unlikely that the soft component could reach ∼10 35 erg s −1 , while the hard component remains at ∼10 33 erg s −1 . For this reason, we prefer the fits with fixed n H on physical grounds.
Further investigation of the fixed-n H fits indicates that the main reason for the poor fits is that the models are above the data at around 1 keV, which leaves apparent residuals that resemble an absorption feature, in both the PN and MOS spectra. We thus tried incorporating a Gaussian absorption line component (gabs) into these models (pow+bbodyrad, comptt+bbodyrad, and cutoffpow+bbodyrad) to compensate for the residuals. We found that the fits were significantly improved (e.g. the pow+bbodyrad fit was improved from χ 2 ν = 1.31 (94) to χ 2 ν = 1.09 (91); χ 2 = 23.95), so the absorption feature might be genuine. The resulting kT bb from each model is slightly higher than in the original model without the gabs component (0.12 ± 0.01 keV versus 0.11 ± 0.01 keV for pow+bbodyrad fit). The absorption line is regularly found between 0.99 and 1.02 keV. As a point of comparison, we also added a gabs component to the corresponding fits in which n H was free, but found no clear improvement (χ 2 ν = 1.00 (93) to χ 2 ν = 0.91 (90); χ 2 = 11.01).
The nsmaxg models intrinsically include redshifted cyclotron lines. For B = 10 11 G, the line is approximately at 0.94 keV (assuming a 1.4 M , 12 km NS; this nsmaxg model is henceforth referred to as #1), which can partially compensate for the residual at 1 keV, so the pow+nsmaxg fit is slightly better than the pow+bbodyrad fit (χ 2 ν = 1.31 (94) versus χ ν = 1.22 (94)). The fit can be further improved by adjusting the line location; for example, one can use a larger radius to reduce the redshift and therefore elevate the line energy. We do find a better fit (χ 2 ν = 1.08 (94)) when we increase R NS from 12 to 14 km. A more physical approach might be using a model with the same mass and radius but a magnetic field slightly above 10 11 G such that the intrinsic line locates exactly at kT = 1.01 keV (as suggested by our gabs * (pow+bbodyrad) fit, corresponding to a magnetic field of 1.07 × 10 11 G; this nsmaxg model is henceforth referred to as #2). This also significantly improved the fit to χ 2 ν = 1.04 (94). However, either approach results in a normalisation greater than unity ( R em /R NS = 6 +6 −3 versus 3.8 +0.2 −0.3 for the former and latter cases, respectively). This model shows residuals at ≈1.36 keV, resembling a second absorption feature. We thus also tried convolving a gabs component to the pow+nsmaxg models. With the introduction of three more free parameters, the fits did not become any better but did yield more reasonable normalisations. For example, χ 2 = 0.95 (91) for model #2, with R em /R NS = 1.3 +2.9 −0.7 . Introducing a second absorption line at a higher energy might be a result of variation of magnetic field due to complex distribution of the magnetic across the NS surface.
We also fit the Chandra spectrum of RX J0812.4-3114 presented by Tsygankov et al. (2017a) to a bbodyrad model. For this purpose, this low-count spectrum was binned to at least 1 count per bin, and we used C-statistics (Cash 1979). Because of the poorer calibration at low energies, all channels below 0.5 keV were ignored during the fit. We fixed n H to the Galactic value, considering the discussion above, and the low-count statistics in this spectrum. The best-fitting model is a blackbody with a low temperature (kT bb = 0.13 +0.03 −0.02 keV) and an unconstrained blackbody radius (R bb ≤ 11.6 km). Adopting the Gaia-estimated distance of 6.76 kpc (Bailer-Jones et al. 2018; Gaia Collaboration 2018), we found an unabsorbed 0.5-10 keV luminosity of 5.5 +5.2 −2.7 × 10 32 (d/6.76 kpc) 2 erg s −1 . We also tried the nsmaxg fits to the Chandra data, using both models #1 and #2, with free normalisations. We noted that there is no sign of absorption feature as in the XMM-Newton spectra. As a result, either model gives an equally fair fit (goodness= 47.4 per cent versus goodness= 40.6 per cent for the #1 and #2, respectively). However, due to low counting statistics, we cannot get proper constraints on the normalisations. To make sure that the absence of a hard component in the Chandra spectrum is not due to low counting statistics, we include a hard power-law component, with the powerlaw index ( = 1.32) from the gabs * (pow+bbodyrad) fit, and fit the normalisation parameter (with kT bb and R bb free) of the power-law. We found an upper limit for the power-law flux to be ≈2 orders of magnitude lower than the flux of the blackbody component (F 2-10 7.94 × 10 −15 erg s −1 cm −2 , or L 2−10 4.34 × 10 31 × (d/6.76 kpc) 2 erg s −1 ), suggesting that the Chandra spectrum is genuinely soft.
In summary, solely based on the fit quality, it seems that the XMM-Newton spectrum can be well fitted by either an absorbed soft model (which could be bbodyrad or nsmaxg) plus a hard component (which could be pow, cutoffpl, or comptt) with increased absorption, or by the same model but with the n H fixed to the Galactic value and modified by an absorption line (at ≈1.0 keV for bbodyrad and ≈1.3 keV for nsmaxg). However, considering the physical implications, the latter model is more favourable (see also Section 3). The inferred kT bb , R bb (or kT and R em inferred from the nsmaxg model), and thermal L X from the XMM-Newton data are consistent with the results from the Chandra data. We summarise all relevant XMM-Newton spectral fitting parameters in Table 1, and the Chandra parameters in Table 2. In Fig. 3, we show the XMM-Newton spectra with the tbabs * gabs * (bbodyrad+pow) and tbabs * gabs * (nsmaxg+pow) models, with fixed n H , overplotted. For comparison, we also show the Chandra data and their best-fitting model (tbabs * bbodyrad) and plot the upper limit of a possible hard component for the Chandra spectrum.

Temporal analyses
The PN camera has sufficiently high timing resolution (73.4 ms in full-frame mode) to search for pulsations in this system. For that purpose, we first applied barycentric correction for the arrival times of photons in the flare-cleaned event lists, and then extracted PN light curves with SAS evtselect task over a soft band that primarily covers the soft excess (0.4-1.0 keV), a band covering the hard spectral component (1-10 keV), and a broad band (0.4-10 keV). These time series were then rebinned to 1 s time bins, which are short enough to search for the expected pulse period (≈31.88 s; see Corbet & Peele 2000).
Timing analyses were performed with tasks from the HEASARC/XRONOS software package. 4 We searched for pulsations by running the powspec task on the rebinned PN light curve in each band using a total of 32 768 frequency bins between 1.52 × 10 −5 and 0.5 Hz. A clear periodicity was revealed at ≈ 0.031 Hz in the hard band power spectrum, represented by a prominent peak (power ≈72.53; see left-hand panels of Fig. 4). The noise in a power spectrum is expected to follow a χ 2 distribution with 2 d.o.f. (Leahy et al. 1983), from which we derived a 5σ significance level given the number of frequency bins in our analysis. To search for the exact period, we used the efsearch tool, which rebins and folds the light curve over a range of period and searches for best period by finding the maximum χ 2 from fitting a constant to the folded light curves (see Fig. 5). We then folded the light curves at the best period with the efold task to create pulse profiles (see right-hand panels of Fig. 4).
We found a best pulse period in the hard-band light curve at 31.908 ± 0.009 s (upper and lower bounds correspond to periods with χ 2 values at half of the maximum). Compared to the pulse period (P 1999 ≈ 31.8856 ± 0.0001 s) found by Corbet & Peele (2000), this indicates a spin-down (Ṗ ) of 3.63 × 10 −11 s s −1 . This spin-down is rather slow, but not particularly unusual in BeXRPs. For example, SAX J0635+0533 was observed to have a long-term spin-down of >3.8 × 10 −13 since its discovery (La Palombara & Mereghetti 2017). However, the difference in spin period could have been affected by the orbital Doppler effect. To estimate the maximal magnitude of the orbital Doppler effect, we assume the system is edge-on, and use the primary mass of 17 M from Reig et al. (2001). The resulting Lorentz factor due to orbital motion (≈4 × 10 −4 ) introduces an uncertainty of 0.013 s in the spin period. This is comparable to the spin difference we have measured, so, without further knowledge of the orbital ephemeris and inclination, the spin-down reported here should be interpreted with caution. Table 1. Best-fitting parameters to the XMM-Newton data. F X is unabsorbed flux over 0.2-10 keV, and L X is the corresponding luminosity. A † indicates parameters that were fixed during the fits, and a * marks parameters for which no valid constraints were found, thus should be taken with care. The strength in the gabs model is related to the σ and the optical depth (τ gabs ) of the line by Strength = √ 2πσ τ gabs . All nsmaxg models assume a 1.4 M , 12 km NS. All kT bb s, R bb s (kTs and R em s for nsmaxg models), and E gabs s are redshifted quantities as observed by distant observers.  We found no clear signs of pulsations in the folded light curve of the soft excess, to which we fitted a constant and found a χ 2 ν = 0.76. However, clear sharp dips are present in the folded hard and broad-band light curves (see Fig. 4), which were previously suggested by Galloway et al. (2001) to be partial eclipses of the emitting region by the channelled accretion column. Table 2. Best-fitting parameters to the 2013 Chandra spectrum. The notations for subscripts and superscripts are the same as in Table 1. F X is the unabsorbed flux over 0.5-10 keV, and L X is the corresponding unabsorbed luminosity. To quantify the light-curve modulation, we calculated the pulsed fraction, which is defined as where C max and C min are the maximum and the minimum count rates, respectively. We found PF = 0.84 ± 0.10 for the hard band, whereas the pulse fraction for the soft band is more uncertain but could be very low (PF = 0.37 ± 0.18). Non-detection of pulsations in the soft excess could be a result of the relatively low counting statistics (≈ 323 counts), or due to the fact that the pulsed fraction is genuinely too low. To test this, we generated a series of simulated light curves that are modulated by sinusoidal functions at the observed pulsed period (31.908 s) but with different pulse fractions (up to 0.99).
Using the observed counts in the soft excess, we generate 100 light curves for each pulsed fraction. For each pulsed fraction, pulsations are considered to be detectable if more than 90 per cent of the realizations result in powers at the expected pulsed period that are 3σ above the corresponding noise levels. With the given counts in the soft excess, we found that the pulsed fraction has to be 31 per cent for a pulsed signal to be detected. In other words, if the pulsation is not detected, the pulsed fraction is then at most 31 per cent.
To check for long-term variability, we rebinned the light curves to 500 s time bins, using the same set of soft (0.4-1 keV), hard (1-10 keV), and broad (0.4-10 keV) bands while defining a hardness ratio with Hardness = log 10 C 0.4−1 C 1−10 . ( We then fitted these rebinned light curves to constants and used the resulting reduced χ 2 s as a measure of variability. Fig. 6 shows the rebinned light curves and the corresponding time series of hardness ratios. The result implies strong variability in the hard band light curve (χ 2 ν ≈ 2.27), while the soft excess shows no sign of variability (χ 2 ν ≈ 1.09). Due to the large error bars, variability in the hardness ratio is hard to determine solely based on the χ 2 ν (≈0.82); however, the best-fitting hardness ratio (≈ −0.32 +0.04 −0.04 ) suggests that the soft excess contributes less than 50 per cent of the total observed flux (F 0.4-1 /F 0.4-10 ≈ 32.3 per cent). To further explore a possible correlation between the soft and hard counts, we calculated a correlation coefficient defined as where Cov(C soft , C hard ) is the covariance between the soft and hard count rates, while σ soft and σ hard are standard deviations in soft and hard rates, respectively. |ρ| = 1 corresponds to linear correlation, and ρ = 0 corresponds to non-correlation. We found ρ = 0.15 ± 0.14 (the uncertainty is propagated from the data), suggesting very weak or no linear correlation between the soft and hard count rates. The soft counts might therefore have a completely distinct origin from the hard counts. A plot of soft count rates against their corresponding hard count rates can be found in Fig. 7. To see if the source returned to quiescence after the active state, we also compared the powspec and efsearch results on the ASM light curve during the active state (between 1997 December 5 and 2000 July 2) with those on the light curve after the active state. We found a clear periodicity only in the former case. The best period was found to be P orb = 80.39 +3.00 −2.18 d for the active epoch with a maximum power of 73.06 (uncertainties in the period are estimated using periods with χ 2 values that are half of the maximum; see Fig. 8). This is consistent with the ∼81.3 d orbital period found by Corbet & Peele (2000). Because the ASM light curves are background-subtracted, some phases contain negative count rates. We approximated a background level by fitting the quiescent light curve to a constant. This gives a best-fitting value at −0.11 counts s −1 with χ 2 ν = 1.77, indicating variability possibly due to some minor source activity. We then applied this background level to the light curves to shift them to unsubtracted levels. The power spectra of the active and quiescent epochs are shown in the left-hand panels of Fig. 9 (the 5σ level is calculated using the same method as for the PN power spectra), while the corresponding light curves folded at the best period are shown in the right-hand panels.

Soft excess
The 2013 Chandra observation revealed a very soft (F 2-10 /F 0.5-10 4.1 × 10 −4 ) spectrum which makes RX J0812.4-3114 the coolest (kT bb ∼ 0.1 keV) among all known quiescent BeXRPs (Tsygankov et al. 2017a) that share similar  spectral shapes. Although not well constrained due to low counting statistics, the fit to the Chandra data does suggest a blackbody spectrum with a large emitting region (R bb ∼ 10 km), which has never been observed in any other quiescent BeXRP.
However, soft excesses have been observed in some luminous (L X 10 37 erg s −1 ) XRPs. Particularly in high-inclination systems, soft X-rays are thought to originate from reprocessing of hard X-rays by the optically thick inner disc region, which leads to a larger effective R bb (e.g. Endo, Nagase & Mihara 2000). Hickox, Narayan & Kallman (2004) have estimated that the corresponding blackbody temperature (T bb ) is related to L X by T bb ∝ L 11/28 X , so in faint XRPs (L X 10 36 erg s −1 ), given that L X s are low, the majority of the reprocessed hard X-rays would instead shift into the EUV regime. Reprocessing of hard X-rays is therefore not likely to be the mechanism at work for the soft excess in J0812. It might seem to be possible that the intrinsic absorption suggested by the n H - free fits might arise from an obscured inner disc region. The high unabsorbed luminosity mentioned in Section 2.1 for these fits might therefore favour this scenario. However, because the soft excess is powered by the hard X-rays, the scenario is valid only when the hard component is brighter or at least comparable to the soft excess  (e.g. Burderi et al. 2000 found a soft excess in Cen X-3 that takes ≈ 58 per cent of the total unabsorbed flux). Moreover, Endo et al. (2000) showed that the cooling time-scale of the irradiated inner disc should be only a fraction of a second, so the soft excess should also be pulsed in accordance with the hard component. Therefore, the above discussion strongly disfavours the n H -free fits with their enhanced absorption.
The companion Be star may partially contribute to the soft X-rays. According to Nazé et al. (2011), we can roughly estimate the expected L X from the companion star adopting the reported B0.2IVe spectral type (Section 1), which suggests that the companion's X-ray luminosity should lie in the range of 10 31-31.5 erg s −1 . This is ∼2 orders of magnitudes lower than the observed luminosity, so unlikely to be a significant contributor to the soft component.
Soft X-ray emission in quiescence has typically been ascribed to thermal emission from the NS, either from small regions of higher temperature -hotspots -or, although not previously detected in a quiescent BeXRP, from the entire NS surface. Hotspots can either be formed externally as channelled accretion columns heat up the polar caps (e.g. pulsed soft excess was noted in the BeXRP RX J1037.5-5647 by La Palombara et al. 2009 with R bb ∼ 128 m), or intrinsically as heat from the core is channelled towards parts of the surface by strong internal magnetic fields (Greenstein, Dolez & Vauclair 1983;Potekhin & Yakovlev 2001;Geppert, Küker & Page 2004). The large inferred R bb from our analyses therefore indicates a rather large hotspot. The former is then not likely the mechanism at work since the spot size in this source is actually predicted to be ∼ 0.1 km (following R pc = (2π R NS /(cP)) 1/2 R NS , e.g. Lyne & Graham-Smith 2006;Forestell et al. 2014). Large hotspots of radius several km have indeed been detected in some NSs (e.g. Gotthelf & Halpern 2009). However, absence of pulsations in the soft excess weakens the hotspot scenario, although it is still possible to have a relatively small pulsed fraction with a particular observer geometry. The large inferred R bb ∼ 10 km in our source suggests that the soft X-ray photons might have primarily originated from the whole NS surface. Even our fits with the smallest R bb (e.g. R bb ≈ 5.3 +2.6 −1.8 km in the bbodyrad+comptt fit) are much larger than predicted hotspot sizes, indicating that the observed spectrum might be comprised of a hotspot plus emission from the overall surface (see e.g. Elshamouty et al. 2016).
Heat thermally radiated during quiescent states is thought to be principally deposited during the previous accretion episodes, especially the bright outbursts. After outbursts, NSs cool via thermally radiating away the heat from the surface, and/or through either slow (e.g. modified Urca) or fast (e.g. direct Urca) neutrino emission processes in the core (Potekhin, Pons & Page 2015, and references therein). Heat is generated both by pycnonuclear reactions in the deep crustal regions ('deep crustal heating'; see Brown, Bildsten & Rutledge 1998), which leaks out over long (∼10 5 yr) time-scales, and by several processes in the outer crust, which leak out of the NS on shorter (months to decades) time-scales (Rutledge et al. 2002;Shternin et al. 2007;Brown & Cumming 2009;Deibel et al. 2015). It is usually assumed that the NS crust and core return back to thermal equilibrium several years after the outbursts (see the review by Wijnands, Degenaar & Page 2017). The heating rate (H) in this case is then simply related to the quiescent luminosity (L q ) and the neutrino cooling rate (L ν ) by Here, H depends on the average mass accretion rate of the system: where Q nuc is the amount of energy generated by pycnonuclear reactions per accreted nucleon (≈1-2 MeV; see Haensel & Zdunik 2008) and m u is the atomic mass unit. We first note that, given the last recorded outburst in 2000 (Fig. 1), and the dates of the Chandra (2013) and XMM-Newton (2018) observations, that shallow crustal heating is unlikely to explain the observed thermal luminosity. We then try to test the deep crustal heating scenario. Since we have some knowledge of the outburst history of our source (Fig. 1), we can make a rough calculation of the average mass accretion rate and compare the inferred L q predicted by the deep crustal heating model with our observations. We extrapolated the Chandra L X down to 0.01 keV for bolometric correction, which gave us an L q of 1.23 +1.66 −0.65 × 10 33 erg s −1 (note that uncertainties in distance are not included in this result).
We estimate the average mass accretion rate of the source based on the observed L X during the outbursts, Using the folded light curve during the active epoch, we obtained a phase-averaged ASM count rate of 0.14 ± 0.06 counts s −1 . We convert this count rate to X-ray flux using the WEBPIMMS tool 5 and account for bolometric correction by assuming a power-law model between 0.1 and 12 keV, and adopting the best-fitting powerlaw index ∼ 1 and the e-folding energy E fold ∼ 12 keV (for the upper limit of bolometric correction) from Reig & Roche (1999).
To account for uncertainty in the bolometric correction due to the unknown spectral shape, we extended the power-law to 30 keV and calculated the combined error. The resulting ASM flux is 6.46 +10.15 −2.54 × 10 −11 erg cm −2 s −1 . If we only account for errors in the flux, and assume a 12 km NS with M NS = 1.4 M , we then have Ṁ active ≈ 3.62 +5.69 −1.42 × 10 −11 M yr −1 . This is only the accretion rate over the active period. To convert it to cover the whole period of ASM observations, we calculate the weighted mean, where we have assumed that Ṁ active Ṁ quiescent . To compare with other quiescent systems, we plot observed quiescent NSs in LMXBs on a L q -Ṁ plot ( Fig. 10) with tracks for possible cooling mechanisms indicated. The theoretical tracks are calculated assuming the BSk24 equation of state (Pearson et al. 2018) with a maximum mass of 2.28 M and a mass threshold for rapid cooling of 1.595 M ; for modified-Urca processes, we included in-medium effect following Shternin, Baldo & Haensel (2018). We used the MSH and BS gap models from Ho et al. (2015) and accounted for effects of triplet superfluidity following Ding et al. (2016). J0812 likely lies above the minimum cooling curves Figure 10. Observations of the quiescent thermal luminosities of low-mass X-ray binaries (qLMXB, black) and RX J0812.4−3114 (red), compared to theoretical predictions of the thermal luminosities (redshifted as seen by a distant observer) produced by deep crustal heating for different timeaveraged accretion rates. Theoretical models for neutron stars of different masses and different heat-blanketing envelope compositions (either iron or accreted helium and carbon), as well as the qLMXB data, are taken from Potekhin et al. (2019). The relatively low-mass neutron stars in qLMXBs (upper curves) undergo the minimal cooling (e.g. Page et al. 2004), whereas the high-mass neutron stars (lower curves) undergo rapid cooling via direct Urca process. The observational accretion rates are scaled from the canonical R = 10 km in Potekhin et al. to the more probable R = 12 km in this paper. RX J0812.4−3114 lies at or above the minimal cooling tracks.
for NSs with iron heat-blanketing envelopes (for a definition of heatblanketing envelope, see Gudmundsson, Pethick & Epstein 1983), but is consistent with low-mass (1.0-1.2 M ) NSs with accreted heat-blanketing envelopes.
J0812's position above some of the minimum cooling tracks may require explanation. We speculate that this NS might be relatively young, so that the NS has not yet lost the internal heat deposited in its supernova explosion. Reference to, e.g. Page et al. (2004) shows that the thermal L X from the supernova should fall below ∼10 33 erg s −1 after 10 5 yr. J0812's companion's B0 spectral type indicates a ∼20 M mass, and thus a 2 Myr lifetime, suggesting a 5 per cent chance of observing a NS in this HMXB before it has lost its supernova heat.
We can attempt to estimate the age since the SN from its height above the Galactic Plane, and proper motion. J0812 would take 10 8 yr, at its Gaia-measured proper motion (Gaia Collaboration 2018), to reach its location ≈1. • 54 above the Galactic Plane. However, considering the much shorter lifetime of the companion star, it was almost certainly born near its current location, and this method cannot constrain this hypothesis.

Hard component: stable low-level accretion
The hard (power-law) component in quiescent accreting NSs is generally thought to originate from low-level accretion (e.g. Wijnands et al. 2015). In high-B systems, the accretion columns are channelled down to the magnetic poles, so the hard component is likely to be pulsed (e.g. in Cep X-4; see McBride et al. 2007). Theoretically, accretion on to the NS poles is only possible when the rotational velocity of the magnetic field lines is lower than the local Keplerian velocity. Otherwise, matter from the accretion disc would be spun away by the centrifugal barrier, making the system enter the propeller regime. Observationally, the onset of the propeller regime is marked by an abrupt drop in accretion luminosity (L acc ) when the luminosity decays below a limiting level (L prop ). The limiting luminosity is estimated by equating the magnetospheric radius (R m ) with the corotation radius (R c ), where the local Keplerian velocity equals the rotational velocity of the NS. One can then derive the following 6 L prop ≈ 4 × 10 37 ξ 7/2 B 2 12 P −7/3 M −2/3 1.4 R 5 6 erg s −1 , where ξ is a parameter that defines the accretion geometry (ξ = 0.5 for accretion from a disc; ξ = 1 for accretion from winds; Ghosh & Lamb 1978); B 12 is the magnetic field strength in units of 10 12 G; P is the spin period of the NS; M 1.4 is the mass of the NS in units of 1.4 M , and R 6 is the NS radius in units of 10 6 cm. We do not know L prop for J0812, its magnetic field is unknown. However, it is plausible to assume that the source was in the propeller regime (L acc L prop ) during the 2013 Chandra observation, based on the fact that the source possessed a very soft spectrum (consistent with emission of stored heat from the NS with no active accretion), and assume that the source in 2018 had left the propeller regime and was accreting (L acc L prop ), based on the fact that a pulsed hard power-law component is present. We estimate the 2018 L acc by using the power-law component ( = 1.49) from the gabs * (pow+nsmaxg) #2 fit and extrapolate it to 0.1-30 keV for bolometric correction. This yields an L 0.1−30 = 1.80 × 10 33 erg s −1 , which is then adopted as an upper limit on the propeller luminosity 6 Notice that the magnetic field here is the dipolar strength at the magnetic poles, which is a factor of 2 higher than that at the equator.
(L prop ). With this constraint on L prop , we can then place an upper limit on the magnetic field strength of the NS using equation (8). We assumed an NS of 1.4 M and R NS = 12 km. J0812 has been quiescent for many years, so it might seem plausible to assume accretion from the winds of the Be star (i.e. ξ = 1). However, a disc may still form in even wind-fed systems (Karino, Nakamura & Taani 2019) or systems fed by a companion star's circumstellar disc (Klus et al. 2014). We therefore calculated for cases of ξ = 0.5 and ξ = 1 and found B 12 (ξ = 1) 0.24 while B 12 (ξ = 0.5) 0.81, corresponding to cyclotron energies 2.81 and 9.44 keV, respectively. A recent study by Campana et al. (2018) on different classes of objects subject to the propeller effect found a best fit with ξ = 0.49 ± 0.05, which results in a different scaling factor of ≈ 3.3 × 10 36 erg s −1 for equation (8). 7 Using this factor, we obtained B 12 0.84 and cyclotron energies 9.78 keV. It should be noted that this is a rough estimate on the magnetic field, subject to errors in distance and fluxes, but the possible absorption lines, either the one at 1 keV indicated from the gabs * (pow+bbodyrad) fit or the one at 1.4 keV suggested by the gabs * (pow+nsmaxg) fit, might then be real, as the inferred upper limit from the propeller argument (in the case of ξ = 1) is very close to what we found in the spectrum. Therefore, the inferred B-field from this line, ∼10 11 G, might be the magnetic field of this system. Verification of this relatively low B field for an HMXB NS could be confirmed by future high-sensitivity observations that cover a broad energy range, e.g. observations with NICER and NuSTAR during outburst.

C O N C L U S I O N
Our 72 ks XMM-Newton observation of the quiescent BeXRP RX J0812.4-3114 revealed a spectrum best described by a soft blackbody component plus a hard power-law, with an absorption component at ≈ 1 keV and/or at ≈ 1.4 keV. The blackbody component implies a large emission region, which we argue likely originates from the whole NS surface. The hard component can be described by a hard power-law ( ∼ 1.3-1.5), likely caused by low-level accretion.
Temporal analyses reveal that the hard component is pulsed at a period of 31.908 ± 0.009 s, slightly longer than the previous studies, from which we estimated a spin-down rate of ≈ 3.63 × 10 −11 s s −1 , which might just be due to orbital Doppler effect. We did not find pulsations in the soft excess (PF 31 per cent), indicating the soft emission has a different origin from the hard emission (consistent with emission from the full NS surface). Long-term light curves reveal variability in the hard component; however, no sign of variability was found in the soft excess. Temporal analysis of the ASM light curves confirms the previously measured orbital period of ≈81.3 d, and that the source returned to a quiescent state after the active epoch between 1997 December 5 and 2000 July 2. Based on the accretion history, we estimated the time-averaged mass accretion rate ( Ṁ ). Assuming the quiescent thermal luminosity is produced by deep crustal heating, we found that J0812 lies above some of the minimum cooling tracks, though with large uncertainty in Ṁ . The NS in J0812 may yet agree with minimum cooling processes. However, it is also possible that the NS in J0812 is too young to have fully cooled after its supernova explosion -a possibility which we estimate to have a 5 per cent chance, given the estimated lifetime of the B0 companion and the time-scale for NS cooling.
J0812 seems to have two distinct X-ray spectral states in quiescence: a soft, or thermally dominated state as observed by Chandra, versus a harder state with possible on-going accretion as observed by our XMM-Newton observation. We suspect the source lies in the propeller regime during the soft state, while it is out of the propeller regime and accreting during the hard state. With these assumptions, we estimate the magnetic field strength of the system to be 8.4 × 10 11 G. Should the ≈1 or ≈1.4 keV absorption feature be real, and represent an electron cyclotron line, then we may further estimate B ∼ 10 11 G, unusually low for BeXRPs, but consistent with the estimate from the propeller argument.

AC K N OW L E D G E M E N T S
COH is supported by NSERC Discovery Grant RGPIN-2016-04602 and a Discovery Accelerator Supplement. COH thanks the organizers and participants of the April 2019 'Investigating Crusts Of Neutron Stars' workshop at the Anton Pannekoek Institute of the University of Amsterdam, supported by JINA-CEE, for discussions. SST was supported by the grant 14.W03.31.0021 of the Ministry of Science and Higher Education of the Russian Federation. WCGH appreciates use of computer facilities at the Kavli Institute for Particle Astrophysics and Cosmology. The work of AYP was supported by RFBR and DFG within the research project 19-52-12013. This work is primarily based on observations obtained with it XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. We also use data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortiu m). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Spectral and temporal analyses in this work made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. This research has also made use of the NASA Astrophysics Data System (ADS) and software provided by the Chandra X-ray Center (CXC) in the application package CIAO.