© The Authors 2023

1. Introduction

Binary sources containing neutron stars (NSs) or black holes (BHs) dominate the Galactic X-ray emission above 2 keV (see, e.g., Grimm et al. 2002). These systems are called X-ray binaries and are usually divided into two major classes, high-mass X-ray binaries and low-mass X-ray binaries, according to the mass of the donor star (mass ≳8 M_⊙ for the former and ≲8 M_⊙ for the latter).

Within the high-mass X-ray binary class, three types of non-transient systems can emit γ-ray radiation (Chernyakova et al. 2019; Chernyakova & Malyshev 2020):

1. Colliding-wind binaries involve the interaction of two massive stars, whose nonrelativistic winds collide and produce γ-ray emission. Prominent examples of such systems are η-Carinae, WR 11, and Apep. In colliding-wind binaries, intense shocks occur in the wind collision region, leading to the formation of a very hot plasma (> 10⁶ K). In addition, these systems have the ability to accelerate relativistic particles (Eichler & Usov 1993; Benaglia & Romero 2003), which classifies them as particle-accelerating colliding-wind binaries (e.g., De Becker & Raucq 2013; del Palacio et al. 2023, and references therein).

2. Gamma-ray binaries are characterized by the presence of a young, magnetized, and rapidly rotating NS that emits a relativistic wind. This wind collides with the nonrelativistic wind of the companion OB star, and through this interaction they emit nonthermal radiation that peaks at high energies (E > 100 MeV). These systems are thought to represent a short phase in the evolution of massive binaries, which comes after the birth of the compact object (CO) and is followed by the X-ray binary phase, in which the CO accretes matter from its companion (e.g., Dubus et al. 2017; Saavedra et al. 2023, and references therein). In the latter phase, the nonthermal emission of the system peaks in X-rays. Examples of γ-ray binaries are LS I+61°303, PSR J1259–63, PSR J2032+4127, LS 5039, 1FGL J1018.6–5856, and HESS J0632+057.

3. Microquasars (MQs) differ from the previous two categories in that their emission of γ-rays does not come from wind collisions. Instead, it comes from the jets ejected by the CO (BH or NS) and their interaction with the environment. Examples of MQs are Cyg X-1, Cyg X-3, and SS 433.

In total, there are nine known non-transient γ-ray-emitting binaries in our Galaxy (e.g., Corbet et al. 2019; Dubus et al. 2017; Chernyakova & Malyshev 2020, and references therein). Five of these systems contain an O-type star: LS 5039, LMC P3, 1FGL J1018.6–5856, HESS J1832–093, and 4FGL J1405.1–6119. The remaining four contain a Be star: HESS J0632+057, LS I+61°303, PSR B1259–63, and PSR J2032+4127. Radio pulsations were observed in three of these systems – PSR J2032+4127 (Camilo et al. 2009), PSR B1259–63 (Moldón et al. 2014), and LS I+61°303 (Weng et al. 2022) – evidence in favor of the presence of a NS. Many more galactic γ-ray binaries are expected to exist. We recall that the nature of many of the γ-ray sources detected by telescopes such as Fermi-LAT is not known, and some of the sources may be γ-ray binaries. In fact, the number of detectable Galactic γ-ray binaries is estimated to be ∼100 (Dubus et al. 2017).

4FGL J1405.1−6119 (also known as 3FGL J1405.4−6119) is a high-mass γ-ray-emitting binary first studied by Corbet et al. (2019). Using Fermi and Swift/XRT data, Corbet et al. (2019) found a strong modulation of 13.713 ± 0.002 days associated with the system’s orbital period. The companion is classified as an O6.5 III star with a mass of about 25 − 35 M_⊙ (Mahy et al. 2015). The absence of partial and total eclipses suggests that this system has a low inclination (60°).

To explain the origin of its γ-ray emission, Corbet et al. (2019) proposed that 4FGL J1405.1−6119 may be a γ-ray binary. This hypothesis draws on analogies with other similar systems, such as 1FGL J1018.6–5856 and LMC P3. Xingxing et al. (2020) modeled the GeV time behavior of 4FGL J1405.1−6119 in the γ-ray binary scenario, assuming a binary consisting of a young pulsar and an O-type main sequence star. Conversely, the radio (5.5 GHz and 9 GHz) and X-ray (0.2−10 keV) luminosities show a positive correlation (Corbet et al. 2019), as expected in MQs (Falcke et al. 2004).

We were able to perform a detailed temporal and spectral study of 4FGL J1405.1−6119 over a broad energy range by analyzing simultaneous XMM-Newton and NuSTAR observations. In this paper we present our results and conclusions about this source. In Sect. 2 we present the X-ray observations and the corresponding tools used for their analysis. Our main results are presented in Sect. 3. A jet model for the source is introduced and discussed in Sect. 4. We discuss our results and present our conclusions in Sect. 5.

2. Observation and data analysis

2.1. XMM-Newton data

The XMM-Newton observatory is equipped with an optical instrument and two X-ray instruments: the Optical Monitor, which is mounted on the mirror support platform and provides coverage between 170 nm and 650 nm of the central 17 arcmin square region; the European Photon Imaging Camera (EPIC); and the Reflecting Grating Spectrometers (RGS). The EPIC instrument comprises three detectors – the pn camera (Strüder et al. 2001) and two MOS cameras (Turner et al. 2001) – which are most sensitive in the 0.3 − 10 keV energy range. The RGS instrument comprises two high-resolution spectrographic detectors sensitive in the energy range 0.3 − 2 keV.

XMM-Newton observed 4FGL J1405.1−6119 on August 17, 2019, with an exposure time of 32 ks (ObsID 0852020101) and on August 24, 2019, with an exposure time of 44 ks (ObsID 0852020201). In both observations, the MOS cameras were in large window mode, and the pn camera was in timing mode.

We reduced the XMM-Newton data using Science Analysis System (SAS) v20.0 and the latest available calibration in early 2022. To process the observation data files, we used the EPPROC and EMPROC tasks. We selected circular regions with radii of 18 arcsec and 36 arcsec for the source and the background, respectively, with the latter away from any source contamination. We then filtered the raw events lists, removing the high-energy, single-pattern particle backgrounds and thus creating cleaned event lists for each camera and observation. The resulting exposure times after background filtering are 19 ks (59% of total) for the first observation, and 44 ks (100% of total) for the second observation. We studied the presence of pile-up with the EPATPLOT task and did not find any deviation of the data from the expected models; thus, no excision radii were applied. We barycentered each cleaned event list with the barycen task in order to perform precise timing studies. EPIC light curves were summed using the LCMATH task, with proper scaling factors for the different source photons collecting areas. We extracted and grouped spectra with a minimum of 25 counts per bin in the 0.5 − 12 keV energy band.

2.2. NuSTAR data

The NuSTAR X-ray observatory was sent into orbit in the year 2012 and is notable for its exceptional sensitivity at hard X-rays. It is equipped with two X-ray grazing incidence telescopes, designated FPMA and FPMB, which are arranged in parallel and contain 2 × 2 solid-state CdZnTe detectors each. NuSTAR can operate in the energy range of 3−79 keV and can achieve an angular resolution of 18 arcsec as reported in Harrison et al. (2013).

NuSTAR observed 4FGL J1405.1−6119 on August 16, 2019 (58711.7031 MJD – ObsID 30502015002), with an exposure of time of ∼61 ks and on August 24, 2019 (58719.3511 MJD – ObsID 30502015004) with an exposure time of ∼86 ks. We processed NuSTAR data using NuSTARDAS-v.2.1.2 from HEASoft v.6.30 and CALDB (v.20211020) calibration files. We took source events that accumulated within a circular region of 85 arcsec around the focal point. The chosen radius encloses ∼90% of the point spread function. We took a circular source-free region with a radius of 160 arcsec to obtain the background events within the same CCD.

We used the nupipeline task to create level 2 data products, with SAA parameters saacalc=1 saamode=strict tentacle=yes to filter for high-energy particle background, obtained from SAA filtering reports¹. We extracted light curves and spectra with the nuproducts task. We obtained the barycenter-corrected light curves using barycorr task with nuCclock20100101v136 clock correction file. We used celestial coordinates α = 211.2472° and δ = +61.3234° for the barycentric correction with JPL-DE200 as the Solar System ephemeris. Finally, we subtracted the background from each detector’s light curve. Then we used the LCMATH task to create FPMA+B light curves. We extracted and re-binned spectra with a minimum of 25 counts per bin in the 3 − 78 keV energy band.

We used the XSPEC v12.12.1 package (Arnaud 1996) to model XMM-Newton and NuSTAR spectra, with parameters uncertainties reported at the 90% confidence level.

3. Results

3.1. Analysis of the light curves

Figure 1 shows the background-corrected light curves obtained from XMM-Newton (0.5−10 keV, top panel) and NuSTAR (3−78 keV, bottom panel) missions, with a bin time of 350 s. The observation conducted on August 17 is labeled as Epoch 1, while the observation on August 24 is labeled as Epoch 2. The long-term exposures of NuSTAR, both on the order of ∼1.5 d, show a hard X-ray flux modulation of ∼1 d seen in both epochs. The shorter but continuous exposures of XMM-Newton do not capture this behavior. Instead, it captures very short (on the order of some ks) changes in soft X-ray flux, as seen in Epoch 1.

Fig. 1. Background-corrected light curves of 4FGL J1405.1−6119 observed by NuSTAR (FPMA+B; red) and XMM-Newton (EPIC pn+MOS; blue) with a binning of 350 s, starting at 58711.705573 MJD. The first section of the light curve corresponds to the observation from August 18 (Epoch 1), while the second light curve is associated with the observation from August 25 (Epoch 2).

Figure 2 shows the orbital flux modulation associated with each observation and mission. Epoch 1 occurred within the orbital phases of 0.93−1.08, while Epoch 2 occurred within the orbital phases of 0.37−0.48. The observed orbital behavior of XMM-Newton and NuSTAR data is very similar to that reported by Corbet et al. (2019) using Swift/XRT data.

Fig. 2. NuSTAR and XMM-Newton folded light curve using 48 phase bins, an orbital period of 13.713 days, and with 56498.7 MJD as the reference epoch (Corbet et al. 2019). A sine function fit is shown in orange (see the main text for details). The observed orbital modulation is similar to that shown by Corbet et al. (2019) using Swift/XRT data.

We used a sinusoidal model with Gaussian measurement errors to visualize the orbital modulation through ∼10⁴ simulations (Buchner 2021). Specifically, the following was used:

$$\begin{array}{c}\hfill y=Asin\left(2\pi [\frac{t}{P}+t_0]\right)+B+ϵ,\end{array}$$(1)

where ϵ ∼ Normal(0, σ), that is, a normal distribution with a mean of zero and a standard deviation of σ. We obtain the following values: A = 0.56 ± 0.3, P = 1.002 ± 0.002, t₀ = 0.623 ± 0.001 and B = 1.13 ± 0.01. The fitted model is shown in Fig. 2.

From phase ∼0.93, the flux starts to decay, and from phase ∼0.37 the source has the maximum local emission. This modulation is anticorrelated with the modulation obtained from the Fermi data in the energy range 200 MeV−500 GeV (Corbet et al. 2019).

We employed spectral timing routines provided by Stingray software (Huppenkothen et al. 2019) to conduct a comprehensive search for any potential pulsation linked to the X-ray source. Light curves from both telescopes across various energy ranges do not show any significant pulsations above noise on the 0.1−100 mHz frequency range.

3.2. Spectral analysis

We simultaneously modeled source and background spectra extracted from all five detectors (pn, MOS1, MOS2, FPMA and FPMB). We introduced calibration constants in our models to account for disparities in effective areas between instruments. The pn constant was fixed to unity, while the remaining calibration constants were permitted to vary: C_MOS1 = 1.10 ± 0.16, C_MOS2 = 1.05 ± 0.15, C_FPMA = 1.40 ± 0.16, and C_FPMB = 1.35 ± 0.16. Each epoch was modeled separately.

The NuSTAR background showed significant activity at energies above 20 keV. As a result, we limited our analysis to the energy range 3−20 keV for both epochs. In the case of the XMM-Newton background, it was significant at energies below 2 keV. Therefore, we focused our analysis on the 2−10 keV energy range for both epochs. Consequently, the total spectrum used for each epoch included the 2−20 keV energy range.

The interstellar absorption was modeled using the Tübingen-Boulder model (tbabs), with solar abundances defined by Wilms et al. (2000) and effective cross sections of Verner et al. (1996). Several continuum models were used to fit the time-averaged spectra, including combinations of power-law variants such as powerlaw, highecut+powerlaw, cutoffpl, and bknpow, and thermal models such as apec, bbodyrad, and diskbb. After several fits, we selected the best fits: a power law (tbabs powerlaw, hereafter Model 1) and a combination of a power law and a blackbody (tbabs(powerlaw+bbodyrad), hereafter Model 2). No emission or absorption lines above the continuum were observed in the spectra of all epochs.

Model 1 yielded a χ² of 221 with 191 degrees of freedom (${\chi }_{\nu }^2=1.16$) for Epoch 1, while Epoch 2 resulted in 184/194 − 0.95. On the other hand, if we apply Model 2, we obtain 216/189 − 1.14 for Epoch 1 and 181/192 − 0.94 for Epoch 2.

To assess the significance level of the blackbody components, we performed spectral simulations using the same observational data using FAKE-IT command from XSPEC. Each fake spectra was constructed from arrays of randomly sampled parameters using the SIMPARS command. By generating a cumulative distribution function of F-values from the simulated data, we determined the minimum significance level of detection by comparing it with the F-value derived from the real data, as described in Hurkett et al. (2008).

Each F-value, for both the real and simulated data, is calculated as $F=({\nu }_0/\delta \nu )\times (\delta {\chi }^2/{\chi }_1^2)$, where the subscripts 0, 1 correspond to the null hypothesis (powerlaw) and the tested hypothesis (powerlaw+bbody). Each hypothesis is associated with a total χ² and ν degrees of freedom. The significance level is determined by computing the corresponding p value, which is the number of simulated spectra with F values greater than the F value derived from fitting the actual data. The uncertainty of this quantity can be calculated using $\sqrt{p(1-p)/Ns}$, where Ns is the total number of simulated spectra.

We ran ∼10⁵ simulations for both epochs and found that the blackbody component is significant at ∼2.6σ level for Epoch 1 and ∼2.8σ for Epoch 2.

From our analysis we conclude that the spectrum was nonthermal dominated during the observed period, possibly of synchrotron origin. This implies the presence of particles with TeV energies for the typical magnetic field strengths in this type of system. In the next section we explore this hypothesis in more detail.

In Fig. 3 we present the time averaged spectra and residuals associated with Model 1 and Model 2 of Epoch 1 (left panel) and Epoch 2 (right panel), while the corresponding best-fit parameters and uncertainties are detailed in Table 1. The powerlaw normalization component is expressed in units of photons keV⁻¹ cm⁻² s⁻¹ at 1 keV. The bbodyrad normalization is equal to $R_{\text{km}}^2/D_{10}^2$, where R_km is the source radius in km and D₁₀ is the distance to the blackbody source in units of 10 kpc. At a distance of 6.7 kpc, the size of the emitting region ranged from 1.5 to 5.4 km during Epoch 1 and from 2.9 to 8.9 km during Epoch 2. Alternatively, at a distance of 8.7 kpc, the size of the emitting region ranged from 2 to 7 km for Epoch 1 and from 4 to 11.6 km for Epoch 2. In Table 1, we report the values assuming a distance of 7.7 kpc.

Fig. 3. Spectral modeling results corresponding to Epoch 1 (left column) and Epoch 2 (right column) derived from simultaneous XMM-Newton and NuSTAR data (top panels). χ2 residuals correspond to an absorbed power-law (middle panel) and absorbed power-law with a black-body component (bottom panel).

To compute the unabsorbed flux, we used the convolution model cflux. Assuming a distance of 7.7 kpc (Corbet et al. 2019), the unabsorbed 2−20 keV luminosity ranges between 7.2 − 7.8 × 10³³ erg s⁻¹ for Epoch 1 and 4.4 − 5 × 10³³ erg s⁻¹ for Epoch 2.

4. Model

4.1. Jet nonthermal radiation

We adopted the hypothesis that a jet is present in the γ-ray-emitting binary 4FGL J1405.1−6119, and tried to evaluate whether it can adequately explain the nonthermal spectral energy distribution (SED) of the system. The radiative jet model used follows that presented in detail in Escobar et al. (2022), which in turn is based on Romero & Vila (2008). In the following, the model setup is outlined.

The scenario consists of a CO accreting material from the companion star with an accretion power L_acc, which can be expressed in terms of the Eddington luminosity as

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=q{L}_{\mathrm{Edd}}\approx q1.3\times {10}^{38}\left(\frac{M}{M_{\odot }}\right)\mathrm{erg}{\mathrm{s}}^{-1},\end{array}$$(2)

where q is a constant that represents the accretion efficiency in terms of the Eddington limit. Coupled with the inner accretion disk we assume the presence of a lepto-hadronic jet of kinetic luminosity, L_jet. This jet power relates to the accretion power through

$$\begin{array}{c}\hfill L_{\mathrm{jet}}=q_{\mathrm{jet}}{L}_{\mathrm{acc}},\end{array}$$(3)

where q_jet is another constant indicating the fraction of the accretion power that is transferred to the jet. The parameters q and q_jet define the accretion regime of the MQ. To compute the SED and the normalizations with the emission power, we chose to use L_jet directly as the parameter instead of a combination of q and q_jet; we defer discussion of the interpretation of the accretion regime to Sect. 5.

The jet propagates with a bulk velocity v_jet, corresponding to a bulk Lorentz factor Γ_jet. A fraction q_rel of this power is converted into relativistic particles by an acceleration mechanism. The relativistic proton and electron luminosities, L_p and L_e, are distributed according to the power ratio a = L_p/L_e.

We use cylindrical coordinates, with the coordinate z along the jet axis and the origin at the CO. The jet is started at a distance z₀ from the CO. The region in which the particles are accelerated extends from z_min to z_max, and its shape is described by

$$\begin{array}{c}\hfill r(z)=r_0{\left(\frac{z}{z_0}\right)}^{\epsilon },\end{array}$$(4)

where r₀ is the radius of the jet at its base and 0 < ε ≤ 1 describes its geometry. We note that the degree of collimation of the jet increases with decreasing values of ε, with the extreme value ε = 1 representing a conical shape. The magnetic field along the jet, B, decreases with z following a power law of index m, namely, B(z) = B₀(z/z₀)^−m, where 1 ≤ m ≤ 2 (e.g., Krolik 1999). The value of B₀ is obtained by assuming equipartition between magnetic and kinetic energy at the jet base.

We parameterize the injection function of energy to relativistic particles as a power law with an exponential cutoff,

$$\begin{array}{c}\hfill Q_{\mathrm{i}}=Q_{\mathrm{i},0}{E}_{\mathrm{i}}^{-\mathrm{p}}exp\left(\frac{E}{E_{\mathrm{i},\max }}\right),\end{array}$$(5)

where i = e, p accounts for electrons and protons, respectively; Q_i, 0 is obtained normalizing the injection function with the total power of each particle population, L_i; and E_i, max is the maximum reachable energy, achieved when the acceleration rate equals that of energy losses.

Relativistic particles are accelerated at a rate of

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{-1}=\frac{\eta ceB}{E_i},\end{array}$$(6)

where η ≤ 1 is the acceleration efficiency, c is the speed of light, and e is the elementary electric charge. On the other hand, relativistic particles lose energy via both radiative and non-radiative mechanisms. The latter are adiabatic and escape losses for both proton and electron populations. Regarding the radiative mechanisms, we consider synchrotron emission for both types of particles. In the case of protons, we also considered photons from the decay of neutral pions, which are the product of inelastic collisions between relativistic and cold protons; the latter consist mainly of protons in the jet bulk and those in the stellar wind of the companion star. In the case of electrons, we also computed the inverse Compton scattering of photons from the stellar radiation.

We assume that the relativistic particle populations reach a steady state. Their spectral densities are obtained by solving the steady state transport equation, taking injection, escape, and continuous losses into account. For a discussion of the general form of the transport equation, we refer the reader to Ginzburg & Syrovatskii (1964).

The complete formulae for calculating of radiation processes in this work can be found in Blumenthal & Gould (1970), Begelman et al. (1990), Atoyan & Dermer (2003), Kafexhiu et al. (2014), and references therein.

4.2. Spectral energy distribution of 4FGL J1405.1–6119

To fit the SED of the source with our model, we used the radio observations of Corbet et al. (2019), the X-ray flux obtained in this work, and the γ-ray flux from the 4FGL-DR2 catalog (Abdollahi et al. 2020; Ballet et al. 2020).

In Fig. 4 we show the acceleration, escape, and cooling rates of the relativistic particles at different locations in the jet. We find that for the electron population, the losses are dominated by synchrotron cooling throughout the emitting region. For the protons, both adiabatic and proton-proton collisions with bulk protons are the dominant mechanisms of energy loss in the region close to the base; the escape rate competes with them in regions further from the base and dominates the losses toward the end of the emitting region.

Fig. 4. Acceleration, escape, and cooling rates of electrons (top row) and protons (bottom row). Each column shows the aforementioned rates calculated at zmin (left column), at the logarithmic midpoint (middle column), and at zmax (right column).

In Fig. 5 we show the derived nonthermal SED of the jet, L_γ, which includes all the above radiative processes, where

$$\begin{array}{c}\hfill L_{\gamma }=E_{\gamma }^2\frac{\mathrm{d}N}{\mathrm{d}{E}_{\gamma }\mathrm{d}t},\end{array}$$(7)

Fig. 5. Nonthermal SED derived from our jet model. We considered the following radiative processes: for protons, synchrotron emission and proton-proton collisions with the cold protons of the jet bulk and the companion’s wind; for electrons, synchrotron emission and inverse Compton scattering off the radiation field of the companion. The figure also shows luminosities derived from XMM-Newton+NuSTAR data (our work) as well as ATCA and Fermi data taken from Corbet et al. (2019) and the 4FGL-DR2 catalog (Abdollahi et al. 2020; Ballet et al. 2020), respectively. All the references are in the figure.

and dN is the number of photons with energies between (E_γ, E_γ + dE_γ) emitted during a time dt. The assumed and derived parameters of the model, with uncertainties reported at the 90% confidence level for the free parameters, are listed in Table 2. The X-ray data points shown in Fig. 5 correspond to Epoch 1. Since there is no significant change in the X-ray flux between the two epochs (see Fig. 1), the same set of parameters also fits the observations of Epoch 2.

To obtain the reported model parameters, we ran a first set of ∼100 simulations² covering a wide range of values, and decided which of these remained fixed, apart from those resulting from observed properties of the system. Then we ran another set of 120 simulations considering the variation of all free parameters (i.e., using 40 simulations for each parameter variation at a time; m, p, ε, and η), and chose the model with the minimum χ²/d.o.f. The set of values of the free parameter space was covered with a uniform grid for p, ε, m, and log η. To estimate the free parameter errors, we computed the chi-squared for each simulation and chose the ${\chi }_{\min }^2+2.706$ contour for each free parameter, which represents the 90% credible region (e.g., Frodesen et al. 1979). The errors are reported in Table 2. In the particular case of m (η) the upper (lower) bound on the error comes from restricting the possible values to those reported in the fourth column of Table 2, while for the case of ε all the values within the bounds fall below the aforementioned contour. We obtained a minimum chi-squared goodness of fit of χ²/d.o.f. = 1.06, with 7 d.o.f.

We show that the broad spectrum of 4FGL J1405.1−6119 can be explained by the nonthermal emission associated with a mildly relativistic (Γ_jet = 1.9), lepto-hadronic model of a MQ jet. In particular, this model represents an approximately parabolic jet (ε ≈ 0.56), with a compact emitting region of size ≈1.0 × 10¹² cm (jet extension could be orders of magnitude larger), and a magnetic induction field at its base of B₀ ≈ 2.8 × 10⁷ G. Relativistic protons and electrons share the total power in the relativistic populations, L_rel = q_relL_jet = 10³⁷ erg s⁻¹, according to an assumed hadron-to-lepton ratio of a ≈ 0.11, and are accelerated through a low-efficiency mechanism for which η ≈ 1.0 × 10⁻⁴. The particle injection function shows a hard spectrum with spectral index p = 1.98.

In all cases, the values of the parameters are consistent with their commonly assumed range of values for MQ jets (see, for example, Romero & Vila 2008; Vila et al. 2012; Pepe et al. 2015; Escobar et al. 2021). The SED is computed assuming a viewing angle of θ = 30° (i.e., the angle between the jet axis and the line of sight). With the Lorentz factor from Table 2, and for lower viewing angles, the same data could be explained with a less powerful jet (∼0.07 times the one reported here), while higher angles would favor a scenario with more powerful jets, up to a factor of ∼50. Instead, maintaining fixed all the values of the parameters but the viewing angle, the model still manages to explain the observations with at most a variation of ≈10°. As we can see, knowing the viewing angle is a crucial factor in determining the accretion regime of the emitter.

5. Discussion and conclusions

In this work we have analyzed simultaneous XMM-Newton (2−10 keV) and NuSTAR (3−20 keV) observations of the γ-ray-emitting binary 4FGL J1405.1−6119. The difference between the two observed epochs is about six days. We found no significant pulsation above statistical noise in the XMM-Newton and NuSTAR light curves. This non-detection leaves open the question of the nature of the CO.

We fitted the time-averaged spectra of both epochs with empirical models, which yielded continuum parameters consistent with γ-ray-emitting binaries (Kretschmar et al. 2019; An et al. 2015; Yoneda et al. 2021) and MQs (Natalucci et al. 2014; Soria et al. 2020; Hirsch et al. 2020; Rodi et al. 2021; Saavedra et al. 2022). For Epoch 1, a power-law model adequately characterized the spectrum. For Epoch 2, however, the addition of a blackbody component significantly improved the fit to the observed data. The power-law index was found to be ∼1.7 for Epoch 1 and ∼1.4 for Epoch 2. Furthermore, the associated blackbody temperature for Epoch 2 was ∼0.8 keV, corresponding to a compact region of radius ≲16 km tentatively associated with the inner region of an accretion disk.

Corbet et al. (2019) present a different explanation for the emission of 4FGL J1405.1−6119. They set out a scenario involving the collision of winds between a rapidly rotating NS and its stellar companion. In addition, Xingxing et al. (2020) modeled this scenario to explain the GeV emission from 4FGL J1405.1−6119. In the colliding-wind scenario, the interaction between the stellar and pulsar winds gives rise to shocked regions that are characterized by a spiral shape due to Coriolis forces (see Molina & Bosch-Ramon 2020). Within this framework, the high-energy emission in a γ-ray binary system can be attributed to the up-scattering of photons from the stellar radiation field due to inverse Compton scattering by the relativistic electrons and positrons present in the shocked fluid. Synchrotron emission contributes to a lesser extent, especially in the energy range around ∼1 MeV (Molina & Bosch-Ramon 2020).

To explain the new X-ray data and archival multiwavelength observations (ATCA and Fermi), we implemented a lepto-hadronic jet model under the hypothesis that 4FGL J1405.1−6119 is a MQ. The parameters associated with the properties of the relativistic particle populations in the jet and the nonthermal emission are consistent with the values commonly assumed for this type of source (see, for example, Romero & Vila 2008; Vila et al. 2012; Pepe et al. 2015; Escobar et al. 2021, 2022). The scenario consists of a parabolic jet with a compact acceleration region, where the high-energy emission is produced by a hard spectrum of relativistic particles driven by a low-efficiency acceleration mechanism.

In this scenario, the γ-ray emission is produced via two processes. First, inverse Compton scattering produces nonthermal radiation in the energy range from about 100 keV to 10 GeV. Second, the decay of neutral pions from proton-proton collisions produces energetic photons at energies above 10 GeV. At lower energies, the synchrotron emission of electrons completely dominates, allowing us to accurately model radio and X-ray data. We thus propose that a lepto-hadronic jet model, which includes both leptonic and hadronic processes, may be sufficient to explain the multiwavelength emission of 4FGL J1405.1−6119.

There are essentially two interpretations of the results we obtained with our jet model. If the accretion regime is sub-Eddington, the CO should be a BH of at least ∼10 M_⊙. This scenario, like the case of Cygnus X-1, would require an X-ray luminosity on the order of ∼10³⁷ erg s⁻¹ in the low-hard state (Di Salvo et al. 2001; Makishima et al. 2008), which is in apparent contradiction with the observations. On the other hand, if the source is in a super-Eddington accretion regime, this would imply that the CO consists of a BH of a few solar masses or a NS with a weak magnetic field. In this case, the X-ray emission from the disk is absorbed by the photosphere of the wind ejected by the same disk (Abaroa et al. 2023). This is more consistent with what is observed in the X-ray emission. The supercritical source may also have an equatorial radio component or equatorial lines with velocities of 10³ − 10⁴ km s⁻¹, as in the case of SS 433, which could be observed in the future (Fabrika et al. 2021; Abaroa et al. 2023).

We note that leptonic jet models can also be adopted to explain observations of other MQs, as in the cases of Cyg X-3 (Zdziarski et al. 2012) and Cyg X-1 (Zdziarski et al. 2014). In our case, however, this type of model would not fully explain the observations for two main reasons. First, there is no clear evidence for disk or coronal emission. The X-ray observations could then be explained by a dominant synchrotron emission, which at these energies hides the radiation from the disk and/or the corona (see, for example, Fig. 11 of Bosch-Ramon et al. 2006). On the other hand, a relativistic proton component seems necessary to explain the observed slope change in the high-energy part of the Fermi data. In the lepto-hadronic picture, this part of the spectrum is dominated by the γ-ray emission, which originates in proton-proton collisions through neutral pion decays.

As shown in Fig. 2, the flux is higher in Epoch 1 than in Epoch 2. This behavior is also consistent with the presence of a thermal component in the spectra of Epoch 2, which can be explained by the contribution of a thermal inner disk and a reduced nonthermal jet component. In the case of a moderate or low viewing angle, some of the emission from the disk can escape from the central funnel in the wind (Abaroa et al. 2023). As for Epoch 1, the contribution to the total flux may be completely dominated by the jet.

The method we implemented for estimating parameters and uncertainties, although not very robust, allowed us to explain the multiwavelength behavior of 4FGL J1405.1−6119, favoring the MQ scenario. A more detailed analysis using a Markov chain Monte Carlo method would improve these estimates. This method would also allow a model comparison test to be included (e.g., to compute the odds ratio between leptonic and lepto-hadronic models, or between the MQ and colliding-wind scenarios). This analysis is beyond the scope of the current work, but we will present it in a future paper.

It is important to collect more observational data to evaluate the pros and cons of the pulsar-star collision wind and MQ scenarios for 4FGL J1405.1−6119. These additional data should be used to focus on timing analysis, especially with respect to the orbital period. By analyzing these data, we can gain a more complete understanding of the nature of this and other γ-ray-emitting binary systems.

Acknowledgments

We thank the anonymous reviewer for their valuable comments on this manuscript. We extend our gratitude to Fiona A. Harrison and Brian W. Grefenstette for their assistance with the instruments’ technical aspects. We thank Sergio Campana for his suggestions that helped improve this work. E.A.S., F.A.F. and J.A.C. acknowledge support by PIP 0113 (CONICET). F.A.F. is fellow of CONICET. J.A.C. is CONICET researchers. J.A.C. is a María Zambrano researcher fellow funded by the European Union – NextGenerationEU – (UJAR02MZ). This work received financial support from PICT-2017-2865 (ANPCyT). J.A.C. was also supported by grant PID2019-105510GB-C32/AEI/10.13039/501100011033 from the Agencia Estatal de Investigación of the Spanish Ministerio de Ciencia, Innovación y Universidades, and by Consejería de Economía, Innovación, Ciencia y Empleo of Junta de Andalucía as research group FQM-322, as well as FEDER funds. G.E.R. acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grants PID2019-105510GB-C31AEI/10.13039/501100011033/ and PID2022-136828NB-C41/AEI/10.13039/501100011033/, and by “ERDF A way of making Europe”, by the “European Union”, and through the “Unit of Excellence María de Maeztu 2020–2023” award to the Institute of Cosmos Sciences (CEX2019-000918-M). Additional support came from PIP 0554 (CONICET). G.J.E. acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017.

References

All Tables

All Figures

	Fig. 1. Background-corrected light curves of 4FGL J1405.1−6119 observed by NuSTAR (FPMA+B; red) and XMM-Newton (EPIC pn+MOS; blue) with a binning of 350 s, starting at 58711.705573 MJD. The first section of the light curve corresponds to the observation from August 18 (Epoch 1), while the second light curve is associated with the observation from August 25 (Epoch 2).
In the text

	Fig. 2. NuSTAR and XMM-Newton folded light curve using 48 phase bins, an orbital period of 13.713 days, and with 56498.7 MJD as the reference epoch (Corbet et al. 2019). A sine function fit is shown in orange (see the main text for details). The observed orbital modulation is similar to that shown by Corbet et al. (2019) using Swift/XRT data.
In the text

	Fig. 3. Spectral modeling results corresponding to Epoch 1 (left column) and Epoch 2 (right column) derived from simultaneous XMM-Newton and NuSTAR data (top panels). χ2 residuals correspond to an absorbed power-law (middle panel) and absorbed power-law with a black-body component (bottom panel).
In the text

	Fig. 4. Acceleration, escape, and cooling rates of electrons (top row) and protons (bottom row). Each column shows the aforementioned rates calculated at zmin (left column), at the logarithmic midpoint (middle column), and at zmax (right column).
In the text

Fig. 5. Nonthermal SED derived from our jet model. We considered the following radiative processes: for protons, synchrotron emission and proton-proton collisions with the cold protons of the jet bulk and the companion’s wind; for electrons, synchrotron emission and inverse Compton scattering off the radiation field of the companion. The figure also shows luminosities derived from XMM-Newton+NuSTAR data (our work) as well as ATCA and Fermi data taken from Corbet et al. (2019) and the 4FGL-DR2 catalog (Abdollahi et al. 2020; Ballet et al. 2020), respectively. All the references are in the figure.

In the text