Radio Polarization of Millisecond Pulsars with Multipolar Magnetic Fields

NICER has observed a few millisecond pulsars where the geometry of the X-ray-emitting hotspots on the neutron star have been analyzed in order to constrain the mass and radius from X-ray light-curve modeling. One example, PSR J0030 + 0451, has been shown to possibly have significant multipolar magnetic fields at the stellar surface. Using force-free simulations of the magnetosphere structure, it has been shown that the radio, X-ray, and γ-ray light curves can be modeled simultaneously with an appropriate field configuration. An even more stringent test is to compare predictions of the force-free magnetosphere model with observations of radio polarization. This paper attempts to reproduce the radio polarization of PSR J0030 + 0451 using a force-free magnetospheric solution. As a result of our modeling, we can reproduce certain features of the polarization well.


INTRODUCTION
A millisecond pulsar (MSP) has an extremely small spin period (P ∼ 1 − 30 ms) and differs from ordinary pulsars in that it has an extremely small spindown rate and exists in a binary system (Lorimer 2008).Their short periods are thought to be caused by accretion of matter from a donor star (Alpar et al. 1982).A pulsar's magnetic pole emits cones of bright radio emission when it rotates rapidly, sweeping around like a lighthouse, and the interval between pulses is as precise as that of an atomic clock.Radiation emitted by pulsars is generally believed to be produced in pair plasma outflow streaming along the dipolar magnetic field lines (Philippov et al. 2020;Melrose et al. 2021;Philippov & Kramer 2022).However, it is still theoretically unclear how pulsars produce such coherent sources of radio waves.Nevertheless, over the past decades, pulsars have proven to be fascinating laboratories for studying fundamental physics.For example: (1) measuring their masses and radii can constrain the equation of state for nuclear matter (Raaijmakers et al. 2019;Bogdanov et al. 2021), (2) measuring their motion can help test general relativity (GR) (Kramer et al. 2006;Kramer & Wex 2009), (3) Pulsar timing can be used to detect gravitational-wave background (Arzoumanian et al. 2018(Arzoumanian et al. , 2020)), and finally (4) they are good probes for studying the local environ-ments and properties of their host galaxies (Coles et al. 2015;Jones et al. 2017).
We can get constraints on the emission physics from the pulsar magnetosphere through multiwavelength analysis.For example, NICER has obtained detailed X-ray observations of hot spots present in the pulsar PSR J0030+0451 (Miller et al. 2019;Riley et al. 2019) which has allowed us to constrain its magnetic field geometry (Bilous et al. 2019).Their results showed that the magnetic field is far from a simple dipole but favors multipolar components at the stellar surface.The shape and location of the hot spots observed in thermal X-rays are the footprints of open magnetic field lines.Here, active pair-production occurs, which leads to particle bombardment of the stellar surface and production of X-rays.The thermal X-ray lightcurves of PSR J0030+0451, and, in addition, radio and γ-rays, produced in the outer magnetosphere, were successfully modeled by Chen et al. (2020); Kalapotharakos et al. (2021).This modeling confirmed that a multipolar magnetic field is key in understanding this system and can be further used to constrain the magnetic inclination angle.
The magnetic field in ordinary pulsars has been shown to evolve through the Hall effect and Ohmic dissipation in the crust mediated by free electrons (Cumming et al. 2004;Pons & Geppert 2007;Gourgouliatos & Cumming 2014;Bransgrove et al. 2018), and ambipolar diffusion driven by binary scattering processes in electron-protonneutron plasma (Goldreich & Reisenegger 1992;Castillo et al. 2017;Bransgrove et al. 2018).Modeling the interior and exterior field consistently is difficult, and in fact, most models for the interior field in the star do not model the magnetosphere and just have some fixed exterior boundary, and vice versa.The dynamical interplay between these regions is crucial to understand some magnetospheric phenomenology.Global magnetohydrodynamics simulations with initial dipolar magnetic field for ordinary pulsars have shown that after birth the system evolves to a state with significant power in the multipolar components (e.g., Sur et al. 2020).Higherorder multipoles near the surface of the star have been proposed to activate pair production that is supposed to power the radio emission (Gil et al. 2006).However, these multipoles would dissipate within a few million years due to Ohmic decay (Gourgouliatos & Cumming 2014) even if they were generated after an ordinary pulsar is recycled into a MSP.The formation process of MSPs through accretion may also change the magnetic field configuration of MSPs, either through burial of magnetic field (Romani 1990;Melatos & Phinney 2001;Payne & Melatos 2004), or due to field migration as the neutron star spins up (Ruderman 1991;Chen & Ruderman 1993;Chen et al. 1998).The magnetic field is therefore challenging to model in MSPs.The dipolar field model has been commonly used in various studies, such as determining the magnetic field strength from the dipole spindown and explaining the inverse relationship between spin period and pulse width.However, the presence of non-dipolar magnetic field configurations has an impact on various aspects of pulsar astrophysics, including birth velocities (Radhakrishnan 1984;Bailes 1989) and the interpretation of multi-wavelength magnetospheric emission.Therefore, it is crucial to have a precise representation of the magnetic field when studying MSP-related phenomena.
One of the most important tools to connect the pulsar magnetic field geometry to observations is the radio light curve and polarization modeling.In the past, the rotating vector model (RVM) (Radhakrishnan & Cooke 1969) was widely used in this respect.The model assumes that the emission comes from near the magnetic pole, and the polarization is determined by the direction of the magnetic field at the emission point.As the line of sight cuts across the magnetic pole, the plane of linear polarization sweeps a characteristic S-shaped swing over a single pulse period.Radio pulsar mean profiles are typically interpreted according to the hollow cone model (Manchester & Taylor 1977).The main assumptions in this model are the following: (1) the emission is generated in the inner magnetosphere (where the magnetic field (B) is thought of as a dipole), (2) the emission travels in a straight line, (3) cyclotron absorption may be ignored, and finally (4) the polarization is determined at the emission point.Analytically, the change of the position angle (PA) with respect to the rotation phase (ϕ) of the pulsar is given by where α is the angle between the magnetic moment axis and the rotation axis while δ is the angle between the rotation axis and the line of sight of the observer.According to the RVM, PA is determined solely by the projection of magnetic field on the sky plane (see figure 1).In this convention, the x-axis is along the projection of the angular velocity Ω on the sky plane.Another modification to the RVM was the model presented by Blaskiewicz et al. (1991) (BCW) which incorporates relativistic effects, such as aberration due to significant corotation component of the plasma velocity and retardation.The BCW model calculates the position angle by finding the direction of the acceleration of radiating charged particles at the emission point.The relativistic effect causes the polarization profile to lag the intensity profile.The magnitude of the lag of the inflection point of the PA profile in the presence of aberrations and retardations can be calculated as follows: where P is the period, r em is the emission radius, and c is the speed of light.Based on the observed lag, the emission radius can be inferred.
Although the RVM with a dipolar field successfully explains pulse profiles of many pulsars qualitatively (e.g., Manchester et al. 1975;Johnston et al. 2023) and even quantitatively (e.g., Desvignes et al. 2019), observations of MSPs often show flat, distorted, or even random PA profiles hinting towards possible non-dipolar configurations (Backer et al. 1976).In addition, propagation of radio waves through magnetospheric pair plasma can further change the polarization properties (Petrova & Lyubarskii 2000;Wang et al. 2010;Beskin & Philippov 2012;Hakobyan et al. 2017;Galishnikova et al. 2020) which are not included in the RVM model.It is important to note that some MSPs have flux densities similar to ordinary pulsars and their radio profiles are only marginally more complex (Kramer et al. 1998).Nonetheless, physical conditions in magnetospheres of MSPs and their surface magnetic field structure could differ considerably owing to its evolutionary history, which can result in changes to their observed appearance (Philippov & Kramer 2022).
In this paper, we are interested in comparing radio polarization with observations for the MSP PSR J0030+0451 taking into account a multipolar magnetic field configuration as obtained in Chen et al. (2020).Developing a better understanding of the multipolar structure and surface field will help to constrain not only radio emission sites, but also the evolution of the magnetic field.As a first step, we use a modified RVM: we assume that the radio emission is produced as plasma normal modes at a radius r em which then propagates along a straight line, where its PA evolves adiabatically following the magnetic field.The PA then freezes at a distance h from the emission point where the magnetospheric plasma density has dropped so that the radio waves follow vacuum propagation afterwards.We take into account the aberration effect self-consistently, but neglect other propagation effects for now.With this approach, we systematically obtain PA sky maps and curves with r em and h as free parameters, which we then compare with the observed PA signatures.The paper is organized as follows: in section 2 we discuss the magnetic field model, in section 3 we derive an expression for PA considering aberration, in section 4 we describe our method, in section 5 we present the results, and finally, in section 6 we discuss conclusions.

THE MAGNETOSPHERIC MODEL
At the polar caps of a plasma-filled magnetosphere, electric current flows on open field lines.It has been reported that PSR J0030+0451 has a spin period of 4.18 ms with hotspots that are beyond a simple dipole, requiring multipolar magnetic field components (Bilous et al. 2019).We use the field configuration first deduced by Chen et al. (2020), who showed that a quadrudipolar magnetic field can be used to reproduce the light curves in X-rays, gamma rays, and radio waves for PSR J0030+0451.The analytical form of this field in vacuum is given by where the dipolar component is with the dipole moment p = p 0 (0, 0.985, 0.174), while the quadrupolar component is with where R is the neutron star radius.The dipole inclination angle is 80 • and the quadrupolar component is centered at (0, 0, −0.4R).These exact set of parameters were used in Chen et al. (2020).A 3-dimensional view of the field is shown in figure 2. The field is symmetric about the y-z plane, and the polar caps are located in the southern hemisphere, one of them circular and the other crescent-shaped.
The structure of the force-free magnetic field is shown in figure 3. Beyond the light cylinder radius (R LC ), which is located at 20R for PSR J0030, the field lines are opened, and opposite polarities are separated by a current sheet (Spitkovsky 2006;Kalapotharakos & Contopoulos 2009;Chen & Beloborodov 2014;Philippov & Spitkovsky 2018;Hakobyan et al. 2019).

EFFECT OF ABERRATION
In addition to our force-free magnetic field solution, we assume that the outflowing plasma has γ = 1/ 1 − v 2 /c 2 ≫ 1 moving along the magnetic field lines.Let us also assume that the emission comes from  a sphere of constant radius (r em ).1 In the case of emission from one magnetic pole, there is a single point on this sphere where particles beam along the observer direction n at a particular time.In ordinary pulsars, the emission direction is solely dependent on the direction of the magnetic field.But since PSR J0030+0451 is an MSP with spin period of 4.18 ms, the aberration angle at the emission point is approximately Ωr em /c (∼ 0.3 at r em = 5R) which cannot be neglected compared to the angular size of the emission cone (1/γ)2 .As inferred from the aberration-retardation modeling, typical values of the emission radius for MSPs are less than 100 km, which translates to ∼ 0.05 − 0.5 R LC (Rankin et al. 2017).On the other hand, ordinary pulsars emit radiation at a height ∼ 300 − 1000 km (Gupta & Gangadhara 2003;Dyks et al. 2004;Johnston et al. 2023) where the magnetic field is dipolar.The aberration angle of an ordinary pulsar, for example, PSR J1808-0813 with period P = 0.876 s, is Ωr em /c ∼ 1×10 −7 , where r em is taken to be ∼ 693 km (Mitra et al. 2023).This is negligibly small compared to the aberration angle of PSR J0030+0451.In this work, we assume 0.1 R LC ≤ r em ≤ 0.5 R LC .
The location of the emission point, r em , is determined by demanding the direction of the particle propagation, along which the radiation is beamed, to be along the observed line of sight, n (Wang et al. 2010;Beskin & Philippov 2012).The speed of the outflowing plasma is close to the speed of light, c, resulting in the following where b = B/B is a unit vector in the direction of the magnetic field, and κ is a dimensionless constant.Its magnitude can be obtained by demanding that the flow is away from the star, and n • n = 1.

Normal modes
Normal modes in plasma physics refer to the oscillations and pattern of waves that a plasma can sustain due to interaction between the electromagnetic forces and charged particles.Assuming a cold pair plasma, the wave dispersion relation in the rest frame is given by (Philippov & Kramer 2022): where ω p is the plasma frequency, ω is the wave frequency, and k || and k ⊥ are the components of the wave vector k parallel and perpendicular to the background magnetic field, respectively.The solution to equation 9 gives three linearly polarized waves: the extraordinary mode (X-mode) having the electric field (E) perpendicular to both the background magnetic field and k-vector; and two ordinary modes (O-modes) with the E-field in the plane of k-vector and background magnetic field.In the pulsar magnetosphere, only radio waves corresponding to plasma normal modes can propagate and escape.

The new Position Angle
Given that we found the location of the emission point, we next derive the polarization as seen by an observer on Earth.Radio waves are emitted as plasma normal modes as discussed before.In the nearly force-free magnetosphere, the plasma is rotating together with the NS so the rotation will affect the plasma normal modes.Here we first obtain the normal modes in the corotating frame, then transform to the laboratory frame, because the normal modes are relatively easy to write down in the corotating frame.In the corotating frame, the background electric field is zero and the polarization of the wave modes is determined by the background magnetic field and the wave vector in the same way as in the plasma rest frame (see, e.g., Wang & Lai 2007).
Transforming back to the laboratory frame, we can obtain the polarization as seen by a distant observer.In what follows, we consider the limit where the refractive index n ≡ ck/ω ≈ 1.Here, the wave vector in the laboratory frame is k and the frequency is ω.The background magnetic field in the laboratory frame is B, and the magnetospheric electric field is E = −β × B, where β = Ω × r/c is the dimensionless rotation velocity.In the corotating frame (we use primed quantities for the corotating frame), the background magnetic field becomes where Γ = 1/ 1 − β 2 is the Lorentz factor corresponding to the rotation velocity.The wave frequency in the corotating frame is and the wave vector becomes Considering an O mode with electric field E ′ w in the plane of the wave vector k ′ and the background magnetic field B ′ , we can write the electric and magnetic field in the wave as the following: Now, transforming back to the laboratory frame, we obtain the wave electric field Plugging equations ( 13) and ( 14) in equation ( 15), and after some algebra3 , we get where we have assumed ω = ck, and ⊥ indicates the component perpendicular to the wave vector k.Therefore, if we choose a coordinate system such that k is Sur et al.
along ẑ, the wave polarization in the laboratory frame is given by Thus, the new position angle including the effect of aberration is given by In order to obtain the PA curves, we first construct sky maps of the polarization.Since the force-free solution reaches a steady state in the corotating frame, we just need one snapshot from the force-free simulation to create the all sky map.As a demonstration of principle, we consider emission coming from a spherical surface with radius r em .This is the point in the pulsar magnetosphere where an outgoing radio wave is assumed to be emitted.Because of the large uncertainties in the emission physics, we consider a large range of emission heights, from close to the surface all the way up to 10 times the radius of the pulsar4 .
Suppose each emission point has a polar angle θ and an azimuth angle ϕ, the emission direction n is determined using equation ( 8).We denote the polar angle of n as θ e and its azimuth angle as ϕ e .After emission, the ray propagates in a straight line, so it will reach an observer with polar angle θ obs = θ e , and the corresponding phase is determined by (e.g., Bai & Spitkovsky 2010) which takes into account the phase delay due to finite light travel times.This step enables the mapping from the emission point (ϕ, θ) to the observer plane (ϕ obs , θ obs ).For each ray, the polarization first evolves following the background magnetic field, then it freezes at a distance h from the emission point.This height is the distance above the r em where the radio wave decouples from the plasma and propagates freely before reaching the observer without changing its polarization.Therefore, we determine the polarization angle at the freezing point r = r em + hn using the local magnetic field, taking into account the rotation of the magnetosphere during the time interval it takes the ray to propagate through a distance h.The polarization angle is then obtained through equation ( 18): for each ray we put the observer along z and the rotation axis of the NS in the x-z plane, then we find the transverse components of the magnetic field B x , B y and the rotation velocity β x , β y at the freezing point.Note that numerically, both the numerator and the denominator in equation ( 18) are zero at the emission point by definition.Therefore, a small but finite h is needed even if we consider the usual RVM limit where the polarization is determined immediately at the emission point.
For our sky map, we select a 200 × 200 grid uniformly spaced in θ obs and ϕ obs .A ray originating from the emission surface at (ϕ i , θ j ) carries the information of the PA and maps to a particular grid point (ϕ k obs , θ l obs ) on the sky map.Once we have our sky map, we can fix the observing angle θ obs , and plot the PA as a function of ϕ obs .

RESULTS
As a first step, we perform a simple test to determine whether we get the usual result for a purely dipolar magnetic field.This is shown in figure 4, where two S-shaped swings from the two magnetic poles are clearly visible.The wider S curve (which gets broken into two parts by the phase beyond 180 • ) comes from the magnetic pole farther from the observer while the narrower one comes from the pole close to the observer.We also see the effect of aberration which makes the curves shift towards right, i.e. slightly towards higher ϕ obs .This effect is because aberration changes the orientation of the Evector in the wave and orientation of the polarization ellipse, which leads to changes in the PA profile.As per the RVM, most polarization observations are used to identify the Ω Ω Ω −µ µ µ meridional plane which is the plane containing the magnetic axis µ µ µ and the pulsar rotation axis Ω Ω Ω.In this scenario, the steepest gradient or the inflection point of the PA curve is contained within the Ω Ω Ω−µ µ µ plane, and coincides with the midpoint of the pulse profile precisely.Nevertheless, in the observer's frame, as Blaskiewicz et al. (1991); Dyks et al. (2004) pointed out, the aberration and retardation effects cause the PA inflexion point to appear delayed relative to the midpoint of the pulse profile.The magnitude of this delay is given by: where C is some constant coefficient (see equation 3).Next, we analyse the force-free magnetic field model for PSR J0030+0451.Due to the radio emission coming from the open field lines, we restrict our analysis to only using open field lines to obtain the sky map corresponding to the two magnetic poles.A view of the polar caps from different heights above the magnetic poles is shown in the left panel of figure 5.As the height increases, the poles become more oval-shaped.The sky   map at r em = 3R is also shown in the right panel of figure 5 as an example.
In figure 6, we plot the PA curves by fixing the viewing angle to θ obs = 54 • as reported in Bilous et al. (2019).Several different effects are explored, including the effect of varying the emission radius, varying the freezing height, and examining the effects of aberration.Our main aim is to reproduce the observed PA curves of PSR J0030+0451 at 430 MHz and at 1.4 GHz, as shown in Gentile et al. (2018).When comparing different PA curves, we focus more on their shape than on their exact values.Let us first discuss the PA curve observed at 430 MHz which is shown in figure 6.It has two distinct parts corresponding to two different poles (Gentile et al. 2018).The one between ϕ obs = (20 • -100 • ) comes from the smaller intensity pulse while the other comes from the higher intensity pulse.The latter demonstrates a flat and distorted PA profile, which hint that non-dipolar magnetic field effects are present.If the emission comes from very close to the surface, for example r em ∼ 3R, we can see only one of the magnetic poles (see figure 5), and a very small part of the PA curve5 .At very far from the surface (r em > 8R), the PA curve becomes S-like due to the field being dipolar6 there (see panel c in figure 6).Varying the freezing height (h) while fixing the emission radius also yields similar results.Hence, we identify a region at r em ∼ 5R where the PA becomes flat.This is shown by the blue dots obtained at r em = 5R and h = 0.4R in panel (a) figure 6.We see that for the higher intensity pulse, our model roughly reproduces the observed PA trend, but there is a significant discrepancy for the low intensity pulse.
In either decreasing or increasing h, an abrupt discontinuity is seen at ϕ ∼ 290 • , as shown by the squares obtained at r em = 5R and h = 1.0R (panel b).This feature is only present when aberration is taken into account in which case the E y component becomes close to zero along the trajectory.Without aberration, the curve reverses its direction as shown by the green diamonds in panel (d).On the other hand, when neglecting the effect of aberration and choosing a higher emission radius/freezing height, for example r em = 7R and h > 4R (panel e), or r em = 8R (panel f), it appears that the observed PA trend for both pulses can be more or less reproduced.However, the aberration effect is particularly important for MSPs and should not be neglected.Our failure to reproduce the PA trend for both pulse peaks at 430 MHz when aberration is included may suggest that either the other propagation effects need to be taken into account, like GR light bending close to the neutron star (Poutanen 2020), and refraction of the radio waves by the plasma (Beskin & Philippov 2012), or the field configuration obtained by Chen et al. (2020) is not the actual field of PSR J0030+0451.A few other possible field geometries exist, as shown by Kalapotharakos et al. (2021).Our modeling indicates that radio polarization is a more stringent test that may be able to break the degeneracy and select a more realistic field configuration.
Next, we examine the PA curve at 1.4 GHz, which again has two parts corresponding to two different intensity peaks in its pulse profile (see figure 1 in Gentile et al. 2018).The PA for the lower intensity pulse, shown in figure 7, resembles the tail of an S-shaped curve resulting from a dipolar magnetic field.To explain this, we explore emission radii/freezing height greater than 5R because the dipolar components of the magnetic field are expected to become dominant over the quadrupolar components at larger distances.As we can see, the PA curve considering aberration and corresponding to r em = 8R, and h = 9.9R exhibit a similar S-shaped swing (panel (a) in figure 7).Decreasing the emission radius (panel e) or the freezing height (panel b) makes the PA curve straight.Furthermore, neglecting the effect of aberration with the same set of parameters causes the PA curve to be concave (panel f) and thus different from the observed feature.A comparison of panels (a), (c), and (d) in figure 7 indicates that a higher freezing height is necessary to recover the observed PA.Furthermore, we cannot go lower than r em = 5R because the  line of sight of the observer does not cross through this magnetic pole.7 In order to explain the PA corresponding to the higher intensity peak which looks non-dipolar, we study emission coming from close to the surface of the NS.The PA is divided into three sections: an initial rise at the beginning, a distorted shape in the middle, and a decreasing tail at the end (figure 8).It turns out that this feature can be explained by varying the freezing height along the pulse.As the polarization is set when the emission decouples from the plasma, i.e. when the density sufficiently drops, decoupling at different distances might occur because of the plasma density distribution across magnetic field lines.The observed features can be recovered at r em = 2R, starting with h = 0.4R for the initial rise, h = 0.1R for the middle, and h = 0.01R for the decreasing tail.The error bars (standard deviation) correspond to a variation of 100% in h and the PA(s) are reported as the mean value.We note that very different emission radius r em and freezing height h are needed at the two poles in order to reproduce the observed PA trend at 1.4 GHz.Although the plasma density at the two poles could be different, leading to different r em and Sur et al.
h, it may be hard to account for such a large difference.This may again indicate that we either need to consider propagation effects, in particular GR light bending and plasma refraction close to the neutron star, or our field configuration is not the exact field possessed by PSR J0030+0451.

CONCLUSIONS AND DISCUSSIONS
Polarization modeling is a powerful tool and a very stringent test to constrain the magnetic field configuration of pulsars.An effective field model must pass all tests, including radio, X-ray, gamma-ray light curves, and radio polarization.In this paper, we compared radio polarization with observations for the MSP PSR J0030+0451 (which has detailed X-ray hotspot modeling) based on a multipolar magnetic field configuration, and we present easy equations for calculating the PA that can be used in parameter search models; if all physics is parameterized using r em and h.
We used the force-free configuration obtained by Chen et al. ( 2020) that can simultaneously produce the radio, X-ray, and gamma-ray light curves.In our models we calculate the polarization of plasma normal modes, taking into account the co-rotational motion of the plasma.We derive an expression of the PA by considering an Omode instead of calculating the properties of the emission of charged particles in vacuum as described in Blaskiewicz et al. (1991).We modified the RVM to see if the multipolar field configuration could produce the observed radio polarization features.To determine the PA curves, we first created sky maps of the polarization taking into account aberration, and then selected a viewing angle.For PSR J0031+0451, the viewing angle was fixed at θ obs = 54 • .We were able to reproduce some of the observed features of the PA swing both in the case of 430 MHz and 1.4 GHz observations by constraining the emission heights in each case at which emission is produced.
When considering aberration, for 430 MHz, the PA swing at the higher intensity pulse could be explained by emissions coming from an effective radius 5R with h = 0.4R.However, for the lower-intensity pulse, the PA swing couldn't be reproduced at any chosen emission radius/height.On the other hand, without taking aberration into account, both pulses can be roughly reproduced simultaneously by emissions from an effective radius greater than 7R.It is not possible to explain simultaneously high and low intensity pulses at 1.4 GHz by any emission radius.Hence, we examined effective emission radii and heights separately for the two pulses to determine whether the observed features were reproducible.Our first finding was that when aberration was not considered, the PA swing for the smaller intensity pulse at any radius did not match the data.Therefore, we considered aberration in all our models for 1.4 GHz.Considering emissions from a radius greater than 8R, it is possible to determine a similar PA swing for the smaller intensity pulse that resembles the tail of an S-curve.The higher intensity pulse with a distorted PA swing could be explained by emission from 2R and three different freezing heights.We note, therefore, that an important physics behind the MSP radio polarization model is including aberration effect, since without it many of the observed characteristics cannot be explained.
All our models, however, ignore propagation effects.Radio waves propagate through dense magnetospheric plasma, where the polarization signature first evolves adiabatically following the magnetic field, before becoming permanent as they reach a larger distance where the plasma density has dropped (e.g., Petrova & Lyubarskii 2000;Wang et al. 2010;Beskin & Philippov 2012).Effects like wave refraction, cyclotron absorption, and transition from geometrical optics to vacuum propagation can all influence the observed polarization.These propagation effects will depend on the properties of the magnetospheric plasma, e.g., the plasma density and its radial dependence, the drift motion of plasma particles, and the distribution function of the outgoing plasma.It has also been found that pulsars whose PA profiles are not fitted with the RVM exhibit a much higher fraction of circular polarization than those with linear polarization (e.g., Johnston et al. 2023).Circular polarization is usually considered to be a result of propagation effects in the magnetospheric plasma (e.g., Melrose & Luo 2004).Furthermore, if the emission is produced close to the neutron star surface, GR light bending also needs to be taken into account (e.g., Beloborodov 2002;Poutanen 2020), which we neglected in this work.A complete framework thus needs to include all the ingredients: a self-consistent magnetic field configuration, GR light bending, and propagation effects in the magnetospheric plasma.We plan to develop such a framework for polarization modeling in the future and combine it with multiwavelength light curve fitting.As a result of high-precision observational measurements, we may be able to constrain not only magnetic field configurations and radio emission sites but also magnetospheric plasma properties, which include density, bulk Lorentz factor, etc.This will give us further insights into the properties of pulsar plasma and radio emission physics.
The multiwavelength light curve and radio polarization modeling of PSR J0030+0451 suggest that a multipolar magnetic field may be important near the neutron star surface.A new analysis of the NICER data for PSR J0030+0451 to constrain the size and location of its Xray hot spots shows that there is a multi-modal structure in the posterior surface (Vinciguerra et al. 2023).Given the uncertainty and degeneracy in X-ray fitting, radio polarization may turn out to be an important constraint to help distinguish between different field configurations.Obtaining a good understanding of the field configuration can have implications for other branches of pulsar physics as well, including the interior magnetic field evolution, formation of gravitational-wave mountains, and damping of r-modes when considering a superconducting core.These will be investigated in our future works.

Figure 1 .
Figure1.The rotating NS is located in the center of XYZ frame with its magnetic moment axis (⃗ p) and rotation axis ( ⃗ Ω). ⃗ Ω lies in the x-z plane.The observer is located on the z-axis, the green dotted line shows the magnetic field line, while the curly black solid line shows a ray propagating towards the observer.

Figure 2 .
Figure 2. The vacuum magnetic field configuration as seen from the meridional view along the x-axis.The left figure shows the closed magnetic field lines with some of the quadrupolar loops while the right figure shows the open dipolar field lines.The rotation axis is along the z-axis and passes through the center of the star.

Figure 3 .
Figure 3.The force-free magnetic field configuration as seen from (left) meridional view (along x-axis), and (right) equatorial view.The green lines represent the closed field lines while the cyan lines represent the open field lines.Also shown in the background is the location of current sheets.This field configuration was obtained in Chen et al. (2020).

Figure 4 .
Figure 4. Left panel: the sky map for a vacuum dipolar magnetic field with an inclination angle of α = 80 • , and emission radius rem = 3R.The color scale represents PA values while the horizontal black solid line represents PSR J0030+0451's viewing angle from earth.Right panel: the corresponding PA curves with and without aberration were observed from a viewing angle of 79 degrees and freezing height h = 0.1R.

Figure 5 .
Figure5.Left panel: the shape of the polar caps at r = 2R (red), r = 3R (blue), r = 5R (green), and r = 7R (violet), as seen by the observer.Close to the surface, one of the polar caps at 2R is crescent-shaped, while the other pole is more circular in nature.With increasing height, the poles become oval-shaped.Right panel: the sky map at rem = 3R with the colorscale representing the PA values.Again, the horizontal black line is the line of sight of the observer set at 54 degrees.

Figure 6 .
Figure 6.A comparison of the different PA curves represented by the scatter plots, obtained from our models by varying the freeing height (panel b), emission radius (panel c), and not including the effect of aberration (bottom-row) for the same parameters.The black dots shows the observation at 430 MHz.

Figure 7 .Figure 8 .
Figure 7.Comparison between our models and observed PA of PSR J0030+0451 at 1.4 GHz corresponding to the smaller intensity peak in its pulse profile.