© S. Barnier et al. 2021

1. Introduction

X-ray binaries (XrBs) are formidable laboratories for the study of the accretion-ejection processes around compact objects. Most of the time in a quiescent state they can suddenly enter an outburst that can last from a few months to a year, increasing their overall luminosity by several orders of magnitude. The X-ray emission is commonly believed to be produced by the inner regions of the accretion flow, whereas the radio emission is thought to originate from relativistic jets. Simultaneously with the X-ray spectral evolution along the outburst (from hard to soft states; e.g., Remillard & McClintock 2006; Done et al. 2007), the radio emission switches from jet-dominated states during the hard X-ray states, at the beginning and the end of the outburst, to jet-quenched states during the soft X-ray states, in the central part of the outburst (e.g., Corbel et al. 2004; Fender & Belloni 2004). These outbursts are usually represented in the so-called hardness-intensity diagram (HID) where they follow a typical “q” shape (see, e.g., Dunn et al. 2010).

A physical understanding has yet to be found to explain the complete behavior of these outbursts, even if a few points have reached consensus. For instance, the start of the outburst is believed to originate from disk instabilities in the outer regions of the accretion flow, driven by the ionization of hydrogen above a critical temperature (e.g., Hameury et al. 1998; Frank et al. 2002). The nature of the soft state X-ray emission, peaking in the soft X-rays, is commonly attributed to the presence of an accretion disk, down to the Innermost Stable Circular Orbit (ISCO), and the standard accretion disk model (SAD, Shakura & Sunyaev 1973) seems to describe the observed radiative output reasonably well. The exact nature of the hard X-ray emitting region, the hot corona, is however less clear. Given its intense luminosity, it is expected to be located close to the black hole where the release of gravitational power is the largest. The variability of the source is consistent with a very compact region (De Marco et al. 2017), but the exact geometry is still a matter of debate. It could be located somewhere above the black hole (the so-called lamppost geometry, e.g., Matt et al. 1991; Martocchia & Matt 1996; Miniutti & Fabian 2004). It can also partly cover the accretion disk (the patchy corona geometry, e.g., Haardt et al. 1997) or it can fill the inner part of the accretion flow, the accretion disk being present in the outer part of the flow (the corona-truncated disk geometry, e.g., Esin et al. 1997). The hot corona is most probably a combination of all these geometries; it could even evolve from one geometry to the other depending on the state of the source.

The dominant radiative process producing the hard X-rays is generally believed to be external Comptonization, meaning Comptonization of the external UV and soft X-ray photons produced by the accretion disk off the hot electrons present in the corona. These electrons are generally supposed to follow a relativistic thermal distribution to explain the presence of a high-energy cutoff, generally observed in the brightest hard states (see, e.g., Fabian et al. 2017 for a recent compilation). While the release of the gravitation power is undoubtedly the source of the corona heating, how this heating is transferred to the particles and how particles reach thermal equilibrium is still not understood. The magnetic field is expected to play a major role (e.g., Merloni & Fabian 2001), but the details of the process are unknown. The correlation between X-ray emission (from the corona) and radio emission (from the jet) (e.g., Gallo et al. 2003, 2012; Corbel et al. 2000, 2003, 2013; Coriat et al. 2011) also indicates a strong link between accretion and ejection, supporting the presence of a magnetic field in the disk.

Models of the hot corona commonly used in the literature are generally oversimplified. The geometry is assumed to have a basic shape, for instance spherical, slab, or even point-like. Its temperature and density are supposed to be uniform across the corona. External Comptonization is usually the unique radiative process taken into account and the spectral emission is often approximated by a cutoff power law (or similar) shape. Most of the time no physically motivated configuration or comparison with numerical simulations of this simplified model is proposed. More importantly, the jet emission and its impact on the accretion system is generally entirely ignored.

Magnetized accretion-ejection solutions that self-consistently treat both the accretion disk and the jets have been in development for more than 20 years (Ferreira & Pelletier 1995; Ferreira 1997) and have since been validated through numerical simulations (e.g., Zanni et al. 2007, Jacquemin-Ide et al. 2021). In these works the accretion disk is assumed to be threaded by a large-scale magnetic field. In the regions where the magnetization μ(r) = P_mag/P_tot is on the order of unity (with P_mag the magnetic pressure and P_tot the total pressure, i.e., the sum of the thermal and radiation pressure), the magnetic hoop stress overcomes both the outflow pressure gradient and the centrifugal forces, and self-confined non-relativistic jets can be produced. In these conditions the accretion disk is called a jet emitting disk (JED). The effect of the jets on the disk structure can be tremendous since the jets’ torque can efficiently extract the disk angular momentum, significantly increasing the accretion speed. Consequently, for a given accretion rate, a JED has a much lower density in comparison to the standard accretion disk (e.g., Ferreira et al. 2006). The parameter space for stationary JED solutions corresponds to magnetization μ in the range [0.1,1], small ejection index p < 0.1 defined by ṁ(r) ∝ r^p (with ṁ the mass accretion rate measured at a given radius r), large sonic mach number m_s = u_r/c_s in the range Arnaud (1996), Begelman & Armitage (2014), Bel et al. (2011) (with u_r the accretion speed and c_s the local speed of sound), and jet power fraction b = P_jets/P_acc between 0.1 and almost 1 for very thin JED (Ferreira 1997), where P_jets is the power feeding the jets and P_acc the total gravitational power released by the accretion flow within the JED. We note that for weak magnetization (μ ≪ 0.1), no collimation occurs and uncollimated winds are produced. The accretion disk structure is then not very different from the standard solution (e.g., Jacquemin-Ide et al. 2019).

Ferreira et al. (2006) proposed a hybrid disk paradigm to address the full accretion-ejection evolution of XrBs in outbursts (see also Petrucci et al. 2008). The accretion disk is threaded by a vertical magnetic field and extends all the way down to the innermost circular orbit. The outer part of the flow has a low magnetization, resulting in an outer SAD. On the contrary, the inner part has an important mid-plane magnetization and the accretion flow has a JED structure. This radial distribution of the magnetization appears to be a natural outcome of the presence of large-scale magnetic fields in the accretion flow, the magnetic flux accumulating toward the center to produce a magnetized disk with a fast accretion timescale (Scepi et al. 2020, Jacquemin-Ide et al. 2021).

Marcel et al. (2018a) developed a two-temperature plasma code to compute the spectral energy distribution (SED) of any JED-SAD configuration. In addition to the parameter μ, p, m_s, and b that characterizes a JED solution, the output SED also depends on the transition radius r_J (in units of gravitational radius R_G = GM/c²) between the JED and the SAD, as well as the mass accretion rate ṁ_in (in units of Eddington accretion mass rate Ṁ_Edd = L_Edd/c²) reaching the ISCO.

In our view, r_J and ṁ_in are expected to be physically linked through the evolution the magnetization across the accretion flow, dividing it in an inner strongly magnetized part (the JED) and an outer weakly magnetized one (the SAD). This link is far from being trivial to estimate and requires global 3D MHD simulations, which is far beyond of the scope of our present modeling, to catch it. In the absence of any physical law that could be used as an input, we consider the parameters r_J and ṁ_in to be independent parameters.

Varying r_J and ṁ_in, Marcel et al. (2018a) showed that the JED-SAD model was able to qualitatively reproduce the spectral evolution of an entire outburst. The JED radiative properties agree with the X-ray emission observed in compact objects (Petrucci et al. 2010; Marcel et al. 2018b,a), meaning that the JED can play the role of the hard X-ray emitting hot corona. At the beginning of the outburst the disk is characterized by a low-luminosity hard component that can be represented by a JED with large radial extend (r_J ≫ 1), while the mass accretion rate ṁ_in is low (ṁ_in < 0.1). During the rising part of the outburst the luminosity increases, thus the mass accretion rate increases, but r_J is still several times the ISCO radius. While transitioning to the soft states, r_J starts to decrease until the complete disappearance of the JED when r_J ∼ r_ISCO. This coincides with the disappearance of the radio emission (no more JED implies no more jets). During the soft states, the disk is dominated by the thermal component of the SAD and r_J stays equal to r_ISCO. Eventually, during the decaying phase of the outburst, a JED reappears when the system transitions back to the hard state, ṁ_in decreases, r_J increases again, and the system returns to the hybrid JED-SAD configuration. With the following decrease in the accretion rate, the XrB then fades to the quiescent state.

Marcel et al. (2019) (hereafter M19) performed the first application of the JED-SAD model to real data by qualitatively reproducing the spectral evolution of GX 339-4 during the 2010 outburst observed by RXTE. A similar study has been recently extended to three other outbursts of GX 339-4 (Marcel et al. 2020, hereafter M20). These authors did not directly fit the data given the large number of observations as well as the lack of a consistent reflection model component. Instead, they produced a large grid of spectra for a set of parameters (r_J, ṁ_in) and fitted each simulated spectrum with a disk plus power law (disk + power law) model. This provided spectral characteristics (disk flux, power law luminosity fraction, X-ray spectral index) that were then compared to the best fit results obtained by fitting the RXTE/PCA data with a similar disk + power law model (Clavel et al. 2016). This procedure allowed us to derive the qualitative evolution of r_J and ṁ_in that reproduces JED-SAD spectra with the closest spectral characteristics to the observed ones.

An important output of the JED-SAD model is the estimate of the jet’s power in a consistent way with the JED structure. The radio signature produced by this jet is less straightforward to estimate, however, since it depends on the detailed treatment of the jet particle emission all along the jet. Following Heinz & Sunyaev (2003) (hereafter HS03), M19 proposed an expression for the radio flux produced at a radio frequency ν_R by a jet launched from a JED characterized by r_J and ṁ_in :

$$\begin{array}{c}\hfill F_R={\stackrel{\sim }{f}}_R{\dot{m}}_{\mathrm{in}}^{17/12}{r}_{\mathrm{isco}}{(r_J-r_{\mathrm{isco}})}^{5/6}\frac{F_{\mathrm{Edd}}}{{\nu }_R}.\end{array}$$(1)

Here ${\stackrel{\sim }{f}}_R$ is a scaling factor and F_Edd = L_Edd/4πd² the Eddington flux. Equation (1) has the same dependency with the accretion rate as in the self-similar approach of HS03, but there is an additional multiplicative term (r_J − r_isco)^5/6 that reflects the necessarily finite radial dimension of the jet due to the finite radial dimension of the JED (see discussion in M19). By trying to simultaneously reproduce the X-ray and radio emission, M19 were able to put constraints on ṁ_in, r_J, and ${\stackrel{\sim }{f}}_R$ in the case of GX 339-4, assuming a constant ${\stackrel{\sim }{f}}_R$ for all the outbursts¹.

The present paper aims to make a step forward in the comparison of the JED-SAD model to real data through a direct fitting procedure of simultaneous radio and X-ray data of an X-ray binary. The improvements compared to the previous works are twofold. First we added a consistent reflection component in the model and second we obtained more reliable and precise constraints on our model parameters (r_J and ṁ_in). To do so, we developed the required tools to apply our JED-SAD model to standard fitting software (e.g., XSPEC, Arnaud 1996).

We focus in this paper on the simultaneous radio–X-ray coverage of the XrB GX 339-4 during the lifetime of the RXTE satellite (1995–2012). This corresponds to four major outbursts starting in 2002, 2004, 2007, and 2010. The data selection is discussed more precisely in Sect. 2. Our fitting procedure and first fit results, using Eq. (1) for the radio emission, are discussed in Sect. 3. These results suggest, however, a different functional dependency of the radio emission with r_J and ṁ_in compared to Eq. (1). A deeper analysis of the radio behavior is then performed in Sect. 4 and supports two different functional behaviors of the radio emission between the beginning and the end of the outburst. The implications of these results are discussed in Sect. 5, and we conclude in Sect. 6.

2. Data selection

To test our JED-SAD paradigm we focused on simultaneous or quasi-simultaneous radio–X-ray observations of GX 339-4. We only used pure hard states (i.e., those at the very right part of the HID), either in the rising or decaying phase, and did not include the transition phases of the outburst even when radio emission was detected (during the so-called hard intermediate state, HIS). The reasons for this choice are twofold. First, the radio flux evolves smoothly during pure hard states, a signature of stationary processes hopefully easier to catch. Conversely, an important radio variability is observed during the transition phases, especially during the hard-to-soft transition. Second, during the transition states, the hard-tail component progressively appears. As this component is not well understood and is not self-consistently included in the JED-SAD model, we did not select the transition states.

We selected X-ray spectra from the RXTE-PCA archive of GX 339-4 during the 2000–2010 decade. In order to have a uniform data analysis, we did not include the data from the RXTE-HEXTE instrument since they were not always usable (e.g., in the case of low flux observations or after March 2010 when it stopped observing). The data processing is detailed in Clavel et al. (2016). Since the instrumental background was generally found to be on the order of, or larger than, the source emission above 25 keV, we limited our spectral analysis to the 3–25 keV energy range of the PCA instrument. We plot the 2000–2010 PCA X-ray light curve of GX 339-4 in Fig. 1. During this period, GX 339-4 underwent four complete outbursts, in 2002, 2004, 2007, and 2010, hereafter outbursts 1, 2, 3, and 4. The hard-only or failed outbursts of 2006 and 2008 were not selected for this study² since they may be intrinsically different from the ones accomplishing an entire HID. We follow Clavel et al. (2016) for the definition of the hard-state periods of each outburst. We list in Table 1 the corresponding starting and ending Modified Julian Dates (MJD) of both the rising and the decaying hard-state phases.

Fig. 1. GX 339-4 X-ray light curve in the 3–200 keV energy band of the 2000–2010 decade obtained with the Clavel et al. (2016) fits. The violet filled region shows the power law unabsorbed flux, while the cyan region represents the disk unabsorbed flux. The selected spectra for this study are highlighted in gray: the rising and decaying hard states of the four outbursts (1, 2, 3, and 4). At the top the red lines represent the date when steady radio fluxes were observed at 9 GHz (from Corbel et al. 2013).

In radio we used the 9 GHz fluxes obtained with the Australia Telescope Compact Array (ATCA) and discussed in Corbel et al. (2013)³. Compared to the X-ray observations, the radio survey is quite sparse (see Fig. 1), so we selected only the radio fluxes close to X-ray pointings by less than one day (which we call quasi-simultaneous radio–X-ray observations).

This selection corresponds to a total of 452 hard X-ray spectra and 57 radio fluxes distributed among the four outbursts. Outburst 4 is the one with the best X-ray and radio coverage, with about 80 X-ray spectra and 24 radio measurements that are evenly distributed along the outburst. Thanks to this large radio coverage we chose to linearly interpolate the radio light curve to estimate the radio fluxes for each of the 80 X-ray spectra of this outburst. This is supported by the smooth evolution of the radio light curve during the pure hard states. The resulting interpolation is plotted in Fig. A.1. This interpolation was not possible for the other outbursts, due to the insufficient number of radio pointings.

3. X-ray and radio fits

3.1. Methodology

Similarly to M19, we assume a distance d = 8 kpc for GX 339-4 (Zdziarski et al. 2004, 2019; Parker et al. 2016). The black hole mass is estimated between 4 and 11 M_⊙ (Parker et al. 2016; Zdziarski et al. 2019), and we assume a black hole mass of 10 M_⊙. The innermost stable circular orbit r_isco was assumed equal to 2 in R_g units. This is equivalent to a black hole spin of 0.94 (Miller et al. 2008; García et al. 2015). Finally, for the Galactic hydrogen column density we used 0.6 × 10²² cm⁻² (Zdziarski et al. 2004; Bel et al. 2011).

In the JED-SAD model the two parameters left free to vary during the fitting procedure are r_J and ṁ_in. All the other parameters of the JED-SAD (see Sect. 1) were set to the same values as in M19: b = 0.3, m_s = 1.5, and p = 0.01⁴.

We created XSPEC model tables for the JED and the SAD components separately with 40 values of ṁ_in and 25 values of r_J logarithmicaly distributed in the range [0.001,10] and [1,300], respectively. We also produced a reflection table. For this purpose, we used the XILLVER reflection model (Garcia et al. 2013). For each couple (r_J, ṁ_in) of the JED table, we fitted the corresponding JED spectrum with a cutoff power law model. This fit provides a spectral index and a high-energy cutoff, which we injected in the XILLVER table to produce different reflection spectra for different values of the disk ionization log(ξ) and iron abundance A(Fe) (in solar units). The disk inclination was set to 30°, an inclination consistent with the value expected for GX 339-4 (Parker et al. 2016). The resulting table thus possesses five different parameters for each spectrum: the two JED-SAD parameters (r_J, ṁ_in), and the three reflection parameters log(ξ), A(Fe), and the reflection normalization.

We then used an automatic fitting procedure using the pyxspec library (a python interface to XSPEC). We fit the X-ray spectra with the following XSPEC model: tbabs * (atable(JEDtable) + atable(SADtable) + kdblur*atable(Refltable)). Here JEDtable, SADtable, and Refltable are the XSPEC tables for the JED, SAD, and reflection spectra, respectively, and KDBLUR is a convolution model of XSPEC to take into account the relativistic effects from the accretion disk around a rotating black hole (according to the original calculations by Laor 1991). The parameters r_J and ṁ_in are linked between the tables and we set the inner radius of KDBLUR to the inner radius of the SAD (i.e., r_J). In KDBLUR, we set the index of the disk emissivity to 3 (its default value), the outer disk radius to 400 R_g, and the inclination to 30°.

3.2. X-ray fits

Using the fitting procedure described above, we obtained the best fit values for r_J and ṁ_in for each X-ray observation in our sample. The iron abundance values are clustered around seven times the solar abundance, in agreement with similar spectral analysis of GX 339-4 (e.g., García et al. 2015; Fürst et al. 2015; Parker et al. 2016; Wang-Ji et al. 2018), and we set it to this value in the following⁵. As examples, we show in Fig. 2 a few of our X-ray best fits obtained for different observations distributed in the hard X-ray states of outburst 4. In the top left panel of this figure we represent the HID as well as the hard states (blue diamond) that we fit. We also highlight the five observations whose spectral fits are presented in the other panels of the figure.

Fig. 2. Best fits of some observations of outburst 4. Top left: hardness intensity diagram of outburst 4. The blue diamond shows the hard state used for this outburst. The five red points are the five observations plotted in the different figures from a to e. a–e: best fit spectra and the data/model ratio for the five observations indicated in red in the HID. The gray region shows the PCA energy range used for the fit. The data are in black and the best fit model in gray; the JED spectrum is in red, the SAD spectrum in green, and the reflection component in blue. The best fit parameters for each observation are listed in Table 2.

During the rising phase (observations a, b, and c), the high-energy cutoff slowly appears in the model with the rise in luminosity⁶. At the same time, the iron line changes shape under the influence of both the evolution of the disk ionization parameter and the black hole gravity as the transition radius r_J decreases (general relativity effects). During the decaying phase (observations d and e), as the luminosity decreases, the standard accretion disk component disappears with the increase in r_J.

The evolutions of r_J and ṁ_in for our entire data sample are shown in Fig. 3; the left panel shows the light curves of r_J and the right panel those of ṁ_in. We divided each panel into two parts, showing the rising phase first and then the decaying phase. The large green region represents the area where the minimization function used by M19 to constrain r_J and ṁ_in varies by less than 10% with respect to its minimum. The blue points with black error bars and connected by dashed lines represent the results of this paper obtained by fitting the X-ray spectra in XSPEC.

Fig. 3. Results of the fitting procedure. Left: transition radius rJ (from risco to 300) between the JED and the SAD. Right: mass accretion rate ṁin. Each side is divided vertically between the four outbursts and horizontally between the rising and decaying phase of each outburst. The green solid line represent the results from M19 and M20, and the green region indicates where the minimization function varies by less than 10% with respect to its minimum. The blue dashed line shows the results of the fitting procedure and the black vertical bar the associated 90% confidence range. The decaying phase of outburst 4 is the subject of Appendix B.

We obtained much tighter constraints compared to M19, especially for r_J during the rising phase. There are two reasons for this. First, M19 did not directly fit the data, their main objectives being to qualitatively reproduce the outburst spectral and flux evolution. Second, M19 did not use a χ² statistic to constrain their parameters because the χ² statistic is not well adapted to their methodology.

Nevertheless, the constraints obtained with our fitting procedure are almost always embedded within the green area obtained by M19, which shows the good agreement between the two approaches. This is noticeably the case for ṁ_in, which is well constrained in the two methods and are in very good agreement with each other. Interestingly, our values for r_J are apparently better constrained in the rising phase of the outbursts, its behavior being more erratic and with larger error bars in the decaying phase. This could be a natural effect of the decrease in the data statistics when the flux decreases, but this trend is not observed for ṁ_in. This instead suggests that our JED-SAD spectra are less dependent on r_J at low accretion rates.

It should be noted that the χ² space does not always follow a Gaussian shape (see Fig. B.2) and the error bars should not be taken as sigma errors, but instead as lower and upper limits with a 90% confidence. Thus, none of the pure hard states are consistent with r_isco, and a JED is always required in the fit.

The evolution of r_J during the decaying phase of outburst 4 is the subject of Appendix B where we detail how we obtained the presented values of r_J using a maximum likelihood method. When a small r_J solution was found in the automatic procedure (r_J < 10), we checked the parameter space for a statistically equivalent solution for a larger value of r_J. Whenever such a solution was found, we selected it (see Appendix B). The motivations for this choice are twofold: higher values of r_J are observed in the decaying phase of the other outbursts (see Fig. 3, outbursts 2 and 3) and the resulting increase in r_J when going to quiescence is consistent with the JED-SAD dynamical picture.

3.3. Taking into account the radio emission

We now reproduce the radio fluxes using Eq. (1) to model the radio emission. Since the radio survey is generally quite sparse, we first concentrate on the rising phase of outburst 4, where the radio coverage is sufficiently dense to interpolate the radio fluxes for all X-ray spectra (see Fig. A.1). We use the results of the X-ray fits (see previous subsection) and set the parameters r_J and ṁ_in to the best fit values. Then we compute the radio flux F_R with Eq. (1) using $\stackrel{\sim }{f}=1.5\times {10}^{-10}$, the value used in M19. In Fig. 4 we plot the ratio of the observed (and interpolated) radio flux F_obs to the expected radio flux F_R from Eq. (1) as function of r_J and ṁ_in. A clear anticorrelation is observed ${\displaystyle \frac{F_{\mathrm{obs}}}{F_R}\propto r_J^{\alpha }}$ with α ∼ −1.25. Similarly, ${\displaystyle \frac{F_{\mathrm{obs}}}{F_R}}$ is correlated with ṁ_in, with a power β ∼ 1.56. In conclusion here our fitting procedure suggests a functional dependency of the radio emission at least on r_J and/or ṁ_in that is not taken into account correctly when using Eq. (1). This is studied in more detail in Sect. 4.

Fig. 4. Ratio of the observed radio fluxes to the results of Eq. (1) (using $\stackrel{\sim }{f}=1.5\times {10}^{10}$) as a function of rJ (left) and ṁin (right) for outburst 4. In the case of X-ray observations without simultaneous radio measurement, the radio flux was interpolated from the radio light curve (see Fig. A.1). The dashed line shows the best fit power law: $\frac{F_{\mathrm{obs}}}{F_R}\propto$$r_J^{-1.25}$ (left) and $\frac{F_{\mathrm{obs}}}{F_R}\propto$ ${\dot{m}}_{\mathrm{in}}^{1.56}$ (right).

4. Functional dependency of the radio emission

In the JED-SAD paradigm, the evolution of the X-ray spectrum (hardness, energy cutoff, and flux) is described through the changes of two parameters, r_J and ṁ_in, controlling the balance between the power released through advection and radiation. In a similar way, we use both of these parameters to describe the radio flux.

We thus assume in this section a more general expression for the radio flux:

$$\begin{array}{c}\hfill F_R={\stackrel{\sim }{f}}^{\ast }{r}_J^{\alpha }{\dot{m}}_{\mathrm{in}}^{\beta }(1-\frac{r_{\mathrm{isco}}}{r_J}{)}^{5/6}\frac{F_{\mathrm{Edd}}}{{\nu }_R}·\end{array}$$(2)

This expression is similar to Eq. (1), but the indexes of the dependency on both ṁ_in and r_J are now free parameters. This new expression allows us to put all the dependency on r_J and ṁ_in in the parameters α and β, ${\stackrel{\sim }{f}}^{\ast }$ acting then as a true constant in this respect. The term (r_J − r_isco)^5/6 linked to the radial extension of the jet is re-expressed to isolate the dominant power dependency with r_J in α. We look for a unique triplet (${\stackrel{\sim }{f}}^{\ast }$, α, β) that could reproduce the whole radio data set.

4.1. Rising phase of the 2010 outburst

We first test Eq. (2) in the rising phase of outburst 4. We fit all 16 radio observations with simultaneous or quasi-simultaneous radio–X-ray data in XSPEC. We set r_J and ṁ_in of each observation to the best quasi-simultaneous X-ray fit values (obtained in Sect. 3). When simultaneously fitting X-ray and radio data, if the JED-SAD and reflections parameters are left free to vary simultaneously to ${\stackrel{\sim }{f}}^{\ast }$, α, and β, the X-ray fit is found to be significantly worse, especially around the iron line, for the benefit of a perfect match of the radio fluxes. By freezing the JED-SAD and reflections parameters to their best fit values obtained by fitting the X-rays, we instead chose to favor the X-ray fit for which we have a fully developed physically motivated spectral model. We implement in XSPEC a model to fit the radio emission following Eq. (2). We impose the same value of ${\stackrel{\sim }{f}}^{\ast }$, α, and β for all the observations. The number of radio fluxes we use can be found in Table 1.

As the errors on the radio are sometimes quite small, we introduce a 10% systematic error on the radio fluxes to account for the non-simultaneity between the radio and X-ray observations and the few percent radio intrinsic variability (Corbel et al. 2000). This is done so that the fit is not driven by only one radio flux, but tries to reproduce all the fluxes within this 10% error margin. Furthermore, the maximum variation observed in the radio light curve is about 20% variation over three days (see Fig. A.1). Thus, within the one day delay between the radio and X-ray observations we do not expect variations exceeding the 10% systematic error we add, justifying the use of not exactly simultaneous X-ray and radio pointings. The effects of adding systematic errors is discussed in Appendix C.

The best fit gives ${\stackrel{\sim }{f}}^{\ast }=7.1_{-5}^{+15}\times {10}^{-8}$, α = −0.66 ± 0.32, and $\beta =1.{00}_{-0.38}^{+0.39}$. The contours α–β are also reported as thin blue solid lines in Fig. 5. The fit reproduces all the radio fluxes within an error lower than 10% (see examples of residuals in Fig. 6) suggesting that Eq. (2) works adequately. The positive value of β is consistent with the observed correlation between the radio emission and the luminosity of the binary system. The negative value of α agrees with a decrease in the inner radius of the SAD when the system reaches bright hard states with stronger radio emission as expected in our JED-SAD approach (and similarly to most of the truncated disk models like Esin et al. 1997).

Fig. 5. Contour plots β–α for the rising (blue) and decaying (red) phases of the outburst of 2004 (outburst 2, dashed line), 2007 (outburst 3, dot-dashed line), and 2010 (outburst 4, thin solid line). The contours in thick solid lines represent the dependency when fitting all rising (blue) or decaying (red) phase radio fluxes simultaneously. Confidence contour levels correspond to 68%, 90%, and 99% (Δχ2 of 2.3, 4.61, and 9.2, respectively). The contours are obtained when fitting only the quasi-simultaneous radio–X-ray observations (not the interpolated radio observations).

Fig. 6. Ratio of the data to the model for the best fit of 5 of the 16 multiwavelength observations (radio–X-ray) of the rising phase of the 2010 outburst (MJD 55217 in blue, 55259 in red, 55271 in black, 55288 in green, and 55292 in violet). Only five ratios are shown for purposes of visualization, but the best fit was obtained by using all the simultaneous or quasi-simultaneous radio–X-ray observation, fixing the JED-SAD parameters to the best fitting values obtained by fitting the X-ray spectra first, and then fitting the radio points with Eq. (2).

In the second step we apply the same procedure to all the interpolated radio fluxes of the rising phase of outburst 4. Following the first step, we set r_J and ṁ_in to their best fit values obtained when fitting the X-ray alone. Then we reproduce the radio using the best fit values of ${\stackrel{\sim }{f}}^{\ast }$, α and β obtained previously to compute the expected radio flux F_R using Eq. (2). The corresponding ratios F_obs/F_R are reported in Fig. 7. There is almost no remaining dependency on r_J or ṁ_in. Compared to Fig. 4, this now shows a much more clustered distribution around 1, with a dispersion of about ±15%.

Fig. 7. Ratio of all radio fluxes of the rising phase of outburst 4 to the results of Eq. (2) as a function of rJ (left) and ṁin (right). The filled blue points represent the 16 quasi-simultaneous radio–X-ray observations, while the empty points represent the interpolated radio fluxes (see Appendix A). The radio observations are well reproduced using the values ${\stackrel{\sim }{f}}^{\ast }=7.1\times {10}^{-8}$, α = −0.66, and β = 1.00.

4.2. Decaying phase of the 2010 outburst

For the decaying phase of outburst 4 we proceeded similarly to the rising phase. We chose all eight observations, with simultaneous or quasi-simultaneous (differences less than 1 day) radio–X-ray observations of the decaying phase of outburst 4. We set the values of r_J and ṁ_in to the best X-ray fit, then we fit the radio fluxes using Eq. (2).

As a first test, we set ${\stackrel{\sim }{f}}^{\ast }$, α, and β to the best fit values obtained in the rising phase. The corresponding data/model ratio is plotted in Fig. 8. The top panel shows that using the value of ${\stackrel{\sim }{f}}^{\ast }$ of the rising phase in the decaying phase induces an error in the radio flux of up to a factor of 5. The bottom panel shows that even if we leave the scaling factor ${\stackrel{\sim }{f}}^{\ast }$ free, converging to the value 2.0 × 10⁻⁷, the radio flux is incorrect by a factor up to 1.8. Thus, the parameters ${\stackrel{\sim }{f}}^{\ast }$, α, and β cannot have the same values as obtained in the rising phase.

Fig. 8. Ratio of the data to the model for the best fit of five of the eight multiwavelength observations of the decaying phase of outburst 4 (MJD 55613 in blue, 55617 in red, 55620 in black, 55630 in green, and 55639 in violet). Top panel: all fits were done simultaneously, fixing the parameters of Eq. (2) to those found in the rising phase: ${\stackrel{\sim }{f}}^{\ast }=7.1\times {10}^{-8}$, α = −0.66, and β = 1.00. Bottom panel: fixing α and β to the values found for the rising phase, ${\stackrel{\sim }{f}}^{\ast }$ is free to vary and converges to the value 2.0 × 10−7.

Following what we present for the rising phase, we now leave ${\stackrel{\sim }{f}}^{\ast }$, α, and β free to vary, but linked between all the observations. The best fit values are ${\stackrel{\sim }{f}}^{\ast }=2.9_{-0.9}^{+1.0}\times {10}^{-8}$, $\alpha =-0.{13}_{-0.16}^{+0.15}$, and β = 1.02 ± 0.23. The corresponding confidence contour α–β is plotted as red thin solid lines in Fig. 5. It is clearly inconsistent with the blue contour obtained in the rising phase.

4.3. Comparison with the other outbursts

We constrain the functional dependency of the radio emission of the other outbursts by repeating a similar analysis. We thus need at least three observations taken in the corresponding rising and decaying phases to constrain the three free parameters α, β, and ${\stackrel{\sim }{f}}^{\ast }$. Only outburst 2 (in 2004) and the decaying phase of outburst 3 (in 2007) have the sufficient number of simultaneous or quasi-simultaneous radio and X-ray observations to apply our procedure. The number of radio fluxes we use for each phase of the outbursts can be found in Table 1. The corresponding contour plots of α–β are overplotted in Fig. 5 as dashed and dot-dashed lines, respectively. Two results are remarkable. First, and similarly to outburst 4, we need different functional dependencies of the radio emission with r_J and ṁ_in between the rising and decaying phase for outburst 2. Even more interestingly, the values obtained for α and β are in quite good agreement between the different outbursts, the contour of the decaying phase of outburst 3 also close to the contours of the decaying phases of outbursts 2 and 4. While this could be surprising given the simple expression used to model the radio emission, we believe that this result reveals intrinsic differences in the jet emission origin (see Sect. 5 for this discussion).

In the last step we use all the quasi-simultaneous radio observations, simultaneously fitting all the rising phase observations together with the same parameters α and β for all outbursts, but with different normalization ${\stackrel{\sim }{f}}^{\ast }$ for each outburst. We did the same for all the decaying phase observations. The resulting α–β contours are plotted in Fig. 5 as thick solid lines.

The process confirms the two different and mutually inconsistent functional dependencies of the radio emission on r_J and ṁ_in between the rising and decaying phases observations. The radio flux observed in the rising phases is nicely reproduced (within about 15%) by the relation

$$\begin{array}{c}\hfill F_R^{\mathrm{rise}}\propto r_J^{{-0.67}_{-0.22}^{+0.21}}{\dot{m}}_{\mathrm{in}}^{{0.94}_{-0.24}^{+0.25}},\end{array}$$(3)

instead, in decaying phases it follows

$$\begin{array}{c}\hfill F_R^{\mathrm{decay}}\propto r_J^{-0.15\pm 0.06}{\dot{m}}_{\mathrm{in}}^{0.9\pm 0.1},\end{array}$$(4)

with a weaker dependency on r_J.

Some variations in ${\stackrel{\sim }{f}}^{\ast }$ are required, however, to significantly improve the radio emission modeling. This can be seen in Fig. 9 where we report the ratio F_obs/F_R using Eq. (3) to compute the radio flux if the observation is in the rising phase and Eq. (4) if in the decaying phase. In the top panel we use the same value ${\stackrel{\sim }{f}}^{\ast }$ for all outbursts. The ratios cluster around 1, although there is some scattering between the different phases of the different outbursts. We report in the bottom panel of Fig. 9 the same ratio but letting ${\stackrel{\sim }{f}}^{\ast }$ free to vary between outbursts and between the rising and decaying phases. The improvement is clear and almost all radio fluxes can be reproduced within a 20 % margin error. The different values of ${\stackrel{\sim }{f}}^{\ast }$ found are reported in Table 3. We observe variation up to a factor of three (e.g., between the rising phase of the 2002 and 2010 outbursts). This could be related to local changes in the radiative efficiency of the radio emission from outburst to outburst.

Fig. 9. Ratios of the observed radio fluxes to the modeled radio fluxes for all the quasi-simultaneous radio–X-ray observations of the four outbursts. Each outburst is represented using a different symbol (see legend). In blue are shown the rising phases and in red the decaying phases. The modeled radio fluxes were obtained using Eq. (2). The parameters are (α = −0.67, β = 0.94) for the rising phases and (α = −0.15, β = 0.9) for the decaying phases. Top panel: ${\stackrel{\sim }{f}}^{\ast }=4.1\times {10}^{-8}$ used for all outbursts. Bottom panel: different ${\stackrel{\sim }{f}}^{\ast }$ used for each phase of the outbursts. All values of ${\stackrel{\sim }{f}}^{\ast }$ used are reported in Table 3. The horizontal dashed lines represent a 20% error margin (ratio of 0.8 and 1.2, respectively).

5. Discussion

We present in this paper the first X-ray spectral fits of an X-ray binary using the JED-SAD model. Compared to previous works, we constructed model tables that enable the use of XSPEC for a direct fit procedure. We also constructed a reflection table, based on RELXILL, using as inputs the photon index and high-energy cutoff that best fit the JED-SAD spectral shapes. We obtained good fits for all the X-ray observations of GX-339-4 during the hard states observed by RXTE on the period 2002–2010 by only varying the accretion rate ṁ_in of our system and the transition radius r_J between the JED and the SAD.

As said in the Introduction, in the absence of a known physical law that would link these two parameters, we left them free to vary independently of each other in the fit procedure. This is the simplest approach to try to understand, and hopefully to physically interpret (see Ferreira et al., in prep.), their behavior.

Then radio emissions simultaneous to the X-rays (or quasi-simultaneous, with a one-day difference) were reproduced using a generic formula only depending on ṁ_in and r_J. One of the main results of this spectral analysis is the necessity of a different functional dependency of the radio emission on r_J and ṁ_in. The two different expressions for the radio emission are reported in Eqs. (3) and (4). We believe that this difference in functional dependency is a “back product” of the true physical link between r_J and ṁ_in.

5.1. Indications of different radiative behaviors between the rising and decaying phases

Observational clues on the jet behavior can be derived from a set of different diagnostics: (i) the radio spectral index α_R, (ii) the measure of the spectral break frequency ν_break, (iii) timing properties, (iv) the correlation L_R − L_X, and (v) linking the radio luminosity L_R to disk properties (ṁ_in,r_J). Items (i)–(iv) are discussed in this section, while item (v) is discussed in Sect. 5.2.

The radio spectral index α_R can be analytically derived under the assumption of a self-absorbed synchrotron emission smoothly distributed along the jet. It depends on the particle distribution function, the jet geometry, and the way the dominant magnetic field varies with the distance (see Eq. (A.8) in the Appendix of Marcel et al. 2018b). There is no a priori reason to assume that these parameters should not vary in time. Observationally, however, there is no clear evidence of differences in the radio spectral index α_R between the rising and decaying phases of GX 339-4 (Espinasse, priv. comm.; see also Koljonen & Russell 2019 for more detailed discussion on this point and Tremou et al. 2020 for the quiescent state case where the radio spectrum is clearly inverted). Although this is already an important piece of information, we note that these α_R are derived within a rather limited radio band and might therefore not be fully representative of the whole jet spectrum (see, e.g., Péault et al. 2019).

The evolution of the spectral break frequency, ν_break, marking the transition from self-absorbed to optically thin jet synchrotron radiation, could however be different in the two phases (rising and decaying). The radio spectral index being flat or inverted in the hard state, the power of the jets is mainly sensitive to the position of the spectral break. Gandhi et al. (2011) measured this break at ∼5 × 10¹³ Hz in a bright hard state during the rise of the 2010–2011 outburst. By comparison, Corbel et al. (2013) constrain the break to be at lower frequency in the decaying phase, suggesting a less powerful jet in this phase. There is also a potential link between the X-ray hardness and the jet spectral break frequency, harder X-ray spectra having a higher ν_break (Russell et al. 2014; Koljonen et al. 2015). Interestingly, GX 339-4 shows on average a softer power law index in the decaying phase compared to the rising phase (see Fig. 10). Given the observed correlation between ν_break and the X-ray hardness, this also suggests a different behavior for ν_break (and consequently of the jet power) between the two phases.

Fig. 10. Evolution of the hard X-ray power law index Γ during the four outbursts of GX339-4. Left: hard X-ray power law index Γ as a function of the 3–9 keV X-ray luminosity of the pure hard state observations during the four GX339-4 outbursts (data from Clavel et al. 2016). In blue the rising phase and in red the decaying phase. Right: histograms of Γ. These distributions are subject to a certain number of observational biases: inclusion of error bars and the number of observations per phase.

There are other indications that the accretion (through X-ray emission) and ejection (through radio emission) processes could behave differently at the beginning and the end of the outburst. At first sight the radio–X-ray correlation followed by GX339-4 agrees with a linear correlation of index ∼0.7 in log-log space (e.g., Corbel et al. 2000, 2003, 2013) even down to very quiescent states (Tremou et al. 2020), but a more careful analysis shows the presence of wiggles along this linear correlation, especially between the high- and low-luminosity states (e.g., Corbel et al. 2013, Fig. 8). When looking more precisely at the rising and decaying phase, two different correlations may even be observed (Islam & Zdziarski 2018).

These differences may be linked to a change in the radiative efficiency of the X-ray corona with luminosity. The low X-ray luminosity states, below 2–20% of the Eddington luminosity, are potentially less radiatively-efficient than the high X-ray luminosity states (Koljonen & Russell 2019; Marcel et al., in prep.). As noticed by Koljonen & Russell (2019), these changes of the accretion flow properties could affect the jet launching, and therefore its radio emission properties.

In the JED-SAD model, the accretion power available in the accretion flow, ${\displaystyle P_{\mathrm{acc}}=\frac{GM\dot{M}}{2R_{\mathrm{isco}}}[1-{\left(\frac{r_{\mathrm{isco}}}{r_J}\right)}^{1-p}]}$ (see Sect. 1 for the definition of p and b), is released in three different forms: advection, radiation, and ejection. The first two occur inside the JED, and their sum is defined as P_JED = (1 − b) P_acc. The ejection power is released in the jets and is defined as P_jets = b P_acc. We also define the ratio η_R = L_R/P_jets and the ratio η_X = L_{3 − 9 keV}/P_JED that can be respectively interpreted as the radiative efficiency in the radio and X-ray bands. We report in Fig. 11 the ratio η_R⁷ as a function of the ratio η_X for the rising phase (blue points) and the decaying phase (red points) of the outbursts. We highlight in Fig. 11 the observations (labeled a to e) presented in Fig. 2 to mark the chronological evolution along an outburst. Figure 11 mostly depends on the well-constrained mass accretion rate obtained with our fits of each X-ray observations.

Fig. 11. Radio emission efficiency (LR = L9 GHz/PJets) vs. X-ray emission efficiency (LX = L3 − 9 keV/PJED) during the outbursts of GX339-4. The blue points are the rising phases. The red points are the decaying phases. The symbols distinguish the different outbursts: diamonds for 2002, squares for 2004, dots for 2007, and triangles for 2010 (filled for the quasi-simultaneous observations and empty for the interpolated radio fluxes). We highlighted the five observations (labeled a to e) presented in Fig. 2 to provide the chronological evolution of an outburst.

The blue points of the rising phases follow a similar trend for all the outbursts with a change in the X-ray and radio radiative efficiency by a factor of ∼4 and ∼2, respectively. In the decaying phase, however, each outburst clusters at the same radio and X-ray radiative efficiency. Interestingly, the radio radiative efficiency changes from outburst to outburst, while the X-ray radiative efficiency stays roughly constant at η_X ∼ 3 − 4 × 10⁻², the lowest values observed in the rising phase. These results indeed suggest a change in the radiative properties of the accretion-ejection structure between the beginning and the end of the outburst, and it is possible that it has some impact on the functional dependency of the radio emission highlighted in this paper. Contrary to the conclusion of Koljonen & Russell (2019), however, the accretion rate does not seem to be the only parameter that controls the evolution of η_R. Looking at outburst 2 and 4 separately, η_R stays roughly constant in the decaying phase of each outburst, whereas ṁ_in varies by at least a factor of 10 (see Fig. 3). In addition, different radio efficiencies are observed between each outburst during the decaying phases even at similar values of ṁ_in. Something else also seems to be at work.

5.2. Possible changes in the dynamical ejection properties

The existence of two functional dependencies F_R(ṁ_in,r_J) raises a profound question. Radiative processes in jets are local and are independent of disk parameters such as ṁ_in and r_J. However, the time evolution F_R(t) can be quite accurately reproduced with a function of (ṁ_in,r_J), which shows that global jet parameters do actually depend on them. These parameters, which constitute the jet dynamics, are for instance the magnetic field strength and geometry, the jet collimation degree, the existence of internal shocks or even jet instabilities. Our findings seem therefore to highlight two different jet dynamics.

5.2.1. A possible threshold in ṁ_in

The hard state data sets of the rising and decaying phases used in this analysis do not overlap in terms of accretion rate. Only two out of the 28 radio observations of the rising phases require a mass accretion rate comparable to those observed during the decaying phases. All the others have a higher mass accretion rate than the decaying phases. This is an observational bias due to the difficulties in catching the source as quickly as possible at the beginning of the outburst. However, the detected difference in the functional dependency of F_R could be due to some threshold in ṁ_in that could, in turn, translate into some difference in the way the radio emission scales with the disk parameters. Above the threshold the radio emission would follow Eq. (3), and below the threshold Eq. (4). Given the insufficient number of low accretion rate hard states in the rising phases, our analysis cannot test this possibility. Clearly, more observations are needed to assess this hypothesis.

We do not favor this interpretation, however. The main reason is that the rising and decaying hard states are temporally disconnected. The source stays in the soft state between these two hard-state phases for several months, and thus do not share the same history. The hard states in the rising phase come from a quiescent, already radio emitting state, while the hard states in the decaying phase come from soft radio silent states. This supports a link with the global jet structure (as proposed in Sect. 5.2.2) rather than a threshold in ṁ_in.

5.2.2. A change in the dominating ejection process

Since the commonly invoked radiative process is synchrotron, the first thing that comes to mind to explain this difference is the magnetic field strength. The only reasonable assumption to make is that this field is proportional to the field anchored at the JED, which writes (see Marcel et al. 2018b for more details)

$$\begin{array}{c}\hfill B_z(r)={(\mu {\mu }_0{P}_{\mathrm{tot}})}^{1/2}\simeq (\mu {\mu }_0{P}_{\ast }\frac{{\dot{m}}_{\mathrm{in}}{r}^{-5/2}}{m_s}{)}^{1/2},\end{array}$$(5)

where μ is the magnetization, m_s the accretion Mach number, P_* = m_in_*c², and $n_{\ast }=\frac{1}{{\sigma }_T{R}_g}$. Assuming constant JED parameters μ = 0.5 and m_s = 1.5 used in our model, we evaluate the magnetic field strength measured in r_isco at around 10⁸ G during the outbursts.

It should be noted that a JED exists within a small interval [μ_min, μ_max] of disk magnetization μ, with μ_min ∼ 0.1 and μ_max ∼ 0.8 (Ferreira 1997). The existence of such an interval has led Petrucci et al. (2008) to propose that the hysteresis observed in XrBs could be a consequence of a JED switch-off with μ = μ_min and switch-on at μ = μ_max. In the spectral analysis shown in the present paper, we suppose a constant μ since, as shown in Marcel et al. (2018b), the JED spectra are poorly affected by μ within the allowed parameter space. However, the possible difference in magnetization between the rising and decaying phase could also have a direct impact on the jet dynamical and radiative properties, explaining the change of the observed radio behavior. According to Eq. (5), a dichotomy of the magnetization μ at a given value of the mass accretion rate ṁ_in entails a dichotomy in the magnetic field strength. Thus, the rising phase, switching-off with μ = μ_min, would present a weaker magnetic field strength compared to the decaying phase, switching-on with μ = μ_max. This difference in the magnetic field strength could play a role in the difference of functional dependency of the radio emission. This could also explain the higher radio efficiencies observed in Fig. 11 during the decaying phases (e.g., Casella & Pe’er 2009).

Another possibility could be suggested by the most recent numerical simulations showing that the vertical magnetic field is carried in and accumulates around the black hole (building up a magnetic flux Φ_bh) until the surrounding disk magnetization reaches a maximum value near unity (see, e.g., Tchekhovskoy et al. 2011; Liska et al. 2020). In our view, the inner disk regions are nothing else than a JED driving a Blandford & Payne (BP, Blandford & Payne 1982) jet, although a Blandford & Znajek (BZ, Blandford & Znajek 1977) spine launched at its midst has attracted more attention in the literature⁸. Then another possible explanation for the existence of two functional dependencies for F_R(ṁ_in,r_J) could be that jets are two-component MHD outflows: a BZ spine, tapping into the rotational energy of the black hole, surrounded by a BP jet, tapping into the accretion energy reservoir of the disk. The jet dynamics and subsequent radio emission then depend on the relative importance of these two flows, which can be roughly measured by the ratio of the magnetic flux associated with each component, namely Φ_bh for the spine and Φ_JED for the outer BP jet. By construction, Φ_bh builds upon Φ_JED and reaches large values, such as ${\stackrel{\sim }{\mathrm{\Phi }}}_{\mathrm{bh}}={\mathrm{\Phi }}_{\mathrm{bh}}/{(<{\dot{M}}_{\mathrm{in}}>r_g^{2}c)}^{1/2}\sim 50$, only if a large magnetic flux is available initially in the disk (Tchekhovskoy et al. 2011; Liska et al. 2020). The functional dependency on r_J that we observe for the radio flux in our fits of the rising phase spectra could thus come from the dependency of ${\stackrel{\sim }{\mathrm{\Phi }}}_{\mathrm{bh}}$ on r_J.

A simple scenario can then be designed and is sketched in Fig. 12. During the rising hard-state phase, r_J is initially large and decreases in time (top right). The system comes from a quiescent state and the presence of a JED over a large radial extent allowed the disk to build up a maximum Φ_bh. The spine is very important and affects the overall jet dynamics, which translates into a radio flux described by Eq. (3). When the disk magnetization becomes too small, the JED transits to a SAD accretion mode (left, top and bottom). The magnetic field diffuses away, thereby decreasing Φ_bh, and no more jets are observed (neither BP nor BZ). As long as the system remains in the soft state, the field keeps on diffusing away until some equilibrium is eventually reached. At some point however the outburst declines, which translates into a decrease in the inner disk pressure, and thus an increase in the disk magnetization. In this decaying phase an inner JED becomes re-ignited inside-out, with its bipolar BP jets, but with a limited magnetic flux available (bottom right). By construction, Φ_bh remains small and the BZ spine has a limited impact on the overall jet dynamics. This would translate into a radio flux described by Eq. (4) until the JED is rebuilt over a large enough radial extent.

Fig. 12. Sketches of the inner regions of an accretion flow around a black hole during the different phases of the outburst. Top right: in the rising hard phase the JED is settled over a large region (RJ is large), leading to efficient magnetic flux accumulation on the black hole. The radio emission arises from a two-component outflow, made of an important BZ spine surrounded by a BP jet. Top left to bottom left: during the soft state there is no more JED, the magnetic field diffuses away, and the BP and BZ jets both disappear. Little or no weak radio emission is expected (jetline). Bottom right: in the decaying phase a JED reappears in the innermost region, and the magnetic field advection becomes efficient again. The magnetic flux on the black hole is still weak, and the BZ spine has little or no impact on the jet dynamics and subsequent radio emission.

There are many uncertainties in our different interpretations, since our JED-SAD modeling has its own simplifications. This last scenario is only an attempt to provide an explanation to our puzzling finding. Quite interestingly, it also provides a means of observationally testing it. It relies on the existence of a BZ spine in the case of GX 339-4, which is a black hole candidate. Around a neutron star the invoked scenario of magnetic flux accumulation into the central object clearly should not work. It would therefore be useful to investigate any changes in the radio properties during the rise and decay phases for neutron star binaries.

This scenario may look similar to the ones proposed by Begelman & Armitage (2014) or Kylafis & Belloni (2015) where the presence of a hot inner corona (an ADAF-like accretion flow in both cases) would help in accumulating or creating the required magnetic field that will eventually produce a jet. However, in these two approaches it is not clear why the process would differ between the rising and decaying phases, and how the functional dependency of the radio emission would depend on the ADAF properties. Clearly, more dedicated work should be done in this respect.

5.3. Effects of the JED-SAD parameters

In the results shown in this work we use the same values of the JED-SAD parameters b, m_s, and p as used in M19 and M20. A detailed study of the JED-SAD parameter space has already been performed (see Marcel et al. 2018b Sect. 4) and converges on these values in the case of GX 339-4. The parameter p has almost no spectral impact in the RXTE/PCA energy range (3–25 keV); this can be seen in Fig. 10 of Marcel et al. (2018b) (p was called ξ at the time). Thus, we do not expect any variation in the fits. Letting p remain free would result in an unconstrained parameter. The main impact of the parameters b and m_s is a variation of the maximum temperature in the JED, and results in a variation of the high-energy cutoff of the hard X-ray emission. However, the high-energy cutoff is not visible within the RXTE/PCA energy range used in our fitting procedure. We are then unable to constrain these parameters from the data and choose to set them to the same values used by M19 and M20. This also allows us to compare the evolution of the main JED-SAD parameters r_J and ṁ_in with the qualitative results obtained by M19 and M20.

However, we study the impact of using other values of m_s or b. The global trend of the parameters is similar, with only slightly different values of r_J and ṁ_in (see Figs. D.1 and D.2 in the Appendix). The biggest changes are observed with the different values of m_s. We thus concentrate on this parameter to see the effect of its value on the α − β contours (see Fig. D.3). The contour plots of the rising and decaying phase vary for different values of m_s; however, they are never consistent between each other, leaving our main conclusions unchanged: we need two different functional dependencies between the rising and decaying phases.

6. Conclusion

We presented in this paper the first direct fit of the X-ray data of an X-ray binary with our JED-SAD model. We constructed fits format tables that can be used in XSPEC. This includes a reflection table based on the XILLVER reflection model (Garcia et al. 2013). We applied our model to the X-ray observations of GX339-4, focusing on the “pure” hard-state phases of the four outbursts observed during the RXTE lifetime. We deduced from the fits the temporal evolution of the main parameters of the JED-SAD configurations: the inner accretion rate ṁ_in and the transition radius r_J between the inner JED and the outer SAD (see Fig. 3). This evolution is in relatively good agreement with the qualitative estimates done by M19 and M20; however, our spectral fit procedure puts much stronger constraints on r_J.

We were also able to put constraints on the functional dependency of the radio emission with r_J and ṁ_in for all outbursts of GX 339-4. Assuming a general radio flux expression $F_R=\stackrel{\sim }{f}{\dot{m}}_{\mathrm{in}}^{\beta }{r}_J^{\alpha }$, we were able to constrain the values of α and β for the different outbursts. These values appear consistent between the different rising phases or the different decaying phases. However, the rising or high mass accretion rate and decaying or low accretion rate phases solutions differ significantly. In the rising phase the radio emission varies as $\sim {\dot{m}}_{\mathrm{in}}^1{r}_J^{-0.6}$, while in the decaying phase the radio emission has a weaker dependency with r_J and varies as $\sim {\dot{m}}_{\mathrm{in}}^{0.9}{r}_J^{-0.2}$. While the exact values of the indexes depend slightly on the JED-SAD parameters, the two solutions obtained for the rising and the decaying phases are always mutually exclusive. A significant improvement in the fit of the radio fluxes is obtained by letting the scaling factor $\stackrel{\sim }{f}$ free to vary between the different phases. The observed variation (up to a factor 2) could correspond to a change in the radiative efficiency of the radio emitting process from outburst to outburst.

We suggest a few explanations for the difference in the functional dependency of the radio emission with r_J and ṁ_in between the rise and decay phases. A clear understanding is challenging given the scarce information we have on crucial jet parameters like the magnetic field strength and geometry, the jet collimation degree, the existence of internal shocks, or even jet instabilities. A possible scenario relies on a change of the relative importance of the Blandford & Znajeck versus Blandford & Payne processes in the radio emitting process due to the expected evolution of the magnetic field strength in the inner part of the accretion flow.

All results of this paper should be tested on other XrBs and on more recent outbursts of GX339-4. It would be also interesting to see how the case of neutron stars (where no BZ is expected) or objects belonging to the so-called outlier population (Coriat et al. 2011) would compare with our present results. Especially since the outliers present a steep power index in the radio–X-ray plane at high luminosity, similar to that observed in neutron star binaries. This will be the focus of a forthcoming paper.

Acknowledgments

The authors acknowledge funding support from CNES and the French PNHE.

References

Appendix A: Outburst 4 radio interpolation

With 80 almost daily X-ray observations and 24 radio flux measurements, the radio–X-ray survey of the hard states of the 2010 outburst of GX339-4 is the best in the available archive of RXTE. Considering the steady jet expected during the hard states of an outburst and the evenly spread radio survey, we opted to linearly interpolate the radio light curve to obtain a radio flux for each X-ray observation. In Fig. A.1 we plot the radio light curve.

Fig. A.1. Radio light curve of outburst 4 during its hard states: between MJD 55208 and 55293, and between MJD 55608 and 55646. The green squares are the observed radio data. The blue dashed line is the linearly interpolated function. The black diamonds are the interpolated radio fluxes at the date of the X-ray observations.

Appendix B: Fitting the decaying phase of outburst 4

In Fig. B.1 we plot the evolutions of r_J and ṁ_in during outburst 4. The blue dashed line represents the final results presented in Fig. 3 used for the study. The red dashed line represents the initial results obtained by our fitting procedure. In this appendix we explain why we rejected the results of this fitting procedure, and how we obtained the final results for the decaying phase of outburst 4 using a maximum likelihood method.

Fig. B.1. Results of the fitting procedure for outburst 4. On the left side is shown the transition radius rJ between the JED and the SAD, on the right side the mass accretion rate ṁin. Each side is divided horizontally between the rising and decaying phase of each outburst. The green solid line represents the results from M19 and M20, and the green region where their minimization function varies by less than 10% with respect to its minimum. The blue dashed line shows the results presented in Fig. 3 and the black vertical bar the associated 90% confidence range. The red dashed line represents the initial fitting results for the decaying phase of outburst 4.

The main disagreement of the fitting procedure with M19 was observed in the decaying phase of outburst 4. The fits yielded quite small values of r_J (∼5R_G) and a decrease in the transition radius with time (Fig. B.1, red), which is not expected in the JED-SAD paradigm and in contradiction with the evolution observed by M19, although their constraints are quite large. We report in Fig B.2 the evolution of the χ² statistic along the parameter space of r_J for a few observations of this decaying phase. These observations show a non-trivial χ² space with multiple local minimums. Most of them presented either a better fit or a statistically equivalent fitting solution at higher r_J. For some reason, Xspec did not find these solutions even during the error calculations. In Fig. B.3 we represent Markov chain Monte Carlo (MCMC) tests for one of these observations (MJD 55630), starting from the higher (red) and lower (black) r_J solutions. In both cases, the procedure converged to a reasonable solution within a few steps (< 5). When starting with the lower r_J solution, the procedure first explores only the lower r_J values. It is only after ∼ 2500 steps that it explores the higher r_J solutions, and in the last 1000 steps, both procedures only explore the higher r_J solution. These higher r_J solutions were further motivated by the fact that high values of r_J naturally appear in outbursts 2 and 3 (see Fig. 3). For these different reasons, we set r_J to this solution for each observation of the decaying phase of outburst 4.

Fig. B.2. Evolution of the $\mathrm{\Delta }{\chi }^2=({\chi }^2(r_J)-{\chi }_{\text{min}}^2)$ with rJ for a few observations from the decaying phase of outburst 4. From top to bottom in chronological order: MJD 55613, 55617, 55620, 55630, 55639. The 90% confidence threshold Δχ2 = 2.71 is shown as a dashed green horizontal line. We indicate with a star the solutions with higher values of rJ (plotted in red in Figs. 3 and B.1).

Fig. B.3. MCMC procedure for observation 55630 (see the χ2 parameter space in Fig. B.2 in green) starting from two different initial priors. In black starting with a lower rJ prior, in red starting with a higher rJ prior. Bottom: histograms of rJ during the first 1000 steps (left) and the last 1000 steps (right). Both procedures explore only the higher rJ solution.

A similar process could have been applied to the decaying phase of the other outbursts; however, the lack of radio flux where similar problems were encountered made this issue non-critical for the rest of the study.

Appendix C: Effects of the systematic error

In this section we show the effects of adding systematic errors to our radio fluxes. The reasoning behind this addition is that we accept a 10% error on the reproduction of the radio fluxes. As the χ² statistic gives more weight to the observations with small errors, without systematic errors added the fit is mostly driven by a few fluxes and does not match our quality criterion. Our goal is to reproduce all the radio fluxes within a 10-20% error margin. Adding a systematic error to all fluxes increases the relative weight of the observations with higher uncertainty compared to the ones with low flux error.

However this has a major impact on the size of the confidence contour β-α that we obtain. These contours depend directly on the χ² plane of our parameter space. In Fig. C.1 we plot the confidence contour for a 90% confidence level (Δχ² = 4.61) for all the rising phase observations (in blue) and for all the decaying phase observations (in red), and for different values of systematic errors added to the radio fluxes (0%, 5%, and 10% depending on the thickness of the line). As expected, the larger the systematic error added, the larger the contour plot, but even with a 10% systematic error the rising phase solution and decaying phase solution are inconsistent.

Fig. C.1. 90% Confidence contours β–α for different values of systematic error added to the radio fluxes (0%, 5%, and 10%). The thicker the line, the larger the systematic error added. Shown are the contours for all the rising phase observations (in blue) and the contours for all the decaying phase observations (in red).

Appendix D: Impact of the other JED-SAD parameters

In this study we set the values of the parameters m_s, p, and b to that used by M19. Marcel et al. (2018a) already performed a study of the parameter space and converged on these values. However, as the main spectral impact of m_s and b is a variation of the high-energy cutoff, which is not visible inside the RXTE/PCA energy range, we are faced with unconstrained parameters. In this paragraph we present the same study reproduced using other set values of m_s and b. In Figs. D.1 and D.2, we plot the fitting results obtained for different values of m_s and b. The global trend of the parameters is similar; however, the values of r_J and ṁ_in are slightly different. The biggest difference is observed with the different values of m_s. In Fig. D.3 we plot the corresponding contour regions α-β for the different values of m_s. For the decaying phase observations, we used the procedure described in Appendix B to obtain the transition radius r_J for eight quasi-simultaneous observations used for the contour plots. The contours obtained for different values of m_s are not always consistent with each other. However, in each case the rising phase and decaying phase contours are inconsistent. This implies that we always have two different behaviors for the radio between the rising and decaying phase regardless of the values of m_s and b that we use in the fits.

Fig. D.1. Results of the fitting procedure for outburst 4. On the left side is shown the transition radius rJ between the JED and the SAD, on the right side the mass accretion rate ṁin. Each side is divided horizontally between the rising and decaying phase of each outburst. The green solid line represents the results from M19, and the green region where the minimization function varies by less than 10% with respect to its minimum. The blue, orange, and purple dashed lines show the results of the automatic fitting procedure using different values of ms (1.0, 1.5, and 2.0, respectively).

Fig. D.2. Results of the fitting procedure for outburst 4. On the left side is shown the transition radius rJ between the JED and the SAD, on the right side the mass accretion rate ṁin. Each side is divided horizontally between the rising and decaying phase of each outburst. The green solid line represents the results from M19 and M20, and the green region where the minimization function varies by less than 10% with respect to its minimum. The blue, orange, and purple dashed lines show the results of the automatic fitting procedure using different values of b (0.1, 0.3, and 0.5, respectively).

Fig. D.3. Confidence contours β-α for different values of ms used in the fitting procedure of outburst 4: 1.0 (solid line), 1.5 (dashed line), and 2.0 (dotted line). In all cases the rising phase contour and decaying phase contour are inconsistent.

All Tables

All Figures

Fig. 1. GX 339-4 X-ray light curve in the 3–200 keV energy band of the 2000–2010 decade obtained with the Clavel et al. (2016) fits. The violet filled region shows the power law unabsorbed flux, while the cyan region represents the disk unabsorbed flux. The selected spectra for this study are highlighted in gray: the rising and decaying hard states of the four outbursts (1, 2, 3, and 4). At the top the red lines represent the date when steady radio fluxes were observed at 9 GHz (from Corbel et al. 2013).

In the text

Fig. 2. Best fits of some observations of outburst 4. Top left: hardness intensity diagram of outburst 4. The blue diamond shows the hard state used for this outburst. The five red points are the five observations plotted in the different figures from a to e. a–e: best fit spectra and the data/model ratio for the five observations indicated in red in the HID. The gray region shows the PCA energy range used for the fit. The data are in black and the best fit model in gray; the JED spectrum is in red, the SAD spectrum in green, and the reflection component in blue. The best fit parameters for each observation are listed in Table 2.

In the text

Fig. 3. Results of the fitting procedure. Left: transition radius rJ (from risco to 300) between the JED and the SAD. Right: mass accretion rate ṁin. Each side is divided vertically between the four outbursts and horizontally between the rising and decaying phase of each outburst. The green solid line represent the results from M19 and M20, and the green region indicates where the minimization function varies by less than 10% with respect to its minimum. The blue dashed line shows the results of the fitting procedure and the black vertical bar the associated 90% confidence range. The decaying phase of outburst 4 is the subject of Appendix B.

In the text

Fig. 4. Ratio of the observed radio fluxes to the results of Eq. (1) (using $\stackrel{\sim }{f}=1.5\times {10}^{10}$) as a function of rJ (left) and ṁin (right) for outburst 4. In the case of X-ray observations without simultaneous radio measurement, the radio flux was interpolated from the radio light curve (see Fig. A.1). The dashed line shows the best fit power law: $\frac{F_{\mathrm{obs}}}{F_R}\propto$$r_J^{-1.25}$ (left) and $\frac{F_{\mathrm{obs}}}{F_R}\propto$ ${\dot{m}}_{\mathrm{in}}^{1.56}$ (right).

In the text

Fig. 5. Contour plots β–α for the rising (blue) and decaying (red) phases of the outburst of 2004 (outburst 2, dashed line), 2007 (outburst 3, dot-dashed line), and 2010 (outburst 4, thin solid line). The contours in thick solid lines represent the dependency when fitting all rising (blue) or decaying (red) phase radio fluxes simultaneously. Confidence contour levels correspond to 68%, 90%, and 99% (Δχ2 of 2.3, 4.61, and 9.2, respectively). The contours are obtained when fitting only the quasi-simultaneous radio–X-ray observations (not the interpolated radio observations).

In the text

Fig. 6. Ratio of the data to the model for the best fit of 5 of the 16 multiwavelength observations (radio–X-ray) of the rising phase of the 2010 outburst (MJD 55217 in blue, 55259 in red, 55271 in black, 55288 in green, and 55292 in violet). Only five ratios are shown for purposes of visualization, but the best fit was obtained by using all the simultaneous or quasi-simultaneous radio–X-ray observation, fixing the JED-SAD parameters to the best fitting values obtained by fitting the X-ray spectra first, and then fitting the radio points with Eq. (2).

In the text

Fig. 7. Ratio of all radio fluxes of the rising phase of outburst 4 to the results of Eq. (2) as a function of rJ (left) and ṁin (right). The filled blue points represent the 16 quasi-simultaneous radio–X-ray observations, while the empty points represent the interpolated radio fluxes (see Appendix A). The radio observations are well reproduced using the values ${\stackrel{\sim }{f}}^{\ast }=7.1\times {10}^{-8}$, α = −0.66, and β = 1.00.

In the text

Fig. 8. Ratio of the data to the model for the best fit of five of the eight multiwavelength observations of the decaying phase of outburst 4 (MJD 55613 in blue, 55617 in red, 55620 in black, 55630 in green, and 55639 in violet). Top panel: all fits were done simultaneously, fixing the parameters of Eq. (2) to those found in the rising phase: ${\stackrel{\sim }{f}}^{\ast }=7.1\times {10}^{-8}$, α = −0.66, and β = 1.00. Bottom panel: fixing α and β to the values found for the rising phase, ${\stackrel{\sim }{f}}^{\ast }$ is free to vary and converges to the value 2.0 × 10−7.

In the text

Fig. 9. Ratios of the observed radio fluxes to the modeled radio fluxes for all the quasi-simultaneous radio–X-ray observations of the four outbursts. Each outburst is represented using a different symbol (see legend). In blue are shown the rising phases and in red the decaying phases. The modeled radio fluxes were obtained using Eq. (2). The parameters are (α = −0.67, β = 0.94) for the rising phases and (α = −0.15, β = 0.9) for the decaying phases. Top panel: ${\stackrel{\sim }{f}}^{\ast }=4.1\times {10}^{-8}$ used for all outbursts. Bottom panel: different ${\stackrel{\sim }{f}}^{\ast }$ used for each phase of the outbursts. All values of ${\stackrel{\sim }{f}}^{\ast }$ used are reported in Table 3. The horizontal dashed lines represent a 20% error margin (ratio of 0.8 and 1.2, respectively).

In the text

Fig. 10. Evolution of the hard X-ray power law index Γ during the four outbursts of GX339-4. Left: hard X-ray power law index Γ as a function of the 3–9 keV X-ray luminosity of the pure hard state observations during the four GX339-4 outbursts (data from Clavel et al. 2016). In blue the rising phase and in red the decaying phase. Right: histograms of Γ. These distributions are subject to a certain number of observational biases: inclusion of error bars and the number of observations per phase.

In the text

Fig. 11. Radio emission efficiency (LR = L9 GHz/PJets) vs. X-ray emission efficiency (LX = L3 − 9 keV/PJED) during the outbursts of GX339-4. The blue points are the rising phases. The red points are the decaying phases. The symbols distinguish the different outbursts: diamonds for 2002, squares for 2004, dots for 2007, and triangles for 2010 (filled for the quasi-simultaneous observations and empty for the interpolated radio fluxes). We highlighted the five observations (labeled a to e) presented in Fig. 2 to provide the chronological evolution of an outburst.

In the text

Fig. 12. Sketches of the inner regions of an accretion flow around a black hole during the different phases of the outburst. Top right: in the rising hard phase the JED is settled over a large region (RJ is large), leading to efficient magnetic flux accumulation on the black hole. The radio emission arises from a two-component outflow, made of an important BZ spine surrounded by a BP jet. Top left to bottom left: during the soft state there is no more JED, the magnetic field diffuses away, and the BP and BZ jets both disappear. Little or no weak radio emission is expected (jetline). Bottom right: in the decaying phase a JED reappears in the innermost region, and the magnetic field advection becomes efficient again. The magnetic flux on the black hole is still weak, and the BZ spine has little or no impact on the jet dynamics and subsequent radio emission.

In the text

	Fig. A.1. Radio light curve of outburst 4 during its hard states: between MJD 55208 and 55293, and between MJD 55608 and 55646. The green squares are the observed radio data. The blue dashed line is the linearly interpolated function. The black diamonds are the interpolated radio fluxes at the date of the X-ray observations.
In the text

Fig. B.1. Results of the fitting procedure for outburst 4. On the left side is shown the transition radius rJ between the JED and the SAD, on the right side the mass accretion rate ṁin. Each side is divided horizontally between the rising and decaying phase of each outburst. The green solid line represents the results from M19 and M20, and the green region where their minimization function varies by less than 10% with respect to its minimum. The blue dashed line shows the results presented in Fig. 3 and the black vertical bar the associated 90% confidence range. The red dashed line represents the initial fitting results for the decaying phase of outburst 4.

In the text

Fig. B.2. Evolution of the $\mathrm{\Delta }{\chi }^2=({\chi }^2(r_J)-{\chi }_{\text{min}}^2)$ with rJ for a few observations from the decaying phase of outburst 4. From top to bottom in chronological order: MJD 55613, 55617, 55620, 55630, 55639. The 90% confidence threshold Δχ2 = 2.71 is shown as a dashed green horizontal line. We indicate with a star the solutions with higher values of rJ (plotted in red in Figs. 3 and B.1).

In the text

	Fig. B.3. MCMC procedure for observation 55630 (see the χ2 parameter space in Fig. B.2 in green) starting from two different initial priors. In black starting with a lower rJ prior, in red starting with a higher rJ prior. Bottom: histograms of rJ during the first 1000 steps (left) and the last 1000 steps (right). Both procedures explore only the higher rJ solution.
In the text

	Fig. C.1. 90% Confidence contours β–α for different values of systematic error added to the radio fluxes (0%, 5%, and 10%). The thicker the line, the larger the systematic error added. Shown are the contours for all the rising phase observations (in blue) and the contours for all the decaying phase observations (in red).
In the text

Fig. D.1. Results of the fitting procedure for outburst 4. On the left side is shown the transition radius rJ between the JED and the SAD, on the right side the mass accretion rate ṁin. Each side is divided horizontally between the rising and decaying phase of each outburst. The green solid line represents the results from M19, and the green region where the minimization function varies by less than 10% with respect to its minimum. The blue, orange, and purple dashed lines show the results of the automatic fitting procedure using different values of ms (1.0, 1.5, and 2.0, respectively).

In the text

Fig. D.2. Results of the fitting procedure for outburst 4. On the left side is shown the transition radius rJ between the JED and the SAD, on the right side the mass accretion rate ṁin. Each side is divided horizontally between the rising and decaying phase of each outburst. The green solid line represents the results from M19 and M20, and the green region where the minimization function varies by less than 10% with respect to its minimum. The blue, orange, and purple dashed lines show the results of the automatic fitting procedure using different values of b (0.1, 0.3, and 0.5, respectively).

In the text

	Fig. D.3. Confidence contours β-α for different values of ms used in the fitting procedure of outburst 4: 1.0 (solid line), 1.5 (dashed line), and 2.0 (dotted line). In all cases the rising phase contour and decaying phase contour are inconsistent.
In the text