Observational Properties of Accreting Neutron Stars with an Optically Thick Envelope

Pulsating Ultra Luminous X-ray sources (PULXs) are thought to be X-ray bright, accreting, magnetized neutron stars. Their measured soft ( < 20 keV) X-ray luminosity can exceed the Eddington luminosity for a neutron star (NS) by a few orders of magnitude. Although several scenarios have been proposed to explain the diﬀerent components observed in the X-ray spectrum and the characteristics of the X-ray pulse proﬁle of these system, detailed quantitative calculations are still missing. In particular, the observed soft X-ray pulse proﬁles are almost sinusoidal and show an increase in the pulsed fraction (from 8% up to even 30%) with increasing energy. In this work, we discuss how emission originating from an optically thick envelope, expected to be formed during super-Eddington accretion, can match the observed spectra and pulse proﬁles.

PULXs are point-like X-ray bright sources with apparent luminosities > 10 39 erg s −1 in the 0.1 − 20 keV energy band (see e.g.Kaaret et al. 2017;Fabrika et al. 2021;Mushtukov & Tsygankov 2022;King et al. 2023), fairly exceeding the Eddington luminosity for a ∼ 1.4 object (∼ 3×10 38 erg s −1 ).The mechanism that produces an apparent super-Eddington luminosity is still unresolved, although there are two main candidate scenarios in the case of a neutron star (NS) accretor.One possibility is to invoke geometric beaming by optically thick outflows from the accretion disc (King & Lasota 2020); in this case, the local accretion rate is anyhow limited by the Eddington luminosity.Another scenario is super-Eddington accretion (possibly through an advection dominated disc) onto a NS with a magnetar-like magnetic field.The reduction in the opacity produced by the ultrastrong magnetic field pushes the Eddington limit to higher values, allowing a larger luminosity to escape (Mushtukov et al. 2015;Brice et al. 2021).
The crucial property of PULXs that identifies their central engine as a NS is the observed sinusoidal pulses, which were first discovered in M82 X-2 (Bachetti et al. 2014), and subsequently observed in the PULXs NGC7793 P13 (Fürst et al. 2016;Israel et al. 2017b) and NGC5907 ULX1 (Israel et al. 2017a).The observed pulse profile of these sources remains sinusoidal throughout the whole energy range, but interestingly the pulsed fraction increases with increasing energy, e.g. from ∼ 12% for < 2.5 keV to ∼ 20% for > 7 keV for NGC5907 ULX1 (Israel et al. 2017a).Mushtukov et al. (2021) showed that pulsed fractions 20% make the interpretation of PULXs (e.g.NGC5907 ULX1, NGC7793 P13, M82 X2) as strongly beamed X-ray sources via optically thick ★ E-mail: nabil.brice.17@ucl.ac.uk outflows unlikely.However, no similar quantitative study has been performed for assessing how the pulsed fraction behaves in the super-Eddington accretion scenario.Mushtukov et al. (2017) investigated some of the observational implications of the super-Eddington accretion scenario, focusing in particular on the absence of cyclotron resonance features (CRFs) in PULXs with super-strong magnetic fields.They found that an optically thick envelope around the NS could form when there is a sufficiently large ( 10 19 g s −1 ) mass accretion rate, resulting in a reprocessing of the accretion column emission and hence the disappearance of the CRF.Moreover, the emission from the accretion column emission is thermal as a result of multiple scatterings when the radiation travels through the optically thick envelope walls, and appears as a multicolour blackbody to an observer.These spectral predictions were tested, and Koliopanos et al. (2017), showed that the observed spectra of ULXs (and some PULXs) are well fitted by a dual thermal emission component model.
In this paper we study the observational properties of radiation emitted by an optically thick envelope in more detail with the specific aim to derive the pulse profile and compare it with the observations.The paper is laid out as follows: in §2, we describe the model of the optically thick envelope and its emission following Mushtukov et al. (2017), and specify the method of calculation for the pulse profiles.In §3, the results from analysis of the model is given, with emphasis on those models for which the pulsed fraction trend matches the overall observed one in the sample of PULXs.In §3.2, we apply the model to two astrophysical sources, which were previously analysed by Brice et al. (2021).We discuss the implications in §4.

ENVELOPE MODEL
In the standard picture of accretion onto a highly magnetised NS (see Ghosh & Lamb 1979), the accretion disc is assumed to be truncated at the magnetospheric radius, given by ≈ 7 × 10 7 Λ 1/7 10/7 6 4/7 d,12 where Λ is a dimensionless parameter that depends on the mode of accretion (typically Λ = 0.5 for accretion via thin disc), d,12 is the surface dipole field strength at the magnetic poles in units of 10 12 G, and 39 is the accretion luminosity in units of 10 39 erg s −1 .In addition, the disc penetrates into the magnetosphere by some length, here indicated as .However, the exact physical process and hence a quantitative expression for are still poorly known.In our model, we extrapolate the expression of for a thin disc to the scenario with a thick disc (see Mushtukov et al. 2015) as: where is the height of the disc at the magnetospheric radius, but we impose an upper limit so that / ≤ 0.2 (see Brice et al. 2021 and references thererin).
Let us now consider a simple description of the discmagnetosphere interaction, in which two sets of magnetic field lines cross the plane of the disc at a distance from the star − and respectively.Inside the region bounded by these field lines, the accreting material is channelled onto the NS surface, see figure 1.
If the accretion rate is 10 19 g s −1 , the accreting plasma forms an optically thick envelope and encloses a cavity devoid of material but filled with radiation.(see again figure 1).Since the envelope is optically thick, the radiation injected into the inner cavity cannot immediately escape, and thus thermalises, attaining a blackbody distribution at temperature in (Mushtukov et al. 2017).
The spectrum of the thermal radiation that emerges from the en-velope, and is observed at infinity, depends on the temperature of its outer boundary, out , which in first approximation is related to in and to the local optical depth of the envelope, , via out = in −1/4 . (3) Since the local optical depth is determined by the dynamics of the accreting plasma, and the velocity field is not necessarily uniform inside the envelope, out will vary along the envelope wall as well, and the spectrum of the escaping radiation is expected to be better approximated by a multicolour blackbody, rather than a Planckian at a single temperature.We relate in to the luminosity by integrating the local flux 4 out (where is the Stefan-Boltzmann constant), over the surface of the envelope.Using equation (3), we obtain where d = d ( , ) is the area element of the envelope in spherical coordinates (with the zenith aligned with the direction of the magnetic moment), and is the latitude at which the same magnetic field line meets the surface of the NS.In equation ( 4), we only account for 1/2 of the full accretion luminosity.This is done because the geometry of the accretion column results in half of the emission first entering the cavity and thermalising before escaping, whereas the other half of the emission escapes directly through the exterior envelope wall, i.e. does not contribute to the cavity temperature.The photons that directly escape through the exterior envelope wall may be reprocessed from multiple scatterings, which would result in spectral changes.

Local optical depth
In order to estimate the optical depth of the envelope, we indicate with the distance of the approximate path on which the radiation diffuses outward between the two magnetospheric lines that border the envelope, at each latitude (see again figure 1).This direction is approximately the normal to the envelope surface.The local optical depth is then approximated by where is the electron scattering opacity, and is the local electron density along (which we approximate as constant).
The local density is given by where is the total mass accretion rate onto the NS, is the velocity of the charges along the field lines (which is assumed to be constant along the photon escape path), and is the cross-sectional surface area perpendicular to the magnetic field lines, i.e. perpendicular to the velocity direction of the charges.In this case, we have assumed the mass transfer rate to each hemisphere to be equally distributed for all magnetic field lines.
Finally, can be calculated from the geometry of the envelope cross-section and for purely dipole magnetic field lines is given by where is the magnetic latitude and = cos .We refer to appendix A for more details.

Dynamics
To obtain an expression for , we consider the dynamics of a charge moving along a dipole field line of a NS, which rotates around a spin axis with angular velocity Ω.The NS magnetic moment is misaligned from the spin-axis by an angle .In the following, we use a orthonormal coordinate system in which the -axis is aligned with the magnetic moment and refer to this as the magnetic reference frame.The radial coordinate of the dipole field lines, in spherical polar coordinates, is given by where and are the latitude and the azimuthal angle in the magnetic axis reference frame, is the maximal distance that the field line reaches from the source.In our model, we use = − 1/2 for all (e.g.see appendix A).
The charge moves under the influence of gravitational and centrifugal forces.The gravitational force contributes a potential field grav ( ) where Ω = ( / 3 ) 1/2 is the Keplerian angular velocity at a distance from the point mass.The centrifugal force is taken into account by introducing a constraint to the Lagrangian: where = ( , , ) is the distance of the moving charge from the spin-axis (see appendix B for a derivation of the full expression, given by equation (B6)).Hence, the Lagrangian of a charge with mass moving on the dipole field line given by equation ( 8), under the influence of gravitational and centrifugal forces is We define the parallel velocity of the charge to be and thus the effective potential is given by eff = Since the effective potential includes the gravitational and rotational potentials, the NS spin affects the velocity along the field lines.Note the effective potential is conservative, so we can integrate the Euler-Lagrange equations to obtain an analytic expression for the parallel velocity along the magnetic field lines, where 0 = cos 0 is the cosine of the initial latitude angle of the charge, 0 is the distance of the charge at = 0 from the spin-axis and 0 is the charge initial velocity.We give the derivation of the initial angle expression for a geometrically thin disc in appendix C. Following Mushtukov et al. (2017), we approximate the initial velocity of the particle at any disc height to be the characteristic thermal velocity of the protons at the inner disc radius (for a geometrically thin disc).This gives 0 ≈ 0.056 where is the speed of light, is the adiabatic index and is the hydrogen mass fraction of the accreting material.

Calculation of the pulsed fraction
Once in and are known, out can be calculated for every point on the outer boundary of the envelope using equation (3).For each , the specific flux along a particular line of sight (LOS) is given by where n is the normal to the surface of the envelope at , ℓ is the LOS unit directional vector, ( ) is the Planck function at temperature , and Δ is the area element of the envelope surface at .Note that n • ℓ changes according to the NS phase of rotation, denoted here by .
The total specific flux along a particular LOS, , is obtained by summing over all the points on the surface of the envelope that are in view, i.e.

= in view
, (ℓ). (17 In particular, and depending on the viewing geometry, the points of the envelope can be obscured by other parts of the envelope and/or by the accretion disc.The points in view, hence, are dependent on , and hence = ( ).The specific flux is then integrated over some energy interval, Δ, to obtain In addition to the flux from the points on the envelope, we include the disc contribution to the total flux.This is done in a similar way to equation ( 16), where the disc temperature at the point is used instead of out , and the disc angular area element is used instead of Δ .We use the standard multi-colour disc temperature profile for a thin disc (Makishima et al. 1986), which is given by For simplicity, we consider all disc points from the inner disc radius, , to infinity.Finally, we compute the pulsed fraction for the energy interval Δ using the definition where Δ,max = max Δ and Δ,min = min Δ denote the maximum and minimum values of Δ during a NS rotational period.

NUMERICAL RESULTS
We now consider a scenario in which the accreting material flows through a thin accretion disc (Shakura & Sunyaev 1973), that is perpendicular to the spin axis, and then through the envelope as described in §2.
For each particular LOS, we calculated the phase-dependent spectrum and the pulsed fraction by ray-tracing the thermally emitting points on the envelope and the accretion disc to an observer at infinity.To this end, we re-adapted the numerical ray-tracing code described in Taverna & Turolla (2017), which was originally developed to calculate the spectra of magnetar bursts originating from a magnetically trapped fireball.We modified the code to allow for magnetic field lines reaching several hundred NS radii and added a planar disc component.The code takes into account the visibility of the envelope and disc surface points by checking for (self-)shadowing by envelope and disc.We considered a set of models, with each one being characterised by five parameters: the accretion luminosity in units of 10 39 erg s −1 , 39 , the (surface dipole) magnetic field strength at the pole in units of 10 12 G, 12 , the tilt of the dipole magnetic moment with respect to the spin axis, , the observer LOS angle (i.e. the inclination) with respect to the spin-axis of the NS, , and the spin period of the NS, .

O R I G
We calculated the synthetic spectrum in the 0.1 − 20 keV energy range and for 60 grid-points in phase.For every model, we found that the phase-resolved synthetic spectrum varies sinusoidally in each energy interval considered.An example of the energy resolved pulsed profile, corresponding to a model for NGC7793 P13 considered in §3.2, is shown in figure 2. We also analysed the phase-averaged synthetic spectrum by using the IDL routine MPCURVEFIT to fit it to a multi-color disc plus blackbody model (mcd+bb).In each case, we obtained a statistically good fit.
In the following, the characteristic temperatures of the components from the best-fit spectral model are indicated by mcd for the mcd component and bb for the bb component, which we refer to collectively as the best-fit component temperatures.In addition, we calculated the maximum temperature on the envelope wall, max , which we refer to as the envelope temperature, and the inner disc temperature, disc , as diagnostic measures that are independent of the ray-tracer and spectral fitting.
We present the results of varying the luminosity and dipole magnetic field strength on the envelope temperature and inner disc temperature in figures 3 and 5 alongside the best-fit component temperatures for the same models.These plots confirm that the best-fit component temperatures in the mcd+bb spectral fit (of the phaseaveraged spectra), mcd and bb , accurately tracks the underlying model variables disc and max .Otherwise, we found that changes in , , and did not significantly alter the best-fit component temperatures nor the underlying model variables.
In figures 4, 6, 7, and 8 we present the results of varying the luminosity, dipole magnetic field strength, , and respectively on the pulsed fraction profile in energy.Otherwise, we observed no change to the pulsed fraction when was varied within its valid range for the models.We found that all models shared a common characteristic in their pulsed fraction profile with energy, namely a change in the gradient of the pulsed fraction profile, which we refer to as the "break" in the pulsed fraction profile.
Below, we discuss the response of the synthetic spectrum and pulsed fraction profile when varying each model parameter in turn.We describe how each parameter relates to the underlying model variables as a result of the choices made in the modelling.
In the first case, the reason that a variation in the accretion luminosity parameter leads to a change in the underlying model variables can be immediately understood from equations ( 19) and (4).Namely, the inner disc temperature depends directly on the accretion rate in equation ( 19) and the envelope temperature is determined by the cavity temperature and minimum optical depth, both of which are related to the accretion luminosity and accretion rate in equations ( 4) and ( 6) respectively.
Our analysis of the synthetic spectra show that the best-fit temperature trends are aligned with the underlying model variable trends.In the primary cases shown in figure 3 where the model parameters were fixed to 12 = 1.0, = 10 • , = 30 • , = 1s, we found the increase in the mcd characteristic temperature to be more rapid than the increase in the bb characteristic temperature.Using a simple linear regression model for the logarithmic luminosity and temperature values, we calculated the mcd temperature is related to the accretion luminosity as mcd ∝ 0.44 39 whereas the bb temperature is related to the accretion luminosity as bb ∝ 0.39 39 .This trend, where the power law index is steeper for the component associated with the disc, is reflected in the relationship of the accretion luminosity with the underlying model variables, i.e. disc ∝ 0.46 39 and max ∝ 0.39 39 .In addition, the power law index for the inner disc temperature and envelope temperature were consistent for the alternate magnetic field strengths shown in figure 3.
Since a variation of the accretion luminosity changes the synthetic spectrum, the energy-dependent pulsed fraction is subsequently affected.In particular, we find that models with a higher accretion luminosity (with all other parameters held constant) exhibit a lower pulsed fraction in all energy intervals, as shown in figure 4 dition, the break from steep increase in pulsed fraction with energy occurs at a higher energy for models with higher accretion luminosity.
The fact that the pulsed fraction is lower for all energy intervals, including for energies after the break, can be attributed to the specifics of the temperature distribution on the envelope.Namely, the highest temperatures lie in a small region close to where the disc meets the magnetosphere, which means the pulsed fraction at high energies is representative of the change in effective area, i.e. the projected area onto a plane perpendicular to the observer's LOS, of these high temperature regions during one full rotation of the system.For increased accretion luminosity, the higher temperature regions make up a larger proportion of the total envelope area, and hence there is a smaller change to the effective area during one full rotation.
The fact that the pulsed fraction break occurs at higher energy can be attributed to a hotter inner disc due to an increased accretion luminosity (and simultaneously smaller magnetospheric radius).It is the temperature of the hottest regions of the disc that determines at what energy the disc flux peaks, after which the envelope emission (with the higher characteristic temperature) dominates.
A change in the dipole magnetic field strength parameter leads to a change in the underlying model variables, which (similarly to the analysis of the accretion luminosity parameter) we can see from the equations of the envelope model outlined in §2.In this case, the inner disc temperature and envelope temperature both depend inversely on the magnetospheric radius, which in turn depends directly on the magnetic field strength, given in equation ( 1).
The best-fit temperatures from the synthetic spectra again show aligned trends to the underlying model variables, with an increase in the dipole magnetic field strength resulting in a decrease of the bestfit component temperatures.For the models shown in figure 5, we fixed model parameters to 39 = 5, = 10 • , = 1s, = 30 • , and varied the magnetic field strength parameter between 12 = 0.1 to 3.2.For these models, we found (using a linear regression model on the logarithmic d and temperature values) the relation between bestfit component temperatures and the dipole field strength parameter to be mcd ∝ −0.42 d,12 and bb ∝ −0.28 d,12 .This aligns with the relations between the underlying model variables and the dipole field strength parameter, where disc ∝ −0.43 d,12 and max ∝ −0.29 d,12 .Similarly to the case of increasing accretion luminosity, the  changes to the synthetic spectra from an increase in the magnetic field strength results in changes to the pulsed fraction profile.In the case of increased magnetospheric radius, the higher temperature regions on the envelope make up a smaller proportion of the total envelope area, and so there is a larger change to the effective area during one full rotation.Hence the pulsed fraction is larger for all energy intervals.The pulsed fraction break occurring at a lower energy also naturally follows from the explanation in the accretion luminosity case.In this instance, the magnetospheric radius is increased with increasing magnetic field strength, which results in the truncation of the disc farther from the compact object.Thus the inner disc temperature is cooler and so the peak of the disc flux is at a lower energy.
In testing the model parameters and , we found no significant change to the best-fit component temperatures nor in the underlying model variables when these parameters were varied within their  respective range of validity, except for parameter values close to the limits of validity, e.g. when the period is close to the propeller regime.However, an examination of the model behaviour in these scenarios is beyond the scope of this work, since the parameter limits are only approximated.Despite no change in the underlying model variables, we found that a variation of , well within the validity limits, produces a substantial difference to the pulsed fraction, which we show in figure 7.In this case, the pulsed fraction is not changed as a result of an altered contribution from the disc and the envelope emission components.Instead, the pulsed fraction varies with increasing tilt because of the altered viewing geometry to the envelope from the observer LOS.The high temperature region on the envelope is located close to the magnetospheric equator so when the magnetic moment is tilted relative to the spin axis, the high temperature region presents an effective area that depends on the phase of the rotation.With increasing , the change in the effective area during one full rotation is larger and hence producing a larger pulsed fraction at all energies without a difference in the energy at which the break occurs.
One notable divergence is in the pulsed fraction profile for the model with = 30 • , which can be seen to flat-line after the break.This change can be understood as a result of the specifics of the viewing geometry to the envelope, in particular that the points of hottest temperature on the envelope are visible throughout the entire rotation.Hence, the envelope flux minimum (even at high energies) is non-vanishing, which results in an upper limit of less than 100% to the pulsed fraction.
In the case of varying the parameter , the underlying model variables are completely independent of the viewing angle so, again, there is no change to the phase-averaged spectrum for most models.However, we found that at high inclinations, 60 • , the best-fit temperature for the mcd component decreased.This can be understood as an effect of a partial obscuration of the inner disc by the envelope.The fraction of the inner disc obscured by the envelope increases as the LOS angle increases so the hotter parts of the disc contribute less to the overall disc emission.Hence, the resulting characteristic mcd temperature obtained from the fit is lowered.
For the same reasons as in the analysis of , the changes of produce a substantial difference to the pulsed fraction, which we show in figure 8.However, In this case, the break of the pulsed fraction profile occurs at a lower energy for models with larger .This is an effect of the inner disc obscuration by the envelope (described above), which results in a lower disc specific flux at energies for which the inner disc is the main contributor.Hence, the energy at which the envelope flux dominates is down-shifted.

Parameter constraints
The envelope model as described in §2 is valid provided that the accretion column emission is thermalised inside the magnetosphere cavity.This condition holds when the envelope is optically thick, i.e. 1.As previously discussed by Mushtukov et al. (2017), the condition of an optically thick envelope translates to a constraint on the model parameters for self-consistent models.In the case of models with = 0 • , the constraint can be approximated by an accretion luminosity bound: 39 1/4 12 .In the case of models with > 0 • i.e. a tilted magnetic axis to the spin axis, the equivalent minimum accretion luminosity bound depends on and must be computed numerically.
Another constraint on model parameters comes from the fact that the radiation contained in the cavity exerts a pressure on the envelope walls.Since the material is channelled along magnetic field lines, the magnetic field must be strong enough to confine the accretion flow against the radiation pressure.As a rough approximation, this means that the magnetic pressure (at the magnetospheric radius) is required to be higher than the radiation pressure for the model to be self-consistent.This translates into an upper bound for the accretion luminosity as a function of the magnetic field strength.The exact maximum luminosity also depends on , due to the effect of the magnetic field tilt on in , and thus must be computed numerically.
We calculated the minimum and maximum accretion luminosity, fixing = 1s and for various values of , as a function of the magnetic field strength.The results are shown in figure 9. We have truncated the curves for the maximum accretion luminosity at the approximate magnetic field strength for which the maximum accretion luminosity becomes lower than the minimum accretion luminosity for the cor- responding models.Beyond this point, the cavity temperature is not physical.
Specifically for the minimum accretion luminosity, the sharp increase reflects the progressive onset of a propeller mechanism.This is the case when centrifugal forces from the rotating magnetosphere (which are included in the velocity expression) become progressively dominant for larger magnetosphere radii, which in turn decrease the velocity along the magnetic field lines and eventually halting accretion altogether.
The self-consistency conditions rule out models with accretion luminosities 1039 erg s −1 or 10 40 erg s −1 , and in general models with tilts > 40 • .Since the accretion luminosity is twice the envelope luminosity, i.e. the luminosity used to calculate the emission as used in equation ( 4), the maximum luminosity constraint would exclude astrophysical sources with observed X-ray luminosities 5 × 10 39 erg s −1 .However, the anisotropy of the emission from the envelope means that the observed luminosity (calculated assuming isotropic emission) can be greater than the envelope luminosity.We found that this luminosity amplification increases from a factor ∼ 1.1 to 2 for increasing viewing angle .Thus, for real astrophysical sources, we can approximate the apparent X-ray luminosity as up to the accretion luminosity.We leave the full investigation of the luminosity amplification in this model for a future work.

Application to sources
We present the application of our the model to two PULX sources: NGC 7793 P13 (Fürst et al. 2016), and NGC 5907 ULX1 (Israel et al. 2017a), below.These sources were chosen because the pulsed fraction data was available for multiple energy bands and because they were the focus of our previous study (Brice et al. 2021).Our aim was to identify whether the optically thick envelope model as described in §2 is capable of explaining the spectra and pulsed fraction data.In particular, we have attempted to find a model that simultaneously agrees with the modelled observational spectrum (as calculated by Koliopanos et al. 2017) and the pulsed fraction profile of each source.
To find a suitable model for each source, firstly we narrowed the possible set of model parameters (specifically accretion luminosity and magnetic field strength) by using the temperatures from the bestfit mcd+bb spectral model from Koliopanos et al. (2017).Secondly we varied the values of and to find a geometric configuration that best reproduces the pulsed fraction profile of the observed X-ray data.These data are taken from Fürst et al. (2016) and Israel et al. (2017a) for the sources NGC7793 P13 and NGC5907 ULX1 respectively.

NGC 7793 P13
The spectrum of the XMM-Newton observation of NGC 7793 P13 (obs.0748390901), in the 0.3 − 10 keV energy band, can be fitted by a combination of two thermal components, one representing the disc (mcd model), and one representing the hotter envelope (bb model), with mcd ≈ 0.6 keV and bb ≈ 1.7 keV, respectively (see Koliopanos et al. 2017).The X-ray luminosity for the XMM-Newton observation is ∼ 7 × 10 39 erg s −1 (assuming isotropic emission and a distance of ≈ 3.6 Mpc), and the spin period of the source is ≈ 0.4 s (Fürst et al. 2016).
We generated the synthetic spectra specifically for models with Brice et al. and = 45 • was able to closely match the pulsed fractions as reported by Fürst et al. (2016) in the energy intervals 0.3−0.8,0.8−1.6,1.6− 3.0, 3.0 − 5.0, 5.0 − 7.0, 7.0 − 9.0 keV and 3.0 − 5.0, 5.0 − 7.0, 7.0 − 10.0, 10.0 − 20.0 keV (for XMM-Newton and NuSTAR data sets respectively).The results are shown in figure 10.For this particular model, the overall trend of the pulsed fraction profile is consistent with the observational data and the higher energy (> 3.0 keV) pulsed fraction values are all within the 1 error margin.

NGC 5907 ULX1
NGC5907 ULX1 is the brightest PULX observed to date, with a luminosity that exceeds ∼ 500 times the NS Eddington luminosity in some observations (Israel et al. 2017a).For our application, we attempted to match to the spectral data from XMM-Newton observation 072956AS1301 (Koliopanos et al. 2017), which suggest an X-ray luminosity of ∼ 8.3 × 10 40 erg s −1 (assuming isotropic emission and a distance of ≈ 17 Mpc).The best-fit mcd+bb model in this case has characteristic temperatures mcd ≈ 0.7 keV and bb ≈ 1.4 keV respectively.We used a spin period of ≈ 1.1 s for our modelling (Israel et al. 2017a).
To match the observed flux of the source would require 39 ∼ 100, which is above the maximum accretion luminosity of our models.Instead, by using an accretion luminosity of 39 = 10 and a dipole magnetic field strength parameter of 12 = 0.8, our models had a similar synthetic spectrum to the observed spectrum, with best-fit temperatures of mcd ≈ 0.7 keV and bb ≈ 1.4 keV.This result suggests that the spectrum can be adequately described by the optically thick envelope model and that there is an additional mechanism responsible for the amplification of the observed flux.
Together with the accretion luminosity and magnetic field strength parameters used above, we found that a model with = 10 • and = 25 • was able to reproduce the pulsed fractions as reported by Israel et al. (2017b) in the energy intervals < 2.5 keV and > 7 keV.The results are shown in figure 11.

DISCUSSION AND CONCLUSIONS
The model presented in this work is applicable in scenarios where the accreting NS has a magnetic field 10 11 G, and the accretion rate is 10 39 erg s −1 , such that an optically thick envelope can develop.If there is little or no beaming and if the dipole magnetic field strength is typical of X-ray pulsars, PULXs fall in the range of applicability.
This study is the first to account for the tilt of the dipole moment from the spin-axis in the calculations of the velocity of the accreting material and the optical depth of the envelope.We found an analytical solution to the velocity of the accreting material, which decreased computation time and enabled comparison of different model parameters.Crucially, a tilt of the dipole moment leads to a non-homogeneous temperature distribution on the envelope, which contributes to a greater pulsed emission as the NS rotates.
Our aim was to understand the observational implications of an optically thick envelope.In particular, we studied the pulsed fraction profile and phase-averaged spectrum, comparing the synthetic model data with observational data from PULXs.
We found that our models reproduced the trend of increasing pulsed fraction with increasing photon energy, which is observed in PULXs (and in fact a large pulsed fraction at high photon energies is a typical feature of X-ray pulsars in general).In particular, we were able to match the pulsed fraction profile of the two PULXs studied for this work: NGC7793 P13 and NGC5907 ULX1.In fact, we were able to obtain a range of pulsed fractions, from < 5% to ∼ 50% (depending on the tilt and viewing angle).Thus, the optically thick envelope scenario is also consistent with low pulsed fractions, such as that measured for the PULX NGC1313 X-2 by Sathyaprakash et al. (2019), provided the viewing angle is close to the spin-axis (as suggested by these authors) and the dipole magnetic field tilt is small.
We also found that the synethic spectra produced by our models are broadly consistent with a double thermal spectral model, which has previously been used to analyse ULX spectra (Koliopanos et al. 2017).In particular, we were able to find model parameters that reproduced the mcd+bb best-fit parameters for NGC7793 P13 and NGC5907 ULX1.
While we were able to find model parameters that reproduced the spectral best-fit parameters of NGC5907 ULX1, we used a lower accretion luminosity parameter in the models than the value of the observed luminosity.This was necessary because the value of the observed luminosity is an order of magnitude larger than the maximum accretion luminosity for the optically thick envelope models (see §3.1), which suggests the need for another mechanism to explain the observed luminosity, e.g.luminosity amplification through geometrical beaming by optically thick outflows from a super-Eddington accretion disc (King & Lasota 2020;Mushtukov et al. 2021;Mushtukov & Portegies Zwart 2023).
A possible mechanism for luminosity amplification in the optically thick envelope model comes from the accretion column emission that is not injected into the cavity.We did not include this part of the emission in our modelling since we focused on the envelope and disc emission.This emission from the accretion column accounts for half of the total accretion luminosity and it is funnelled along a narrow cone aligned with the magnetic axis (since the envelope exterior is optically thick).Hence, this emission is highly an-isotropic and would contribute significantly to a luminosity amplification effect for observers that are within the emission cone.The contribution of the accretion column to the spectrum requires a full radiative transfer calculation, which is beyond the scope of this work.However, we speculate that this emission would be dominated by hard X-rays since it comes from a relatively small region.In fact, Walton et al. (2018) find a hard X-ray power-law excess in the pulsed part of the broadband spectra of NGC7793 P13 and NGC5907 ULX1 after fitting the lower energy spectra with an mcd+bb model, which could indicate that this component of the accretion column emission does contribute in these sources.
The spectral model best-fit parameters for both of the PULXs studied here suggest that the dipole magnetic field strength is < 10 12 G at the NS surface.However, the locally super-Eddington luminosity in the accretion column requires a magnetar-like surface magnetic field strength, i.e. 10 14 G (Mushtukov et al. 2015).This contrast can be explained by the presence of significant multipole components in the magnetic field, i.e. a more complex magnetic field topology than a pure dipole (Brice et al. 2021).
As part of our modelling, we made several simplifying assumptions that may change the specifics of the phase-resolved synethic spectra.One of these assumptions is in our use of the Shakura & Sunyaev (1973) thin accretion disc model and its spectral profile.This could be inaccurate at super-critical accretion rates (locally super-Eddington for the disc) because the disc becomes geometrically thick due to radiation pressure and the disc may be advection dominated (Poutanen et al. 2007;Chashkina et al. 2019).
A more suitable disc model for this regime, (e.g.see Chashkina et al. 2017Chashkina et al. , 2019)), would take into account the effect of the disc thickness on the disc emission and on the envelope emission.The spectral profile of a thick disc could differ substantially from that of a thin disc due to changes in the density and viscosity of the accreting material.In addition, a geometrically thick disc changes the spectral profile of the envelope by blocking a portion of the emitting area of the envelope, which would result in a change of the temperature distribution on the envelope wall.
In our modelling of the two PULXs, the accretion disc is in the super-critical regime and so is expected to produce outflows (Shakura & Sunyaev 1973).Material outflows would change the total mass accretion rate through the envelope (Mushtukov et al. 2019) and hence change the spectral profile from both the disc and envelope.In addition, certain LOS would be affected by the collimation of emission by optically thick outflowing material, in a similar scenario to that presented by King & Lasota (2020).In this work, we assumed the accretion luminosity to be equal to the mass accretion rate without any mass loss in outflows for simplicity.In principle, a ray-tracing computation could include the effect of collimation by outflows but we leave this for a future work.
In conclusion, we found that the observed pulsed fraction from NGC7793 P13 and NGC5907 ULX1 can be reproduced in our modelling of an optically thick envelope.In contrast, Mushtukov et al. (2021) found that the pulsed fraction values from these sources would be difficult to explain in the scenario with highly geometrically beamed emission.

Figure 1 .
Figure 1.A diagram of the two bounding dipole field lines (blue solid curves) reaching out to radii and − .The purple line shows, for a particular latitude, the shortest path along which the radiation diffuses outward and is labelled by .The NS is at centre (black circle).

Figure 2 .
Figure 2. The normalised specific flux in phase.The phase-resolved specific flux is normalised with the phase-averaged specific flux, .The pulsed profile was taken from the phase-resolved synthetic spectrum of a model with 39 = 10, 12 = 0.55, = 0.4 s, = 10 • , = 45 • .

Figure 3 .
Figure 3. Inner disc temperature and envelope temperature, disc and max respectively, and the characteristic temperatures from the best-fit mcd+bb spectral model, mcd and bb respectively, as a function of the accretion luminosity.We used fixed model parameters: 12 = 1.0, = 10 • , = 1s, and = 30 • .In addition, max is shown for models with various magnetic field strengths for comparison.

Figure 5 .
Figure 5. Same quantities as in Figure 3, here shown as a function of the (dipole) magnetic field strength, dip .We used fixed model parameters: 39 = 5, = 10 • , = 1s, and = 30 • .Another line of max for models with 39 = 10 is shown for comparison.

Figure 6 .
Figure 6.Pulsed fraction as a function of energy for models with different values of the magnetic field strength parameter, shown with each line.We used fixed model parameters: 39 = 5, = 10 • , = 1s and = 30 • .

OFigure 9 .
Figure 9.The accretion luminosity bounds for an optically thick envelope, > 1, as a function of the dipole field strength.All models are calculated with = 1s.The black solid, dashed, dot-dashed, and dotted lines indicate the minimum accretion luminosity for models with = 0 • , 20 • , 30 • , 40 • respectively.The minimum accretion luminosity line for models with = 10 • overlaps with the minimum accretion luminosity line for models with = 0 • so are omitted.The progressively lower solid purple lines indicate the maximum accretion luminosity for models with = 0 • , 10 • , 20 • , 30 • , 40 • .These lines have been truncated at the last d for which the maximum luminosity is greater than the minimum luminosity for the corresponding models.The shaded red region indicates the set of parameters for which the disc is super-critical.

Figure 10 .
Figure 10.Pulsed fraction as a function of energy for the source NGC7793 P13.The red and blue points with error bars are the observed XMM-Newton and NuSTAR data, respectively, reported by Fürst et al. (2016).The pulsed fraction profile shown is from a model with 39 = 10, 12 = 0.55, = 0.4 s, = 10 • .The blue shaded region indicates the pulsed fraction profiles obtainable from this model by varying the inclination (from = 10 • for the region bounding line at the bottom to = 70 • for the region bounding line at the top).The dashed line in the blue region indicates the pulsed fraction for a model with = 45 • , which we found to most closely match the pulsed fraction data.

Figure 11 .
Figure 11.Pulsed fraction as a function of energy for the source NGC5907 ULX1.The red bars are pulsed fraction values given in specific energy intervals (0 − 2.5 keV and 7 − 12 keV) as reported by Israel et al. (2017a).The pulsed fraction profile shown is from a model with 39 = 10, 12 = 0.8, = 1.1 s, = 10 • .The blue shaded region indicates the pulsed fraction profiles obtainable from this model by varying the inclination (from = 10 • for the region bounding line below to = 70 • for the region bounding line above).The dashed line in the blue region indicates the pulsed fraction profile for the model with = 25 • , which we found to most closely match the pulsed fraction data.