Constraining the magnetic field geometry of the millisecond pulsar PSR~J0030+0451 from joint radio, thermal X-ray and $\gamma$-ray emission

With the advent of multi-wavelength electromagnetic observations of neutron stars, spanning many decades in photon energies, from radio wavelengths up to X-rays and $\gamma$-rays, it becomes possible to significantly constrain the geometry and the location of the associated emission regions. In this work, we use results from the modelling of thermal X-ray observations of PSR~J0030+0451 from the NICER mission and phase-aligned radio and $\gamma$-ray pulse profiles to constrain the geometry of an off-centred dipole able to reproduce the light-curves in these respective bands simultaneously. To this aim, we deduce a configuration with a simple dipole off-centred from the location of the centre of the thermal X-ray hot spots and show that the geometry is compatible with independent constraints from radio and $\gamma$-ray pulsations only, leading to a fixed magnetic obliquity of $\alpha \approx 75\deg$ and a line of sight inclination angle of $\zeta \approx 54\deg$. We demonstrate that an off-centred dipole cannot be rejected by accounting for the thermal X-ray pulse profiles. Moreover, the crescent shape of one spot is interpreted as the consequence of a small scale surface dipole on top of the large scale off-centred dipole.


Introduction
Millisecond pulsars (MSPs) are old neutron stars whose spin periods (∼ 1-30 ms) find their origin in the recycling scenario during which angular momentum is transferred to the neutron star via the accretion of matter from a companion star.They are also known to exhibit complex radio pulse profiles associated with strong variations in their polarization position angle (PPA) (Yan et al. 2011).This must be contrasted with young pulsars showing much more regular variation in the PPA, well described by the rotating vector model (RVM) (Radhakrishnan & Cooke 1969) because of the emission region occurring at higher altitude, in an almost dipole magnetic field geometry (Mitra 2017).The irregular behaviour of MSPs is probably related to the small size of the light cylinder and therefore to strong deviation from the dipole field and to the low altitude of emission where the multipolar components still give a contribution which is not negligible compared to the dipole component (Pétri 2019).
Nevertheless, multi-wavelength observations are able to put some constraints on their magnetic configuration.For instance, a joint modelling of the radio and γ-ray pulsed emission pins down the geometry of a centred dipole for young pulsars (Pétri & Mitra 2021), moreover supported by good quality radio polarization data when available.For MSPs, the situation is less constraining because the RVM cannot explain the PPA likely due to small scale magnetic field loops close to the surface, affecting the radio pulse profile and polarization angle.However, as shown by Benli et al. (2021), phase aligned γ-ray pulse profile fitting helps to extract the magnetic obliquity α of a pulsar as well as the line of sight inclination angle ζ.
As radio and γ-ray photons originate from well above the stellar surface, around several tens or hundreds of stellar radii, we do not expect these methods to allow us to disentangle the surface magnetic field structure.For example, for the pulsar studied in this paper, PSR J0030+0451, rotating at a period of P = 4.87 ms, the light-cylinder radius is about r L ≈ 232 km corresponding to more than 19 times the neutron star radius of R ≈ 12 km.To explore the surface magnetic field, thermal X-ray pulse profile modelling looks much more promising.In the last years, the launch of the Neutron Star Interior Composition Explorer (NICER) mission led to a breakthrough in the modelling of the thermal X-ray pulsations.One of the first targets of this mission was the MSP PSR J0030+0451, showing two prominent X-ray pulses separated by almost half a period.A careful and detailed analysis of this signal allowed to constrain the mass and radius of the neutron star (Riley et al. 2019;Miller et al. 2019) by fitting the thermal X-ray pulses with a model of stellar surface hot spots.This work also presented several possible shapes for the two hot spots responsible for the emission.One of the spots deviates significantly from a circular shape and resembles more to a crescent shape implying a strong non-dipolar surface magnetic field.It should be noted that other good fits were found independently by another group (Miller et al. 2019) implying possibly three hot spots.These spots are thought to be heated by particles accelerated to relativistic energies in the pulsar magnetosphere before impacting the region where the magnetic field lines meet the neutron star surface (Sznajder & Geppert 2020, see also Bauböck et al. 2019;Salmi et al. 2020) and should therefore be related to the magnetic field structure and plasma flow within the neutron star magnetosphere.Before the NICER era, Ruderman (1991) and more recently Bogdanov et al. (2008) already suggested the presence of an off-centred dipole to fit the X-ray pulse profile.However this seems insufficient to fully account for the recent NICER results because some quadrupolar components seemed to be also required.A comparison between the expectations from the NICER observations about the geometry of the magnetic moment and the inclination of the line of sight with older works like Johnson et al. (2014) is discussed in Bilous et al. (2019).In the outer gap model of Johnson et al. (2014) the line of sight inclination angle was found to be about ζ = 68 • ± 1 • whereas Riley et al. (2019) From a theoretical point of view, several groups attempted to connect these observations to their results about neutron star magnetosphere simulations.For instance, Kalapotharakos et al. (2021) used an off-centred dipole+quadrupole field structure to reproduce the polar cap shapes as well as the radio and γ-ray light curves.Following the same idea, Carrasco et al. (2023) performed general-relativistic force-free simulations of a centred dipole magnetic field only, introducing non standard emission regions based on the current density hitting the surface in order to model the thermal X-ray radiation for four NICER MSPs (PSR J0437−4715, PSR J1231−1411, PSR J2124−3358, and PSR J0030+0451).In this model, no off-centred dipole or multipoles are required.This contrasts also with Chen et al. (2020) who modelled a magnetic field geometry including an off-centred dipole+quadrupole to reproduce the polar cap shape found by the NICER collaboration.However in order to produce efficiently pair cascading close to the surface, it is known that small curvature radius are required, one to two orders of magnitude smaller than for the dipole.These smaller radii can be achieved by adding a small scale dipole anchored in the neutron star crust (Gil et al. 2002).Solving for the magnetic field topol-ogy is certainly one of the most difficult and central problems in neutron star physics remaining to be solved (Pétri 2019).
In this paper we show that a large scale off-centred magnetic dipole configuration is compatible with PSR J0030+0451 multiwavelength observations.We show that the joint radio and γray pulse profile modelling leads independently to the same conclusion as the thermal X-ray pulse profile fitting.The dominant magnetic field structure is consistent with an off-centred dipole and the peculiar hot spot crescent shape can be attributed to a small scale dipole localized in the vicinity of one pole.Section 2 summarizes the multi-wavelength data used in the present study and the analysis technique.Section 3 describes the method to deduce the magnetic field structure from the hot spot location.A discussion of the results is given in Section 4 before concluding in Section 5.

NRT
The Nançay Radio Telescope (NRT) is a meridian telescope equivalent to a 94-m parabolic dish located at the Nançay Radio Observatory near Orléans, France.Since the early 2000s, PSR J0030+0451 has been routinely observed with the NRT for timing purposes, with an average cadence of about 10 observations per month and with the bulk of the observations conducted at a central frequency near 1.4 GHz.
As mentioned in the subsequent Section 2.3, the γ-ray pulse profile for PSR J0030+0451 used in this work was taken from the Second Fermi Large Area Telescope Catalog of γ-ray Pulsars, hereafter 2PC (Abdo et al. 2013).For the γ-ray analysis of PSR J0030+0451, a pulsar timing solution for the pulsar constructed using NRT data was used.The timing solution was built by analysing data recorded with the Berkeley-Orléans-Nançay (BON) pulsar backend, between October 2004 and August 2011.The corresponding 1.4 GHz BON pulse profile shown in Fig. 1 and the timing solution can be retrieved from the 2PC archive of auxiliary files1 .
The NICER X-ray data analysed as part of this work were recorded more recently, as further described in Section 2.2.We therefore used NRT data to build another timing solution for PSR J0030+0451, contemporaneous with the NICER data.The radio data considered in this new analysis were recorded with the NUPPI backend (a version of the Green Bank Ultimate Pulsar Processing Instrument designed for the NRT), which replaced the BON backend as the principal pulsar timing instrument in operation at Nançay in August 2011.We used pulsar timing data recorded between August 2011 and April 2022.Although the two pulsar timing solutions described in this Section were obtained by analysing data from the same telescope, the two timing solutions were obtained from two separate analyses; therefore, the two ephemerides did not necessarily assume a common reference phase for the radio pulse profiles.We compared the BON and NUPPI reference profiles to determine the offset between the two reference phases.This phase offset of ∼ 0.056 was added to the NICER photon phases, in order for the radio, X-ray and γ-ray pulse profiles to share a common phase reference.

LOFAR
We also observed J0030+0451 using the international LOFAR station SE607 in Onsala, Sweden.LOFAR is fully described in Stappers et al. (2011) andvan Haarlem et al. (2013).LOFAR stations have two different frequency bands: We used the HBA band, for which the station consists of 96 antenna tiles.Each of these tiles is made up of 16 dual-polarization antenna elements.For each tile, an analog sum is performed.The tile signals are then digitized, and added numerically to form a coherent beam.For this project, we recorded data between 110.0 to 187.9 MHz.We re-folded the LOFAR observations using the NUPPI timing solution, except for the dispersion measure, for which we directly calculated one DM value per LOFAR observation.We combined 168 observations recorded between 2019 and 2023, correcting for time-dependent DM variations.

NenuFAR
We also observed J0030+0451 using the new low frequency radio telescope NenuFAR (New extension in Nançay upgrading LOFAR, see Zarka et al. 2015Zarka et al. , 2020) ) NenuFAR is a compact phased array and interferometer formed of hexagonal groups of 19 dual-polarization antennas called mini-arrays.It is located in Nançay, and operates in the 10-85 MHz frequency range .We used the dedicated pulsar instrumentation UnDysPuTeD and associated software LUPPI.LUPPI is based on NUPPI, and is designed to dedisperse and fold pulsar observations in real time, with up to four bands of 37.5 MHz simultaneously (Bondonneau et al. 2021).NenuFAR is an array of antennas, rolled out in several construction phases since the start of its observations in 2019.Observations in 2019 -2022 used 56 mini-arrays and observations in 2022 -2023 used 80 mini-arrays, out of which a minimum of 75 are active at any point, i.e. ≥1425 antennas in total, allowing for sensitive observations (Pétri et al. 2023).We recorded data between 20.3 and 84.3 MHz; in the final analysis, we only kept data above 42 MHz as scattering started to considerably affect the profile shape at the lowest frequencies.We re-folded the NenuFAR observations using the NUPPI timing solution, except for the DM value, for which the value determined directly from the NenuFAR observations is considerably more precise owing to the low-frequency coverage of NenuFAR.We combined 95 observations, correcting for time-dependent DM variations.We then added a constant phase offset to the Nen-uFAR observations to correct for different instrumental delays between NenuFAR and NUPPI.
J0030+0451 is among the twelve MSPs detected by Nenu-FAR, see Bondonneau et al. (2021) and Bondonneau et al. (in prep.).The NenuFAR profile of the pulsar differs from that obtained by NRT, but it shows a comparable width.Furthermore, our observations are compatible with the radio emission detected by NenuFAR and NRT being emitted at the same rotational phase of the pulsar.

Thermal X-ray pulse profile with NICER data
The thermal X-ray pulse profile was obtained from NICER observations of PSR J0030+0451.We used all available data from 2017-07-24 20:36:00 to 2022-02-12 15:10:00, corresponding to ObsIDs 1060020101 to 4060020619.To process the data, we employed the standard recipes with nicerl2 (from NICERDAS v9 distributed with HEASOFT 6.30 and NICER calibration file  ).Furthermore, we used the NICERsoft2 suite to perform additional filtering, with the following criteria beyond the default ones.We excluded the detector # 34, known to be particularly noisy; we excluded portions of the NICER orbit during which the Earth magnetic cut-off rigidity was lower than 1.5 GeV/c at the satellite's location.We also restricted the data to observing time when the space weather index KP was lower than 5, and when undershoot were lower than 200 c/s (see the NICER Data analysis threads3 for detailed description of these filtering criteria).
The NICERsoft package was also used to merge the event files into a single event file.Finally, we also performed a count rate cut on the merged event file, first when detector #14 (also known to be noisy on occasion) had a count rate above 1 c/s (with 8-sec time bins), and second when the total count rate (all remaining detectors) was larger than 6 c/s to remove generally noisy time intervals remaining after the filtering described above.The phases of all events were calculated with the photonphase task of the PINT package4 , and using the NUPPI ephemeris, including a phase offset of ∼ 0.056, as described in Section 2.1 Given the relative faintness of thermal X-ray emission of MSPs compared to the background level of NICER observations, we also calculated the optimal energy range that maximises the detection of the pulsations (see Guillot et al. (2019) for the description of this optimisation).This energy range will vary from pulsar to pulsar, and depends on the spectral shape of the pulsar, the energy range where the NICER effective area is the largest, and on the overall background spectral shape of the observations (which can vary depending on the strength of the various background components).For PSR J0030+0451, we find an optimal energy range of 0.28-1.46keV, giving a detection significance of 223 sigma.The pulse profile is shown in Fig. 1.

γ-ray pulse profile
Similar to the BON radio pulse profile for PSR J0030+0451 displayed in Fig. 1, the γ-ray pulse profile was taken from the 2PC auxiliary files archive.The profile, also displayed in Fig. 1, was constructed as part of the preparation of 2PC by analysing γ-ray photons recorded by the Fermi Large Area Telescope (LAT) be-tween 2008 August 4 and 2011 August 4. γ-ray photons with reconstructed energies from 0.1 to 100 GeV were selected and phase-folded with the BON ephemeris for PSR J0030+0451 as described in Section 2.1.

Pulsar magnetic field determination
Based on the radio and γ-ray data, the geometry of the dipole magnetic field of PSR J0030+0451 has already been investigated with the striped wind model in Pétri (2011) relying on the split monopole configuration introduced by Bogovalov (1999).More recently Benli et al. (2021) used the force-free magnetosphere solution for the special class of millisecond pulsars.In the latter work, the best fit parameters were found to be around α = 70 • and ζ = 60 • where α is the magnetic obliquity and ζ the line-ofsight inclination angle relative to the spin axis.In the next section, we re-explore these results, showing several possible configurations for the joint radio and γ-ray light curve fitting.Then we constrain the off-centred position of the dipole from the thermal X-ray hot spot location and shape.

Dipole geometry from joint radio and γ-rays modelling
The light-curve fitting procedure has been explained in depth in Benli et al. (2021) as well as in Pétri & Mitra (2021).We simply recall that the two main characteristics to be adjusted are the γray peak separation ∆ and the phase lag between the first γ-ray peak and the peak of the radio profile δ.To a very good approximation, with an error of less than 1%, Pétri (2011) showed that this relation is ( To find a good fit to the data, we minimize a kind of reduced χ 2 ν value defined by where N is the number of data points, ν = N − 2 the degree of freedom, I obs i the observed γ-ray flux and I mod i the predicted γ-ray flux.The number two parameters to adjust, the magnetic obliquity and the inclination of the line of sight.The most important characteristics of our model is the γ-ray peak separation ∆ and the γ-ray lag δ.Our model does not produce off-pulse emission.Therefore we select only the phase intervals containing the two peaks, as shown in the gray coloured boxes in Fig. 2  and agree reasonably well between radio, γ-rays and the geometry constrained from thermal X-rays.We found α = 75 • and ζ = 54 • for the upper panel and α = 70 • and ζ = 60 • for the lower panel.These results rely on a centred dipole, thus assuming that the radio and γ-ray emission emanate from high altitude, where the perturbation of surface fields and the off-centering of the large scale dipole become imperceptible due to the faster decrease of these components with distance.We next show how to determine the close region field structure, where the off-centred position of the large scale dipole is crucial for the thermal surface emission.

Off-centred dipole geometry from thermal X-rays
We constrain the geometry of the off-centred large scale dipole from both the location and size of the two polar caps that are deduced from the NICER observations of Riley et al. (2019) and Miller et al. (2019) and the half-opening angles of the cones subtending their rims.Fig. 3 highlights the geometry and shows the relevant angles and axes.Let us assume that the centre of the north polar cap is located at the point N with spherical coordinates (R, θ n , ϕ n ) and the position vector given by R n = R n and similarly for the south polar cap at point S with coordinates (R, θ s , ϕ s ) and a position vector R s = R s where R is the neutron star radius.n and s are therefore unit vectors pointing towards the centre of the hot spots.Let us also denote the unit vector joining both polar cap centres by The obliquity of the magnetic dipole is therefore given by projection onto the rotation axis e z according to We choose the solution corresponding to α ∈ [0, π/2].The minimal distance of the line joining N and S to the centre of the neutron star is given by This point on the segment NS is denoted by M. Equation ( 5) gives a lower bound to the distance d from the magnetic dipole moment location to the star centre.Would the dipole be located at any other point along the segment joining the point N to the point S , the distance d between the dipole position M and the centre of the star O would be larger.The angle between the line joining the dipole moment to the star centre and the rotation axis is where The angle β is the angle between the vector −−→ OM projected onto the plane (xOy) written as − − → OP and the vector joining S and N thus n − s or the unit vector t The relative angular size difference between the north and the south hot spot region helps to constrain the location of the magnetic dipole along the line joining the two poles N and S .Indeed, knowing the north and south polar cap radius respectively denoted as ξ n and ξ s , their ratio ξ n /ξ s is related to the distance d between the location of the dipole moment and the centre of the star.For instance, for an aligned rotator, if the dipole is located at the stellar centre, this ratio is equal to unity because both spots are identical.Moreover, using a spherical coordinate system (r, θ, ϕ), the field lines are given by r = λ r L sin 2 θ with r L = c/Ω the light-cylinder radius, c the speed of light, Ω the stellar spin frequency and λ a parameter labelling each field line.However, already when the dipole is shifted along the rotation axis, the geometry becomes asymmetrical and one spot grows whereas the second shrinks.More quantitatively, the halfopening angle ξ (north or south spot) is solution to the transcendental equation with a = R/r L and λ = ϵ a sin θ 0 cos(ϕ − ϕ 0 ) ± 1 − ϵ 2 a 2 sin 2 θ 0 sin 2 (ϕ − ϕ 0 ), ( 9) θ 0 and ϕ 0 the spherical polar colatitude and azimuth of the magnetic moment position.For an arbitrary location of this magnetic moment, even for an aligned rotator, the polar cap shapes deviate from a circular rim and the axial symmetry is broken.Therefore we need to solve for the radius at each azimuth ϕ to deduce the opening angle ξ for given ϵ and a.We need to find the roots ξ(ϕ) of Eq.( 8).The sign of λ is arbitrary because the solution always shows a symmetry between the north pole and the south pole.For a centred dipole ϵ = 0 and λ = 1, thus sin ξ = √ a = √ R/r L as is well known.A similar computation could be done with an orthogonal static rotator.The polar cap shape is then given by Pétri (2018) and depends on the angle ϕ.We do not enter into such complications because accurate determination of the rims would require a full force-free numerical simulation of the magnetosphere which is too computationally expensive in regard to the number of free parameters to explore.However, the aligned case already gives a good insight into the impact of off-centring onto the polar cap respective sizes.
In the particular case of the millisecond pulsar PSR J0030+0451, the ratio of stellar radius to light-cylinder radius is a = R/r L ≈ 0.055 according to Riley et al. (2019) who found an estimate of its radius R = 12.71 +1.14 −1.19 km and its mass M = 1.34 +0.15  −0.16 M ⊙ .The offset value ϵ = d/R required to adjust the relative spot angular sizes for different values of α and ϕ according to Eq.( 8) is shown in Fig. 4. θ 0 and ϕ 0 in the legend are given in units of π/4 for several different values spanning the range [0, π] in steps of k π/4 with k ∈ [0..4].Because some curves overlap, they are not all shown.A ratio ξ s /ξ n larger that 4 requires an offset as high as ϵ ≳ 0.9.A ratio ξ s /ξ n equal to unity means that the two spots are of equal size.In Fig. 4 we show the ratio ξ s /ξ n in log scale.We notice that this plot is symmetric with respect to the x-axis.Translated into the ratio ξ s /ξ n , it means that for each value of displacement ϵ there exist a value r 1 of the ratio ξ s /ξ n = r 1 with a position of the magnetic moment at (θ 1 , ϕ 1 ) associated to the inverse ratio ξ s /ξ n = 1/r 1 with a new position of the magnetic moment given by the angles (θ 2 , ϕ 2 ) uniquely related to (θ 1 , ϕ 1 ).Such symmetry is expected because of the symmetric nature of the magnetic dipole with respect to the inversion between north pole and south pole.For a more realistic force-free field these values could differ but we expect also large offsets for large differences in the spots sizes.
The fit parameters of Riley et al. (2019) are summarized in Table 1 that shows some of the models they used.Among others, these are the single temperature+protruding single temperature (ST+PST), the single temperature+eccentric single temperature (ST+EST), the single temperature+concentric single temperature (ST+CST), the single temperature unshared parameters (ST-U), the concentric dual-temperature unshared parameters (CDT-U).The geometry of the off-centred dipole with (α, β, δ, ϵ) is then deduced and shown on the right columns in the same  Independent constraints from a second group (Miller et al. 2019) are reported in Table 2 for a 2 spot and a 3 spot model5 .The line "3 spot1" considers only the first and second spot and the line "3 spot2" the first and third spot to find the geometry of the off-centred dipole.The line "3 mean" interprets the last 2 hot spots as emanating from a same polar cap where the centre is approximately located in the middle of the 2 spot centres (Θ p + Θ s )/2 and (ϕ p + ϕ s )/2 (which is an approximation of the exact middle joining both centres on the sphere).With this assumption we obtain other parameters that better agree with the 2 spot model referred as ST+PST in Table 1.

X-ray and γ-ray time lag
Our light-curve fitting procedure relies heavily on the phase aligned pulse profiles.The most important characteristics of these pulse profiles are the phase of the peaks in γ-ray and Xray.For the radio pulse profile, we prefer to use the centre of the radio pulse as the reference phase.We think indeed that the centre is more robust because it relies on a geometrical estimate not related to the flux of the radio emission that can strongly vary between components.We define the centre of the pulse to be at a phase half way between the phases of 10% of the maximum flux of the corresponding pulse.The three light-curves for radio, thermal X-ray and γ-ray are shown in Fig. 1 with vertical dashed lines depicting the phase of the mid-point at 10% maximum values for each pulse, see the values in Table 3.
Fig. 5. Geometry explaining the time lag between the radio peaks and the X-ray peaks and depicted by the angles ϕ 1 and ϕ 2 .
The first X-ray pulse peak arrives slightly before the pulse in the radio pulse profile whereas the second X-ray peak slightly lags the second radio pulse by a phase difference of about ∆ϕ ≈ ±0.025.This ordering is compatible with the off-centred dipole geometry because the order of appearance of radio and X-rays is inverted between north and south pole for simple circular shape polar caps.Even for PSR J0030+0451, for which one pole differs significantly from a circular shape, according to Riley et al. (2019) and Miller et al. (2019), the inversion in the phase lag between both wavelengths holds.
The phase lag between the radio peak and the X-ray peak is produced as follows: on one side, for a circular hot spot on the stellar surface, assuming constant and isotropic emissivity, the maximum X-ray flux is observed when the line of sight, the normal to the centre of the hot spot and the rotation axis are all coplanar.On the other side, the radio flux peaks whenever the tangent to the magnetic field line at the emission point lies in the plane defined by the rotation axis and the line of sight.Because the field lines are usually curved and not pointing in the radial direction, we expect a shift in phase between X-rays and radio.What matter is the projection of the magnetic moment and the position vectors of the hot spot centres onto the equatorial plane, the points A and B as shown in Fig. 5.The red straight line corresponds to the projection of the magnetic moment.The two green lines depict the projection of the vector position s and n.The star is rotating counter-clockwise.The first radio peak r 1 leads the first X-ray peak X 1 by a phase shift ϕ 1 .However the second radio peak r 2 trails the second X-ray peak X 2 , separated by a phase shift ϕ 2 .The first and leading radio peak has switched to a second and trailing radio peak.This configuration always holds except when the projected magnetic axis passes by the origin O.In such a case the time lag vanishes and the pulse peaks are aligned.A last important point concerns the phase lag between radio and X-rays which is not taken into account by the thermal X-ray fitting only.The peaks in X-ray are almost aligned Table 3. Phase location of the first and second peak in radio, X-ray and γ-ray.The phase difference ∆ϕ between both peaks is also shown (restricted to the interval [0,0.5]).The w 10 line indicates the phase interval where the radio pulses are detected above 10% of maximal flux.Errors in phase are of the order of the size of several bins thus ≈ ±0.01 in Xrays and γ-rays and ±0.001 in radio.
Fig. 6.Time lag as expected from the NICER observations.The dashed colour lines correspond to the direction of the magnetic axis whereas the solid colour lines to the direction of the normal for each hot spot, all directions being projected onto the equatorial plane.
with the centre of the radio pulse (Table 3) with a time lag of only 0.02 − 0.03 in phase.Fig. 6 shows the different locations of the hot spots projected onto the equatorial plane as deduced from the Table 1.They all give very similar directions of radio and X-ray emission phase of ϕ 1 ≈ 6 • and ϕ 2 ≈ −10 • .Only the ST+PST model produces larger time lags of ϕ 1 ≈ 8 • and ϕ 2 ≈ −26 • .The precise values are summarized in Table 4.Note the change in sign between ϕ 1 and ϕ 2 because of the switch from lagging to leading peak.The phase lags in units of the rotation period are of the order ϕ 2 ≈ −0.03 and ϕ 1 ≈ 0.02 in good agreement with the time lag in Fig. 1.
Accurately estimating the expected phase lag requires a more quantitative geometric study of the location of the peak X-ray emission within the hot spot as well as a careful analysis of aberration, retardation and magnetic sweep-back effects in the different wavelengths (Phillips 1992).Close to the neutron star, especially for the thermal X-ray emission, light-bending and the  4. Time lag ϕ 1 and ϕ 2 in degree and in units of the period, as defined in Fig. 5 and corresponding to the plot in Fig. 6.
Shapiro delay should be included too.This is however out of the scope of this work.
Our multi-wavelength analysis as well as NICER thermal Xray pulse profile fitting independently lead to a line of sight inclination angle of ζ = 54 • .As a consequence we take this value as a robust estimate of ζ.Therefore from relation 1 we found α ≈ 78 • in accordance with our γ-ray fit of 75 • .We then take α ≈ 75 • as another robust result.

Small scale surface dipole
The two independent groups working on the thermal X-ray pulse profile modelling fitted the hot spot emission with different patches, circular shapes for Riley et al. (2019) or ovals for Miller et al. (2019) who then also find possibly three hot spots.The rim of these polar caps must somehow be related to the surrounding magnetosphere, connected to the relativistic plasma flow hitting and heating the surface.It is well known that radio emission physics requires strong curvature field lines to produce sufficient electron-positron pairs and the associated radiation.This can only be achieved by small scale surface magnetic fields, dominating locally the global dipole field, possibly off-centred as discussed for instance by Szary (2013).In this section, we show that a small amplitude surface magnetic dipole located slightly inside the neutron star crust (at a distance 0.05 R below the surface) on top of a large scale off-centred dipole field can lead to polar cap shapes similar to those expected from the NICER observations.
There are several free parameters to specify each dipole configuration making it impossible to explore the full space parameter with high resolution and realistic magnetospheres filled with a pair plasma.Instead, here we only give a "proof of concept" of this idea by using a static vacuum dipole geometry, neglecting the displacement current.The components of the total magnetic field (off-centred dipole and small scale surface dipole) being known, we compute the field line structure by straightforward integration for a sample of representative lines.Each field line crossing the light-cylinder is reputed to be part of the open region where the plasma flows.Their foot points are connected to the hot spots.Fig. 7 shows an example of a two spot configuration, both poles being located in the southern hemisphere, one with a small, almost circular shape and the other much more elongated and thin, a crescent shape as reported in Riley et al. (2019).To obtain these shapes, the parameters used in the figure are as follows: a large scale magnetic moment inclination angle of α 0 = 75 • and a position vector (x 0 , y 0 , z 0 ) = (0, 0, −0.8) R, a ratio between the small scale dipole B 1 and the large scale dipole B 0 given by B 1 /B 0 = 0.3, a position vector of the small dipole (x 1 , y 1 , z 1 ) = (−0.4,−0.6, −0.6) R and an inclination angle of the small dipole amounting to α 1 = 150 • .Note however that the quantitative pictures do not agree as we are not attempting to model accurately the hot spots because of our very restrictive assumptions of an empty and static magnetosphere.

Discussion
On one side, joint radio and γ-ray fitting constrains the two important angles, namely the line of sight and magnetic dipole moment inclination.On the other side, thermal X-ray pulse profiles from the surface determine the location of the hot spots and their shape.Combining both approaches in a unified framework would better constrain the priors used in the procedure of thermal X-ray fitting.Indeed, from a theoretical point of view, the heating of the stellar surface seen as hot spots must be connected to the surrounding magnetospheric activity related to pair cascades and particle bombardments.From an observational point of view, the knowledge gained from the γ-ray fitting would pin down the line of sight inclination angle and the large scale dipole obliquity, narrowing the allowed parameter space to search in the thermal X-ray pulse profile fitting.
However, a severe limitation to this approach is the number of free parameters to be added to describe the magnetic field topology.We can model the field with for instance two dipoles or a dipole plus some quadrupolar components but the number of free parameters to adjust would be similar.The thermal X-ray hot spot geometry supports the idea of an off-centred dipole locally disturbed by a small scale dipole on the surface.However, we are far from a precise quantitative picture because this would require simulating a large sample of force-free magnetospheres with dipole+dipole fields to span a reasonable range of the space parameter.In the current stage, this is computationally not feasible except if the region to search for is already small enough with only a dozen of runs to perform.
Oblateness is also neglected in our study, we assumed a perfect spherical body.According to Silva et al. (2021), the perturbation induced by the rotation of a 200 Hz neutron star like PSR J0030+0451 is weak, the ratio between the equatorial and polar radius is expected to be between 1 and 1.04 so not too large for our qualitative study performed in this work.See for instance the vacuum and force-free oblate and prolate solutions found by Pétri (2022a,b).
Due to the presence of an interpulse component in its pulse profile, PSR J0030+0451 is assumed to be an approximately orthogonal rotator with a line of sight close to the equatorial plane, i.e. ζ ∼ 90 • .But the parameters found by NICER collaboration and by our radio/γ-ray fit shows that the line of sight is closer to ζ = 54 • .This means that the radio beam cone opening angle must be large.Assuming a dipolar magnetic field structure at the radio emission site, we can estimate the emission height.Indeed the opening angle θ em of the emission cone along the magnetic field lines at the surface is (Gangadhara & Gupta 2001): The north pole is visible only which is just on the lower limit for ζ.For the south pole to become visible we moreover require that which is clearly not the case.In order to reconcile the geometry with observations, the photon emission site for radio has to be shifted to higher altitudes around r em ≈ 6 R. In this case, the new cone opening angle becomes θ * em = 3/2 arcsin( just on the edge to observe the south pole.This altitude would correspond to a significant fraction of the light-cylinder radius of about 6 R/r L ≈ 1/3.Nevertheless, the presence of a small scale surface dipole could alter the direction of radio emission of the crescent shaped pole if photons are produced relatively close to the surface, possibly shifting the emission height to lower altitudes and then seen at different phases.This could also change the radio peak separation in our model where we used a simple centred dipole in Fig. 2. Looking at the width of both radio pulses, they are about 110-120 • , thus a half-opening angle of the emission cone of about θ em ∼ 55 − 60 • .This agrees with the above estimate to see both pulses when ζ ≈ 54 • .Consequently, there is a convergence towards the fact that the radio emission emanates from high altitude, close to the light-cylinder in an almost dipolar region. Inspecting the NenuFAR radio pulse profile in Fig. 1, we observe that the emission in both radio peaks are almost phasedaligned with their NRT counterpart (see Table 3).Moreover, the second radio peak pulse width is comparable to the NRT width and located at the same phases.However, the first peak seems to depart slightly slightly from its NRT counterpart.Its width is larger and it leads the high frequency counterpart by 3% in phase.This suggests that both radio bands are produced at approximately the same altitude within the magnetosphere.The same conclusion holds for the LOFAR observations, the width of the pulses are similar to those seen by NenuFAR although their shape differ, especially for the weaker pulse around phase 0.5.Actually because of the small size of the light-cylinder, only r L ≈ 18.3 R, there is indeed not much space left to span a vast range in emission heights.Moreover the first peak morphology changes at low frequency, highlighting a possible change in the local magnetic field geometry at the altitude where these photons are produced and due to the small scale surface dipole.Indeed, Riley et al. (2019) points out that the first X-ray pulse profile belongs to the crescent shape hot spot.In our picture this is related to a small scale surface dipole that then also impacts the radio emission properties between the NenuFAR and NRT frequencies.

Conclusions
With the increase in sensitivity of telescopes at all wavelengths, a simultaneous fitting of radio, X-ray and γ-ray light-curves of pulsars becomes a very powerful tool to constrain the location of the emission regions as well as the magnetic dipole obliquity and observer line of sight inclination angles.For PSR J0030+0451, we showed that on one side, thermal X-ray pulse profile fitting and on the other side a joint radio and γ-ray light-curve fitting gave very similar expectations for those angles.This supports the idea that radio and γ-ray emission models for young pulsars apply equally well to recycled millisecond pulsars even if the presence of multipolar fields could strongly impact the surface emission and to a lesser extent the radio and γ-rays radiation.We found that α ≈ 75 • and ζ ≈ 54 • are very robust results for the multi-wavelength light-curves of PSR J0030+0451.
The above study checked the consistency between both approaches without digging into a precise description of the magnetic configuration.To go further would require a more quantitative analysis, performing numerical simulations of off-centred force-free dipolar magnetospheres on top of small scale surface dipoles.Unfortunately, the total number of parameters makes it intractable computationally because of the huge number of configurations to explore.It is therefore compulsory to narrow down the region of interest in the parameter space for each individual pulsar as done in the present work.
The concomitance between the two independent results in Xray and radio/γ-ray may just be a coincidence.But with the ongoing NICER campaign, the number of pulsars with confident hot spot modelling will increase and such chance occurrence could be rejected with ever increasing confidence.Therefore, our next target is PSR J0740+6620 (Miller et al. 2021;Riley et al. 2021) that shows very similar behaviour to PSR J0030+0451.Results will be discussed in an upcoming paper.
Last but not least, millisecond pulsars might emit thermal X-rays that are polarized.The recent discussion by Suleimanov et al. (2023), based on a self-consistent atmosphere model in local thermodynamics equilibrium, shows that for an unmagnetized neutron star with a hydrogen or a carbon atmosphere, the maximum polarization degree can reach 25% and 40% respectively.Moreover, it is possible for the X-ray emission to become substantially polarized if multipolar components are present, with magnetic field strength exceeding B ∼ 10 6 T (see for example spectra and polarization of magnetized neutron stars in Lloyd 2003).Including polarization degree and polarization angle in the analysis could provide two additional dimensions to help to break the degeneracy in the parameter space and disentangle the surface magnetic field structure.This kind of study would be particularly relevant in the context of future X-ray polarimetry missions with sensitivity around 0.1-1 keV 6 .

Fig. 1 .
Fig. 1.Time lags between peaks in the different wavelengths.Radio pulse profiles are in red for NRT, in dashed blue for LOFAR, in dotted magenta for NenuFAR, thermal X-ray are in green, and Fermi/LAT are in black.The dashed coloured vertical lines show the phase of the peak flux of each pulse. v20210707 1) For PSR J0030+0451, the second Fermi pulsar catalogue (Abdo et al. 2013) reports a value of the phase separation between the two γ-ray peaks of ∆ = 0.450 as well as a γ-ray lag of the first γ-ray peak with respect to the peak of the radio profile of δ = 0.146.Additionally the NICER collaboration reports a value of ζ = 54 • (Riley et al. 2019).This would imply an obliquity of α = 78 • not far from the independent γ-ray light curve fitting done by Benli et al. (2021) who found α = 70 • .Using the input from the thermal X-ray detailed in Section 3.2, we were also looking for solutions of the γ-ray light-curves by constraining the line of sight inclination to ζ = 54 • or to ζ = 50 • as given by the NICER collaboration.
. All light-curves are plotted in arbitrary units, applying a normalization such that the peak value in each energy band is unity.The best values for the parameters α, ζ and the shift in phase ϕ s are those minimizing the χ 2 ν .Two possible fits are shown in Fig.2

Fig. 2 .
Fig. 2. Some fits for the radio and γ-ray light curves of PSR J0030+0451.The parameters used are α = 75 • and ζ = 54 • on the top panel and α = 70 • and ζ = 60 • on the bottom panel.The gray coloured boxes show the phase intervals used for the γ-ray fit.

Fig. 3 .
Fig. 3. Geometry of the off-centred dipole showing the location of the hot spots N and S and the angles defined in the text (α, δ, β).

Fig. 4 .
Fig. 4. Ratio ξ s /ξ n of the polar cap half-opening angles of the south and north pole respectively ξ s and ξ n as a function of the fractional offset from the centre ϵ = d/R.The colours correspond to different values of the angles (θ 0 , ϕ 0 ) in units of π/4.

Fig. 7 .
Fig. 7.A 3D view of the two hot spots from the south pole.They are produced by a large scale off-centred dipole and a small scale surface dipole located close to one magnetic pole.

Table 1 .
table.Off-centred dipole geometry deduced from the polar cap location after Riley et al. (2019).