Super-Eddington Emission from Accreting, Highly Magnetised Neutron Stars with a Multipolar Magnetic Field

Pulsating ultra-luminous X-ray sources (PULXs) are characterised by an extremely large luminosity ($>10^{40} \text{erg s}^{-1}$). While there is a general consensus that they host an accreting, magnetized neutron star (NS), the problem of how to produce luminosities $>100$ times the Eddington limit, $L_E$, of a solar mass object is still debated. A promising explanation relies on the reduction of the opacities in the presence of a strong magnetic field, which allows for the local flux to be much larger than the Eddington flux. However, avoiding the onset of the propeller effect may be a serious problem. Here, we reconsider the problem of column accretion onto a highly magnetized NS, extending previously published calculations by relaxing the assumption of a pure dipolar field and allowing for more complex magnetic field topologies. We find that the maximum luminosity is determined primarily by the magnetic field strength near the NS surface. We also investigate other factors determining the accretion column geometry and the emergent luminosity, such as the assumptions on the parameters governing the accretion flow at the disk-magnetosphere boundary. We conclude that a strongly magnetized NS with a dipole component of $\sim 10^{13} \text{G}$, octupole component of $\sim10^{14} \text{G}$ and spin period $\sim1 \text{s}$ can produce a luminosity of $\sim 10^{41} \text{erg s}^{-1}$ while avoiding the propeller regime. We apply our model to two PULXs, NGC 5907 ULX-1 and NGC 7793 P13, and discuss how their luminosity and spin period rate can be explained in terms of different configurations, either with or without multipolar magnetic components.


INTRODUCTION
Ultra-luminous X-ray sources (ULXs) are X-ray bright compact objects inside or near the optical extent of a galaxy but off-nucleus. Their observed X-ray luminosity exceeds the Eddington limit for a stellar mass, ∼ 1 − 10 , compact object ( > 10 39 erg s −1 ). The nature of the compact object and the exact mechanism which powers the observed luminosity remains debated to this day (see Kaaret et al. 2017 for the most recent review).
For a long time, ULXs were thought to be either intermediate mass black holes (Colbert & Mushotzky 1999) accreting at sub-Eddington rates or stellar mass black holes accreting at super-Eddington rates (King et al. 2001). The discovery of a pulsating ULX (PULX) M82 X-2 in 2014 by Bachetti et al. (2014) revolutionised our understanding of these sources. For the first time, firm evidence was presented that a ULX could host an accreting, neutron star (NS). Since then, further discoveries of PULXs have been made (Fürst et al. 2016;Israel et al. 2017b,a;Rodríguez Castillo et al. 2019;Sathyaprakash et al. 2019). At present, there is no evidence that PULXs differ from non-pulsating ULX basing on their spectral ★ E-mail: nabil.brice.17@ucl.ac.uk properties alone, which hints at the possibility that there are may be many more NS-powered ULXs than previously thought (e.g. King & Lasota 2016). This sparked a new theoretical effort aimed at investigating and modelling the physics of ultra magnetized accreting NSs, which is crucial to our understanding of ULXs.
The very first investigations of whether a strongly magnetized NS may be capable of emitting above its Eddington limit was presented several decades ago by Basko & Sunyaev (1976). In their model, it is assumed that the accretion disk is truncated far from the NS, due to the interaction between the disk and the star magnetic field at the magnetospheric radius. The accreting matter is then channelled along magnetic field lines onto the polar caps. However, this particular model was primarily concerned with how the maximum luminosity can be affected by the presence of a funnelled accretion column and a radiative shock occurring above the star surface. The model neglected the effects induced by the strong magnetic field on the plasma opacities and on the radiation field and consequently, Basko & Sunyaev (1976) found that the maximum luminosity could only be increased by a factor of a few above the Eddington limit. This alone is insufficient to explain the super-Eddington luminosity observed in PULXs, e.g. in M82 X-2. Lyubarskii & Syunyaev (1988) expanded the previous model by calculating the structure of the slow sinking region below the shock in two dimensions. This provided the basis for the most recent model by Mushtukov et al. (2015), which also includes the opacity reduction effect due to the NS's strong magnetic field. Mushtukov et al. (2015) concluded that a maximum luminosity of up to 10 40 erg s −1 can be sustained by the accretion column, depending on the strength of the magnetic field. However, problems arose in trying to explain the observed luminosity and spin period of NGC 5907 ULX-1 due to the assumption of a pure dipole field. Namely, using the model of Mushtukov et al. (2015), the required dipole magnetic field strength of NGC 5907 ULX-1 would place the source in the propeller regime (Israel et al. 2017a). In such a regime, the propeller effect is thought to halt the accretion process due to the transfer of angular momentum at the magnetopsheric radius from the NS to the accretion disk (see Illarionov & Sunyaev 1975).
In order to overcome this issue, Israel et al. (2017a) proposed that the pure dipole field assumption should be relaxed and that higher order multipole moments of the magnetic field may be dominant close to the surface of the NS. Such a magnetic field configuration is reminiscent of the one suggested for magnetars  and is supported by observational data from the magnetar SGR 0418+5729 (Tiengo et al. 2013), the isolated NSs RX J0720.4-3125 and J1308.6+2127 (Borghese et al. 2015(Borghese et al. , 2017, and the millisecond pulsar PSR J0030+0451 (Bilous et al. 2019). A multipolar magnetic field configuration would avoid the problem of the propeller effect induced by a super-strong dipolar component while also allowing for sufficient opacity reduction and in turn the release of a substantial super-Eddington luminosity.
In this paper, we construct a model of column accretion onto a NS that allows for a multipolar magnetic field. As a basis for our calculation, we use the model described by Mushtukov et al. (2015). The paper is laid out as follows: in §2, we give the basic equations and assumptions used in the model calculations. In §3.1 we present results of our numerical computations and compare with the previous model of Mushtukov et al. (2015). In §3.2 we present models based on different degrees of X-mode polarization, while in §3.3 we discuss how our results are affected by different assumptions on disk parameters. The robustness of some of our model assumptions and their regime of validity is discussed in §3.5. The results of an application to two astrophysical sources, which were suggested to have multipolar magnetic fields, is given in §3.6. Finally, we discuss our results and present our conclusions in §4 and §5.

ACCRETION COLUMN MODEL
Through the whole paper, we assumed a neutron star mass and radius of = 1.4 and = 10 6 cm, respectively. As a basis for our calculation, we consider the accretion column model originally developed by Basko & Sunyaev (1976), according to which the free falling plasma is efficiently decelerated in a radiative shock above the surface of the neutron star. Below the shock, the plasma slowly sinks toward the surface and liberates its gravitational potential energy in the form of X-ray radiation. This region is referred to as the sinking region.

Basic Equations
For completeness and to introduce our notation, we summarise the radiative hydrodynamical equations that describe the plasma flow Below it, the plasma sediments towards the surface of the NS. The dashed blue lines represent magnetic field lines confining the accretion column. The maximum shock height is labelled and half the column base width is labelled as 0 /2. For a given , the coordinates of the shock boundary are ( , ). Alternatively, for a given ℎ, the coordinates of the shock boundary are ( ℎ /2, ℎ).
in the region below the radiative shock, using the same assumptions as in Mushtukov et al. (2015).
Since the accretion column is localized at the magnetic poles and the height of the sinking region is less than the star radius, we neglect the curvature of the magnetic field lines and adopt an orthonormal coordinate system ( , ℎ), with the ℎ-axis along the magnetic field lines. We define = 0 to be the centre of the accretion column and ℎ = 0 to be at the surface of the NS.
We indicate with the maximum height of the radiative shock, and with 0 the width of the accretion column base.
is the height of the shock at width along the base and ℎ is the width of the sinking region at height ℎ above the surface. These geometrical quantities are illustrated in figure 1, which shows a vertical crosssection of the sinking region. The footprint of the accretion column on the surface is an annulus with arc length 0 and width 0 . The area of the accretion column base is given by = 0 0 . The calculation of 0 and 0 is detailed in §2.3.
Inside the sinking region, we consider a steady state flow, with velocity directed along the magnetic field lines only. A number of additional assumptions are made in order to simplify the equations further. First, we assume that the radiation pressure, rad , dominates the gas pressure, gas . Following Mushtukov et al. (2015), we assume that the density and velocity profiles are independent of , and coincide with the profiles at the center of the column (see §2.2). Finally, we assume the energy flux to be dominated by the radiative flux, which is a the sum of two components, one directed vertically along the field lines, , and one perpendicular to them, ⊥ .
Accordingly, the equations of continuity, momentum, and en-ergy can be written respectively as where is the plasma density, is the plasma velocity, is the accretion rate, and are the mass and the radius of the neutron star, respectively. The energy fluxes and ⊥ can be expressed as where is the radiation energy density, and and ⊥ are the angle and energy averaged Rosseland mean opacities in the direction parallel and perpendicular to the magnetic field lines, respectively (see §2.4).
In the sinking region below the shock, the flow is decelerated to a velocity much less then the free fall velocity, ℎ 2 /2 − /( + ℎ) 2 , so hydrostatic equilibrium can be assumed in the vertical direction, and thus equation (2) is simply Moreover, since the energy flux is dominated by the radiative flux, equations (4) and (5) can be written as Equations (7) and (8) can be integrated to calculate the radiation pressure distribution within the sinking region, and hence the whole structure of the latter. First, the parallel flux, , can be obtained by coupling equation (7) with the momentum conservation equation (6). This yields which is the local Eddington flux, Edd ( , ℎ). Then, the perpendicular flux, ⊥ , can be obtained integrating the energy conservation equation (3) in and by assuming ℎ ≈ constant in . This yields where ⊥,esc (ℎ) denotes the perpendicular flux escaping from the sinking region at height ℎ, and we used the boundary conditions Equations (7) and (8) can now be integrated, yielding and rad,⊥ ( , ℎ) = where rad, and rad,⊥ are the radiation pressure obtained by integrating the PDE in ℎ and respectively. Here, we used the boundary conditions

Density Profile
The expressions for the radiation pressure in the sinking region, (12) and (13), depend on the plasma density, . The mass continuity equation (1) can be used to obtain the density profile in the sinking region, once the velocity profile is known. In principle, we should solve the radiative hydrodynamical equations including the velocity terms to obtain a fully self-consistent velocity profile. Instead, following Mushtukov et al. (2015), we approximate the velocity profile by a power-law ∝ ℎ , taking a fiducial value of = 1. In §3, we will discuss the effects of varying the value of on the accretion column properties.
As previously mentioned in §2.1, we consider for simplicity the case in which the velocity profile is unchanged along , so that ( , ℎ) = (ℎ). As one of our boundary conditions, we assume that the plasma is in free fall with velocity ff above the shock, and is decelerated to ff ( )/7 below the shock (therefore losing all but ∼ 1/50 of its kinetic energy), where ff ( ) is the free-fall velocity at height . We refer to Becker (1998) for a one dimensional treatment of an adiabatic flow around the shock point that motivates the aforementioned velocity jump. The second boundary condition is given by the velocity vanishing at the surface of the NS. Hence, at the height ℎ above the NS surface, the velocity is given by By the mass continuity equation (1), we obtain the density profile where we introduced the accretion luminosity acc = / . Note that, as already pointed out by Mushtukov et al. (2015), formally the density diverges at ℎ = 0. Thus, the model assumptions become inadequate very close to the surface of the NS. In particular, at some point the gas pressure will start to dominate over the radiation pressure. In order to avoid this, and following again Mushtukov et al. (2015), we truncate the numerical calculation slightly above the surface, where rad ≈ gas .

Geometry of the accretion column
The expressions for the radiation pressure in the sinking region, (12) and (13), explicitly contain terms related to the geometry, and ℎ . In addition, the density profile equation (17) contains the sinking region area . These geometrical quantities depend on the accretion column base geometry, (ℎ = 0) = 0 0 , which in turn is determined by the specifics of the disk-magnetosphere interaction. The simplest and most commonly used model is the one proposed by Ghosh & Lamb (1978), according to which the disc is not sharply truncated at the magnetospheric radius. The result is a boundary region with a finite width, although much smaller than the magnetospheric radius (see figure 2). The two crucial quantities that appear in this simple disk-magnetosphere interaction model are the magnetospheric radius, , and the penetration depth, , which determines the boundary region width.
The magnetospheric radius is given by ≈ 7 × 10 7 Λ 1/7 10/7 6 4/7 d,12 where Λ is a dimensionless parameter depending on the accretion mode, 6 is the radius of the NS in units of 10 6 cm, d,12 is the surface strength of the dipole component of the magnetic field in units of 10 12 G, and 39 is the accretion luminosity in units of 10 39 erg s −1 . In this work we use the canonical value of Λ = 0.5 (Ghosh & Lamb 1978) for disk-fed accretion, and we investigate the effects of using different values of Λ on the luminosity in §3.3. Besides, for the magnetic field configurations considered in this work, the effects of including higher order multipole components on the magnetospheric radius can be neglected, as they decay much faster than the dipole component with increasing radius, and typically ∼ 100 . The penetration depth, , is expected to be of the order of the disk height for a geometrically thin disk (Ghosh & Lamb 1978). However, as it will be discussed later on, in certain models considered in this work, the thin disk approximation is not valid, and the parameters are more consistent with a geometrically thick disk. In these cases, by taking a prescription that relates to the disk height would result in a penetration depth greater than the magnetospheric radius, which is physically unreasonable. To account for this, we introduce an upper limit to the boundary region width with respect to the magnetospheric radius. We introduce the penetration parameter ≡ / , and we assume ≤ max , where max is a maximum penetration parameter. In the calculations, we use the prescription that the penetration depth is equal to the disk height, following Mushtukov et al. (2015), and separate from this assumption only when would be larger than the preset maximum (see §3.3 for more details). For our models, we use a fiducial value of max = 0.2 (see Li & Wang 1999 for a study of the boundary region width). However, in order to account for the uncertainty in the specifics of the disk-magnetosphere interaction, we also made some calculations by leaving max as a free parameter, and in §3.3 we discuss the sensitivity of the accretion column base geometry to changes in max .
Once the penetration depth is set, the shape of the magnetic field lines constrains the accretion flow, and fully determines the accretion column base width, 0 , and length, 0 . Since, for general magnetic field topologies, the equation of the magnetic field lines cannot be expressed analytically, we compute their shape numerically. This is done by finding the vector expressions for the multipolar magnetic field in polar coordinates and proceeding to integrate in the direction parallel to the magnetic field lines until reaching the surface of the NS. A graphic of the result of one such calculation, involving a multipolar magnetic field, is shown in figure  2.
The footprint of the magnetic field lines which also pass through the disk-magnetospheric boundary region forms an annulus centred on the magnetic axis with width 0 . The accretion column base length, 0 , is then given by the mean of the inner circumference and the outer circumference.
In addition to the column base geometry, the sinking region geometry depends on the shape of the magnetic field lines close to the surface (0 ≤ ℎ ). This must be taken into account when , there is no significant departure from the shape of a pure dipole magnetic field. The outer field lines (blue) are drawn in correspondence to the magnetospheric radius, while the inner ones (purple) to the inner radius of the boundary region. The two red segments show the part of the disk which enters the magnetosphere. calculating the sinking region width above the surface, ℎ , as well as the accretion column area above the surface, . In the case of a pure dipole magnetic field, it is where and are the accretion column width and length, respectively, at a height ℎ above the surface. As stated above, in the case of a multipolar magnetic field, a numerical computation is required. However, we find that the accretion column dimensions can still be written in the form with the area thus being given by In figure 3, we show the variation of the accretion column dimensions, for a selection of multipolar magnetic field configurations. The values of and in the multipolar magnetic field case can differ considerably with respect to those of the pure dipolar case.

Scattering opacity
The expressions for the radiation pressure in the sinking region, (12) and (13), contain the Rosseland mean opacity terms and ⊥ , which are the angle and frequency averaged opacities in the directions parallel and perpendicular to the magnetic field lines, respectively.
In a strongly magnetised plasma, and assuming large Faraday depolarization, radiation propagates in two normal modes, the ordinary (O) and the extraordinary (X) mode, with different polarization and opacity properties (see e.g. Meszaros 1992;Harding & Lai 2006). In this paper, we consider a pure scattering medium and calculate the electron scattering opacities in the (magnetic) Thomson limit, neglecting both the ion and vacuum contributions. The accreting plasma in the sinking region is assumed to be cold ( B 2 and | − cyc | 2 B / 2 1/2 | cos( )|, where is the mass and is the temperature of the electrons). We use the expression for the frequency and angle dependent electron scattering opacity of the two modes as discussed in Kaminker et al. (1982) (see also Zane et al. 2000), and we assume a fully ionized solar mix plasma with mean molecular weight = 1.17. The Rosseland mean opacity parallel to the magnetic field lines is given by and the Rosseland mean opacity perpendicular to the magnetic field lines is given by Here ( , ) is the electron scattering opacity integrated over all possible outgoing photon directions (see the appendix of Zane et al. 2000 for further details), the index denotes the polarization mode, where = 1 is the X-mode and = 2 is the O-mode, ( ) is the Planck function, is the photon energy, is the cosine of the angle between the photon propagation direction and the magnetic field lines, and ⊥ is the cosine of the angle between the photon propagation direction and the direction perpendicular to the magnetic field lines. ⊥ is related to by = √︃ 1 − 2 ⊥ cos , where is the azimuthal angle relative to .
Following Mushtukov et al. (2015), we estimate the effective opacity for mixed polarization modes by where is the fraction of radiation in the X-mode. To make comparisons with the purely dipolar model presented by Mushtukov et al. (2015), we adopt the same approach of considering an accretion column with X-mode photons only ( = 1). We discuss the variation of accretion column properties with a change in the X-mode fraction in §3.2. Note that and ⊥ depend on the temperature, , of the plasma. Since the sinking region is optically thick, we can approximate the plasma to be in thermal equilibrium with the radiation field. We relate the radiation pressure rad to the temperature locally using the Eddington approximation and the blackbody approximation, where is the radiation constant. Thus, the Rosseland mean opacities are calculated once the radiation pressure is known (see §2.6 for an outline of the computation scheme).

Scattering opacity for a hot plasma
The cold plasma approximation is valid only when the thermal motions of the electrons are negligible compared with the phase velocity of the wave (Harding & Lai 2006). While this is a sound assumption for many of the photon energies and plasma temperatures encountered in the models studied here, there are a few cases for which the cold plasma approximation no longer holds. To estimate the effects of a hot plasma on the scattering opacity, we average the cold plasma scattering opacity with the thermal motions of the electrons and introduce a line broadening of the cyclotron resonance. This treatment has the advantage of providing an approximation without resorting to a full computation of the magnetic Compton cross-section (which is beyond the purpose of this investigation). In our treatment of the thermal motions of the electrons, we consider the electron velocity distribution to be a one dimensional relativistic Maxwellian, given by where is the electron momentum along the magnetic field lines, and the proportionality constant is given by imposing the normalization condition, The averaged scattering cross-section is given by where rf is the energy of the photon in the rest frame of the electron and rf is the incident angle that the photon makes with the magnetic field lines in the rest frame of the electron. = / is the dimensionless electron velocity. For line broadening, we increase the resonance width by adding to the resonance damping term. The full calculation is beyond the scope of this paper (see Meszaros 1992). Following these estimations, the hot plasma scattering opacity is frequency-and angle-averaged to obtain the Rosseland mean opacities, as described in §2.4.

Model Estimates
An estimate of the radiation pressure, shape, and luminosity of an accretion column can be made for a sinking region with constant density profile (ℎ) = , and constant parallel and perpendicular opacity , ⊥ . In this case, rad, and rad,⊥ can be expressed analytically as where only the functions and ℎ are left to be determined. Note that is approximately the vertical optical depth of the sinking region, while ⊥ ℎ /2 is the horizontal optical depth of the sinking region at ℎ.
The expression for the normalised escaping flux can be obtained by equating the radiation pressure values obtained through equations (30) and (31) Following Mushtukov et al. (2015), we now use rad, ( , ℎ = 0) = rad,⊥ ( , ℎ = 0), and assume / 1, to obtain the relation Hence, in this simplified case the shape of the shock is approximately quadratic near the base of the column and becoming less so near the top. The luminosity of the column is obtained by integrating the escaping flux over the surface of the column. Doing so yields which is not integrable analytically due to the dependence of ℎ on ℎ. However, we can obtain a lower bound for the luminosity by setting ℎ ≈ , i.e. approximating the horizontal optical depth of the accretion column with its maximum value. Since the column is optically thick, i.e. ⊥ ℎ 1 and 1, then the luminosity is approximately given by as also obtained by Mushtukov et al. (2015), where However, note that the above equations are valid for any magnetic field configuration, provided = , which for instance turns out to be the case when the magnetic field multipoles are aligned (as it is assumed in this work). Equation (35) gives an approximate relation between the luminosity and the accretion column base geometry.
As previously discussed by Mushtukov et al. (2015), we can also see that ( / ) grows only logarithmically for large / , and this sets a natural scale for the maximum luminosity at / = 1, since the luminosity increases only marginally for higher .

Computational scheme
The procedure for computing rad from the radiation pressure equations, (12) and (13), is non-linear owing to the dependence of the opacity terms, and ⊥ , on the plasma temperature, which itself depends on rad through equation (25). For this reason we use an iterative method, and again we follow the scheme as given by Mushtukov et al. (2015), and we refer to this paper for all details.
We assume a magnetic topology of either purely dipolar or made up of a dipole plus an octupole component. The model parameters are the NS mass , radius , the accretion luminosity acc = / , the velocity power-law index , the polarization fraction , the penetration parameter upper bound max , and the strength of each of the two magnetic field components at the NS surface. The total surface magnetic field strength at the poles, , is the sum of the strength of each component at the poles, i.e. = dip + oct . Subsequently, the magnetospheric radius is calculated according to equation (18). The column base arc length 0 and the column base width 0 are calculated from acc and by following the magnetic field lines from the magnetospheric radius to the surface of the NS (described in §2.3).
For every set of parameters, the calculation is a double iterative process. First, we assume a trial value for the maximum shock height , and we calculate iteratively the radiation pressure profile and the shape of the sinking region. This inner loop is repeated until the calculated accretion column luminosity is within 1% from its value in the previous iteration. At convergence, the luminosity corresponding to this trial value of is calculated by integrating the escaping flux over the surface of the accretion column: A second iteration is then started, where is then adjusted until the luminosity of the accretion column matches the accretion luminosity to within 1%.

Effects of the magnetic field strength and topology
In this section, we investigate how changes in the magnetic field strength and topology affect the accretion column properties, due to the changes they produce in the radiative opacities. We consider a magnetic configuration made of two components: a dipole and a higher order multipole. The dipole component dominates the behaviour of the field at large distances, i.e. at the magnetospheric radius, while the higher order multipole regulates the behaviour of the field near the NS surface. In the following calculations, we take the multipole moment to be the octupole. We chose the octupole moment over the quadrupole moment to better localize the effects of the surface magnetic field and avoid potential problems with null points of the magnetic field above the magnetic poles. In principle, other multipoles can be used. We start calculating a series of models by varying the strength of the octopolar component. In order to separate the effects of the change in opacity from the changes in the column geometry (base and thickness) which are, in principle, also introduced by the magnetic field, we keep the accretion column base variables, 0 and 0 , fixed to the values that they assume in the pure dipolar case. For these models, we also neglect the curvature of the magnetic field lines. Numerical results are presented in black in figures 4 and 5 for a particular accretion luminosity and a velocity profile with power law index = 1. We present two set of models, for two different values of the accretion luminosity and magnetic dipole strength ( 39 = 1.0 and dip = 3 × 10 12 G, on the left, and 39 = 10 and dip = 5 × 10 13 G, on the right). In particular, in figure 4 we show the vertical cross-section of the sinking region (only half of the region is shown, due to symmetry), to investigate the changes in the shock shape. Figures 5 shows the profiles of the central internal temperature, of the effective temperature, and of the perpendicular mean opacity in the accretion column. These properties allow an easy comparison with the model presented by Mushtukov et al. (2015), in which the field was assumed to be a simple dipole.
For a more self-consistent treatment of the effects of a multipolar magnetic field, 0 and 0 must also be allowed to vary. As discussed in § 2.3, in a general case the base size and depth of the accretion column differs from those of a pure dipole magnetic field, so that 0 and 0 can both be reduced by a factor of several. In turn, this will affect the plasma density, the internal temperature, the escaping flux, and the maximum shock height. We therefore computed a set of models, by accounting for this effect and using the approach outlined in § 2.3 to calculate numerically 0 and 0 in correspondence of every assumed magnetic topology. For these models, the curvature of the magnetic field lines was taken into account. Numerical results are shown in red in figures 4 and 5, alongside the results for models with fixed column base geometry.
When using our numerical scheme we find that the maximum height of the shock, , is slightly smaller than in the models presented by Mushtukov et al. (2015), which may be due to a difference in the calculation of the opacity (see §2.4). However, the internal temperature profile, effective temperature profile, and opacity profile of our models (including those with an octupolar component) are qualitatively similar with those of Mushtukov et al. (2015). This indicates that the different field topology changes the quantitative details of the models but not the overall trend of the column properties.
The first thing to note from figure 4 is that, as also pointed by Mushtukov et al. (2015), the shape of the shock is not quadratic in , but instead the accretion column is quite narrow and the height of the shock above the surface drops practically to zero at a certain width˜ 0 /2 (we will refer to˜as the "sinking region width"). The sinking region width is the width at which the radiation pressure at the base becomes equal to the Eddington flux pressure, i.e. rad (˜, ℎ 0 ) = 2 3 Edd (ℎ 0 )/ . For >˜, the radiation pressure at the base is smaller than the Eddington flux pressure and no shock height can be supported.
The radiation pressure at ( , ℎ 0 ), i.e. at some distance along the column base, determines the height of the shock at , , by equation (12). Thus, the radiation pressure profile along the column base determines the shape of the shock. rad ( , ℎ 0 ) is determined by equation (13), in which we have assumed a linear function for the perpendicular flux, ⊥ , in our model. A different choice of function for ⊥ will give a different shock shape.
In figure 6 we present a series of models calculated varying the velocity index. As expected, a shallow velocity profile yields a shape more similar to that found for the analytical model, which predicts a quadratic column shape and is based on the simplified assumption of constant density. On the other hand, when a velocity index > 1 is used, the sinking region becomes narrower. This is a consequence of the greater deceleration of the particles in the lower layers of the accretion column. The upshot is that the radiation energy released by the particles is concentrated in the lower layers, which results in a lower shock height. Figure 4 shows that the maximum shock height, , decreases for an increasing surface magnetic field strength. By reversing the argument, for a fixed maximum shock height, a higher luminosity can be obtained by increasing the strength of the multipolar components. Increasing the magnetic field produces a larger opacity reduction in the X mode, therefore radiation escapes more readily from the sides of the sinking region. As a consequence, a smaller maximum shock height is sustained from the vertical radiation pressure.
The internal central radiation temperature profile also shows an anti-correlation with the magnetic field strength (see the top plots of figure 5). This is because models with a stronger magnetic field have a lower , for reasons outlined previously. In models with a stronger magnetic field, particles start to be significantly decelerated by the shock at a point nearer to the NS surface. Hence by continuity, the density in the sinking region is lower, and in turn the internal radiation temperature is lower.
The effective temperature, eff , obtained from the escaping flux using ⊥, = 4 eff , is also shown in figure 5. As already noticed by Mushtukov et al. (2015), this quantity does not have a profile that simply reflects that of radiation temperature. Instead eff increases with increasing height above the NS surface and then drops near the top of the accretion column. At the bottom of the accretion column, both the density, , and the geometrical thickness of the sinking region, ℎ , are large, which results in a large horizontal optical depth and a smaller escaping flux in that direction. Higher up, the accretion column becomes smaller in size, and the horizontal layers have a lower optical depth. The reduction in the optical depth is greater than that in the central temperature, and this is why the effective temperature generally increases with column height in the deeper regions. In fact, the peak of the effective temperature profile identifies the altitude at which the escaping flux is the greatest, and in turn this depends on the assumed density profile. For the accretion columns with velocity index = 1, the effective temperature peak is close to the maximum shock height, where the optical depth is lower. For < 1, the peak in eff is at a lower altitude than when = 1. The perpendicular Rosseland mean opacity,˜⊥ (shown in figure 5), is calculated at the central plane of the sinking region, = 0. This quantity depends on both the total magnetic field strength in the accretion column and the temperature. In general, a higher magnetic field strength or lower temperature reduce the perpendicular Rosseland mean opacity. However, the perpendicular mean opacity is not a monotonic increasing function of temperature or magnetic field strength. In fact, it is largest when the photons in the sinking region have energies close to the electron cyclotron resonance energy.
Comparing the models with fixed column base geometry (black lines in figures 4 and 5) and models including the curvature calculation (red lines in figures 4 and 5), it is immediate to note that there are no simple trends that explain the changes from the black curves to the red ones. This is because the change in accretion column base geometry (which decreases when the strength of the multipolar component is increased) results in a squeezing of the accretion column into a smaller area while the change in the curvature of the magnetic field lines results in an increase of the accretion column area moving higher up the accretion column. For each model, the overall outcome of these competing effects is different. However, it is worth noting that the change in the accretion column properties is modest and does not affect the qualitative behaviour discussed so far.
The effect of including the column geometry calculation is easier to explain in the models with multipolar magnetic fields shown in figures 4(b) and 5(b), which are low enough in height that the curvature does not make a substantial difference and the main effect is the reduction of the base size. For these models (in red), the internal temperature at a given height in the sinking region is increased compared with the fixed base geometry models (in black), since the same amount of energy is produced in a smaller area, and the perpendicular mean opacity rises following this increase in internal temperature. In addition, the density of the sinking region is increased, which makes it more difficult for radiation to escape from the sinking region, hence decreasing the effective temperature.
On the other hand, the effect of the curvature alone can be understood by comparing the purely dipole models with dip = 5 × 10 13 G (figures 4(b) and 5(b)). In this case, the accretion column base area is unchanged, since there are no multipolar magnetic fields, and the density is sufficiently high that the change in area due to the curvature is a consequential factor. For this model (in red), the density of the accreting plasma is lower near the top compared with the fixed base geometry model (in black), and this subsequently decreases the perpendicular mean opacity, allowing for more radiation to escape. Hence, a lower maximum shock height is sustained. At lower fields (as for dip = 3 × 10 12 G, figures 4(a) and 5(a)), the density inside the column is not high enough to make this effect substantial.

Mixed Polarization
So far, we have considered only models with X-mode fraction = 1. However, a more complete description of the radiation field in the sinking region includes a mix of X-mode and O-mode photons as well as scattering between polarization modes. In particular, one may be worried that when a substantial fraction of O-mode photons are present, the resulting decrease in the total opacity induced by the magnetic field is not sufficiently large to allow for super-Eddington emission. Solving this problem self-consistently requires a complete angle and frequency dependent solution of the radiative transfer problem, which is not the purpose of this paper. Instead, we investigate the issue by varying the value of in our calculation, i.e. we build up solutions by assuming that the radiation field consists of a given, fixed fraction of X and O-mode photons throughout the entire accretion column.
In general and as expected, we find that when a fixed fraction of O-mode photons is included, for an accretion column with a given magnetic field strength and accretion luminosity, the total opacity increases (as can be seen in the bottom plot of figure 8). However, even with a significant fraction of O-mode photons ( = 0.3), the average perpendicular mean opacity is still well below the Thomson scattering opacity in the cases considered here. Solutions with a luminosity well above the Eddington limit are still possible, but since the opacity is larger, a higher shock height is sustained. This results in a higher internal temperature and lower effective temperature. A set of fixed X-mode fractions, = 1, 0.7, 0.3, was used to illustrate the effect of changing the polarization degree. In figure 7, we show the vertical cross-section of the sinking region for an accretion column with magnetic field dip = 3×10 12 G, oct = 3×10 13 G and accretion luminosity = 10 39 erg s −1 . Figures 8 show the internal temperature, effective temperature, and perpendicular effective opacity profiles respectively for an accretion column computed by using the same parameters as in figure 7.
It is worth noticing that for low magnetic field strengths ( < 10 13 G), the opacity decreases when a fraction of O-mode photons is included. This occurs when a large portion of the photons have energy close to or higher than the electron cyclotron resonance energy, ≈ 11.6 12 keV. However, we are primarily interested in modelling sources with high magnetic field strengths ( > 10 13 G), in which case most photons in the accretion column will have energies below .
(a) (b) Figure 5. Plots of properties of the accretion column sinking region, with (a) and (b) corresponding to the models shown in Fig. 4. From the top to the bottom, the different panels show: the internal ( = 0) temperature profile of the accretion column, the effective temperature of the emitted radiation, and the perpendicular mean opacity, ⊥ (see text for details).

Disk-magnetosphere interaction
In our modelling, calculation of the magnetospheric radius from equation (18) requires an input parameter Λ, while calculation of the penetration depth requires taking a particular prescription for (see §2.3 for details). However, both the exact value of Λ and the expression for are poorly known. Thus, to test the robustness of the model results, we studied the response of the accretion column base geometry variables, 0 and 0 , to changes in Λ and respectively. For the models presented previously, we used the canonical disk accretion value of Λ = 0.5 (Ghosh & Lamb 1978). However, the exact value of Λ depends on the extent to which the NS magnetic field threads the accretion disk (Wang 1996), and for instance Dall'Osso et al. (2016) suggest Λ in the range 0.3 − 1 as a conservative estimate of the possible values. Repeating our calculations for 0 with various Λ, we find that 0 changes less than an order of magnitude when Λ is varied from 0.3 to 1. Since ∝ 0 and ∝ −1 , the overall accretion column properties are not very sensitive to changes in Λ.
With regards to , we assumed a penetration depth proportional to the disk height at the magnetospheric radius, as done by Mush-tukov et al. (2015). According to this prescription, the penetration parameter is given by However, since many of the models studied in this paper have a large accretion luminosity (with 39 ∼ 10) and low dipole magnetic field strength (with d,12 ∼ 1), the penetration parameter in equation (38) can be close to or in excess of = 1. Such values for correspond to the disk penetrating through the entire magnetosphere to the surface of the NS, which is a physically unlikely scenario. Thus, we introduced a maximum penetration parameter by hand (accordingly constraining 0 ), although the behaviour of remains unchanged from equation (38) unless the disk becomes geometrically thick. Alternatively, a self-consistent approach would be to introduce a new prescription for the penetration depth based on some set of physical principles, such as was done by Li & Wang (1999). However, this would require extending previous disk-magnetosphere interaction models to the case of a geometrically thick disk, which is not the purpose of this paper. Hence, as a substitute to considering many different prescriptions, we tested the response of 0 to changes (a)  in in general, without assuming a particular disk-magnetosphere interaction model. To do this, we repeated the calculation for 0 (see §2.3) while varying ∈ (0, 1). In each case, the accretion luminosity and magnetic field configuration are fixed. We considered = 10 39 erg s −1 and several magnetic field configurations with dip = 3 × 10 12 G. The results are reported in figure 9. Our calculations show that across the domain of , the change in 0 can be an order of magnitude or more. Since ⊥,esc ∝ −1 0 (see §2.5), the luminosity of the accretion column is sensitive to changes in . Hence, different prescriptions for the penetration depth can lead to dissimilar results for the accretion column properties, in particular the maximum luminosity.

Maximum Luminosity
One of the central aims of this work is to investigate the maximum possible luminosity from a highly magnetized, accreting NS, given some set of assumptions (such as the magnetic field configuration). In order to calculate the maximum luminosity that can be sustained by the NS accretion column, we compute the maximum acc for each set of model parameters by fixing the maximum shock height at = . In fact, at higher accretion column heights the luminosity only grows more slowly (see §2.5) and also the curvature of the magnetic field lines affects the vertical pressure balance equation (2), making our approximation unsuitable. We repeated the calculation of the maximum luminosity for several magnetic field configurations, namely a pure dipole field, a field with oct = 3 dip , and a field with oct = 10 dip , while the other model parameters have been fixed at = 1, = 1, max = 0.2. In agreement with the findings of Mushtukov et al. (2015), we find that, for total surface field strengths of < 10 13 G, the value of the maximum luminosity is mainly dictated by the accretion column geometry. This is because, at internal temperatures typical of the accretion column and for these low magnetic field strengths, most of the photons have > cycl and therefore are not subject to the reduction in opacity induced by the magnetic field. For higher total surface field strengths ( > 10 13 G), the scattering opacity of the X-mode is instead significantly reduced (by a factor of several to several orders of magnitude from T ) such that this becomes the determining factor in constraining the maximum luminosity. This can be seen by the change in slope of the maximum luminosity line in figure 10.
The decreasing trend for the maximum luminosity for magnetic field strengths up to 10 13 G can be explained by an increase in the temperature of the accretion column. Since the accretion column becomes thinner for higher magnetic field strengths (due to the choice of a maximum penetration paramater of = 0.1), the temperature increases, which also increases the overall scattering opacity.
As expected, when multipolar magnetic field configurations are Figure 11. Same plot as in Fig. 10 for models with different polarization fraction. The circles, triangles, and diamonds in green are computed using = 0.7, and assuming a pure dipole, a oct = 3 dip , and a oct = 10 dip surface magnetic field configuration, respectively. The solid, dashed, and dot-dashed black lines shows the same configurations but computed with = 1.
accounted for, we find that the maximum luminosity is increased when a stronger octupole component is present. However this is simply due to the fact that the magnetic field increases in strength: in fact, the maximum luminosity of a multipolar magnetic field corresponding to a total surface strength matches closely with the maximum luminosity for a purely dipolar magnetic field at the same . Thus, the change in the column geometry due to the presence of higher order multipoles does not affect the maximum luminosity in a significant way. Instead, as we will discuss in the next subsections, the maximum luminosity is more sensitive to the accretion column geometry and hence to the prescription for the penetration depth into the magnetosphere.

Maximum luminosity with mixed polarization
Since the introduction of a mixed polarization radiation field changes the accretion column properties (see §3.2), we also calculated the maximum luminosity for a fixed X-mode fraction of = 0.7 (as might be expected in a more realistic case, for a scattering dominated model). Results are presented in figure 11. As the maximum shock height is typically increased, the maximum luminosity is lower than in the pure X-mode case. However, this trend is reversed for magnetic field strengths below ∼ 10 13 G. This is due to a lowering of the average Rosseland mean opacity when a fraction of O-mode photons is included, which occurs when a large portion of the photons have energy close to or higher than the electron cyclotron resonance energy cycl ∼ 11.6 12 keV.
Otherwise, from figure 11, the relation between the surface dipole field strength and the maximum luminosity is seen to follow a shallower slope, as generally predicted. However, the deviation is not too significant, less than a factor of 2 in all cases we have examined. This is due to the fact that most of the flux that supports the accretion column, in the diffusion approximation, is due to Xmode photons. Thus, the effective mixed mode opacity, calculated in Rosseland approximation, is dominated by the X-mode opacity.  Fig. 10 for models with a different disk model. The orange circles show the computed maximum luminosity for a pure dipole surface magnetic field configuration, and with max ∼ 1. The black line is the maximum luminosity according to the relation given by Mushtukov et al. (2015).

Comparisons with previous models
Compared with the model of Mushtukov et al. (2015), we have used a different method of calculating the scattering opacity as well as a different disk-magnetosphere interaction model (as described in §3.3). In figure 12, we show the curve that represents the maximum luminosity obtained using our opacity files and assuming the same disk model as Mushtukov et al. (2015), namely by using max ∼ 1. In this case, the maximum luminosity differs only by a factor of a few with respect to the calculation presented by these authors, indicating a good agreement between the two codes.

Constraints on the parameter space
Before applying our accretion column model to observed PULXs, there are several considerations that we need to take into account. First, the region of the parameter space in which the model holds self-consistently is bound by the model assumptions, and in primis by the fact that we assumed a geometrically thin accretion disk at the magnetospheric boundary. This means where , the disk height at the magnetospheric radius, depends on the assumed disk model. By using the standard thin accretion disk model of Shakura & Sunyaev (1973), we find that, for the strong accretion luminosities in which we are interested ( > 10 39 erg s −1 ), the magnetospheric radius is always within the radiation pressure dominated zone of the accretion disk. In this zone, the disk height expression is given by which is independent of the radius. Using equations (40) and (18), the condition < is equivalent to a lower bound on the dipole field strength, above which our model is valid, which is given by For magnetic configurations with smaller dipole component, the thickness of the disk becomes large at the magnetospheric boundary, causing it to envelop the magnetosphere. In this case, our estimates of the accretion column geometry and our assumptions about the distribution of infalling plasma are no longer applicable. A proper analysis of this scenario requires a new disk-magnetosphere interaction model, which is beyond the purpose of this work. Second, since PULX are rotating NSs, the strength of the dipole component must also be sufficiently small so the propeller effect is avoided. This means the magnetospheric radius must be smaller than the Keplerian corotation radius, so The Keplerian corotation radius is given by where is the NS spin period. Equation (42) For magnetic configurations with a larger dipole component, the propeller effect prevents accretion onto the poles (Illarionov & Sunyaev 1975). Third, we will assume the spin period derivative to be dominated by the accretion torque. A simple accretion torque model is used to estimate the minimum average accretion rate that can give rise to the measured secular spin period derivative. In this model, we assume the angular momentum of the accreting matter is transferred to the NS at the corotation radius, which is the largest distance at which accretion can still occur. This produces the largest torque on the NS for a given amount of accreted material (and thus the lowest accretion luminosity). Thus, where ≈ 10 45 g cm 2 is the moment of inertia of the NS, and are the spin period and its derivative. Equation (45) can then be rearranged to give a lower bound on the accretion luminosity, acc = > 0.66 −10 −7/3 10 39 erg s −1 , where −10 = 10 −10 , and is in seconds. Values of the accretion luminosity lower then this limit would be insufficient in explaining the observed , according to our simple model.

Applications
Working with the parameter space restrictions derived in §3.5, we can diagnose the necessity of higher order multipole magnetic field components in observed astrophysical sources. We apply the model to two PULXs, for which the face value application of the model by Mushtukov et al. (2015) has led to the suggestion of the presence of multipolar magnetic fields, namely NGC 5907 ULX1 (Israel et al. 2017a) and NGC 7793 P13 (Israel et al. 2017b). The luminosity of both of these sources show a variation by a factor of ∼ 8, which is large but still more likely to be due to a variation in the accretion rate rather than a transition to the propeller effect. Hence, we apply the dipole magnetic field strength upper bound condition given by equation (44) for the entire luminosity range exhibited.
By discussing the source in the context of the Mushtukov et al. (2015) model, Israel et al. (2017a suggested the need of multipolar magnetic components. In order to test this argument, we plot again a figure analogous to Fig. 3 of Israel et al. (2017a) but using our computed maximum luminosity. The parameter space constraints in thedip plane, and few example configurations for NGC 5907 ULX-1 are shown in figure 13.
In agreement with previous findings, we find that in order to explain the whole range of observed luminosities up to peak ∼ 10 41 erg s −1 while assuming a pure dipole magnetic field, a strong magnetic field strength is necessary, i.e.
> 10 14 G (from our model) or even > 10 15 G (from the model of Mushtukov et al. 2015). Specifically, our model suggests that a pure dipole field of surface strength dip ≈ 3.1 × 10 14 G can give rise to a luminosity of 2.3 × 10 41 erg s −1 . Although the strength of this magnetic field is an order of magnitude lower than inferred from the Mushtukov et al. (2015) model, the source would still be in the propeller regime (see figure 13).
For the PULX to be emitting with peak without entering the propeller regime, a multipolar magnetic field is required, in particular a dipole component surface strength of dip ≈ 5.5 × 10 13 G and octupole component surface strength of oct 5.5 × 10 14 G (see figure 13). However, for all observed values of the luminosity, this configuration falls within the thick disk regime, i.e. the dipole magnetic field strength is lower than the bound given in equation (41). While it can not be excluded that this contradiction may be resolved once a different and less simplified disk model is adopted, to answer this question would require a thorough treatment of disk accretion and disk-magnetosphere interaction, which is beyond the scope of this paper.
As it can be seen in figure 13, both the super-Eddington disk accretion rate and the propeller regime can be avoided by introducing a moderate beaming factor of 0.15. In this case, a model based on a multipolar magnetic field configuration with dipole component surface strength of dip ≈ 2.8 × 10 13 G and a slightly larger octupole surface strength of oct 8.4 × 10 13 G can explain the entire observed range of luminosities, up to peak .
Stronger beaming factors 0.02 allow for a reduced accretion luminosity and therefore again make possible a pure dipole field configuration. However, in this case, the average accretion luminosity falls below the threshold required to give the secular spin period derivative ≈ 8 × 10 −10 s s −1 , i.e. the luminosity is lower than the bound given in equation (46). Hence, according to our model, the most favourable configuration for NGC 5907 ULX-1 includes a moderate beaming factor and crucially a multipolar magnetic field.

NGC 7793 P13
The PULX NGC 7793 P13 was observed to have a peak luminosity peak = 1.6 × 10 40 erg s −1 and luminosity variation between ∼ 2.0 × 10 39 erg s −1 and = 1.6 × 10 40 erg s −1 . A spin period of ∼ 0.42s was measured and a secular spin period derivative ∼ Figure 13. The parameter space plot of the magnetic field dipole component strength and the accretion luminosity for the source NGC 5907 ULX-1. The light red and darker red shaded area indicate the region for which exceeds the thick disk and NS Eddington luminosity at the magnetospheric boundary, respectively. The blue shaded area indicate region for which the source is in the propeller regime. The green shaded area shows the region for which the accretion rate is too low to provide sufficient secular spin period derivative = −8 × 10 −10 s s −1 , calculated based on equation (46). The solid, dashed, and dash-dotted black lines show the maximum luminosity in the case of a pure dipole field, a field with oct = 3 dip , and a field with oct = 10 dip respectively. The solid gray line shows the maximum luminosity according to the relation given by Mushtukov et al. (2015). Several example configurations for the source are shown by the black, red, and orange dots. In each case, the vertical lines represent the observed luminosity range when the source was in a high luminosity state (Israel et al. 2017b). In particular, the black dot shows the configuration required to remain under the maximum luminosity when using the relation given by Mushtukov et al. (2015). The red circle, diamond, and triangle show possible configurations (pure dipole, oct = 10 dip , oct = 3 dip respectively) with beaming factors of = 1.0, 1.0, 0.1475 respectively. The orange circle shows a configuration with high radiation collimation ( < 0.02).
−4.0 × 10 −11 s s −1 inferred from two observations one year apart (Israel et al. 2017b). Our constraints on the possible values of and dip are shown, together with some example configurations for this source, in figure 14.
In this case, the whole range of observed luminosities, up to the peak value of 1.6 × 10 40 erg s −1 , can be achieved by a configuration with a multipolar magnetic field of dipole component surface strength dip ≈ 7.3 × 10 12 G and a much stronger octupole component with surface strength oct > 7.3 × 10 13 G (see the red triangular point in figure 14). Under these conditions, the source is not in the propeller regime and no beaming is required to avoid the super-Eddington disk regime. This particular configuration has the advantage of being comfortably above the lower luminosity bound required to also explain the observed spin period derivative. However, the largest observed flux levels are not compatible with the assumption of geometrically thin disk, which again may demonstrate that disk model is over simplified.
When a mild beaming factor 0.25 is introduced, the effective luminosity peak becomes small enough that it can be explained by a configuration with a pure dipole field of surface strength dip ∼ 1.4 × 10 12 G. This is at variance with respect to the conclusions based on the calculation of Mushtukov et al. (2015). The only downside is that the lowest observed luminosities fall below the bound required to explain the observed spin period derivative. On the other hand, the secular spin period derivative may be the cumulative result of alternating accretion phases and may have been accumulated during epochs of larger mass transfer. Thus, this particular configuration is not severely disfavoured, according to the accretion torque model used here.

DISCUSSION
Motivated by the recent discovery of pulsating ULXs (Bachetti et al. 2014;Fürst et al. 2016;Israel et al. 2017b,a;Carpano et al. 2018;Rodríguez Castillo et al. 2019;Sathyaprakash et al. 2019) and their proposed interpretation in terms of accreting magnetars (see Tong & Wang 2019), we have reconsidered the problem of columnated accretion onto a highly magnetized NS. The main aim was to find model configurations capable of producing a high, super-Eddington luminosity while not in the propeller regime (for the values of the spin period typical of PULXs, ∼ 1s).
We worked in a scenario similar to the one recently discussed by Mushtukov et al. (2015) but we relaxed the assumption of a purely dipolar magnetic field. Instead, we considered combinations of dipolar and octupolar components. This combination was chosen due to the fall off of the octupole component strength with distance from the surface, thereby providing a magnetic field more concentrated close to the NS than a quadrupole component of similar strength.
We computed a series of models, characterized by either a sole low dipolar field (at 3 × 10 12 G) or a low dipolar field in addition to a stronger octupolar component (∼ 3, 10 larger). We first investigated the solutions by assuming that radiation is dominated by the lower opacity X-mode photons, and we indeed found a super-Eddington solution with ∼ 10 39 erg s −1 and ∼ 10 40 erg s −1 is always possible.
With respect to models based on a pure dipole field, we find that when an octopolar component is accounted for, the accretion column height is lower for a given luminosity. We find that in order to compensate the decrease in height, the effective temperature of the accretion column sinking region, eff , is larger in models with octupolar components. Typically, eff is in the range ∼ 3 − 15 keV, with the peak temperatures at the higher end of the sinking region. This thermal component is not observed directly as the radiation escaping the sinking region must pass through a region of freefalling material inside the accretion column. It may yet then be reprocessed in a thick accretion curtain (Mushtukov et al. 2017). This accretion column model relies on the Ghosh & Lamb (1978) disk model, specifically in order to derive the truncation radius of the disk as well as the penetration depth of the flow into the magnetosphere. One can question to what extent the standard disk model parameters of Ghosh & Lamb (1978) affect the accretion column properties as well as more generally the existence of a solution for a given set of model parameters. We investigated this aspect, and found that an accurate value of the truncation radius is almost irrelevant, while the assumption regarding the penetration depth is more crucial (see §3.3). In particular, a variation of 30% in the depth can decrease the base linear size of the column by an order of magnitude.
Another limitation concerns the way in which we treat the curvature of the field lines that constitute the boundary of the accretion column. A more detailed calculation, which is beyond the purpose of this particular paper, may account for the influence of the geometry changing with ℎ on the basic hydrodynamical equations (see Canalle et al. 2005).
The accretion column properties also depend on the assumed fraction of X-mode photons present in the radiation field. We investigated the impact of this, by introducing a fixed fraction of O-mode photons (see §3.2). Although we found variations in the X-mode fraction introduce differences in the accretion column properties for most of the magnetic field strengths considered here, this does not lead to a significant variation in the maximum luminosity curve.
In general, the opacity local to the NS surface is a decisive factor in allowing for the super-Eddington luminosities observed in ULXs ( > 10 39 erg s −1 ). However, the geometry of the accretion column footprint on the surface of the NS is also of crucial importance in determining the column properties. In particular, we found that the thickness of the accretion column has a significant effect on the luminosity (see 3.3). Indeed, we found a greater maximum luminosity for a given magnetic field strength compared with the one calculated by Mushtukov et al. (2015) in part due to our different approach to calculating the accretion column thickness.
Accretion columns with luminosity ∼ 10 41 erg s −1 are in principle obtainable with our modelling. However, in order to avoid the propeller regime for a source with pulse period of the order of ∼ 1s, the dipole field strength must be sufficiently low (∼ 10 13 G). A low strength dipole component ( dip ∼ 10 12 G) together with high accretion luminosity results in a thick accretion disk, as already noticed by Israel et al. (2017a). This is in contradiction with the assumption of a thin accretion disk of our model. Consequently, the problem of an upper limit to the luminosity related to the strength of the dipole component remains. We applied our model to the two sources NGC 5907 ULX-1 and NGC 7793 P13 (see §3.6). The necessity of a multipolar magnetic field configuration is different for each source, once beaming is taken into account.
For NGC 7793 P13, the observed luminosity is ≈ 1.6 × 10 40 erg s −1 . When taken face value, this luminosity level is too large to be compatible with the calculation of Mushtukov et al. (2015) since it would require a magnetic field so high that the source would be deep in the propeller regime. On the other hand, according to our model, the lowest observed flux levels are compatible with a purely dipolar configuration of strength dip ≈ 7.3 × 10 12 G, and the addition of a stronger octupole component with surface strength of oct > 7.3 × 10 13 G can explain the whole range of observed luminosities, up to the peak value of 1.6 × 10 40 erg s −1 . This latter particular configuration does not conflict with the propeller effect, nor with the super-Eddington disk accretion. Furthermore, it is compatible with the interpretation of the source spin period derivative according to a simple treatment (see §3.5). The only problem is that the largest observed luminosities are not compatible with the assumption of a geometrically thin disk, which in principle may be indicative of an overly simplistic disk model. Other possibilities include a moderate beaming, in which case the observed flux levels can be reached even for a purely dipolar magnetic field, with dip ∼ 1.4 × 10 12 G. However, this configuration requires that the observed spin period derivative is due to a secular torque, most of which is produced by the accumulation of material at the diskmagnetosphere interface during epochs of high flux level.
The PULX NGC 5907 ULX-1 has a much larger peak luminosity of 2.3 × 10 41 erg s −1 . In this case, both a super-Eddington disk accretion regime and the propeller regime can be avoided by invoking a moderate beaming factor of 0.15 ( figure 13). If the source has a dipolar magnetic field dip ≈ 3.2×10 13 G and a slightly larger octupole surface strength oct 9.6 × 10 13 G, then the entire observed range of luminosities can be reached. Even for this source, stronger beaming factors 0.02 allow for a reduced accretion luminosity and therefore make possible to explain the observed flux levels with a pure dipole field configuration, provided that, at the same time, it is assumed that most of the torque that gives rise to the spin period derivative is accumulated during phases of larger accretion rate.
There is still a number of open issues that needs to be addressed, before a self consistent explanation of PULXs can be reached. As already mentioned, the presence of multipole magnetic field components can change the properties of the accretion column significantly. The maximum shock height, , is reduced in comparison to a pure dipole magnetic field case. This in turn results in a higher effective temperature, which may manifest in the spectral data. In principle, a super strong magnetic field would be able to lower the maximum shock height very close to the surface so that . However, since our model assumes the radiation primarily escapes perpendicular to the sinking region, its validity would become more dubious as → 0. In addition, as the maximum shock height of the column is lowered, the temperature of the sinking region may exceed 100keV, whereupon we expect electron-positron pair creation and annihilation to play an increasingly important role in limiting the temperature of the accretion column while also increasing the gas pressure (see Mushtukov et al. 2019). A calculation including the gas pressure as well as pair creation and annihilation will be necessary for a more accurate description of the accretion column properties.
Several other simplifying assumptions were made in the model presented in this paper. First, we assumed that the radiation pressure dominates over the gas pressure in the sinking region of the accretion column. Through numerical calculation of several models, we found this assumption breaks down at the lower layers of the sinking region. Since we used a power-law ansatz for the velocity profile of the accreting plasma, the model is not expected to give an accurate picture of the lower layers of the sinking region, where the plasma flow becomes stagnant and hence the density becomes infi-nite. However, the contribution to the luminosity from these lower layers is negligible compared to higher up in the column, where the radiation pressure does indeed dominate over the gas pressure.
Second, in our calculation of the scattering opacity, we neglected the contribution from ions in the plasma and the contribution from vacuum polarization effects, which both become significant exactly in the strong magnetic field regime we consider here ( 10 13 G). Additionally, we assumed that a fixed fraction of Xmode photons made up the radiation field throughout the column and approximated the opacity as an effective scattering opacity (see §2.4). A more physically realistic treatment consists of a careful treatment of the scattering between the polarization modes as well as including mode switching due to resonant scattering. This will be the focus of our future work in development of the accretion column model. Finally, we did not take into account the role of energy advection by the accreting plasma and cooling via neutrino emission. These processes were studied by Mushtukov et al. (2018) and are expected to be relevant in the case of very luminous sources ∼ 10 41 erg s −1 . In §3.6, we opted instead to assume the was beamed by some mechanism. The exact details of admissible beaming factors for each of these sources is beyond the scope of this paper.

SUMMARY
We developed a simplified model of the accretion column for strongly magnetized NSs, building on and altering the model of Mushtukov et al. (2015). Crucially, we relaxed the assumption of a purely dipolar magnetic field, which we found allows for a larger maximum luminosity.
We found that when a magnetic field configuration with a significantly strong multipolar component is assumed, the luminosity released in the accretion column is limited only by the accretion rate from the disk. This, in turn, calls for more refined models of disk accretion and disk-magnetospheric interaction at the near-Eddington regime.
We applied the model to two PULXs, NGC 5907 ULX-1 and NGC 7793 P13, and discussed how their observed properties (luminosity and spin period derivative) can be explained in terms of different configurations, either with or without multipolar magnetic components. Generally speaking, the latter scenario is more favorable in case the emission is assumed to be highly beamed. Although at this level it may be difficult to differentiate further, we notice that strong multipole components may manifest in the spectra or polarization signal, an issue that we plan to investigate further in following work.

DATA AVAILABILITY
Data available on request.