Modeling X-Ray and Gamma-Ray Emission from Redback Pulsar Binaries

We investigated the multiband emission from the pulsar binaries XSS J12270−4859, PSR J2039−5617, and PSR J2339−0533, which exhibit orbital modulation in the X-ray and gamma-ray bands. We constructed the sources’ broadband spectral energy distributions and multiband orbital light curves by supplementing our X-ray measurements with published gamma-ray results, and we modeled the data using intrabinary shock (IBS) scenarios. While the X-ray data were well explained by synchrotron emission from electrons/positrons in the IBS, the gamma-ray data were difficult to explain with the IBS components alone. Therefore, we explored other scenarios that had been suggested for gamma-ray emission from pulsar binaries: (1) inverse-Compton emission in the upstream unshocked wind zone and (2) synchrotron radiation from electrons/positrons interacting with the kilogauss magnetic field of the companion. Scenario (1) requires that the bulk motion of the wind substantially decelerates to ∼1000 km s−1 before reaching the IBS for increased residence time, in which case the formation of a strong shock is untenable, inconsistent with the X-ray phenomenology. Scenario (2) can explain the data if we assume the presence of electrons/positrons with a Lorentz factor ∼ 108 (∼0.1 PeV) that pass through the IBS and tap a substantial portion of the pulsar voltage drop. These findings raise the possibility that the orbitally modulating gamma-ray signals from pulsar binaries can provide insights into the flow structure and energy conversion within pulsar winds and particle acceleration nearing PeV energies in pulsars. These signals may also yield greater understanding of kilogauss magnetic fields potentially hosted by the low-mass stars in these systems.


INTRODUCTION
Pulsar binaries, composed of a millisecond pulsar and a low-mass (< M ⊙ ) companion (Fruchter et al. 1988), are thought to form when past accretion from the companion had spun up the pulsar over a timescale of ∼Gyr (Alpar et al. 1982).In pulsar binary systems, the pulsar is in a tight orbit with a companion whose spin is tidally locked with the < day orbit.The pulsar heats the pulsar-facing side of the companion, driving a strong wind from it.Additionally, the companion may have a strong magnetic field (B) of ≥ kG (e.g., Sanchez & Romani 2017;Kansabanik et al. 2021).The wind-wind or wind-B interaction (e.g., Harding & Gaisser 1990;Wadiasingh et al. 2018) be-Corresponding author: Hongjun An hjan@cbnu.ac.kr tween the pulsar (wind) and the companion (wind or B) can lead to shocks at their interface.The shocks wrap around the pulsar if the companion's pressure/wind is stronger than the pulsar's, or vice versa.
The hallmark properties of emission from pulsar binaries are multiband orbital modulations.The companion's blackbody (BB) emission orbitally modulates due to pulsar heating and the ellipsoidal deformation of the star, exhibiting characteristic day-night cycles (e.g., Breton et al. 2013;Strader et al. 2019).Modeling these day-night cycles has proven useful for determining system parameters such as orbital inclination and has facilitated estimations of the masses of the neutron stars in pulsar binaries (e.g., van Kerkwijk et al. 2011;Linares 2020).Nonthermal X-ray emission, originating from the intra-binary shock (IBS), also shows orbital modulation with a single-or double-peak structure (e.g., Huang et al. 2012) caused by Doppler beaming (relativistic aberration) by the IBS bulk flow.This modulation, when combined with orbital parameters, aids in probing particle acceleration and flow in the shock (e.g., Kandel et al. 2019;van der Merwe et al. 2020;Cortés & Sironi 2022).
Gamma-ray emission from pulsar binaries is dominated by pulsed magnetospheric radiation from the pulsar.Other than a sharp dip in the light curves (LCs) of eclipsing systems (Corbet et al. 2022;Clark et al. 2023), the gamma-ray emission has been thought to be orbitally constant.However, the Fermi Large Area Telescope (LAT; Atwood et al. 2009) has discovered orbital modulation in the GeV band in a few pulsar binaries (e.g., Ng et al. 2018;An et al. 2018;Clark et al. 2021), suggesting that there should be other physical mechanisms for gamma-ray production in these sources.The GeV modulation observed in the 'black widow' (BW; < 0.1M ⊙ companion; see Swihart et al. 2022, for a list of BWs) PSR J1311−3430 was interpreted as synchrotron radiation from the IBS particles (An et al. 2017;van der Merwe et al. 2020) because its gamma-ray LC peak appears at the phase where the IBS tail should be in the direction of the observer's line of sight (LoS).However, this synchrotron interpretation remains inconclusive due to the absence of evidence for IBS X-rays (i.e., an orbitally-modulated signal) in this source.
Orbitally-modulating X-ray emission with a broad maximum at the inferior conjunction of the pulsar (INFC; pulsar between the companion and observer) has been well detected in bright 'redbacks' (RBs; > 0.1M ⊙ companion) (e.g., Roberts 2013;Wadiasingh et al. 2017), and XSS J12270−4859, PSR J2039−5617, and PSR J2339−0533 (J1227, J2039, and J2339 hereafter) are no exception.These RBs are particularly intriguing because they exhibit orbital modulation in the GeV band as well (Ng et al. 2018;An et al. 2020;Clark et al. 2021;An 2022).Their LAT LCs have a maximum at the superior conjunction of the pulsar (SUPC).At this phase, the sources' X-ray emission is at a minimum level, indicating that the IBS tail is in the opposite direction of the LoS (Kandel et al. 2019).Hence, the scenario that ascribes the gamma rays to synchrotron radiation from the IBS electrons (van der Merwe et al. 2020) seems implausible for these RBs.Instead, it was suggested that inverse-Compton (IC) scattering off of the companion's BB photons by upstream (unshocked) wind particles may generate the gamma-ray modulation due to orbital variation of the scattering geometry.However, in this case, energy injection from the pulsar seemed insufficient (An et al. 2020;Clark et al. 2021).Alternatively, van der Merwe et al. (2020) suggested that synchrotron radiation from the upstream particles, which pass through the IBS and interact with the compan- b XMM13 and XMM16 for the earlier and later observation, respectively (Figure 1).
ion's B, may be responsible for the gamma-ray emission (see also Clark et al. 2021).This scenario may explain the emission strength and phasing of the gamma rays, but it is yet unclear whether the observed LC shapes can also be reproduced.
In this paper, we explore scenarios for the gammaray modulations discovered in the three RBs, J1227, J2039, and J2339, by modeling their multiband emission.We construct X-ray spectral energy distributions (SEDs) and orbital LCs of the targets, and supplement them with published LAT measurements (Section 3).We then apply an IBS model to the data and investigate mechanisms for the modulating GeV emissions from the sources in Section 4. Finally, we discuss the results and present our conclusions in Section 5.

X-RAY DATA ANALYSIS
We analyze the X-ray observations of the targets (Table 1), taking into account the potential BB emission from the pulsar.This consideration is crucial because the BB emission may significantly impact the measurements of the nonthermal IBS X-ray spectrum and LC.
Although there are archival X-ray observations of J1227 taken just after the source's state transition in late 2012 (∼MJD 56250; Bassa et al. 2014), we do not include them in our analysis.This is due to the potential effects of a residual accretion disk, as indicated by strong variability in the source's flux and LC shape observed between 2013and 2015(de Martino et al. 2015)).Consequently, we acquired new NuSTAR data in 2021, ≥8 yrs after the transition.The XMM observation of J2039 was analyzed by Romani (2015) and Salvetti et al. (2015); however, we reanalyze the data to accurately measure the 'nonthermal' emission from the source.This reanalysis involves a careful examination of contamination from the pulsar's BB emission.The data for J2339 were analyzed by Kandel et al. (2019), with a focus on phaseresolved spectroscopy and LCs.We reanalyze the data to construct a phase-averaged SED.

Data reduction
We processed the XMM data using the emproc and epproc tools integrated in SAS 2023412 1735.We further cleaned the data to minimize contamination by particle flares.Note that we did not use XMM timing-mode data due to the low signal-to-noise ratio.The Chandra observation was reprocessed with the chandra repro tool of CIAO 4.12, and the NuSTAR data were processed with the nupipeline tool in HEASOFT v6.31.1 using the saa mode=optimized flag.For the analyses below, we employed circular source regions with radii of R = 3 ′′ , R = 16 ′′ and R = 30 ′′ for Chandra, XMM, and NuSTAR, respectively.Backgrounds were extracted from source-free regions with radii of R = 6 ′′ for Chandra, R = 32 ′′ (or R = 16 ′′ for small-window data) for XMM, and R = 45 ′′ for NuSTAR.

X-ray light curves
To generate orbital LCs of the sources for use in our modeling (Section 4), we barycenter-corrected the photon arrival times and folded them using the orbital period (P B ) and the time of the ascending node (T ASC ) measured for each source (Roy et al. 2015;Clark et al. 2021;An et al. 2020).We should note that the exposures of the observations are not integer multiples of the orbital periods, resulting in uneven coverage of orbital phases.Furthermore, there are observational gaps caused by flare-removal (XMM) and Earth occultation (NuSTAR), which also introduce nonuniformity in the phase coverage.
We investigated the effects of the nonuniform phase coverage and found that the observational gaps present in the XMM and NuSTAR data were randomly distributed in phase, resulting in spiky features in the LCs.However, these random variations did not significantly impact the overall flux measurements, which remained within a range of < ∼ 1-2%.In the cases of the Chandra LC of J2339 and the XMM LC of J2039, systematic trends were observed.The Chandra observation had 2× more exposure for a phase interval when the source appeared bright, while the XMM data had ∼50-60% longer exposure near the minimum phase.Due to the significant distortion of the LCs caused by both random and systematic exposure variations, we corrected the LCs for the unequal exposure.We computed the source exposures using the good time intervals of the observations, folded the exposures on the orbital periods, and divided the LCs by the folded exposure to account for the exposure variations.The exposure-corrected LCs are presented in Figure 1, where they exhibit a single-or double-peak structure.To minimize contamination from the orbitally-constant BB emission (Section 2.3) while ensuring high signalto-noise ratios, we utilized the 2-10 keV and 3-20 keV bands for the Chandra/XMM and NuSTAR LCs, respectively.The LC of J1227, measured with the 2021 NuSTAR data, displays a broad single bump (Figure 1 left) similar to the 2015 NuSTAR LC (de Martino et al. 2020).This similarity suggests that the LC of J1227 has remained stable since 2015.It is worth noting that the LC of this source exhibited significant morphological changes between 2013 and 2015, immediately after the state transition (Bogdanov et al. 2014;de Martino et al. 2015).The multi-epoch LCs of J2339 appear almost the same (Figure 1 right; see also Kandel et al. 2019), indicating that the source has been stable for ∼2600 days.

Spectral analysis
As mentioned previously, the random variations in exposure do not pose a concern for spectral analyses.However, the systematic excesses in exposure in the Chandra (J2339) and the XMM (J2039) data affected the spectral measurements.Hence, we removed the time inter- vals corresponding to the excess exposure from the data, which resulted in a ∼10% reduction in the livetime.
We generated the X-ray spectra of the targets and created corresponding response files using the standard procedure suitable for each observatory.The spectra were grouped to contain at least 20 events per bin and were then fit with an absorbed BB or power-law (PL) model using the χ 2 statistic.For Galactic absorption, we employed the tbabs model with the wilm abundances (Wilms et al. 2000) and the vern cross section (Verner et al. 1996).The hydrogen column densities (N H ) for these high-latitude sources were inferred to be low, but their values were not well constrained by the fits (see also de Martino et al. 2020;Romani 2015;Kandel et al. 2019).Therefore, we held N H fixed at the values estimated based on the HI4PI map (HI4PI Collaboration et al. 2016). 1  Below, we present results of our spectral analysis for each source.
J1227: The new 2021 NuSTAR spectra of J1227 were adequately described by a simple PL (Figure 2 and Table 2), and the results are very similar to the previous ones obtained from the 2015 NuSTAR data (see de Martino et al. 2020).However, the source flux was ∼20-30% lower in 2021 compared to 2015.Although this flux difference is not highly significant ( < ∼ 3σ), we use the new results for our modeling (Section 4).We thoroughly examined the archival XMM data acquired in 2013-2014 (after the transition) to identify a potential BB component, but no significant evidence for it was found, as previously noted by de Martino et al. (2015).
J2039: For this source, the PL model was acceptable, whereas the BB model was statistically ruled out.Our PL results are consistent with the previous measurements (Romani 2015;Salvetti et al. 2015).However, we noticed a trend in the fit residuals.To address this, 1 https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.plNote.Numbers in parenthesis represent the 1σ uncertainty.
we fit the spectrum with a BB+PL model (Figure 2).We then conducted an f test to compare the PL and BB+PL models, and found that the BB+PL model was favored over the PL model with an f -test probability of 5 × 10 −6 ; we verified this result using simulations.We report the best-fit parameters of the BB+PL model in Table 2. Salvetti et al. (2015) also favored the BB+PL model, but they did not report the best-fit parameter values.J2339: We jointly fit the Chandra, XMM, and NuS-TAR spectra of J2339, allowing a cross normalization for each spectrum to vary.Both the BB and PL models were statistically ruled out with the χ 2 probabilities of < 10 −5 , and we observed a discernible trend in the residual, similar to the case of J2039.Therefore, we fit the data with a BB+PL model, and this model provided an acceptable fit with the χ 2 probability of 0.18.This finding significantly strengthens the previous suggestion of BB emission from the source (f -test probability of 10 −3 ; Yatsu et al. 2015).

Broadband SED and LC data
We construct the broadband SEDs and LCs of our target RBs (Table 1).Our analysis of the X-ray data pro- vided nonthermal X-ray SEDs and LCs (Figure 1 and Table 2).For the modeling, we converted the count units of the X-ray LCs into flux units by comparing the phase-averaged flux to the observed counts for each source.For the LAT data, we adopted the published results (Ng et al. 2018;An et al. 2020;Clark et al. 2021;An 2022).These previous analyses used different energy bands: 60 MeV-1 GeV for J1227, 100 MeV-100 GeV for J2039, and 100 MeV-600 MeV for J2339.To ensure consistency, we converted the count units of the LCs into flux units using the LAT models provided in the Fermi-DR3 catalog (Abdollahi et al. 2022).Additionally, we subtracted constant levels from the gamma-ray LCs, assuming that the constant emission originates from the pulsar's magnetosphere rather than the IBS or upstream wind (see also Clark et al. 2021).
Since the spectrum of the orbitally-modulated LAT signals has not been well measured, we present in Figure 3 flux levels (horizontal lines) estimated by scaling down the pulsar SEDs according to the modulation fractions of the LAT LCs (typically ∼30%).The actual IBS spectra might be softer than those of the pulsars as the observed modulations of the targets were more pronounced at low energies (e.g., Ng et al. 2018;An et al. 2020;An 2022).Therefore, the gamma-ray SEDs of the modulating signals are likely to peak at ≤ GeV energies (see also Figure 4 of Clark et al. 2021).The spectra of the optical companions, which provide seeds for IC emission, were obtained from the literature (Rivera Sandoval et al. 2018;Kandel et al. 2019;Clark et al. 2021) and the VizieR photometry database.2These spectra represent the observed emission near the optical-maximum phase.The broadband SEDs and gamma-ray LCs of the targets are displayed in Figure 3.

Emission scenarios
In IBS scenarios, a relativistic electron/positron plasma (advected in the MHD pulsar wind) originating from a pulsar is injected into an IBS formed by wind-wind or wind-B interaction.The electrons (electrons+positrons) are accelerated at the shock and flow along the IBS.These IBS electrons emit synchrotron radiation in the X-ray band, which is Doppler-boosted along the bulk flow direction (Figure 4).The companion provides seeds for IC scattering to the electrons in the IBS and in the pulsarwind region.While this scenario has been successful in modeling the X-ray SEDs and LCs of pulsar binaries (Romani & Sanchez 2016;Wadiasingh et al. 2017;Kandel et al. 2019;van der Merwe et al. 2020;Cortés & Sironi 2022), it was suggested that the scenario cannot explain the recently-discovered GeV modulations from our RB targets (An et al. 2020;Clark et al. 2021, see Section 3.2.1).We check this basic scenario (Scenario 1) to confirm the previous suggestion, and we adjust the parameters within this scenario to offer a phenomenological explanation for the data.Then, we explore an alternative scenario to explain the LAT measurements.Note that these scenarios share the same mechanism for the X-ray emission (synchrotron radiation from IBS electrons), and therefore, our descriptions of the scenarios concentrate on the gamma-ray emission mechanisms.
• Scenario 1: This is the basic IBS scenario where electrons in the cold wind and in the IBS ICupscatter the companion's BB photons to produce gamma-ray emission (Figure 4 left).
• To summarize, we consider three emission zones as listed below.
• Wind zone (cyan in Figure 4): This is an emission zone between the pulsar's light cylinder and the IBS.Electrons in this zone are assumed to be cold (i.e., δ distribution) and relativistic (but see Section 4.1.1 for Scenario 1a).In our phenomenological study, we consider only the IC emission from the electrons, assuming that their synchrotron emission is weak (but see Section 5).
• IBS zone (pink in Figure 4): The electrons in the wind zone are injected into this IBS zone, and thus the number and energy of the electrons in this zone are connected to those in the wind zone (Equations ( 12)-( 14) below).In this IBS zone, shock-accelerated electrons flow along the IBS surface, and they produce both synchrotron and IC emissions.We do not consider synchrotron-self-Compton emission from this zone, as its flux has been assessed to be negligibly small (van der Merwe et al. 2020).
• Companion zone (green in Figure 4 right): This emission zone is used only for Scenario 2. Most of the upstream pairs interact with the IBS zone, but a sufficiently energetic component of them with large gyro radii are assumed to pass through the IBS unaffected and reach this zone.These electrons can produce both synchrotron and IC emission with the former dominating over the latter due to the strong B of the companion.Therefore, we consider only the synchrotron emission.
We investigate the fundamental properties of the emission scenarios using simplified calculations in Sections 3.2.1-3.2.2, where, for simplicity, we assume that the energy distributions of the particles and their radiation are δ functions and that the IC scattering occurs in the Thomson regime.We also neglect Doppler boosting caused by the bulk motion of the IBS flow (e.g., Wadiasingh et al. 2017) and carry out our analytic investigation in the flow rest frame (equivalent to the observer frame in this section), as Doppler factors of IBS flows in pulsar binaries have been inferred to be small (e.g., Romani & Sanchez 2016;Kandel et al. 2019, see also Table 4).These analytic calculations, made utilizing a one-zone approach and mono-energetic distributions, provide rough estimates of the model parameters that will be used as inputs for detailed multi-zone IBS computations performed without the aforementioned assumptions (Section 4).

Scenario 1: Basic IBS scenario
The energy distribution of electrons in IBSs is often assumed to be a power law dN dγedt ∝ γ −p1 e .Since a pulsar supplies energy to the IBS, we require where η s represents the energy conversion efficiency of the IBS, ĖSD is the pulsar's spin-down power, and f Ω (assumed to be 1 in this section) is the fraction of the solid angle subtended by the IBS.A fraction (η γ ) of ĖSD is converted to the pulsar's gamma-ray radiation (Table 3; see also Abdo et al. 2013), while the remaining energy is eventually converted into the particles' energy (fraction η s ) and the magnetic energy (fraction η B ).We assume that the magnetic energy and radiative energy loss are negligibly small within the IBS (e.g., Kennel & Coroniti 1984).For a mono-energetic electron distribution dN dγedt = Ṅs δ(γ e − γ s ), Equation (1) can be rewritten as where t s is the residence time in the emission zone.It is given by l s /v IBS , where l s is the length of the IBS (Figure 4) and v IBS represents the bulk-flow speed within the IBS.The synchrotron-emission frequency of these electrons is given by and the observed flux would be where σ T is the Thomson scattering cross section, 8π is the magnetic energy density, and d is the distance between the observer and the source.For computation of the flux in a certain energy band, particle cooling needs to be considered.t s is typically shorter than the cooling timescale in IBSs of pulsar binaries (Table 4; see also van der Merwe et al. 2020).So we assume the emission timescale in the IBS (τ s ) to be approximately the residence time.By combining Equations ( 2)-( 4), we obtain in the IBS, where d kpc is d in units of kpc, F X,−13 is the X-ray flux in units of 10 −13 erg s −1 cm −2 , and ĖSD,35 is ĖSD in units of 10 35 erg s −1 .
In IBS scenarios, gamma-ray emission is assumed to arise from IC processes involving electrons in the IBS and wind zones (Figure 4).The emission power generated by IC scattering between an electron with a Lorentz factor γ e (≫ 1) and seed photons with energy density u * and frequency ν * is given by (e.g., Dubus 2013), where β e = 1 − 1/γ 2 e ≈ 1 and µ is the cosine of the scattering angle (µ = −1 for head-on collisions).µ is one of the key parameters that determines the shape and phasing of the gamma-ray LCs.The observed frequency of the IC-upscattered photons is given by The flux of the IC emission from the IBS electrons can be determined by independent of the number of electrons (N s ) in the IBS, as the same electrons produce both synchrotron and IC emission.
On the contrary, to estimate the IC flux from the cold 'wind' particles (see Equation ( 15) below), we need to determine their number.Assuming a δ distribution for the cold wind particles, we have where γ w represents the Lorentz factor of the upstream particles, Ṅw is the number of particles injected (per second) by the pulsar into the wind zone (cyan region in Figure 4), and η w (< 1) is an efficiency factor that accounts for the conversion of the pulsar's energy output into particles within the wind zone.The total number of particles in the wind zone is then given by where t w is the residence time, l w is the size of the zone (cyan in Figure 4), and v wind is the bulk speed of the upstream wind with v wind = c 1 − 1/γ 2 w in this 'cold-wind' case.Because the upstream wind pairs are subsequently injected into the IBS and the B energy is further converted to particle energy in the IBS (Kennel & Coroniti 1984;Sironi & Spitkovsky 2011), we require Ṅs = Ṅw , and as magnetic energy may not be fully dissipated in the wind zone, and radiative energy losses in the IBS are negligible.These equations imply We should note that Equation ( 14) is applicable exclusively to mono-energetic distributions for a representative spatial zone.For arbitrary phase space distributions, one should substitute γ w and γ s with their spatial and momenta averages in the volume of interest (Section 4).For homogeneous one-zone models, as considered here, Equation ( 14) involves calculating the averages using the expression: γ dN dγdt dγ/ dN dγdt dγ.By combining Equations ( 6), ( 7), ( 10) and ( 11), the IC flux of the upstream particles, e.g., in the case of head-on collisions, can be computed as where τ w is the emission timescale in the wind zone.In the case of the cold wind with v wind ≈ c in this scenario (Scenario 1), the residence time is t w ≈1 s, which is shorter than the cooling timescale of electrons with γ w < ∼ 10 8 .Since γ w is expected to be ≈ 10 4 in this scenario (see below), we assume τ w ≈ t w .
These computations can be compared with the observed X-ray and LAT data (Figure 3).For an IBS that extends to the size of the orbit (a orb ≈ l s ), we find τ s ≈ a orb vIBS > a orb c ≈ 10 s.From this, we can infer B ≈ 1 G (Equation ( 5)) and γ s ≈ 10 6 (Equation (3)) from the observed X-ray SEDs with F SY > ∼ 10 −13 erg s −1 cm −2 at hν SY ≈ 10 keV (e.g., Figure 3).The optical seeds provided by the companion have hν * ≈ 1 eV and u * ≈ 0.1 erg cm −3 at the position of the IBS (≈ a orb ).Then, the IC flux of the IBS particles would be F IC,IBS ≈ 10 −13 erg s −1 cm −2 (Equation 8), which may explain (part of) the observed LAT fluxes of the targets (Figure 3).However, the peak of this IC-in-IBS emission is expected to be in the ∼TeV band (Equation (7) and Figure 3; see also van der Merwe et al. 2020).Therefore, IC-in-IBS emission cannot explain the LAT measurements.
In this scenario, additional gamma-ray emission arises from IC scattering by upstream particles.We can adjust γ w of the 'cold and relativistic' upstream electrons such that their IC emission peaks at < ∼ GeV, which requires γ w ≈ 10 4 (Equation ( 7)) and η w ≈ 0.01η s ≈ 0.01 (Equation ( 14)).This requires that the wind zone is Poynting-flux dominated (e.g., Cortés & Sironi 2022).The wind zone, extending from the pulsar's light cylinder to the IBS (cyan in Figure 4), is < a orb , and thus τ w < a orb c < 10 s.Consequently, F IC,wind is typically < 10 −14 erg s −1 cm −2 (Equation ( 15)), which is orders of magnitude lower than the observed GeV fluxes of the targets (Figure 3).This issue can be alleviated if τ w (≈ t w ) becomes longer, e.g., by deceleration of the bulk speed of the upstream flow (e.g., t w ≈ a orb c vs t w ≈ lw v wind with v wind ≪ c).This means that the upstream wind is not cold any more and is thermalized (see Section 4.1.1).For electrons with γ w ∼ 10 4 , the cooling timescale (> 10 4 s) is longer than the residence time if v wind ≥ 10 km s −1 .Therefore, we assume τ w ≈ t w = l w /v wind and proceed to estimate the required v wind to match the GeV flux.
Then Equations ( 2), ( 4) and ( 15) give Assuming v IBS ∼ c, l s ∼ l w ∼ a orb , u * ∼ 0.1 erg cm −3 , B ∼0.1 G, and γ w /γ s ∼ 0.01, we estimate v wind to be ∼ 10 −3 c ≈ 300 km s −1 if F IC /F SY ≈ 1 (Figure 3).We can accommodate various values of F IC in this scenario by adjusting v wind and/or l w (n.b. this is also equivalent to artificially reducing the spatial diffusion coefficient to an unphysically small value).Furthermore, this scenario can explain the phasing of the LAT LC since the modelpredicted gamma-ray flux would be maximum at SUPC because of the favorable scattering geometry (Figure 4).However, our order-of-magnitude estimate above hinges on a considerably lower v wind than upstream precursor regions presented in PIC simulations of pulsar wind shocks (e.g., Sironi & Spitkovsky 2011).Furthermore, the formation of strong shocks and particle acceleration (as demanded by the X-ray phenomenology) appears implausible if the speed of the pulsar wind is indeed decelerated to such a low value.Notwithstanding these discrepancies, we conduct more detailed computations based on this scenario (termed 'Scenario 1a' in Figure 3 and below) in Section 4, where we adjust v wind arbitrarily for a comprehensive analysis.Electrons with such a high Lorentz factor (primaries) can be accelerated by the total potential drop available to the pulsar (e.g., van der Merwe et al. 2020) reaching ∼ 0.1 PeV energies.Such high Lorentz factors proximate to the pulsar light cylinder are also required for GeV emission in the primary curvature radiation scenario (e.g., Harding & Kalapotharakos 2015;Kalapotharakos et al. 2019;Harding et al. 2021;Kalapotharakos et al. 2023).A large fraction of the primary pairs spawns numerous secondaries with low energies, and the remaining primaries can easily pass through the IBS unaffected since their gyro radius is very large ( > ∼ a orb ).Thus, in this scenario, there are two populations of electrons in the upstream wind zone: one ∼ 0.1 PeV component with a high Lorentz factor (primaries, γ p ≈ 10 8 ), and the other with a lower Lorentz factor (secondaries with γ w ≈ 10 4 ; referred to as 'wind' in Scenario 1).The relationships between these two populations are given as follows.A fraction (η p ) of the pulsar's spin-down power is converted to the energy of the primary electrons: where Ṅp represents the number of primary electrons injected (per second) by the pulsar.Assuming a fraction (ζ) of these electrons penetrates the shock, the number of secondary electrons in the wind zone can be calculated using where M stands for the pair multiplicity.Based on energy conservation it follows that: and As these secondary electrons are subsequently injected into the IBS, Equations ( 12)-( 14) still hold.In this scenario, the upstream electrons are assumed to be cold.These high-energy electrons can pass through the IBS and emit synchrotron radiation in the companion's magnetosphere.In this case, the observed GeV flux is primarily contributed by electrons traveling along the LoS.These electrons can interact with a strong B (e.g., ∼kG) when in close proximity to the companion, particularly during the SUPC phase.The combination of high B and γ p results in a very short synchrotron cooling time (≪1 s), given by This cooling time is much shorter than the residence time (t comp = l comp /c) for any reasonable emission-zone size l comp (e.g., ∼ a orb ).So the emission timescale τ comp can be approximated to be t cool during orbital phases around SUPC.Then, the synchrotron flux arising from the companion's magnetosphere can be estimated (e.g., Equation ( 4)) to be This scenario can explain the LAT fluxes of our targets if ζη p > ∼ 0.1.Notice that B does not appear in Equation ( 23).This omission is a result of utilizing t cool (∝ B −2 ) for the emission timescale, under the assumption that it is significantly shorter than t comp .This assumption is valid specifically during orbital phases near SUPC.However, at other phases, it is more appropriate to employ t comp instead of t cool , and the usual B 2 dependence of the flux is reinstated (see below).
While it might seem that this scenario does not predict orbital modulation of the gamma-ray flux (Equation ( 23)), changes in B depending on the distance r * (B ∝ r −3 * for dipole B of the companion) between the companion and the emission zone can induce gamma-ray modulation for two reasons.First, the frequency of the synchrotron emission varies proportionally to B for a given γ p .The observed flux will be high if this peak frequency falls within the observational band, e.g., during the SUPC phase (Figure 5).Second, low B at certain orbital phases (e.g., large r * ; Figure 4 right) increases t cool , potentially making it longer than t comp when B is sufficiently low.In such cases, the kinetic energy of the electrons does not fully convert into radiation within t comp , leading to a decrease in the emission flux (e.g., φ = 0.55 in Figure 5).These two processes can result in a variety of LC shapes (Figure 6 and Section 4.1.3).

MODELING OF THE MULTIBAND DATA
The analytic exploration, performed with a one-zone approach and mono-energetic distributions, in the previous section provides general properties of the emissions from the RBs and establishes initial values.In this section, we leverage these findings to conduct more detailed and precise investigations of the emissions through our numerical model, utilizing a multi-zone approach and non-mono-energetic distributions.

The computational methods
In this section, we describe the emission zones and computational methods used in our numerical model.See Kim et al. (2022) for a comprehensive understanding of the model components, parameters, and their covariance.

Pulsar wind zone
We assume that the pulsar wind zone (blue region in Figure 7) starts from the light cylinder at a distance of r p = R LC = cP 2π from the pulsar, where P denotes the spin period of the pulsar.This zone extends to the The wind zone, situated between the pulsar's light cylinder and the IBS (depicted in blue), contrasts with the green companion zone.The blue and green dashed arrows represent the regions within the wind and companion zones, respectively, where electrons travel parallel to the LoS.Only these electrons contribute to the observable emission due to strong Lorentz boosting.The pulsar-centric coordinates (rp, θp, and φp) describe emission regions, while companion-centric coordinates (r * , θ * , and φ * ) detail seed photon density, direction and the companion's B. R 0 denotes the distance between the pulsar and the shock nose, and s is the distance along the IBS from the shock nose.
IBS surface at r p = r IBS marked by the pink region in Figure 7.
We assume that relativistic mono-energetic electrons are injected into this zone by the pulsar (Equations ( 9) and ( 10)).As these electrons move ballistically in the radial direction within this zone at highly relativistic speeds, their emission is strongly beamed along the direction of motion.Therefore, we compute IC emission only from the electrons propagating along the LoS, denoted by blue dashed arrow in Figure 7.Given that both the density of the seed photons (from the companion's BB) and the IC scattering geometry vary based on the location and flow direction of the emitting electrons with respect to the companion, we divide the emission zone (essentially a line) into 100 segments.
In each segment characterized by a length of dr p , we assign dr p c Ṅw (24) electrons with the Lorentz factor of γ w and compute their IC emission (Section 4.2.2).Note that the factor for the solid-angle fraction subtended by the observer (1/4πd 2 ), which is necessary because of the approxima- tion that the scattered photons are in the same direction as the scattering electrons, is not included in the above equation because it is accounted for in the emission formula (Equation ( 44)).As previously mentioned in Section 3.1, the cooling timescale within this zone exceeds the flow timescale, and thus we opt to neglect the IC cooling of the electrons.The emission frequency and flux are determined by γ w , η w , and the size of the wind zone, given the parameters of the pulsar, companion, and orbit (Table 3).The IC scattering geometry, the angle between the electron's motion (blue arrow in Figure 7) and the incident seed photons from the companion, drives orbital modulation in the GeV LC.This effect is incorporated into the IC formulas (Section 4.2.2).The size of the wind zone changes based on the IBS parameters (see below), but this variation has little influence on the flux.Additionally, in this study, we adopt the maximum possible value for η w explore the limiting case and reduce the number of adjustable parameters.It is worth noting that smaller values of η w (magnetically-dominated wind; see Cortés & Sironi 2022) can also explain the data (Section 4.3).We optimize γ w to match the GeV flux (Figure 8).
Given that this basic scenario (Scenario 1) fails to yield sufficient GeV flux (Figure 3), we modify this scenario by arbitrarily reducing the flow speed (v wind ) in the wind zone (increasing their residence time by decreasing their effective spatial diffusion coefficient).In this decelerated wind case (Scenario 1a; Section 3.2.1),we assume a departure from the relativistic cold wind scenario described above (Scenario 1), i.e., decelerated and thermalized, impacting both the electron distribution, and consequently their emission.Alongside bulk deceleration, we posit that electrons become isotropized in space and follow a relativistic Maxwellian distribution (instead of a δ distribution) in the electron rest frame: where β e = 1 − 1/γ 2 e , Θ shapes the distribution, and K 2 is the modified Bessel function of the second kind.We adjust Θ to ensure that the distribution peaks at γ w (we use γ w instead of Θ as our model parameter).Due to the isotropization, electrons traveling along the LoS exist at every point within the wind zone.To account for this, we integrate over the entire wind zone, dividing it into 10×80×180 radial, polar, and azimuthal (r p , θ p , and φ p in Figures 7) regions.The decelerated bulk speed in the wind zone is modeled as a constant function: The number of electrons in each of the 10×80×180 regions is given by Ṅw where f ∆Ω represents the solid-angle fraction of the emission region.We calculate the IC emission from these electrons, accounting for Doppler beaming and anisotropic IC scattering (Section 4.2).We should note that the energy conservation for these "non-monoenergetic" electrons γ e m e c 2 dN dγ e dt = η w ĖSD ( 28) requires (c.f., Equation ( 10)) modification by a factor of 1/ 1 − (v wind /c) 2 due to the bulk motion.However, we neglect this effect as it is less than 1% for the estimated v wind ≈ 1000 km s −1 .Additionally, we ignore the IC cooling of the electrons, as its impact is insignificant (Section 3.1).

IBS zone
We assume that an IBS is formed by the interaction between two isotropic winds from the pulsar and companion.In this case, the IBS's shape is dictated by the momentum flux ratio where Ṁ * v * denotes the momentum flux of the companion wind (Canto et al. 1996).Given specific values of a orb and β values (Tables 3 and 4), the curve  4 serve as the baseline in this example.
r IBS (θ p ) (pink in Figure 7) that defines the IBS in the cross-section plane is computed using the formulas from Canto et al. (1996): To construct the IBS surface, we rotate this r IBS (θ p ) curve around the line of centers.Although ĖSD was well-measured for the pulsars in our target RBs, the companion's Ṁ * v * remains unknown.Therefore, we adjust β to explain the observed SEDs and LCs with the model.The IBS opening angle increases with increasing β, leading to a larger phase separation between the LC peaks (Figure 9).Another parameter that needs to be specified is the length of the IBS (l s ).Energetic electrons, responsible for generating high-energy radiation, are predominantly situated in regions close to the shock nose (e.g., Dubus et al. 2015).In this work, we adopt l s = 4a orb .
However, we note that different values of l s can also explain the data well.The resulting LCs and SEDs do not alter much, and any difference caused by changes in l s can be compensated for by adjusting other parameters.
For the computation of the IBS emission, we discretize the IBS surface into 80 × 180 emission regions along the θ p and φ p directions.This is necessary because both B and the seed photon density (u * ) within the IBS vary over its surface.The latter can be computed using the orbit and companion parameters (Table 3).The geometry of the IC scattering, which produces anisotropy in the emission, is computed at each of the 80 × 180 regions.B is modeled according to (e.g., Romani & Sanchez 2016): where B s and R 0 denote B and r p at the IBS nose (s = 0), respectively.The parameter B s influences the flux and is optimized during our modeling process.Moreover, we consider the bulk motion of the electrons in the tangent direction of the IBS as follows.We assume that the bulk motion of the electrons undergoes acceleration as they flow along the shock (e.g., Bogovalov et al. 2008), such that the bulk Lorentz factor of the flow is given by Γ The energy distribution of the electrons within the IBS is modeled as a power law with an index p 1 between γ e,min and γ e,max in the flow rest frame.In reality, the energy distribution is anticipated to vary across regions at different distances (s; Figure 7) from the shock nose due to various effects, including radiative cooling and bulk acceleration.For simplicity, we do not consider the spectral changes between emission regions caused by these effects; hence, the same spectral shape is applied across the IBS.However, the energy budget is taken into consideration by comparing the particle energy flowing out of the IBS (at its tail) with the pulsar's injection, as follows.
As electrons flow along the IBS in the θ p direction, continuously replenished by the pulsar, the number distribution along the IBS (integrated over φ p ) is described by continuity (e.g., Canto et al. 1996) as Assuming a (normalized) power-law energy distribution of the electrons in the flow rest frame the particle energy in the observer frame is computed as The number of particles exiting the emission zone per second is Ṅ (θ p,max ), and thus the energy budget of the pulsar is governed by where θ p,max corresponds to θ p at l s,max (IBS tail).Optimization of γ e,min , γ e,max , and p 1 performed to explain the observational (primarily X-ray) data while satisfying Equation ( 37).
The bulk flow induces emission anisotropy and orbital modulation in both synchrotron and IC emissions.The strength of the emissions relies on the Doppler factor of the flow (e.g., Dermer et al. 2009;An & Romani 2017): where v IBS = c 1 − 1/Γ 2 IBS and θ V is the angle between the flow direction and the LoS.Orbital changes in the LoS relative to the IBS (and the binary system) and thus θ V induce the orbital modulation of the LCs.For a given i (Table 3) and the direction of flow (IBS tangent) at an orbital phase, δ IBS in each of the 80 × 180 emission regions is determined by only one parameter, Γ D .This parameter is optimized to explain the SEDs and LCs.This optimization of Γ D impacts the sharpness of the LC peaks (e.g., see Figure 9) and the number of electrons in the IBS (Equations ( 36) and ( 37)).
The emission from the IBS remains essentially the same among the scenarios discussed in this paper, with minor variation attributable to differences in the numbers of electrons injected from the wind zone, induced by shock penetration (e.g., Equations ( 17)-( 21)).We should note that the IBS parameters are linked to the wind parameters (e.g., Equation ( 34)); these parameters are consistently adjusted.

Companion zone
This zone starts at the IBS location (r p = r IBS ; Equation (30)) and extends outward from the pulsar (green region in Figure 7).Given that the electrons in this zone are energetic and cold pulsar-wind particles that have penetrated the IBS (Section 3.2), they move ballistically away from the pulsar.Similar to the wind zone (Section 4.1.1),we consider emission solely from electrons traveling along the LoS (green arrow in Figure 7).Hence, the emission zone essentially takes the form of a line.We split this line into radial segments to account  for an emission region within the companion zone in J2339 (Table 4) during the SUPC phase.The emission region, characterized by a linear profile, initiates at R IBS ≈ 0.3a orb , where electrons are energetic.In this phase, electrons traveling along the LoS approach the companion at a orb .Owing to the strong B, these electrons efficiently lose energy in the inner regions, and the cooling timescale increases with rp.The analysis involves 7000 radial bins to comprehensively characterize the emission region at this phase.
for varying conditions such as the companion's B, assumed to have a dipole structure: The strong dependence of B on r * implies a rapid decline in emission from these electrons with increasing r * .We consider this zone only in Scenario 2. In this scenario, we assume that the electrons injected by the pulsar follow a δ distribution with a very high Lorentz factor (γ p > ∼ 10 8 ): While most of these unshocked primary electrons in the wind zone transform into lower-energy electrons (secondaries), constituting the "wind zone" (Section 4.1.1),a small fraction (ζ) of the primaries passes through the IBS and emits in this companion zone.The energy budget of the pulsar is governed by Equations ( 17)-( 21).Due to our assumption of a very large value for γ p (≥ 10 8 ), the IC emission from these shock-penetrating electrons is highly suppressed, whereas their synchrotron emission under the strong B of the companion is very intense.The synchrotron cooling timescale (Equation ( 22)) of the electrons is significantly shorter than their residence time, especially near the SUPC phase where the electrons closely approach the companion (Figure 4 right).In the computation of the emission from this zone, we account for synchrotron cooling.This necessitates the use of different lengths for the radial segments (dr p ; Figure 10), ensuring that each segment is considerably shorter than ct cool (Equation 22): dr p = 10 −3 ct cool .
In the first segment (closest to the IBS), we inject electrons with a Lorentz factor of γ p and evolve them while accounting for their cooling (middle panel in Figure 10).Similar to the approach for the wind zone, the solid-angle fraction 1/4πd 2 is taken into account in the emission formula (Equation ( 42)).We integrate the emission up to a distance where the Lorentz factor of the electrons drops to 10 −3 γ p .If this length proves to be excessively long, we terminate the integration at l comp = 4a orb .We verified that the emission from electrons beyond the integration region is negligibly small.The emission SED and LC are primarily determined by ζ, γ p , and B c (Figure 6).The combination of γ p and B plays a crucial role in determining the frequency of the synchrotron SED peak (e.g., Equation (3)).Achieving an appropriate balance in Bγ 2 p ensures that the peak emission occurs in the GeV band, resulting in a high GeV flux.Excessively high values of γ p or B (e.g., near the SUPC phase) can push the SED peak beyond the GeV band, leading to a reduction in the ∼GeV flux and causing a dip in the LC at some phases (φ = 0.25 in Figure 6).On the contrary, if B is low (away from the SUPC phase), the peak emission occurs below the GeV band, resulting in a model-predicted GeV emission that is low (Figure 5).Moreover, electrons do not radiate efficiently when B is low.We optimize these parameters to explain the observed GeV fluxes and LCs (Figure 6).

Computation of the synchrotron and IC emissions
We calculate the synchrotron and IC emissions from electrons in each segment of the aforementioned zones using the formulas detailed in Finke et al. (2008) and Dermer et al. (2009).

The synchrotron emission
The synchrotron SED from isotropically-distributed electrons (in the flow rest frame) under the influence of a randomly-oriented B is given in Equation ( 18) of Finke et al. (2008): where primed quantities are defined in the flow rest frame.In this formula, δ D represents the Doppler factor of the bulk flow (e.g., Equation ( 38)), h is the Planck constant, and q e denotes the charge of an electron.The observed and emitted photon energies are expressed in units of m e c 2 as ǫ (≡ hν/m e c 2 ) and ǫ ′ (≡ hν ′ /m e c 2 ), respectively.
dN ′ e (γ ′ e ) dγ ′ e represents the energy distribution of the emitting electrons (Equations ( 35) and ( 40)), and where x ≡ 4πǫ ′ m 2 e c 3 3qeBhγ ′2 e , and K 5/3 is the modified Bessel function of the third kind.For computations of R(x), we use approximate formulas provided by Finke et al. (2008).We carry out the computation of Equation ( 42) using 700 bins for ǫ and 200 bins for γ e .
While the assumptions underlying Equation ( 42) are valid for the synchrotron emission from the IBS, it is important to note that electrons in the companion zone are not isotropically distributed, and B of the companion is not randomly oriented.Since the orientation of the companion's B is unknown (we show an aligned case in Figures 4 and 7 only for illustrative purposes), it is not possible to account for pitch angle scattering.In this study, due to the lack of information about the specific orientation, we take an average over the pitch angle, as reflected in Equation ( 43).The pitch angles experienced by shock-penetrating particles, and consequently, the amplitude and shape of the GeV LC, will be influenced by the orientation of the companion's B. This aspect could be explored in a future paper.

The IC emission
We use the IC emission formula given in Dermer et al. (2009).This formula for the emission SED incorporates various parameters, including Doppler boosting, anisotropic scattering, and the seed photon spectrum, and is expressed as (e.g., Equation (34) of Dermer et al. 2009): where γ e = δ D γ ′ e .As in the case of the synchrotron formula, the photon energies before (ǫ * ) and after (ǫ s ) an IC scattering are normalized by m e c 2 .The parameters φ * and µ * (= cosθ * ) define the direction of the seed photon, with the energy density of u * (ǫ * , Ω * ) in the emission zone.In this context, we consider only the companion's BB spectrum, characterized by R * and T * (Table 3), as the source of the seed photons.Ξ is defined by where y ≡ 1 − ǫs γe and ǭ represents the invariant collision energy ǭ ≡ γ e ǫ * (1 − 1 − 1/γ 2 e cosψ) with ψ denoting the scattering angle.
The integration limits in Equation ( 44) are determined by the scattering kinematics: We calculate the IC SED using 700 bins for ǫ s and ǫ * , and 200 bins for γ e .In this study, we assume that the companion is a point source, simplifying the φ * and θ * integrations in Equation ( 44).
In each emission region (segment) within the emission zones, as described in Section 4.1, we compute the necessary quantities for emissions, such as δ D , B(s), B(r * ), u * (ǫ * , Ω * ), and ψ, taking into account the orbit (e.g., direction of the LoS) and the geometry of the emission region (e.g., r * , r p and the flow direction).We combine these with the particle distribution dN/dγ e within the region to compute the emission SEDs at 50 orbital phases.Note that the energy distributions of the monoenergetic electrons in the wind (Equation ( 9)) and companion zones (Equation ( 40)) are given in the observer rest frame, and thus their emissions are not Doppler boosted (i.e., δ D = 1).Conversely, the distributions of electrons in the decelerated wind (Equation ( 25)) and IBS (Equation ( 35)) zones were expressed in the flow rest frame (these electrons have bulk flow speeds), and their emissions are Doppler boosted.The synchrotron X-rays exhibit orbital modulation due to this Doppler boosting, while the GeV modulation results from changes in either the IC scattering geometry (e.g., ψ; scenarios 1 and 1a) the cooling timescale determined by the companion's B(r * ) (Scenario 2).

Results of modeling
We initiated our model with the rough estimations provided in Section 3 as the starting point for the parameters.Subsequently, we fine-tuned these parameters  4. Due to the large number of parameters, a formal data fitting process was challenging.Instead, we employed a qualitative 'fit-by-eye' approach.Therefore, we do not report uncertainties in the estimated parameter values; we defer this to future work.
Figures 1 and 3 display the computed SEDs and LCs corresponding to the three cases described earlier (see Section 3.2): the basic scenario (red; Scenario 1), its modification (blue; Scenario 1a), and Scenario 2 (green).The predicted GeV fluxes and LCs exhibit significant variations among the scenarios.Similar to Scenario 1, the GeV LCs in Scenario 1a (Figure 3) exhibit orbital modulation due to the changing IC scattering geometry.However, they are slightly broader than the LCs predicted by Scenario 1 due to the isotropization and thermalization of the electrons.Scenario 1 is unable to explain the LAT data, whereas Scenario 2 can readily explain the multiband data.
The parameter values derived from the basic model (Scenario 1) are in accord with those obtained using the analytic approach (Section 3.2).To match the X-ray SED, the IBS electrons should follow a power-law distribution with an index of p 1 ≈1.3-1.6 between γ s,min and γ s,max .These two parameters were adjusted to explain the X-ray data and not to overpredict the LAT flux for the given η s value (Table 3) of each target.The highest-energy electrons in the IBS upscatter the stellar BB emission to TeV energies.Therefore, we relied on the IC emission from the wind zone to generate the peak of the IC SED at < ∼ GeV.As a result, the basic model (red curves in the top row of Figure 3) exhibits two peaks in the >MeV SED.The low-energy peak, centered around γ 2 w hν * ∼100 MeV, corresponds to the IC emission from the wind zone, while the high-energy peak, around γ 2 s,max hν * ∼TeV, corresponds to the IBS IC emission (see van der Merwe et al. 2020;Wadiasingh et al. 2022, for further discussion on this component).While the computed LC shapes of Scenario 1 resemble the observed ones, the predicted gamma-ray fluxes were found to be insufficient, as mentioned in Section 3.2.
The lack of GeV fluxes in this scenario was addressed by increasing the residence time of the electrons within the wind zone, achieved through flow deceleration (Scenario 1a).Despite our analytic investigation suggesting a very low value of v wind (Equation ( 16)), we sought a validation using our numerical model.The model with a constant speed profile resulted in wind speeds (500 km s −1 -5000 km s −1 ) higher than our rough estimations in Section 3.2.1,yet still too low to induce shock formation.Exploring different speed profiles, such as linear or exponential decreases in v wind with increasing distance from the light cylinder (toward the IBS), cannot provide a remedy as these alternatives would require even lower values of v wind at the position of the windwind interaction region (i.e., IBS).B energy in the wind zone may be substantial (Cortés & Sironi 2022), leading to a smaller η w than the assumed η w = η s (Table 4).In this case, further reductions in v wind compound the challenges associated with this scenario.
Scenario 2 can reproduce the multiband data with a companion's surface B c of > ∼ kG, which falls within the range suggested based on physical arguments (Wadiasingh et al. 2018;Conrad-Burton et al. 2023).The gamma-ray SEDs predicted by this scenario are broad due to the effects of energy loss and variation in B. The predicted GeV LCs exhibit a flat top with a rapid decline at the phases corresponding to the IBS tangent (see Figure 5 and Section 4.1).The parameter values for the IBS are very similar to those used for Scenario 1, as ζ is not very high.Due to the loss of particles, the IBS synchrotron emission (X-ray) was slightly lower.This was addressed by using a stronger B s (Table 4), but it could also be achieved by a smaller γ w equally well.
It is important to note that the model parameters are highly covariant, meaning that the values presented in Table 4 are not unique (see Section 5).In this work, we used the maximum possible values for η p (i.e., = 1−η γ ) and η w for Scenario 2. A lower value of η p (e.g., see Cortés & Sironi 2022) would weaken the GeV emis- in Figure 3.The orange and green curves correspond to models with ηp =0.2 and 0.5, respectively.In this example, the decrease in ηp (compared to our baseline value in Table 4) was compensated for by an increase in ζ.
sion, and this could be compensated for by increasing ζ and/or B c (Figure 11).

DISCUSSION AND SUMMARY
We analyzed the X-ray observations of three RB pulsars from which orbitally-modulating GeV signals were detected.We then constructed their multiband SEDs and LCs and investigated potential scenarios for the gamma-ray modulations using a phenomenological IBS model.
Based on our modeling, we found that Scenario 1 is unable to explain the measured GeV fluxes of the RB targets, as previously noted (An et al. 2020;Clark et al. 2021).It is worth noting that our computations might underestimate the GeV flux since electrons in regions with higher inclinations, such as those near the orbital plane, can see stronger BB emission (from the heated surface of the companion) than what we assumed (i.e., observed at an inclination angle i < 90 • ).On average, however, electrons spread over extended emission zones (Figure 4) do not preferentially encounter the most intense photon field, and so the increase in the GeV flux due to this effect would be modest.A further increase in the IC flux can be achieved if the distances to the sources were smaller and the pulsar's ĖSD is larger (Equation ( 10)).The latter may be possible since neutron stars in pulsar binaries may be more massive (e.g., > ∼ 2M ⊙ ; Schroeder & Halpern 2014;Romani et al. 2022) than 1.4M ⊙ used for the ĖSD estimation.These increases are not very large and thus would be still insufficient to explain the orders-of-magnitude discrepancy in the GeV band (red in Figure 3).As demonstrated in the previous sections, addressing this issue involves a bulk deceleration of the unshocked wind to a low speed.However, this approach lacks physical support.If the wind zone does not account for gamma-ray production, there must be an alternative physical process capable of retaining the γ w ≈ 10 4 electrons within the system for an extended period for this scenario to be plausible.
Alternatively, the IBS model constructed based on Scenario 2 could easily accommodate the broadband SEDs and multiband LCs of the three RBs (Figure 3).However, we should note that the parameter values reported in Table 4 are not unique due to parameter degeneracies (see Kim et al. 2022), and so MeV and TeV flux predictions here should be taken as one potential realization among a landscape of possibilities.Some of the degeneracy can be broken by high-quality optical, X-ray, MeV and TeV data or limits.i and η s values can be inferred from optical LC modeling (e.g., Breton et al. 2013) and the LAT measurement of the pulsar flux, respectively.The X-ray data help constrain v IBS (by widths and amplitudes of the LC peaks; Kim et al. 2022), B s within the IBS (Equation ( 5)), and γ s,max (synchrotron cut-off).Current constraints on these parameters are not stringent because the quality of the X-ray LC measurements is poor and the SEDs do not fully cover the energy range of the synchrotron spectrum.For the B s values in Table 4, our model predicts a synchrotron spectral cut-off at ≥100 keV (Figure 3), which can be confirmed by future hard X-ray and/or soft gamma-ray observations.
The observed GeV LCs appear broader than those predicted by Scenario 2 (green in Figure 3 bottom).This might be caused by the observational effect that the probability-weighted LAT LCs show reduced modulation compared to the sources' actual variability (e.g., Corbet et al. 2016), meaning that the intrinsic GeV LCs of the targets may be narrower than the observed ones.Moreover, the LC shape may also alter depending on the flat levels (constant emission) that are subtracted from the LCs.In addition, Scenario 2 predicts broader SEDs in the GeV band, compared, for example, with the J2039 SED obtained from differencing LAT SEDs of the orbital maximum and minimum intervals (Clark et al. 2021).A more direct measurement of the SED of the orbitally-modulating GeV emission, facilitated by the LAT and/or other GeV observatories (e.g., AMEGO-X; Fleischhack & Amego X Team 2022), will contribute to advancing our studies based on Scenario 2.
It was suggested that the 'pulsed flux' is orbitally modulated in J2039 and J2339 (Clark et al. 2021;An et al. 2020), implying that the pulsar's pulsations should be preserved in the emission regions of the modulating gamma-ray signals.To maintain the pulse structure, the flow in the emission zone should have a stripedwind structure, and the emission timescale (cooling timescale) should be considerably shorter than the pulse periods of the pulsars.The latter condition is satisfied by Scenario 2, in which the cooling timescale is shorter than a millisecond (Equation ( 22)).If the flow structure in the companion zone preserves the striped structure, Scenario 2 can potentially explain the orbital modulation of the pulsed signals as well.However, the striped wind might be destroyed in the wind zone or in the IBS; this requires further theoretical studies.
As demonstrated earlier, Scenario 2 exhibits favorable aspects in explaining the observed data for the RB targets.While further theoretical studies with physically-motivated models for pulsars and their winds (e.g., Sironi & Spitkovsky 2011;Cortés & Sironi 2022;Kalapotharakos et al. 2023) and high-quality measurements are needed to confirm the scenario, it is almost certain that the unshocked wind plays an important role in pulsar binaries.Scenario 2 requires high-energy particles attaining ∼PeV, close to the voltage drop available to millisecond pulsars.This would mean that pulsar magnetospheres are very efficient particle accelerators, with implications for pulsed (magnetospheric) TeV emission from pulsars (Harding et al. 2021;H. E. S. S. Collaboration et al. 2023).This is also in accord with some acceleration and emission scenarios for pulsed GeV emission from pulsars (i.e. the curvature radiation scenario).
An important factor that we did not consider in this study is the potential synchrotron emission which may originate from the 'current sheets' within the striped wind.In regions near the pulsar's light cylinder where B is expected to be strong, the synchrotron emission would fall within the GeV band and could add to the pulsed flux of the pulsar (Pétri 2012).On the other hand, in regions near the IBS, we speculate that B is low (e.g., ≈ B s /3 < 1 G; Kennel & Coroniti 1984) due to the dissipation of magnetic energy in the current sheets (Sironi & Spitkovsky 2011).In this case, the synchrotron emission frequency would be ≤ 1 eV (Equation (3)).An important question is "how much of the magnetic energy is converted to particles and to radiation in the wind zone?"This may have a significant impact on the structure of and emission from the pulsar wind.Accurate measurements of the broadband SEDs and LCs of pulsar binaries, and further PIC simulations for pulsar winds (including radiation) are warranted.
We adopted the simplicity of the IBS shape for a two isotropic wind interaction (Section 4.2.1;Canto et al. 1996).The single parameter (β) allows for considering a range of opening angles, but it may not capture the complexities of anisotropic winds or wind-magnetosphere interactions.In these more intricate cases, the IBS shape can exhibit greater complexity, influenced by system parameters such as the pulsar's spin axis orientation or the companion's magnetic field (Wadiasingh et al. 2018;Kandel et al. 2019).While the isotropic wind model reproduces the general features of the IBS in the other cases, subtle discrepancies may arise in the predicted LCs.The wind-B interaction model offers the additional benefit of determining the B orientation encountered by shock-penetrating electrons.Further investigations based on the wind-B interaction model hold the potential to unlock new insights into the gamma-ray emission mechanisms of RBs.
Although our simplified IBS model with phenomenological prescriptions for the wind flows may not encompass all the important physics of the flows, we demonstrated that Scenario 2 offers a plausible explanation for the orbitally-modulating GeV signals detected from a few RBs.This scenario, along with our IBS model, can be further tested with high-quality X-ray and gamma-ray data, and provide an opportunity to comprehend the energy conversion and particle flow within the pulsar wind zone.Future X-ray and gamma-ray observatories such as AXIS (Reynolds et al. 2023), HEX-P (Madsen et al. 2023), and AMEGO-X (Fleischhack & Amego X Team 2022), have potential to furnish high-quality data for more samples, thereby facilitating better understanding of the IBS and pulsar wind.
This work used data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA.We made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA).Z.W. thanks Jeremy Hare for interesting discussions.Z.W. acknowledges support by NASA under award number 80GSFC21M0002.This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean Governmaent (MSIT) (NRF-2023R1A2C1002718).This paper employs a list of Chandra datasets, obtained by the Chandra X-ray Observatory, contained in DOI:10.25574.We thank the referee for the insightful comments that helped improve the clarity of the paper.XSPEC (v12.13.0c;Arnaud 1996)

Figure 1 .
Figure1.Exposure-corrected LCs of J1227 (left), J2039 (middle), and J2339 (right).The left and right ordinates show the count rates and fluxes, respectively, where the latter were estimated by comparing the phase-averaged fluxes to the observed counts.We used the 2-10 keV band for the Chandra and XMM LCs, and the 3-20 keV band for the NuSTAR LCs.The solid curves are LCs computed with our IBS model for Scenarios 1 and 2 (see Sections 3.2 and 4.3).

Figure 2 .
Figure 2. X-ray spectra of J1227 (left), J2039 (middle), and J2339 (right).The solid lines represent the best-fit XSPEC models, and the bottom panels display the fit residuals.

Figure 4 .
Figure 4. Cartoons that depict emission scenarios employed in this work.(Left): Emission components common to Scenarios 1 and 2. Purple, red, and blue arrows denote synchrotron, IC-in-IBS, and IC-in-wind emissions, respectively.Note that the synchrotron emission is stronger in the flow direction due to relativistic aberration.For each emission component, three representative directions (toward INFC, SUPC, and in between) are displayed, where thicker arrows mean stronger emission.(Right): Synchrotron radiation under the companion's B (green region) by the wind particles that penetrate the IBS.For Scenario 2, we add this emission component to those of Scenario 1.
Scenario 2: It was suggested that a sufficiently energetic component of the upstream pairs passes through the IBS unaffected and emits synchrotron radiation under the influence of the companion's B (van der Merwe et al. 2020), producing GeV gamma rays (Figure 4 right).This component is added to the gamma-ray flux of Scenario 1.

Figure 5 .
Figure 5. Gamma-ray SEDs of J2339 at several phases computed by our numerical model (Section 4) based on Scenario 2 (Section 4).

3. 2 . 2 .
Scenario 2: synchrotron radiation under the companion's B Another way to increase the gamma-ray flux of RBs is to utilize the efficient synchrotron process, as was suggested by van der Merwe et al. (2020); some of the upstream electrons may pass through the IBS and emit gamma rays under the strong B of the companion (Figure 4 right).If the companion has a surface B of approximately kG(Sanchez & Romani 2017;Wadiasingh et al. 2018), electrons with a Lorentz factor γ e ≈ 10 7 can produce synchrotron photons at ∼GeV energies (Equation (3)).

Figure 6 .
Figure 6.GeV LCs generated through our numerical model (Section 4) based on Scenario 2, employing optimized parameters specific to J2339 (Table 4).Changes in the LCs are attributed to different values of γp (left), Bc (middle), and ζ (right), reflecting the diverse parameter space under consideration.

Figure 7 .
Figure7.Illustration showing emission zones in the vertical cross section of a system at the SUPC phase.The companion star and pulsar are symbolized by large red and black circles, respectively (not to scale).The thick pink curve displays the IBS.The wind zone, situated between the pulsar's light cylinder and the IBS (depicted in blue), contrasts with the green companion zone.The blue and green dashed arrows represent the regions within the wind and companion zones, respectively, where electrons travel parallel to the LoS.Only these electrons contribute to the observable emission due to strong Lorentz boosting.The pulsar-centric coordinates (rp, θp, and φp) describe emission regions, while companion-centric coordinates (r * , θ * , and φ * ) detail seed photon density, direction and the companion's B. R 0 denotes the distance between the pulsar and the shock nose, and s is the distance along the IBS from the shock nose.

Figure 8 .
Figure 8. IC SEDs resulting from electrons interacting with the companion's BB radiation (Equation (44)).These emissions originate specifically from the wind zone during the SUPC phase.The displayed SEDs represent a range of values for γw.All other parameters are held fixed at their optimized values for J2339 (Table 4).

Figure 9 .
Figure 9. X-ray LCs illustrating emissions from the IBS zone, computed with a range of values for β (top) and Γ D (bottom).The parameters for J2339 reported in Table4serve as the baseline in this example.

Table 2 .
Spectral fit results

Table 4 .
Parameters used for models displayed in Figures 1 and 3 to match the observed SEDs and LCs of the RB targets with the model computations.The optimized parameter values are outlined in Table