OPEN ACCESS

<figure><img src="image_1.png" /></figure>

CrossMark

# Orbital Phase-resolved Analysis of X-Ray and Gamma-Ray Observations of the High-mass Gamma-Ray Binary 4FGL J1405.1–6119

Alexander Lange¹,²,³ ID, Robin H.D. Corbet²,³,⁴ ID, Joel B. Coley⁴,⁵,⁶ ID, Guillaume Dubus⁷ ID, Jeremy Hare²,³,⁸ ID,
Nazma Islam²,³ ID, and Jonathan Barnes⁵

¹Department of Physics, The George Washington University, 725 21st Street NW, Washington, DC 20052, USA

²X-ray Astrophysics Laboratory, Code 662, NASA Goddard Space Flight Center, Greenbelt Road, MD 20771, USA

³Center for Space Sciences and Technology, University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA

⁴CRESST II, USA

⁵Department of Physics and Astronomy, Howard University, Washington, DC 20059, USA

⁶Astroparticle Physics Laboratory, Code 661, NASA Goddard Space Flight Center, Greenbelt Road, MD 20771, USA

⁷Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

⁸The Catholic University of America, 620 Michigan Avenue NE, Washington, DC 20064, USA

Received 2025 January 21; revised 2025 May 13; accepted 2025 May 17; published 2025 June 17

## Abstract

We present the results of multiwavelength observations of the high-mass gamma-ray binary (HMGB) 4FGL J1405.1–6119. A pair of joint XMM-Newton and NuSTAR observations taken in 2019 (sampling the gamma-ray maximum and X-ray maximum) characterize the emission of soft and hard X-rays. We find variability of the hydrogen column density along our line of sight, $N_H$, and photon index, $\Gamma$, and find no evidence of pulsations in X-rays. We also refine a new best-fit orbital period to $P = 13.7157 \pm 0.0014$ days, the first orbital phase-resolved analysis based on nearly 16 yr of Fermi-LAT observations of 4FGL J1405.1–6119 and the evolution of the spectral shape as a function of orbital phase. Finally, the X-ray and gamma-ray spectra for the phases sampled in the new X-ray observations can be interpreted in the framework of the intrabinary shock model, previously applied to HMGBs such as LS 5039.

*Unified Astronomy Thesaurus concepts: Neutron stars (1108); X-ray binary stars (1811); High mass x-ray binary stars (733)*

## 1. Introduction

High-mass gamma-ray binaries (HMGBs) are a rare class of gamma-ray sources in the Galaxy, consisting of a compact object, most likely a neutron star, but potentially a black hole, and a massive companion star (typically, O or B/Be type), that produce a characteristic spectral peak above ~1 MeV. The emission mechanisms of the X-rays and gamma rays from the small population of known HMGBs (≈10) are debated, as well as the nature of the compact objects (G. Dubus et al. 2017; M. Chernyakova & D. Malyshev 2020). As such, HMGBs provide excellent opportunities to test models of high-energy emission and particle acceleration. Currently, there are two prevalent models for the high-energy emission: an intrabinary shock (IBS) model with a neutron star as the compact object, and a jet microquasar model with a neutron star or black hole (for reviews, see I. F. Mirabel 2012; G. Dubus 2013).

The IBS model proposes that the gamma-ray emission is produced by the interactions between the relativistic outflow of a rotation-powered pulsar (i.e., pulsar wind) and the wind of its companion star via the inverse Compton scattering of seed photons by the accelerated charged particles. In contrast, the microquasar model describes the production of gamma rays from particles accelerated in a relativistic jet. The formation of this jet is related to accretion onto the compact object (neutron star or black hole).

Among HMGBs, the search for (and subsequent detection of) radio pulsations has identified the compact objects of PSR B1259–63 (S. Johnston et al. 1992), PSR J2032+4127 (F. Camilo et al. 2009), and LS I +61°303 (S.-S. Weng et al. 2022) as pulsars, with respective spin periods of 47, 143, and 269 ms, supporting the IBS model. However, the companion's stellar wind can suppress radio emission from a putative pulsar due to free-free absorption such as shown by S.-S. Weng et al. (2022), making it hard to confirm the nature of the other HMGBs' compact objects with only radio-based observations.

The gamma rays produced may occur in tandem with multiwavelength flaring and can correspond to certain orbital phases (e.g., X-ray/GeV/TeV, PSR B1259-63; or radio/GeV, Cygnus X-3; G. Piano et al. 2012; M. Chernyakova et al. 2015, 2020), but these correlations are not consistent throughout all HMGBs. Beyond the high-energy (HE) or very high-energy (VHE) gamma-ray regime, HMGBs exhibit modulated emission of X-rays. This places limits on the broadband spectral energy distribution (SED), directly probes the relativistic lepton population, and allows for physical modeling of the IBS and high-energy emission (G. Dubus 2013). Observed X-ray modulation can be attributed to the changes in geometry of the binary and of the IBS (e.g., PSR B1259–63; M. Chernyakova et al. 2020). The modulation of emission and absorption lines of the companion's spectra can further constrain the orbital solution (see, e.g., 1FGL J1018.6–5856; The Fermi LAT Collaboration 2012; J. Strader et al. 2015; B. van Soelen et al. 2022). X-ray spectral fitting can reveal spectral breaks and provide estimates of the magnetic field strength at the source of the high-energy emission. Lastly, multiwavelength observations provide an opportunity to test these emission models and

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
determine the compact object's nature (e.g., HESS J0632+057; J. Kim et al. 2022).

The HMGB 4FGL J1405.1–6119 (hereafter referred to as J1405) was discovered by R. H. D. Corbet et al. (2019) in a blind periodicity search of Fermi-LAT fluxes with an orbital period of 13.7157 ± 0.0014 days. Fermi-LAT probability-weighted aperture photometry data showed an orbital modulation of the gamma rays with two maxima per orbital cycle. From the maximum gamma-ray flux and first peak (defined at $\phi = 0.0$), there is a sharp drop in flux of the binary, followed by a small increase to the second peak ($\phi \approx 0.4$) with a slow decay to a minimum. Both peaks are roughly one day in width (corresponding to ~1/14 in phase), where the peak at $\phi = 0.0$ has an observed flux of roughly twice that of the peak at $\phi \approx 0.4$. There was no spectral analysis beyond fitting the normalization with spectral parameters fixed to the values in the 4FGL catalog (S. Abdollahi et al. 2020), from which J1405 was found to have a spectral peak around 500 MeV.

Follow-up radio measurements from the Australia Telescope Compact Array at 5.5 and 9 GHz and soft (0.2–10 keV) X-ray observations from the X-ray telescope aboard the Neil Gehrels Swift Observatory (Swift-XRT) show a coincident source with modulation on the same orbital period. The radio and X-ray observations feature a single peak at $\phi \approx 0.3$ that is slightly offset from the second gamma-ray peak at $\phi \approx 0.59$. It should be noted that the modulation of the 5.5 GHz flux is weaker than at 9 GHz. There is no X-ray peak at either of the gamma-ray peaks. Lastly, near-infrared spectral measurements by FLAMINGOS-2 classify J1405's companion star as an O6.5 III star (R. H. D. Corbet et al. 2019). The Galactic column density of neutral hydrogen and the photon index were reported as $N_H = 6.9 \times 10^{23} \text{cm}^{-2}$ and $\Gamma_X = 1.5 \pm 0.4$ (R. H. D. Corbet et al. 2019), respectively.

We report on the results of the follow-up XMM-Newton and NuSTAR observations, and provide an update to the Fermi-LAT analysis to include a phase-resolved analysis of 16 yr of data. We fit the X-ray and gamma-ray data with the IBS model and report the results. The remainder of this paper is structured as follows: the details of each observation and their reduction are reported in Section 2, and the spectral and timing analyses of the X-ray data are reported in Section 3. Section 4 provides the results of the timing and orbital phase-resolved spectral analysis of ~16 yr of Fermi-LAT data. Section 5 provides the details of the theoretical models used and the modeling of the phase-resolved broadband SEDs. The implications of these results and how they compare to previous studies are discussed in Section 6, followed by the conclusion in Section 7.

## 2. X-Ray and Gamma-Ray Observations

Joint XMM-Newton/NuSTAR observations of J1405 were conducted around the Chandra source found at the position (R.A., decl.) = (14°05′m14.43, −61°18′ 28′′36) in 2019 August to capture the X-ray behavior at the gamma-ray maximum and minimum as reported by R. H. D. Corbet et al. (2019). We hereafter refer to the observations performed at the gamma-ray maximum and minimum (where the X-ray flux is at a maximum) as “Gamma-ray Maximum” (0.964 < $\phi_{\text{Gamma-ray Max.}}$ < 0.0357) and “X-ray Maximum” (0.393 < $\phi_{\text{X-ray Max.}}$ < 0.464). All analyses use a refined binary orbital period of $P = 13.7157$ days (see Section 4). A reconstructed X-ray phase-folded lightcurve is shown in Figure 1 to compare to Figures 3 and 4 from R. H. D. Corbet et al. (2019). The details of the 2019 NuSTAR

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Gamma-ray and X-ray lightcurves folded on the orbital period from R. H. D. Corbet et al. (2019). Top: gamma-ray orbital phase-folded lightcurve from Fermi-LAT. Data are taken from the probability-weighted aperture photometry data. Bottom: X-ray orbital phase-resolved lightcurves from Swift-XRT, NuSTAR, and XMM-Newton. Swift archival fluxes are shown in black. NuSTAR and XMM-Newton observations are shown as blue circles and squares, respectively. A twin-peaked structure is shown at the Gamma-ray Maximum and the X-ray Maximum in the Fermi-LAT lightcurve. These phases are the focus of the 2019 NuSTAR/XMM-Newton observations.</figcaption></figure>

and XMM-Newton observations are provided in Table 1. Fermi-LAT photon events spanning 16 yr are used in the Fermi-LAT timing and spectral analyses.

### 2.1. XMM-Newton

The XMM-Newton observatory (F. Jansen et al. 2001) consists of three sets of coaligned X-ray telescopes: two EPIC MOS CCDs and an EPIC PN CCD. EPIC MOS provides coverage between 0.2 and 10.0 keV, while EPIC PN extends the range to 0.15–15 keV (L. Strüder et al. 2001; M. J. L. Turner et al. 2001).

Two XMM-Newton observations (50 and 34 ks, respectively) of J1405 were taken in 2019 August. EPIC MOS1/MOS2 cameras observed J1405 in full window mode, while the EPIC PN camera observed in small window mode. EPIC MOS has a timing resolution of 0.3 s and is not used for timing analysis. Small window mode allows for a PN timing resolution of 5.7 ms (L. Strüder et al. 2001) and is used in Section 3.2.
<table><caption>**Table 1**
4FGL J1405.1–6119 NuSTAR and XMM-Newton Observation IDs, Exposures, and Count Rates</caption><thead><tr><th>Telescope</th><th>Phase Range<br>(φ)</th><th>Start Date<br>(UTC)</th><th>Observation ID</th><th>Exposure<br>(ks)</th><th>Count Rate<br>(×10<sup>-2</sup> counts s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>NuSTAR</td><td>0.91–0.05</td><td>2019-08-24 08:24:56</td><td>30502015004</td><td>86.8</td><td>1.04 ± 0.06</td></tr><tr><td>XMM-Newton</td><td>0.93–0.96</td><td>2019-08-24 15:01:01</td><td>0852020201</td><td>49.8</td><td>4.00 ± 0.17</td></tr><tr><td>NuSTAR</td><td>0.35–0.46</td><td>2019-08-16 16:52:41</td><td>30502015002</td><td>61.7</td><td>2.12 ± 0.09</td></tr><tr><td>XMM-Newton</td><td>0.42–0.45</td><td>2019-08-17 16:34:15</td><td>0852020101</td><td>33.9</td><td>3.86 ± 0.32</td></tr></tbody></table>

The XMM-Newton observations were reduced using the standard Scientific Analysis System (SAS) tools version 20.0.0. Both MOS and PN observations are filtered to remove times of high particle background flaring. This filtering reduces the exposures for the Gamma-ray Maximum and X-ray Maximum to 49.8 and 33.9 ks. To correct for the difference in effective area between EPIC and NuSTAR and to align the PN spectra to the NuSTAR spectra, we used the `applyabsfluxcorr=yes` keyword in `arfgen` as per the latest recommendation by the XMM-SAS team.<fnref n="9" />

Circular source and background regions were selected from the same chip around J1405 and a source-free background region. The source and background regions are of radius $r = 25''$ for PN and $r = 35''$ for MOS1/MOS2. Events filtered from source regions were barycentered, using `barycen`, and source and background spectra are extracted from their respective regions. Spectral data were grouped into a minimum of 25 counts per bin.

## 2.2. NuSTAR

The Nuclear Spectroscopic Telescope Array (NuSTAR; F. A. Harrison et al. 2013) observatory uses its two modules, FPMA and FPMB, for broad spectral coverage from 3 to 79 keV and a timing resolution of 65 µs (M. Bachetti et al. 2021).

NuSTAR performed two observations of J1405 (87 and 62 ks) contemporaneously with XMM-Newton (see Table 1). The NuSTAR observations were reduced using HEASoft version 6.31.1 (Heasarc 2014). We used the standard NuSTARDAS pipeline, `nupipeline` and `nuproducts` with keywords `saacalc=2` and `tentacle=Yes` to account for the South Atlantic Anomaly and any space weather, and generate science products, reducing scientific exposure to 74 and 51 ks, respectively.

Source and background regions were chosen from an $r = 25''$ circular region and an $r = 50''$ circular region around J1405 and around a region away from J1405, respectively, on the same detector for both modules. Variations between the source and background regions are corrected for using `nubackscale`. We correct the photon arrival times for barycentric motion and apply the clock correction file `"nuCclock20100101v125.fits"` from the standard CALDB data set to search for pulsations at high frequencies detailed in Section 3.2. The spectral data are binned to contain a minimum of 25 counts per energy bin.

## 2.3. Fermi-LAT

The Fermi Large Area Telescope (LAT) instrument is sensitive to gamma rays at energies from 20 MeV to >300 GeV (W. B. Atwood et al. 2009; W. Atwood et al. 2013; S. Abdollahi et al. 2020). We include all events within a 10° circular region of interest (ROI) from J1405 (R.A., decl.) = (211°2982, -61°3321), ranging from MJD 54682 to MJD 60492 (from 2008 August 4 to 2024 July 1). We chose the energy range between 200 MeV and 500 GeV to avoid the large point-spread function (PSF) associated with Fermi-LAT's lower energies and to match the original analysis of R. H. D. Corbet et al. (2019). The events were filtered with standard data quality selections `"(DATA QUAL>0)&&(LAT_CONFIG==1),"` and excluding zenith angles above 90°. We chose to include both front and back converting events using `"evtype=3"` and `"evclass=128."`

The Fermi-LAT events were processed using `Fermi-Tools` (version 2.2, W. B. Atwood et al. 2009) and the FermiPy analysis package (version 1.2.0, M. Wood et al. 2017) for Python. All Fermi-LAT analyses use Pass 8 data (W. Atwood et al. 2013; P. Bruel et al. 2018), and we adopt the 4FGL-DR4 catalog (S. Abdollahi et al. 2022) source model and include all sources within a larger 15° ROI to account for Fermi-LAT's PSF. We derive a new orbital period of $P = 13.7157 \pm 0.0014$ days (see Section 4) according to the procedure listed by R. H. D. Corbet et al. (2019). This new period is in agreement with the original of $13.7135 \pm 0.0019$ days, refined with more than 5 yr of additional photon events. Close to 16 yr of photon events are folded into 14 phase bins according to this updated orbital period using the same reference epoch, $t_0$, MJD 56498.7. Each phase selection was further divided into 16 logarithmically spaced energy bins (while taking account of energy dispersion), for spectral analysis, using the recommended FermiTools prescription.<fnref n="10" />

## 3. X-Ray Analysis

Due to high background, all X-ray data are heavily dominated by background at soft and hard energies. At 20 keV, NuSTAR FPMA/FPMB events are dominated by background, and so they are filtered to include only events between 3 and 20 keV. XMM-Newton data have evidence of increased count rates due to soft-proton flaring (M. C. Weisskopf et al. 2002; D. Lo & J. Srour 2003; V. Fioretti et al. 2024), as well as heavy background at the soft and hard energies, and so the events are filtered to include events between 1 and 9 keV. These energy ranges are used for all subsequent analyses.

<footnote n="9">⁹ XMM-SAS Calibration note: https://xmmweb.esac.esa.int/docs/documents/CAL-TN-0230-1-3.pdf.</footnote>

<footnote n="10">¹⁰ Further Fermi-LAT data documentation and analysis details can be found at https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.</footnote>
<table><caption>**Table 2**
X-Ray Spectral Parameters</caption><thead><tr><th>Model Parameter</th><th>Gamma-Ray Maximum</th><th>X-Ray Maximum</th></tr></thead><tbody><tr><td>Flux<sup>a</sup></td><td>0.57<sup>+0.02</sup><sub>-0.02</sub></td><td>1.11<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>C<sub>FPMA</sub><sup>b</sup></td><td>1</td><td>1</td></tr><tr><td>C<sub>FPMB</sub><sup>b</sup></td><td>0.99 ± 0.01</td><td>0.99 ± 0.01</td></tr><tr><td>C<sub>PN</sub><sup>b</sup></td><td>0.90 ± 0.07</td><td>0.74 ± 0.05</td></tr><tr><td>C<sub>MOS</sub><sup>b</sup></td><td>1.02 ± 0.08</td><td>0.80 ± 0.06</td></tr><tr><td><i>tbabs</i> N<sub>H</sub></td><td>7.98 ± 0.63</td><td>5.69 ± 0.53</td></tr><tr><td>Γ<sub>X</sub></td><td>1.57 ± 0.10</td><td>1.42 ± 0.08</td></tr><tr><td>Normalization<sup>c</sup></td><td>1.12 ± 0.22</td><td>1.58 ± 0.27</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>137.03/133</td><td>151.33/132</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Fluxes (×10<sup>-12</sup> erg cm<sup>-2</sup> s<sup>-1</sup>) integrated from 1 to 20 keV.

<sup>b</sup> Cross-instrument normalization. C<sub>FPMA</sub> is fixed to 1 at both the Gamma-ray Maximum and the X-ray Maximum.

<sup>c</sup> N<sub>H</sub> measured in units of 10<sup>22</sup> cm<sup>-2</sup>.

<sup>d</sup> Model normalization (×10<sup>-4</sup> photons keV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>) at 1 keV.

## 3.1. Spectral Analysis

Spectral fitting was done using `XSPEC` (version 12.13.0c; K. A. Arnaud 1996). We fit the observations at the Gamma-ray Maximum and the X-ray Maximum with a power law (PL) and an intrinsic Galactic absorption component (`tbabs`, using photoionization cross sections from D. A. Verner et al. 1996 and `wilm` abundances from J. Wilms et al. 2000). The spectra from NuSTAR FPMA/FPMB and XMM-Newton MOS/PN were fitted simultaneously and tied with a constant in the model to account for the cross-normalization between the different instruments, where the cross-calibration constants were normalized to C<sub>FPMA</sub>. Fluxes from each instrument were calculated by including the multiplicative model component `cflux`. The best-fit results are shown in Table 2 along with chi-squared values and degrees of freedom (dof). Spectral fits for the PL model and chi-squared residuals are shown in Figure 2.

A comparison of joint spectral fits between the Gamma-ray Maximum and the X-ray Maximum shows a significant improvement when N<sub>H</sub> for each phased observation is fit independently (χ<sup>2</sup>/dof = 289.79/266), as opposed to when N<sub>H</sub> is tied across phases (χ<sup>2</sup>/dof = 297.08/267).

Due to XMM-Newton's sensitivity in the soft X-ray range, we can investigate line features. A cyclotron resonant scattering feature (CRSF) was reported in the Gamma-ray Maximum and the X-ray Maximum XMM-Newton observations by S.-C. Chiu & K.-L. Li (2024), and we tested for the presence of this feature. In addition to a PL spectral model, we added a Gaussian absorption feature at 2 keV and simulated 10,000 spectra from our best-fit values using the `XSPEC` tool `simftest` to compare the χ<sup>2</sup> fit statistics. We find that adding the Gaussian improves the fit in less than ~68% of trials. Lastly, we do find evidence of a spectral break at ~3.72 ± 0.37 keV during the Gamma-ray Maximum (χ<sup>2</sup>/dof = 134.14/131), but this does not significantly improve upon the PL fit (χ<sup>2</sup>/dof = 137.03/133) and is not considered in any subsequent analysis.

## 3.2. X-Ray Timing Analysis

We fold the observation-integrated (for easy comparison to Swift-XRT and Fermi-LAT fluxes) XMM-Newton and NuSTAR fluxes on the binary's updated orbital period and overlay the 3–10 keV fluxes (calculated using `cflux`) onto

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Joint SED fits of NuSTAR and XMM-Newton observations of J1405. The fits for each data set are shown as solid lines overlaid by each instrument's data. Note that we combined MOS1 and MOS2 for better statistics. NuSTAR FPMA/FPMB are shown in blue and gold, and XMM's combined MOS and PN are shown in green and red, respectively. The top panel shows (a) the Gamma-ray Maximum data and the best-fit model, and (b) its χ-squared contribution to the fit statistic. The bottom panel shows the same (best-fit model and χ<sup>2</sup>) for the X-ray Maximum observation.</figcaption></figure>

archival Swift-XRT orbital phase-folded fluxes. We recalculate the XRT fluxes to ensure each observation had an associated 3σ flux error. The fluxes were extracted using the Swift-XRT data products generator `webtool` (P. A. Evans et al. 2007, 2009) in Figure 1.

We do not fit the X-ray lightcurves and refer to E. A. Saav-edra et al. (2023), who find agreement with the initial results of R. H. D. Corbet et al. (2019) and long-term variability in the X-ray lightcurve throughout the XMM-Newton and NuSTAR observations. Instead, we focus on pulsation searches within the XMM-Newton and NuSTAR observations and modeling of the average flux state.

We searched for pulsations using the Z<sub>n</sub><sup>2</sup> test (R. Buccheri et al. 1983) for the XMM-PN and NuSTAR FPMA/FPMB event data; we elected to use n = 1 and n = 2 given the low signal-to-noise ratio for J1405. Barycentered NuSTAR FPMA/FPMB events were combined and analyzed separately from XMM-PN events. XMM-Newton MOS1/MOS2 were not considered due to their low time resolution. The Nyquist frequencies for XMM-Newton PN and NuSTAR are 87.72 and 7692.31 Hz. We thus chose to search frequency ranges from 5 × 10<sup>-6</sup> Hz (roughly
<figure><img src="image_6.png" /><figcaption>**Figure 3.** An example of two $Z_2^2$ analyses of both XMM-Newton (a), Gamma-ray Maximum) and NuSTAR (b), X-ray Maximum). Only very low-frequency ($1 \times 10^{-4}$ to $1 \times 10^{-5}$ Hz) signals in panels (a) and (b) are significant besides the orbital period (and harmonics) of NuSTAR (red) in panel (b). Panel (a) shows features near 10 Hz, but they fail to surpass a false-alarm probability threshold (shown by the blue horizontal line at the $3\sigma$ level).</figcaption></figure>

twice the observation length) to $1 \times 10^3$ Hz for NuSTAR and from $5 \times 10^{-5}$ Hz to 87 Hz for XMM-Newton PN. We oversample the $Z_n^2$ tests by a factor of 5. Example $Z_2^2$ periodicity searches can be found in Figures 3(a) and (b). We detect no signal (i.e., $>3\sigma$), but features at the $1\sigma-2\sigma$ significance level were found throughout the $Z_{n=1,2}^2$ analyses. Low-frequency signals near $\nu = 1 \times 10^{-5}$ Hz and below are attributed to red noise. We do provide upper limits on pulse fractions of $\approx 16\%$ and $\approx 17\%$ for the Gamma-ray Maximum and the X-ray Maximum, respectively, each calculated from the XMM-Newton $Z_2^2$ analyses.

# 4. Gamma-Ray Analysis

## 4.1. Timing Analysis

Fermi-LAT events from 200 MeV to 500 GeV are selected from within a $3^\circ$ aperture of J1405 and were probability-

weighted using an updated source model from the 4FGL-DR4 source catalog (see Section 4.2). Weights were assigned with the fixed spectral models and the FermiTool `gtsrcprob` (M. Kerr 2011). We choose a $3^\circ$ aperture as it roughly matches the PSF at 200 MeV. Other data selections were chosen according to the prescription detailed by R. H. D. Corbet et al. (2019) to provide an update with 16 yr of events (roughly five more years of data). We conduct a power series search on the weighted aperture photometry data to refine the orbital period. We find an orbital period of $P = 13.7157 \pm 0.0014$ days, fully within the errors of the previously reported period (R. H. D. Corbet et al. 2019).

With an updated orbital period, we investigated the evolution of this orbital modulation (and its harmonics) on a dynamic timescale. We generate 750 day lightcurves from the Fermi-LAT events using FermiTools. We use a dynamic Lomb-Scargle periodogram (LSP; e.g., SMC X-1; W. I. Clarkson et al. 2003), with a 100 day sliding window function, to search for modulation between 0.06 and $0.24 \text{ day}^{-1}$. The resulting LSP, along with the relative strengths of the second and third orbital harmonics, is shown in Figure 4. An inspection of the LSP shows significant strengthening and weakening of the orbital modulation, as well as a decreased significance of the second harmonic. We expect these changes to account for the decreasing significance of the second gamma-ray peak near the X-ray Maximum between Figure 5 and Figure 2 of R. H. D. Corbet et al. (2019).

We made additional time selections to search for changes in the gamma-ray modulation based on the periods of relative strengths of the orbital period and its second harmonic, shown in the lower panel of Figure 4. We also fold photons on the updated binary orbital period (see Section 2.3) from the following ranges: (a) the full interval, (b) MJD 54682–55500, (c) MJD 55500–56500, (d) MJD 56500–57333, and (e) MJD 57333–60500; their resulting time-resolved orbital phase-folded lightcurves are shown in Figure 5. The phase-folded lightcurve for the entire time selection in panel (a) shows similar fluxes to those seen in Figure 2 of R. H. D. Corbet et al. (2019). The twin-peaked structure (from the orbital harmonic) of J1405 is visible in panel (b). Within Figure 5, panels (c) and (e) show that the second peak structure disappears and is less pronounced when compared to panels (b) and (d). The variable strength of the orbital modulation may indicate a superorbital period (e.g., as seen by P. C. Gregory 2002 for LS I +61 303); however, we do not find evidence for one in the power spectrum.

## 4.2. Spectral Analysis

The original analysis of J1405 relied on aperture photometry, and while a likelihood spectral analysis was done, the spectral parameters were frozen to the 4FGL catalog values, allowing only the normalization to be fit (R. H. D. Corbet et al. 2019). With $\sim 5$ yr of additional data and an increased signal-to-noise ratio, we conducted a phase-resolved likelihood analysis.

Event files extracted using the FermiTools were analyzed using FermiPy (version 1.2.0). Our source model consists of the Galactic diffuse and isotropic diffuse models: “gll_iem_v07.fits” and “iso_P8R3_SOURCE_V3_v1.txt.” A phase-integrated fit of all cataloged sources within the $15^\circ$ ROI was first performed, where only each source’s normalization is freed from the 4FGL-DR4’s values and fit to obtain a
<figure><img src="image_7.png" /><figcaption>**Figure 4.** A dynamic Lomb–Scargle periodogram of probability-weighted aperture photometry Fermi-LAT data. Data were mean-subtracted and weighted by the exposure. Observations are binned into 750 day bins, and the relative power is calculated on a sliding window of 100 day width. The fundamental and the second and third harmonics of the orbital frequency are indicated by black, red, and blue arrows, respectively. The relative heights of the fundamental (black), first harmonic (red), and second harmonic (blue) can be seen in the bottom panel. A relative power spectrum of the full data is shown in the right panel. Dashed blue and green lines represent the 99% and 99.9% significance levels, respectively.</figcaption></figure>

“baseline” model. All sources were then frozen to this phase-integrated “baseline” model (including the Isotropic and Galactic Diffuse model values), except J1405. J1405’s spectral model, a LogParabola

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_b} \right)^{-(\alpha + \beta \ln(E/E_b))} \qquad (1) $$

is then fit, with all spectral parameters freed ($E_b$ is fixed to the 4FGL-DR4 value), using a binned likelihood analysis for each phase bin. We do not consider other spectral models for J1405. Each fit required a minimum fit quality of 3 with the MINUIT optimizer.

The evolution of J1405’s best-fit spectral curvature parameters ($\alpha$ and $\beta$) as well as the test statistic (TS; a measure of the goodness of the fit) is shown in Figure 6 as a function of orbital phase. SEDs generated from the best-fit models were extracted from the Gamma-ray Maximum (0.964 < $\phi_{\text{Gamma-ray Max.}}$ < 0.0357) and the X-ray Maximum (0.393 < $\phi_{\text{X-ray Max.}}$ < 0.464), which overlap in phase with the XMM-Newton and NuSTAR observations, and are shown in Figure 7. Here, we approximate the shape of the SED at each energy bin as a PL and allow the local spectral index to vary with `use_local_index=True`. Additionally, spectral fluxes with a TS value < 9 are considered below a $3\sigma$ detection threshold, and upper limits are calculated and used instead. A phase-resolved likelihood analysis of the two epochs MJD 56500–57333 and MJD 57333–60500 (where there is a significant difference between the strengths of the fundamental and the second harmonic) found no significant

changes in the spectral parameters associated with the changes in modulation.

# 5. Broadband SED Analysis

## 5.1. Theoretical Constraints

We made the following assumptions to model the observed high-energy SED. The near-IR spectroscopy of J1405 indicates an O6.5 III type companion star (R. H. D. Corbet et al. 2019), so we assume its mass, temperature, and luminosity to be ~35 $M_\odot$, 3.8 × 10$^4$ K, and 2.75 × 10$^5$ $L_\odot$ (M. M. Hanson et al. 2005; L. Mahy et al. 2015). Using the newly derived orbital period of 13.7157 days and this mass estimate for the O6.5 III companion, we calculated, from Kepler’s third law, a binary separation of $a = 0.37$ au, with the assumption that the binary orbit is circular and the compact object is a 1.4 $M_\odot$ neutron star.

We assume the broadband SED is radiated by a population of electrons emitting synchrotron radiation and inverse Compton-scattering stellar photons in a single zone close to the compact object. The stellar radiation field at the location of the compact object is then

$$ U_{\text{seed}} \approx \frac{R^2 \sigma T^4}{a^2 c} \approx 94 \text{ erg cm}^{-3}. \qquad (2) $$

This is much higher than the magnetic energy density $U = B^2/4\pi$ (if $B \le 34$ G), and so, inverse Compton scattering is expected to dominate the particle energy losses. Electrons upscatter stellar photons with an effective energy of
<figure><img src="image_8.png" /><figcaption>Figure 5. Probability-weighted aperture photometry from 200 MeV–500 GeV phase-folded lightcurves of selections between MJD 54682 and 60500. Phase-resolved lightcurves are shown from (a) MJD 54682–60500, (b) MJD 54682–55500, (c) MJD 55500–56500, (d) MJD 56500–57333, and (e) MJD 57333–60500.</figcaption></figure>

$\epsilon \approx 2.8kT \approx 9$ eV (R. Moderski et al. 2005). The interaction occurs in the Klein–Nishina regime if the energy of the electron is such that $4\gamma\epsilon > m_e c^2$, where $\gamma$ is the Lorentz factor of the electron, and so electrons with energy $\gtrsim 7$ GeV ($\gamma \gtrsim 15,000$) interact in the Klein–Nishina regime. The X-ray emission is thought to be due to synchrotron emission from electrons with this typical energy since their typical synchrotron energy is $\approx 0.1$ keV (with $B = 34$ G).

The case of continuous injection of a PL distribution of electrons cooling due to inverse Compton scattering off an external photon field and synchrotron radiation produces a specific break associated with the transition from dominant inverse Compton cooling to dominant synchrotron cooling for the electron population (R. Moderski et al. 2005; G. Dubus 2006). In the synchrotron spectrum, the characteristic synchrotron energy of the electrons at $E_b$ provides an estimate of the magnetic field that can be written as (G. Dubus 2013)

$$B \approx \left( \frac{750 \text{ keV}}{h\nu_{\text{sync}}} \right) \left( \frac{0.1 \text{ au}}{a} \right)^2 \text{G.} \quad (3)$$

Since we see no significant spectral break with XMM-Newton and NuSTAR, setting $E_b = h\nu_{\text{sync}} \gtrsim 20$ keV implies $B \lesssim 2.7$ G, which is well within the Klein–Nishina regime.

Assuming a PL injection with $E^{-2}$ (with E as the electron energy), the steady-state particle spectrum will be $E^{-3}$ for energies above the break energy, $E_b$, and harder below $E_b$ due to inefficient losses in the Klein–Nishina regime. The electron

<figure><img src="image_9.png" /><figcaption>Figure 6. Phase-resolved Fermi-LAT lightcurve of J1405 using the likelihood analysis and phase evolution of spectral parameters. Each bin represents 1/14th of the binary orbit (~1 day). (a) Phase-resolved lightcurve showing fluxes integrated from 200 MeV to 500 GeV. (b) Spectral index $\alpha$. (c) Spectral index $\beta$. (d) The test statistic associated with J1405.</figcaption></figure>

distribution can harden to $E^{-1.5}$ depending on how strongly inverse Compton dominates and the extent of the particle spectrum (see Figure 2 of R. Moderski et al. 2005). We chose to model this injection population with an exponential cutoff broken power law (ECBPL; Equation (5)).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Overlaid SEDs from 200 MeV–500 GeV are shown for both the Gamma-ray Maximum (red) and the X-ray Maximum (blue) along with their respective upper limits. The spectral models are shown as dashed lines, and 95% confidence intervals are demarcated by the shaded regions.</figcaption></figure>

The maximum energy of the particle distribution is capped by synchrotron burnoff (O. C. de Jager et al. 1996), when the timescale to radiate the energy of the particle is equal to its gyration timescale in the magnetic field (G. Dubus et al. 2015):

$$E_{\text{max}} \approx 60\xi^{-1/2} \left(\frac{1}{B}\right)^{1/2} \text{TeV}, \quad (4)$$

where $\xi \ge 1$ is the acceleration efficiency, and so in the following we set $\xi = 1$. The maximum particle energy is thus about 36 TeV (for $B = 2.7$ G).

The synchrotron spectrum is capped to a maximum around 30 MeV with this assumption on the maximum particle energy, which is reasonable for Fermi acceleration at shocks (O. C. de Jager et al. 1996). This falls short of the GeV emission if one interprets it as synchrotron emission from the same population of electrons that radiate X-rays. Rather than assuming an arbitrary (unphysical) large maximum particle energy to enable the synchrotron emission to reach the GeV regime, we attribute the GeV component to a Maxwellian population of electrons originating from the IBS, where a fraction of the pulsar wind particles are randomized instead of being accelerated to a power law (G. Dubus et al. 2015).

## 5.2. Broadband Fitting

With these theoretical assumptions of the geometry of the binary and the contributions of synchrotron and inverse Compton emission, the broadband modeling of the X-ray and gamma-ray data was conducted using the Naima Python package (version 0.10 V. Zabalza 2015),<fnref n="11" /> a Markov Chain Monte Carlo (MCMC) fitting package. X-ray data from XMM-Newton and NuSTAR were rebinned to contain a minimum of 100 counts per bin and were corrected for absorption by taking the ratio of the absorbed and unabsorbed models and applying it to each data point while taking care to include the cross-constant normalizations. The Fermi-LAT data from the Gamma-ray Maximum and the X-ray Maximum were generated from the SEDs with FermiPy, where upper limits

<table><caption>**Table 3**
Broadband Best-fit Parameters</caption><thead><tr><th>Model Parameters</th><th>Gamma-Ray Maximum</th><th>X-Ray Maximum</th></tr></thead><tbody><tr><td>ECBPL<sup>a</sup></td><td></td><td></td></tr><tr><td>Log<sub>10</sub>(A<sup>b</sup> (eV<sup>-1</sup>))</td><td>2.98 ± 0.03</td><td>3.82 ± -0.04</td></tr><tr><td>α<sub>1</sub></td><td>1.87 ± 0.02</td><td>1.49<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>E<sub>b</sub> (×10 TeV)</td><td>0.26 ± 0.04</td><td>0.30 ± 0.01</td></tr><tr><td>Maxwellian<sup>c</sup></td><td></td><td></td></tr><tr><td>Log<sub>10</sub>(K<sup>d</sup> (eV<sup>-3</sup>))</td><td>-5.77 ± 0.04</td><td>-5.59<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>E<sub>char</sub> (GeV)</td><td>0.77 ± 0.01</td><td>0.60 ± 0.01</td></tr><tr><td>Log-likelihood</td><td>-21.95</td><td>-16.85</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>43.91/35</td><td>33.69/34</td></tr></tbody></table>

**Notes.**

a Exponential cutoff broken power law (see Equation (5)).

b Model normalization eV<sup>-1</sup>.

c Maxwellian distribution (see Equation (6)).

d Model normalization eV<sup>-3</sup>.

were used if the TS of the energy bin was less than 9 (corresponding to a $3\sigma$ detection from noise).

Within the Naima fitting package, as discussed above, we consider the two populations of electrons: an ECBPL

$$\frac{dN}{dE} = \exp(-E/E_{\text{max}}) \times \begin{cases} A(E/E_0)^{-\alpha_1}, & E < E_b \\ A(E_b/E_0)^{\alpha_2-\alpha_1}(E/E_0)^{-\alpha_2}, & E > E_b \end{cases} \quad (5)$$

and a Maxwellian distribution

$$\frac{dN}{dE} \equiv KE^2 \exp(-E/E_{\text{char}}). \quad (6)$$

We set $\alpha_2 = 3$, $B = 2.7$ G, $E_{\text{min}} = 511$ keV (the minimum energy for a gamma ray, i.e., $\gamma = 1$), $E_{\text{max}} = 36$ TeV, and $E_b$ is set to vary around 1 TeV as per the discussion above. $E_0$ is set to 10 TeV. We allowed $\alpha_1$, the spectral index below the break, and $E_{\text{char}}$, the characteristic energy of the particles at the IBS, as well as the normalizations from both electron distributions, to vary and fit using Naima. The best-fit values are given in Table 3 along with their log-likelihood values (and the related $\chi^2$ values). Figure 8 shows the broadband SEDs at the Gamma-ray Maximum and the X-ray Maximum extending from optical to TeV energies, respectively. Using the IBS model to constrain the eccentricity, $e$, of J1405, we fitted the SED with an additional parameter, but did not find a statistically significant estimate for the eccentricity.

# 6. Discussion

## 6.1. X-Ray Timing and $N_H$

We find a $>2\sigma$ detection during XMM-Newton observation 0852020101 at a frequency of $\nu \approx 22$ Hz on the $Z_2^2$ search, but this detection was not consistent throughout the different instruments and observations. Therefore, we report no detection of pulsations between frequencies of $5 \times 10^{-6}$ Hz and $1 \times 10^3$ Hz above a 1% false-alarm probability (as seen in Figures 3(a) and (b)) besides the orbital period of NuSTAR and its harmonics and low-frequency signals attributed to red noise. Since most HMGBs do not exhibit pulsations in soft

<footnote n="11">¹¹ See technical documentation and examples at https://naima.readthedocs.io/en/latest/examples.html.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Broadband SEDs of the Gamma-ray Maximum and the X-ray Maximum. Data include XMM-Newton MOS1/2 and PN, NuSTAR FPMA/FPMB, and Fermi-LAT fluxes. Fermi-LAT upper limits represent a 95% confidence interval and are denoted by one-sided error bars. The cumulative spectral model is shown as a black solid line. The inverse Compton components (shown as dotted lines) represent both the ECBPL (blue) and Maxwellian (red) particle distributions. Similarly, the synchrotron components are shown as dashed lines.</figcaption></figure>

($E < 10$ keV) or hard ($E > 10$ keV) X-rays (M. Chernyakova et al. 2023), this is not enough evidence to reject a rotation-powered neutron star as the compact object. The distance ($\sim 7.7$ kpc) to J1405 and strong $N_H$ absorption contribute to the low signal-to-noise ratio and prevent pulsation searches outside the 1–20 keV energy range in the XMM-Newton and NuSTAR observations. Although our energy range is limited, there is little precedence of finding pulsations at the soft/hard X-ray energies (e.g., PSR B1259−63 or LS I +61°303; S. Johnston et al. 1992; S.-S. Weng et al. 2022).

The simultaneous spectral fits using best-fit models of XMM-Newton and NuSTAR observations at the Gamma-ray Maximum and the X-ray Maximum show two statistically different $N_H$ column densities. At the Gamma-ray Maximum, we detected a slightly greater column density of about $(8.0 \pm 0.6) \times 10^{22} \text{cm}^{-2}$ when compared to $(5.7 \pm 0.5) \times 10^{22} \text{cm}^{-2}$ for the X-ray Maximum. However, a joint fit between observations does not show an improvement when allowing the Gamma-ray Maximum and the X-ray Maximum $N_H$ values to fit independently over a model that fixes $N_H$ to the phase-averaged value. As such, we do not find a statistically significant change in the geometry of the IBS of the binary. We also find a slight hardening of the X-ray photon index $\Gamma_X$ from $\sim 1.6 \pm 0.1$ to $\sim 1.4 \pm 0.1$ from the Gamma-ray Maximum to the X-ray Maximum, which is in agreement with the Swift-XRT measurements reported by R. H. D. Corbet et al. (2019) and the independent NuSTAR/XMM-Newton analysis of E. A. Saavedra et al. (2023). We use the lack of a spectral break below 20 keV as a constraint in the broadband modeling for the magnetic field ($B \approx 2.7$ G) at the IBS and for the synchrotron radiation emission, and a more concrete limit would have implications for broadband modeling.

Previously, the presence of a blackbody component at $\sim 1$ keV was claimed in the X-ray spectrum observed on 2019 August 24 (Gamma-ray Maximum), making J1405 the only HMGB exhibiting a blackbody component (E. A. Saavedra et al. 2023). Moreover, a CRSF was reported in both XMM-Newton observations (S.-C. Chiu & K.-L. Li 2024), which could indicate that the compact object is a slowly rotating neutron star or a magnetar. This would classify J1405 as a source similar to LS 5039 (H. Yoneda et al. 2020), although this has been refuted by I. Volkov et al. (2021) and O. Kargaltsev et al. (2024). We find that there was not a significant number of 10,000 trials ($\gtrsim 68\%$) that show an improvement to the model for either the Gamma-ray Maximum or the X-ray Maximum. CRSFs are typically found in harder X-ray spectra (R. Staubert et al. 2019).

Furthermore, a comparison of spectral fits of the best-fit model between the full energy range of XMM-Newton and an energy range that excludes data below 2 keV shows no significant change in the fit quality or values of fit parameters. Due to the intrinsic high absorption, we attribute any reported feature excess below 2 keV to statistical noise. Both features are found below 3 keV in the heavily absorbed X-ray spectrum, where the background dominates and should be treated with caution. Therefore, it is difficult to argue for the presence of either cyclotron feature in the X-ray spectra.

## 6.2. Gamma-Ray Timing and Spectral Fitting

Using 16 yr of Fermi-LAT data, we constrain the variation of the parameters $\alpha$ and $\beta$ in the LogParabola spectral model. Our phase-average best-fit value of $\alpha = 2.85 \pm 0.06$ is in good agreement with the 4FGL-DR4 catalog; however, the value of $\beta = 0.38 \pm 0.01$ is not—the DR4 catalog value is $0.28 \pm 0.04$. This difference may be attributed to the use of likelihood weights in the DR4 catalog.

There is a slight modulation present for the spectral parameter $\beta$, while $\alpha$ seems to remain constant. $\beta$ appears to be anticorrelated with the gamma-ray flux as shown in Figure 6. One bin, $\phi = 0.46–0.54$, seems to have a significantly different $\alpha = 2.41 \pm 0.11$ and $\beta = 0.06 \pm 0.07$, and may be indicative of a physical change in the binary.
However, due to the nature of the large error bars for α and β, it is otherwise difficult to make meaningful comments on the evolution of the SED as a function of orbital phase.

Lastly, we provide different results for the gamma-ray orbital modulation to the previous study. R. H. D. Corbet et al. (2019) found a double-peaked structure in the gamma-ray region using probability-weighted photometry. However, when using a likelihood analysis, this structure at the X-ray Maximum becomes less significant during different time periods, as seen in Figure 5. Upon revisiting the probability-weighted aperture photometry lightcurve, the structure is variable. Within the top and bottom panels of Figure 6, there is no significant change in the observed flux from J1405 from R. H. D. Corbet et al.'s (2019) discovery analysis, besides the decreased feature at the X-ray Maximum ($\phi \approx 0.4$). Additionally, the dynamic power spectrum of the aperture photometry in Figure 5 shows that the significance of the orbital period harmonics (the second peak) varies. This suggests periods of gamma-ray activity that may be caused by changes such as the structure of the IBS or clumps in the stellar wind (S. Johnston et al. 1992; X. Paredes-Fortuny et al. 2015; E. Kefala & V. Bosch-Ramon 2023). Additionally, some supergiant high-mass X-ray binaries (sgHMXBs) have been observed to display similar variation in strength of the superorbital period harmonics (e.g., N. Islam et al. 2023). Many of these sgHMXBs are thought to have a changing wind geometry such as a corotating interaction region, possibly linking the two classes of binaries.

### 6.3. Comparison between IBS and Microquasar Broadband Models

By setting the spectral break in the X-ray spectrum above 20 keV, we approximate the magnetic field to be ~2.7 G and fit a two-component IBS model similar to LS 5039 (G. Dubus et al. 2015). Within the ECBPL distribution, upon fixing $\alpha_2 = 3$, we find that the electron distribution indeed hardens to nearly 1.5 for the X-ray Maximum but remains softer for the Gamma-ray Maximum. Additionally, the contribution of the ECBPL distribution accounts for the increase in X-ray flux. In both cases, synchrotron emission from the ECBPL sufficiently models the X-ray spectra. As shown in Section 4, there is little modulation in the gamma-ray SEDs, and so it is sensible that the Maxwellian distribution shows little change in K and $E_{\text{char}}$. The synchrotron emission from the Maxwellian distribution is in the optical band, where the companion's emission dominates the spectrum. The IC emission accounts for the MeV–GeV emission well. However, the upper limits of the GeV emission are best fit with the inclusion of the IC emission from the ECBPL for both phases. Additionally, preliminary observations from the High Energy Stereoscopic System (H.E. S.S., T. Martinez 2023)<fnref n="12" /> indicate VHE fluxes consistent with Fermi-LAT's upper limits and our modeling of the GeV/TeV emission.

The microquasar and IBS models offer different mechanisms to explain the high-energy emission in HMGBs. In the microquasar model, the gamma-ray emission arises from relativistic jets from the compact object as it accretes matter from its massive companion. Conversely, the IBS model focuses on the collision between the stellar winds of the binary

and the pulsar wind nebula of the neutron star, leading to the acceleration of particles and the generation of gamma rays.

The IBS model has an advantage over the microquasar jet model in that it requires significantly fewer fitting parameters. Given the quality of the X-ray and gamma-ray data, fitting a model with many free parameters would likely lead to degenerate solutions and large uncertainties, especially without direct observational evidence for components such as disk or coronal emissions.

With a few simplifications, we are able to decrease the number of free parameters to five, fitting just the normalizations of particle distributions, the ECBPL slope α, the break energy of the IBS $E_b$, and the cutoff energy for the GeV emission $E_{\text{char}}$. Additionally, because we used Naima and an MCMC methodology in our fitting, we find a low correlation and variance between parameters (see Figures A1 and A2 in the Appendix). We find statistically satisfying fits of $\chi^2_\nu = 1.25$ and $\chi^2_\nu = 0.99$ (for the Gamma-ray Maximum and X-ray Maximum, respectively) for the broadband data with physical parameters used in previous IBS modeling of HMGBs (G. Dubus et al. 2015). In this interpretation, orbit-to-orbit variability could be due to fluctuations in the IBS position or properties due to variability or clumping of the stellar wind.

## 7. Conclusion

We present detailed orbital phase-resolved multiwavelength observations of the HMGB 4FGL J1405.1–6119, utilizing XMM-Newton and NuSTAR observations at the Gamma-ray Maximum and the X-ray Maximum (at the gamma-ray maximum and minimum, respectively), supplementing nearly 16 yr of Fermi-LAT data. We found no significant pulsations in the X-ray observations after searching a frequency range of roughly $1 \times 10^{-5}–1 \times 10^3$ Hz.

X-ray spectral analysis between 1.0 and 20 keV reveals no evidence of either a thermal component below 2 keV or a CRSF component at ~2 keV. There is a slight hardening of the X-ray spectrum from the Gamma-ray Maximum to the X-ray Maximum, without evidence of a spectral break during either observation. Remarkably, the X-ray spectrum does not exhibit softening up to at least 20 keV. The lack of a significant detection of a spectral break below 20 keV allows for a rough estimate of the magnetic field at the IBS to be at most 2.7 G, which is used for broadband analysis.

Orbital phase-resolved analysis of Fermi-LAT data (200 MeV–500 GeV) provides the evolution of the spectral shape as a function of orbital phase. We find a slight evolution of spectral parameter β, but the SED experiences very little modulation in other parameters, and errors are large. The previously reported second peak at $\phi \approx 0.4$ seems to fade in significance for our binned likelihood analysis and is variable during our probability-weighted aperture photometry analysis for certain time ranges.

We lastly provide an alternative broadband modeling of the X-ray and gamma-ray data to the microquasar model with an IBS model to explain the high-energy emission. Our IBS model provides a good fit with fewer free parameters and results in $\chi^2_{\text{red}}$ values of ~1.25 and ~0.99 for both the Gamma-ray Maximum and the X-ray Maximum.

Due to the high absorption of the soft X-rays and discovery of the weakening gamma-ray modulation, we recommend further (and more sensitive) X-ray observations and GeV and

<footnote n="12">¹² These H.E.S.S. observations are from a Master’s thesis that can be found at: http://astro.uni-tuebingen.de/publications/diplom.shtml</footnote>
TeV gamma-ray observations to better understand the nature
and high-energy emission of 4FGL J1405.1–6119.

Acknowledgments

The authors thank the referee for useful comments and
feedback. The Fermi-LAT Collaboration acknowledges gener-
ous ongoing support from a number of agencies and institutes
that have supported both the development and the operation of
the LAT as well as scientific data analysis. These include the
National Aeronautics and Space Administration and the
Department of Energy in the United States, the Commissariat
à l'Energie Atomique and the Centre National de la Recherche
Scientifique/Institut National de Physique Nucléaire et de
Physique des Particules in France, the Agenzia Spaziale
Italiana and the Istituto Nazionale di Fisica Nucleare in Italy,
the Ministry of Education, Culture, Sports, Science and
Technology (MEXT), High Energy Accelerator Research
Organization (KEK) and Japan Aerospace Exploration Agency
(JAXA) in Japan, and the K. A. Wallenberg Foundation, the
Swedish Research Council and the Swedish National Space
Board in Sweden.

This work was also partially supported by NASA grant 80NSSC20K1304 and also under NASA award number 80GSFC21M0006. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This research has made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. Data analysis was performed using the NuSTAR Data Analysis Software (NuSTARDAS), jointly developed by the ASI Science Data Center (SSDC, Italy) and the California Institute of Technology (USA). This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

Appendix
Broadband Modeling: MCMC Corner Plots

From the broadband modeling with the Naima Python
module, we generate corner plots (Figures A1 and A2) to show
the variance and correlations between the fit parameters from
the ECBPL and Maxwellian distributions listed in Table 3.
These corner plots show the two-dimensional distributions of
posteriors and the maximum likelihoods for each set of
parameters. On the diagonal, one-dimensional distributions are
also shown for each parameter. For both the Gamma-ray
Maximum and the X-ray Maximum, we find no strong
correlations and no bimodal distributions of parameters.

<figure><img src="image_12.png" /><figcaption>**Figure A1.** A corner plot of the sample density in the two-dimensional parameter space of all parameter pairs of the broadband fit of the Gamma-ray Maximum. Confidence intervals of $1\sigma$, $2\sigma$, and $3\sigma$ for a two-dimensional Gaussian distribution are overlaid for each panel, as well as the one-dimensional 68% confidence interval along the diagonal. Dashed horizontal and vertical black lines indicate the one-dimensional 68% confidence interval for each parameter. The maximum likelihood value for each parameter is shown in blue.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure A2. As in Figure A1, but for the X-ray Maximum fit.</figcaption></figure>

ORCID iDs

Alexander Lange https://orcid.org/0000-0003-3540-2870
Robin H.D. Corbet https://orcid.org/0000-0002-3396-651X
Joel B. Coley https://orcid.org/0000-0001-7532-8359
Guillaume Dubus https://orcid.org/0000-0002-5130-2514
Jeremy Hare https://orcid.org/0000-0002-8548-482X
Nazma Islam https://orcid.org/0000-0002-2413-9301
References

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514v1

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bachetti, M., Markwardt, C. B., Grefenstette, B. W., et al. 2021, ApJ, 908, 184

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv:1810.11394

Buccheri, R., Bennett, K., Bignami, G. F., et al. 1983, A&A, 128, 245

Camilo, F., Ray, P. S., Ransom, S. M., et al. 2009, ApJ, 705, 1

Chernyakova, M., & Malyshev, D. 2020, in Multifrequency Behavior of High Energy Cosmic Sources - XIII (Trieste: SISSA), 45

Chernyakova, M., Malyshev, D., Mc Keague, S., et al. 2020, MNRAS, 497, 648

Chernyakova, M., Malyshev, D., Neronov, A., & Savchenko, D. 2023, MNRAS, 525, 2202

Chernyakova, M., Neronov, A., van Soelen, B., et al. 2015, MNRAS, 454, 1358

Chiu, S.-C., & Li, K.-L. 2024, RNAAS, 8, 89

Clarkson, W. I., Charles, P. A., Coe, M. J., et al. 2003, MNRAS, 339, 447

Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2019, ApJ, 884, 93

de Jager, O. C., Harding, A. K., Michelson, P. F., et al. 1996, ApJ, 457, 253

Dubus, G. 2006, A&A, 456, 801

Dubus, G. 2013, A&ARv, 21, 64

Dubus, G., Guillard, N., Petrucci, P.-O., & Martin, P. 2017, A&A, 608, A59

Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177

Fioretti, V., Mineo, T., Lotti, S., et al. 2024, A&A, 691, A229

Gregory, P. C. 2002, ApJ, 575, 427

Hanson, M. M., Kudritzki, R.-P., Kenworthy, M. A., Puls, J., & Tokunaga, A. T. 2005, ApJS, 161, 154

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

Heasarc, 2014 HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Islam, N., Corbet, R. H. D., Coley, J. B., Pottschmidt, K., & Fuerst, F. 2023, ApJ, 948, 45

Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1

Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJL, 387, L37

Kargaltsev, O., Hare, J., & Lange, A. 2024, ApJ, 975, 199

Kefala, E., & Bosch-Ramon, V. 2023, A&A, 669, A21

Kerr, M. 2011, ApJ, 732, 38

Kim, J., An, H., & Mori, K. 2022, ApJ, 936, 32

Lo, D., & Srour, J. 2003, ITNS, 50, 2018

Mahy, L., Rauw, G., De Becker, M., Eenens, P., & Flores, C. A. 2015, A&A, 577, A23

Martinez, T. 2023, Master's Thesis, http://astro.uni-tuebingen.de/publications/diplom/martinez-master.pdf

Mirabel, I. F. 2012, Sci, 335, 175

Moderski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005, MNRAS, 363, 954

Paredes-Fortuny, X., Bosch-Ramon, V., Perucho, M., & Ribó, M. 2015, A&A, 574, A77

Piano, G., Tavani, M., Vittorini, V., et al. 2012, A&A, 545, A110

Saavedra, E. A., Fogantini, F. A., Escobar, G. J., et al. 2023, A&A, 680, A88

Staubert, R., Trümper, J., Kendziorra, E., et al. 2019, A&A, 622, A61

Strader, J., Chomiuk, L., Cheung, C. C., Salinas, R., & Peacock, M. 2015, ApJL, 813, L26

Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18

The Fermi LAT Collaboration 2012, Sci, 335, 189

Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27

van Soelen, B., Mc Keague, S., Malyshev, D., et al. 2022, MNRAS, 515, 1078

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Volkov, I., Kargaltsev, O., Younes, G., Hare, J., & Pavlov, G. 2021, ApJ, 915, 61

Weisskopf, M. C., Brinkman, B., Canizares, C., et al. 2002, PASP, 114, 1

Weng, S.-S., Qian, L., Wang, B.-J., et al. 2022, NatAs, 6, 698

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Wood, M., Caputo, R., Charles, E., et al. 2017, ICRC (Busan), 35, 824

Yoneda, H., Makishima, K., Enoto, T., et al. 2020, PhRvL, 125, 111103

Zabalza, V. 2015, ICRC (The Hague), 34, 922