https://orcid.org/0000-0001-5674-4664
ABSTRACT
Multiwavelength spectral energy distributions of low-mass X-ray binaries (LMXBs) in the hard state are determined by the emission from a jet, for frequencies up to mid-infrared, and emission from the accretion flow in the optical to X-ray range. In the last years, the flat radio-to-mid-IR spectra of black hole (BH) X-ray binaries was described using the internal shocks model, which assumes that the fluctuations in the velocity of the ejecta along the jet are driven by the fluctuations in the accretion flow, described by the X-ray power density spectrum (PDS). In this work, we attempt to apply this model for the first time to a neutron star (NS) LMXB, i.e. 4U 0614 + 091. We used the multiwavelength data set obtained in 2006, comprising data from radio to X-ray, and applied a model that includes an irradiated disc model for the accretion flow and an updated version of the internal shocks code for the ejection. The new version of the code allows to change the geometry of the jet for the case of non-conical jets. Only two alternative scenarios provide a satisfactory description of the data: using the X-ray PDS but in a non-conical geometry for the jet, or either using a conical geometry but with a ‘flicker-noise’ PDS. Both scenarios would imply some differences with the results obtained with similar models on BH X-ray binaries, shedding light on the possibility that jets in NS and BH binaries might somehow have a different geometry or a different coupling with the accretion flow.
1 INTRODUCTION
The ejection of collimated outflows of matter in the form of jets is an ubiquitous phenomenon in astrophysics, and it has been associated with a wide range of celestial objects, from active galactic nuclei (AGNs) to young stellar objects and stellar-mass X-ray binaries (XRBs). Jets have been studied in detail for decades but many points about them, in particular concerning their launching mechanisms and their coupling with the accretion flow, are still debated (see for a review, Belloni 2010).
In XRBs, jets are usually not observed as extended structures, but their presence is witnessed by a radio-to-mid-IR flat spectrum (in this case, jets are referenced to as ‘compact’, see e.g. Corbel et al. 2000; Fender 2001; Corbel & Fender 2002). The spectrum of a compact jet is characterized by a jet break, corresponding to the base of the jet, where emission goes from optically thick to optically thin (see e.g. Gandhi et al. 2011; Russell et al. 2013; Koljonen et al. 2015) and at lower frequencies by a continuum that is the result of the superposition of the self-absorbed synchrotron spectra emitted from the different regions of the jet (Condon & Dressel 1973; de Bruyn 1976; Marscher 1977; Konigl 1981; Ghisellini, Maraschi & Treves 1985). Moving along the jet, the magnetic field decays and the particles lose energy as the jet expands: the peaks of the local single synchrotron spectra are expected to decrease in both intensity and frequency, leading to an observed inverted radio spectrum (Marscher 1980). However, compact jet spectra are usually almost flat, and this has been explained in the past invoking the existence of some continuous energy replenishments mechanisms that would compensate the adiabatic losses due to the expansion of the jet (Blandford & Königl 1979). Alternatively, sufficiently collimated jets would emit flat spectra even without taking into account dissipation mechanisms (Kaiser 2006), but this would require such a fine-tuning of the jet geometry that it is unlikely to explain the majority of the observed cases.
The nature of the dissipation mechanism is still unclear, but several models have been proposed in the past and involve, e.g. magnetic reconnection (Sobacchi & Lyubarsky 2020), acceleration in relativistic shear flows (Rieger 2019), or internal shocks. The latter models take into account the conversion from kinetic into internal energy that arises when two shells in the jet, ejected at different velocity, catch up and collide. Internal shocks models have been applied in the past to γ-ray bursts (Rees & Meszaros 1994; Daigne & Mochkovitch 1998), AGNs (Rees 1978; Spada et al. 2001; Boettcher 2010), and BH XRBs (Kaiser, Sunyaev & Spruit 2000; Jamil, Fender & Kaiser 2010; Malzac 2013). In the last decade, Malzac (2013), Malzac (2014) showed indeed that internal shocks can also explain the flat Spectral Energy Distributions (SEDs) of XRBs if we assume that the fluctuations in the velocity of the ejecta are matched by the fluctuations in the accretion flow. Under this hypothesis, which underlies a profound accretion/ejection connection for jets in XRBs, one can use the observed X-ray power spectral density as a proxy for the fluctuations of the accretion flow. The model, dubbed internal shocks model (ISHEM), has been successfully applied to a number of XRBs hosting black holes (BHs) as the primary star in the past (Drappeau et al. 2015; Baglio et al. 2018; Malzac et al. 2018; Péault et al. 2019; Bassi et al. 2020), but never to a neutron star (NS) XRB.
Jets are observed in XRBs hosting NSs as well. They are very common in low-mass X-ray binaries (LMXBs) systems, which host NSs with usually weak (B ≲ 109 G) magnetic fields, but recently the phenomenon has been also associated with one highly magnetized NS in a high-mass XRB (van den Eijnden et al. 2018). Unfortunately, the wealth of studies of jet phenomenology in BH XRBs is unmatched when it comes to NS XRBs mainly because they tend to have weaker radio emission (from hundreds to tens of μJy), sometimes below the observational capabilities of the most sensitive interferometers on the Earth. In addition, NS LMXBs show faster state transition time-scales (Muñoz-Darias et al. 2014; Marino et al. 2019), which make it harder to schedule coordinated space and ground observations. As mentioned before, jets in NSs are less radio-loud than in BHs, i.e. by a factor of ∼30 at similar X-ray flux levels (Fender & Kuulkers 2001; Migliari et al. 2003; Tudor et al. 2017; Gallo, Degenaar & van den Eijnden 2018). Furthermore, in BH binaries jet emission is always suppressed when the source is in the soft state, while jets are never entirely quenched in NS XRBs (see e.g. Migliari et al. 2004), with just a few exceptions (Miller-Jones et al. 2010; Gusinskaia et al. 2017). Accretion–ejection coupling in binaries has been traditionally studied with radio–X-ray luminosity diagrams LR:LX in both BH and NS XRBs and yet again a clearer picture seems to emerge for the former class of systems. In such diagrams,1 BH XRBs tend to populate two branches |$L_{\rm R}\propto L^{\beta }_{\rm X}$|, i.e. with β ∼ 0.6 for the ‘radio-loud’ systems and β ∼ 1.4 for the ‘radio-quiet’ systems (see e.g. Corbel et al. 2013), and this behaviour has been proposed to originate in different physical properties in the accretion flow (Coriat et al. 2011) or in the jet (Espinasse & Fender 2018) over the two branches. A similar dichotomy can be found in NS XRBs, but the distribution appears more scattered and harder to interpret (Tetarenko et al. 2016). Alternatively, Gallo et al. (2018) proposed a single-track population for both BH and NS XRBs but with different values of β, i.e. the former with β ∼0.7 and the latter with β ∼0.4. It is also noteworthy that even the jet launching mechanism could not be the same in BH and NS LMXBs since of the two traditionally proposed jet launching mechanisms, i.e. Blandford & Znajek (1977) and Blandford & Payne (1982), only the second could be at work in both classes (but see also Migliari, Miller-Jones & Russell 2011, for a discussion on the possibility of having spin-powered jets also on NS LMXBs.) In summary, the emerging picture seems to imply that the nature of the compact object (and, as in the case of an NS, the presence of a magnetic field) could play some role in determining the characteristics and the origin of the jet.
1.1 4U 0614 + 091
Discovered by the Uhuru survey in the 70s (Forman et al. 1978), 4U 0614 + 091 was later identified as an LMXB hosting an NS by the detection of type-I X-ray bursts (Swank et al. 1978). From the study of the bursts, a measure of the distance was obtained, i.e. around 3.2 kpc with a 15 per cent uncertainty (Kuulkers et al. 2010). Due to its short orbital period of around 50 min (Shahbaz et al. 2008; Baglio et al. 2014), the system has been classified as an ultracompact XRB or UCXB, implying a likely degenerate-helium dwarf or white dwarf nature for the companion star (see e.g. Kuulkers et al. 2010, for an extensive discussion on the companion star nature). The source is classified as a persistent atoll source2 and is expected to spend in the hard (island) spectral state almost 90 per cent of its time (van Straaten et al. 2000), with a constantly high X-ray spectral variability level (higher than 5 per cent) and only episodical transitions to softer states (Muñoz-Darias et al. 2014). In the past, a few authors gave evidence for the presence of a so-called hard tail in the spectrum, reaching energies beyond 100 keV, which was modelled with non-thermal models (Piraino et al. 1999; Migliari et al. 2010) or thermal Comptonization from a very hot corona (Ford et al. 1997; Piraino et al. 1999; Fiocchi et al. 2008). A reflection component has been commonly used to describe the X-ray spectral emission too, although the Fe–K line was usually found to be absent or weak. The apparent absence (or weakness) of this feature may be related to an underabundance of Fe (Madej et al. 2014; Ludlam et al. 2019) in the secondary star with respect to solar abundances, and it is compatible with the hypothesis of an out-of-main-sequence companion.
The flat radio to mid-IR spectrum, reported by Migliari et al. (2010) in the first complete multiwavelengths spectral study of the source, witnesses the presence of a compact jet in the system, as confirmed by the polarimetric study by Baglio et al. (2014).
In this paper, we report on the application of the internal shocks model to the spectral energy density of 4U 0614 + 091, the first attempt ever to describe the entire SED of an NS LMXB with a model including both the jet and the accretion flow emission. Furthermore, ISHEM has been used so far only with systems hosting a BH and never for an NS LMXB, as in this work.
2 DATA
The ISHEM model depends on three essential ingredients: (1) a multiwavelength SED, (2) a PDS that is used as tracer of the accretion flow variability, and (3) a synthetic spectrum, simulated on the basis of the PDS to be compared with the real SED. In the previous applications of ISHEM (see Drappeau et al. 2015; Malzac et al. 2018; Péault et al. 2019; Bassi et al. 2020), it was proven that the X-ray PDS quasi-simultaneous to the SED can be satisfactorily used as ingredient (2). The methods to obtain these three ingredients are described in the following sections. In particular, in this section we describe the data set and the timing analysis while the model used, the main parameters adopted and the spectral fitting procedure will be described in Section 3.
This work takes advantage of the multiwavelength observational campaign performed on the source within 5 d, between 2006 October 30 and November 4, from radio to X-ray. For the radio-to-IR realm, we used radio observations collected by the Karl G. Jansky Very Large Array (VLA), mid-IR/IR observations taken by the Infrared Array Camera (IRAC) onboard Spitzer Space Telescope. We used near-IR/optical data by the ground-based small and moderate aperture research telescope system (SMARTS), while the optical/UV observations were taken with UVOT onboard the Neil Gehrels Swift Observatory (Swift in the following). Finally, X-ray data were obtained with the proportional counter array (PCA) and the high-energy X-ray timing experiment (HEXTE) onboard the Rossi X-ray timing explorer (RXTE) and XRT onboard Swift. This vast data set was already used for a comprehensive spectral analysis by Migliari et al. (2010). We refer to this paper for the details concerning the data reduction and analysis performed on the observations, except for the Swift/UVOT data, which were re-extracted in this work.
2.1 SMARTS–UVOT data treatment
As pointed out by Migliari et al. (2010), the data in the optical-ultraviolet (covered by SMARTS–UVOT) region of the electromagnetic spectrum show indeed an unexpected shape, which cannot be ascribed to irradiation of the outer disc or to the blackbody emission from the (likely very faint, Nelemans et al. 2004; Shahbaz et al. 2008) companion star. In order to check if a different treatment of the data might improve their results, we re-analysed UVOT data (ObsID 00030812001) with heasoft v. 6.26 following the standard procedure3 and we used the calibration files updated to the latest available version (CALDB 2017-09-22). This observation was carried out using all the six UVOT filters. Using the task UVOTDETECT, we clearly detected 4U 0614+091 in each image. We defined a source region with a 5 arcsec radius and several different background regions for each filter around the source. Finally, the photometry of 4U 0614 + 091 has been performed with the task UVOTSOURCE. The de-reddening applied by Migliari et al. (2010) on the SMARTS data was removed using equation (1) and the values reported in table 3 (for the bands V, I, and J) of Cardelli, Clayton & Mathis (1989) and considering AV = 2, as reported in Migliari et al. (2010) in order to get A(λ), i.e. the absorption at wavelength λ. Considering then Fλ = Fλ, 0 × e−A(λ)/1.086, with Fλ and Fλ,0 the absorbed and unabsorbed fluxes at wavelength λ, respectively, we obtained the required values for the uncorrected SMARTS data. The new UVOT data and the uncorrected SMARTS data have been then de-reddened by us via xspec, using a proper model (see SubSection 4.2).
With the new treatment of the optical-UV data, the odd IR-UV spectral shape reported by Migliari et al. (2010) has now disappeared. In Fig. 1, we compare the ‘old’ and the new data sets de-reddened, in order to investigate the nature of the previously reported tricky result. We therefore retain this discrepancy arises from the different extraction methods, in particular on the choice and sizes of source and background regions used for the photometric measures.
Comparison between the unabsorbed UVOT data in Migliari et al. (2010; the yellow points) and the unabsorbed UVOT data used in this paper (the green points). In particular, for the new UVOT data set we corrected for an absorption coefficient of AV ∼1.5, found using the model redden on Xspec (see subsections 4.2 and 4.3 for further details).
2.2 X-ray timing analysis
For the timing analysis of the observation 92411-01-06-07 (2006 October 30), we used RXTE PCA data in the event mode configuration with a time resolution of ∼125 μs, allowing to obtain PDS up to a Nyquist frequency of 4096 Hz. We averaged multiple PDS data calculated over 128 s subintervals covering a total data set of 2048 s, using fast fourier transform techniques. No deadtime corrections nor background subtraction were performed before creating the PDS. We subtracted the Poisson noise power, derived from the PDS in the frequency range between 1536 and 2048 Hz, following Zhang et al. (1995). Fig. 2 shows the traditional ν, Pν representation where we applied the Leahy normalization (Leahy et al. 1983) before converting the PDS to squared fractional rms. The resulting PDS was fitted with a model consisting of the sum of two Lorentzians, one broad Lorentzian to fit the low-frequency noise and one narrow to fit the Quasi-Periodic Oscillation (QPO) in the range ∼500–700 Hz. The best-fitting parameters found were used as input for ISHEM and are listed in Table 1.
PDS in the normalized power (Pν) times frequency (ν) representation, with the best-fitting multi-Lorentzian model.
Fit results of the PDS described with a sum of two Lorentzians, each of them given by P(ν) = r2/π[Δ2 + (ν−ν0)2], with r the integrated rms over the full range of frequencies −∞ to +∞, Δ the full width at half-maximum of the Lorentzian and ν0 its central frequency. Values in round parentheses were kept frozen during the fit. Quoted errors reflect 68 per cent confidence levels.
| Lorentzian | ν0 | Δ | r |
|---|---|---|---|
| Component | |||
| 1 | (0) | 62.0 ± 5.0 | |$0.106^{+0.003}_{-0.006}$| |
| 2 | |$650^{+30}_{-24}$| | <250 | |$0.030^{+0.010}_{-0.008}$| |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(126) | |||
| Lorentzian | ν0 | Δ | r |
|---|---|---|---|
| Component | |||
| 1 | (0) | 62.0 ± 5.0 | |$0.106^{+0.003}_{-0.006}$| |
| 2 | |$650^{+30}_{-24}$| | <250 | |$0.030^{+0.010}_{-0.008}$| |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(126) | |||
Fit results of the PDS described with a sum of two Lorentzians, each of them given by P(ν) = r2/π[Δ2 + (ν−ν0)2], with r the integrated rms over the full range of frequencies −∞ to +∞, Δ the full width at half-maximum of the Lorentzian and ν0 its central frequency. Values in round parentheses were kept frozen during the fit. Quoted errors reflect 68 per cent confidence levels.
| Lorentzian | ν0 | Δ | r |
|---|---|---|---|
| Component | |||
| 1 | (0) | 62.0 ± 5.0 | |$0.106^{+0.003}_{-0.006}$| |
| 2 | |$650^{+30}_{-24}$| | <250 | |$0.030^{+0.010}_{-0.008}$| |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(126) | |||
| Lorentzian | ν0 | Δ | r |
|---|---|---|---|
| Component | |||
| 1 | (0) | 62.0 ± 5.0 | |$0.106^{+0.003}_{-0.006}$| |
| 2 | |$650^{+30}_{-24}$| | <250 | |$0.030^{+0.010}_{-0.008}$| |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(126) | |||
3 THE JET MODEL
3.1 The spectral shape
The ishem model aims to describe the spectral energy distribution of jets based on how the energy is dissipated along the jet.
As mentioned before, the emission of jets is ascribed to the superposition of the self-absorbed synchrotron spectra emitted locally from the different regions of the jet, peaking at decreasing energies as we move away from the base of the jet. However, the almost ‘flat’ final spectrum requires that in roughly each region the energy lost in the expansion of the jet is somehow gained back by the particles. If we imagine the jet as the result of the periodic ejection of discretized shells of matter and if these shells are ejected with variable velocity or Lorentz factor Γ, then we expect a fraction of the energy of these shells to be released for each collision between ejecta travelling with different Γ. Internal shocks turn out then to be a viable mechanism to replenish the energy lost and thereby to flatten the spectra. However, shocks have to occur homogeneously all over the jet axis. The required dissipation pattern is ultimately determined by two factors: how fast the flow expand and the particles lose energy, which is in turn determined by the geometry of the jet (1) and how the velocity of the ejected shells fluctuates over time (2). Indeed, a fast variability, i.e. over short time-scales, mainly produces collisions close to the base of the jet, while on the other hand shells subject to slow variability would tend to produce shocks at higher distances. In ishem, these two factors are regulated by a geometry parameter ζ and the input Power Density Spectrum (PDS) P(ν; with ν frequency) of the Lorentz factor Γ fluctuations. A third important ingredient is p, which determines the slope α = (p−1)/2 of the spectrum at high energies, in the optically thin part of the jet spectrum. The final SED shape is determined by a combination of these three ingredients. In the following, we will give more details on the effects of the two aforementioned factors.
3.1.1 PDS
In order to have a homogeneous dissipation pattern, variability in the ejecta must be almost homogeneously distributed over a large range of time-scales. For example, in a conical geometry, the exact compensation of the energy losses requires that the PDS of the jet Lorentz factor fluctuations corresponds to a ‘flicker noise’, i.e. P(f) ∝ 1/f over a broad range of Fourier frequencies (Malzac 2013). Such flicker noises occur if variability in the jet bulk Lorentz factor Γ is dominated by the so-called ‘flicker noise’ (see e.g. Press 1978), which occurs in processes of different nature, e.g. biological, economic, and physical, but also in astrophysics, especially in the X-ray variability of XRBs (Gilfanov 2010). Interestingly, the Fourier PDS extracted from the X-ray light curves of BH XRBs in the hard spectral state is, at low frequencies, very similar to a flicker-noise PDS (see Malzac 2014, and references therein). Furthermore, in a disc–jet coupling scenario, shells are ejected from the disc and it is reasonable to expect that the fluctuations in the accretion flow might be transmitted to the ejecta. According to these clues, the X-ray PDS can serve as the required P(ν) in ishem, as in, e.g. Drappeau et al. (2015).
3.1.2 Geometry
Usually jets are assumed to be conical, i.e. the radius of the jet r at a height z follows a simple linear relation. However, internal or external agents, e.g. a toroidal component of the magnetic field in the jet or the pressure exerted by the interstellar medium (see Section 5 for a discussion on the collimation agents) may collimate the jet and change its shape from conical to parabolic. More specifically, we can describe the geometry of the jet with a parameter ζ such as r ∝ zζ, where ζ = 1 for the ‘standard’ conical geometry, while ζ < 1 holds for a parabolic jet (see Fig. 3) and the latter is likely a physically more realistic description of the jet structure. The case of an ‘overpressured jet’ with ζ > 1 is plausible, but it would be a highly unstable structure that will tend to evolve spontaneously to a situation where ζ ≤ 1 (Kaiser 2006) .4 In the previous applications of ishem, the geometry was assumed to be conical and ζ was fixed to 1 by default. This is the first time that the dependence of the results on ζ is tested. In order to correctly take into account a non-conical geometry, the code used in this work has been updated with respect to the previous versions used by, e.g. Drappeau et al. (2015) and Péault et al. (2019). The new version of the model also includes some improvement of the treatment of radiation transfer, detailed in Appendix A. Some examples of simulations illustrating the effects of a non-conical geometry are shown in Fig. 4. As shown in the figure, reducing ζ allows for a more collimated jet, where the energy losses are reduced and the spectrum is naturally flatter. Furthermore, for a fixed length of the jet, the range of wavelengths emitted by the different regions will shrink as a result of the contained energy losses. This leads to the appearence of a low-frequency turn-over that marks a transition from the flat partially absorbed region of the SED to optically thick emission Fν ∝ ν5/2 at lower frequencies. This optically thick emission arises from the terminal part of the jet at the largest scale. As can be seen in Fig. 4, the low frequency termination turn-over gradually moves towards higher frequencies at lower ζ and might be observable if the jets are strongly confined. Lowering ζ results also in an increase of the overall flux emitted, as both the amount of energy lost and the frequency range over which jet power is distributed diminish. The frequency of the termination break in the SED of the source also depends on the size of the jet, which in our model corresponds to the distance to which the shells have been able to propagate during the time of the ishem simulation, i.e. ∼ctsimu ≃ 3 × 1015 cm in the spectra shown in Fig. 4.5 As shown in Fig. 5, in the case of strongly parabolic jets, increasing the simulation time pushes the spectral turnover towards lower frequencies, without affecting the shape of the spectrum above the turn-over frequency (see Fig. 6 for a distinction between the effect of reducing ζ and increasing tsimu). In contrast, in the conical jet model the SED is barely affected by the duration of the simulation.
Schematic view of a jet in both conical (black) and non-conical geometry (red).
Simulated SEDS with ISHEM using several values of ζ, from 1.0 to 0.5, without rescaling or shift factors applied, showing the impact of geometry on the flux and the shape of the emitted jet spectrum. In order to simulate these spectra, the X-ray PDS was used as input.
Simulated SEDS with ISHEM using several values for tsimu in both non-conical (left) and conical (right) geometry, showing how for strongly non-conical jets the low frequency turnover is dependent on the choice for tsimu. For both plots, tsimu was fixed to 105 s.
Simulated SEDS with ISHEM using several values for tsimu and ζ in order to show the separate effect of reducing ζ and increasing tsimu.
3.2 ISHEM parameters
A wide set of other physical parameters is taken into account by the code as well, but they can only shift in frequency or scale the normalization of the whole SED without modifying its shape. These parameters describe the system, the jet, and the distribution of the radiating particles. The main parameters are the distance (D) to the source, the inclination of the jet axis with respect to the line of sight (θ), the mass of the compact object (MNS in this case), the jet power (PJ), the jet opening angle (ϕ), the radius at the base of the jet (Rb), the average Lorentz factor (Γav) of the ejecta, the volume filling factor fV (Malzac 2013), and the maximum/minimum energy limits (γmax, γmin) of the electron distribution. As mentioned in subsection 1.1, the distance of the source is well constrained to be around 3.2 kpc (Kuulkers et al. 2010), while no constraints have been ever reported to our knowledge on the mass of the NS, which in the following will be fixed to 1.5 M⊙,6 or to the inclination of the jet θ, or to the inclination i of the system itself. We therefore tentatively fix θ to 60°. The jet power PJ is not known, so we fixed it to be of the same order of magnitude of the X-ray luminosity of the source, i.e. 0.01 LEdd (Migliari et al. 2010). The effects of PJ and θ will be explored in more detail in Section 5.
We chose fV = 0.7 (Malzac 2014) and Rb equal to 10 RG, which is plausible for XRBs, although the impact of these parameters on the overall results is negligible. For Γav, which for XRBs is expected to vary between 1 and 10 (see e.g. Casella et al. 2010; Saikia et al. 2019), we started with a value of 2 (Gallo, Fender & Pooley 2003; Heinz 2004). We adopted a value of ϕ of 2°, as opening angles are expected to be ≤10° (Miller-Jones, Fender & Nakar 2006) and a value of 2° has been used in the past7 (see e.g. Stirling et al. 2001). For the lower and upper limits of the electron distribution, we started with some standard values for XRBs, i.e. 10 and 106, respectively (Gandhi et al. 2011; Malzac 2014; Drappeau et al. 2015). A study of how the resulting simulated SED depends on these parameters is presented in Péault et al. (2019), fig. 2.
As already discussed in Section 3.1, in the case of non-conical jets (ζ < 1), the choice for tsimu becomes crucial because it determines both the size of the jet and the location of the low-frequency jet termination turn-over. We choose to set this parameter to 105 s, which corresponds to a final jet extension of ∼3 × 1015 cm. As the observed radio spectrum in 4U 0614+091 is rather flat reproducing the data with a strongly non-conical model will require to have the spectral turn-over well below 10 GHz. So, we want the jet to be as large as possibly allowed by the observational constraints. Regarding 4U 0614+091 the constraints on the jet extension are very poor. Namely, the jets from 4U 0614 + 091 should not be significantly bigger than ∼1017 cm, otherwise they would have been resolved with the VLA. In general, however, the observations of compact jets suggest much smaller dimensions for the radio-emitting region. In the case of the resolved jet of Cyg X-1, the extension of the radio jet at 8.4 GHz indicates scales of the order of 1014−1015 cm (see e.g. Stirling et al. 2001). Our choice for tsimu is therefore the most favourable for non-conical jet models while being still roughly compatible with the expected scale of the radio jet.
A summary of the parameters used in the simulation is reported in Table 2.
Parameters used in the ishem code that were kept fixed in all the simulations run in this paper. a: Simulation running time; b: Initial radius of the ejecta, imposed of approximately the same order of magnitude of the inner radius of the accretion disc; c: effective adiabatic index of the flow, (Malzac 2013).
| Simulation parameters | |
|---|---|
| MNS (M⊙) | 1.5 |
| t|$_{\rm simu}^a$| (s) | 105 |
| r|$_{\rm dyn}^b$| (RG) | 10 |
| ϕ (°) | 2.0 |
| fvol | 0.7 |
| |$\gamma _a^d$| | 4/3 |
| γaverage | 2.0 |
| Ejecta scheme | Constant shell mass |
| Pjet (LEdd) | 0.01 |
| γmin | 10 |
| γmax | 106 |
| θ (°) | 60 |
| Simulation parameters | |
|---|---|
| MNS (M⊙) | 1.5 |
| t|$_{\rm simu}^a$| (s) | 105 |
| r|$_{\rm dyn}^b$| (RG) | 10 |
| ϕ (°) | 2.0 |
| fvol | 0.7 |
| |$\gamma _a^d$| | 4/3 |
| γaverage | 2.0 |
| Ejecta scheme | Constant shell mass |
| Pjet (LEdd) | 0.01 |
| γmin | 10 |
| γmax | 106 |
| θ (°) | 60 |
Parameters used in the ishem code that were kept fixed in all the simulations run in this paper. a: Simulation running time; b: Initial radius of the ejecta, imposed of approximately the same order of magnitude of the inner radius of the accretion disc; c: effective adiabatic index of the flow, (Malzac 2013).
| Simulation parameters | |
|---|---|
| MNS (M⊙) | 1.5 |
| t|$_{\rm simu}^a$| (s) | 105 |
| r|$_{\rm dyn}^b$| (RG) | 10 |
| ϕ (°) | 2.0 |
| fvol | 0.7 |
| |$\gamma _a^d$| | 4/3 |
| γaverage | 2.0 |
| Ejecta scheme | Constant shell mass |
| Pjet (LEdd) | 0.01 |
| γmin | 10 |
| γmax | 106 |
| θ (°) | 60 |
| Simulation parameters | |
|---|---|
| MNS (M⊙) | 1.5 |
| t|$_{\rm simu}^a$| (s) | 105 |
| r|$_{\rm dyn}^b$| (RG) | 10 |
| ϕ (°) | 2.0 |
| fvol | 0.7 |
| |$\gamma _a^d$| | 4/3 |
| γaverage | 2.0 |
| Ejecta scheme | Constant shell mass |
| Pjet (LEdd) | 0.01 |
| γmin | 10 |
| γmax | 106 |
| θ (°) | 60 |
3.3 Simulations
4 SPECTRAL ANALYSIS
In NS LMXBs, the jet emission is expected to dominate only the radio-to-IR wavelengths, while the emission from optical to X-ray should be mainly ascribed to the accretion disc (since usually the radiation emitted by the faint companion is negligible). Therefore, we began with a separate analysis for the radio-to-IR data, fitted with ish. Then the optical-to-X-ray data were described mainly with discir, an irradiated disc plus Comptonization model (Gierliński, Done & Page 2008). However, although X-ray reprocessing from the outer disc or even direct emission from the outer disc is expected to dominate the near-infrared (NIR) – optical region in NS LMXBs (Russell et al. 2006; Russell, Fender & Jonker 2007), some level of contribution from the jet emission might still be present (several examples can be found in, e.g. Lewis et al. 2010; Harrison et al. 2011; Baglio et al. 2016, 2019). We therefore performed a fit of the whole data set, in order to check if accretion and ejection do dominate over two separate frequency ranges or either if there is a border territory, i.e. the NIR-optical region, where these phenomena cannot be easily singled out and have to be taken into account together.
4.1 From radio to IR: the jet emission
We first tried to reproduce the observed SED with a standard set of parameters, which are listed in Table 3.2, with p = 2.0, the X-ray PDS, and adopting the usual conical geometry (ζ = 1). We chose p = 2.0 according to the fit to the optically thin part of the jet spectrum by Migliari et al. (2010) and also to the standard diffusive shocks acceleration theory. We then tested the ‘synthetic’ spectrum on xspec, using the ish model built on it to fit the data. We also checked if using a one Lorentzian model instead of a double Lorentzian model, i.e. ignoring the component for the high-frequency QPO (see subsection 2.2) that could probably be related to the orbital motion of the system (Stella & Vietri 1998), could influence the results of the fit. We found that both models lead to the same results therefore in the following we will refer only to the results obtained including the QPO.
Even if IR fluxes are expected to be only slightly affected by interstellar reddening, we included the model redden, which estimates the extinction in the optical band, E(B−V). The latter was frozen to 0.5 (see Section 4.2) since it was left unconstrained by the fit. The outcome of the fit is quite poor, as witnessed by the resulting |$\chi ^2_\nu$| (d.o.f.) of 3.96 (5). Furthermore, using equations (1)–(2) to explore the parameters needed to improve the simulation, we found that, in order to have reasonable values for the jet power of the expected order of magnitude, i.e. 0.01 LEdd, one has to invoke oddly high opening angles (see Section 5). As mentioned in Section 3 the shape of the SED can be affected by only three elements: the shape of the electron distribution (which modifies the slope of the optically thin part of the spectrum), the geometry chosen and the dissipation pattern of the ejecta, in this case based on the X-ray PDS. Since the IR data are well fitted by the optically thin region of the synthetic SED, the chosen value of p seems to be correct, as expected. In the following, we will then try to change the geometry first and the PDS then in order to see if, with a different choice for these ingredients, we can still find a good model for the data. In order to check if a different value for ζ might improve the results of the fit, we run again the simulations with different ζ between 0.5 and 1.0, and we repeated the whole procedure. The resulting best fits as a function of ζ is shown in Fig. 7, (a) and (b). The results of each fit is reported in Table 3. The values of ζ for which we have the lowest |$\chi ^2_{\nu }$| are 0.57 and 0.6 (1.08 and 1.13, respectively, both with 5 degrees of freedom). A paraboloidal jet with ζ in this range of values therefore represents an acceptable fitting scenario, contrary to the typical conical geometry.
Simulated SEDS with ISHEM using several values of ζ, from 1.0 to 0.6 (panel a) and in the critical region 0.5 to 0.6 (panel b), normalized and shifted in order to fit the data. In all these models, the X-ray PDS was used as an input for ISHEM. Panel b shows the interesting evolution of the SED for values of ζ spanning in the crucial region between ζ = 0.6 and ζ = 0.5: when the jet becomes too collimated, the contribution from the external regions of the jet (the lower frequencies) becomes dominant, and it leads again to an inverted spectrum.
Results of the VLA-Spitzer/IRAC fits with ish with different ζ values. In all the fits, the number of degrees of freedom is equal to 5.
| Fit results for different jet geometries | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| ζ | |||||||||
| 0.50 | 0.53 | 0.55 | 0.57 | 0.60 | 0.70 | 0.80 | 0.90 | 1.00 | |
| |${\bf \chi ^2_{\nu }}$| | 5.95 | 2.34 | 1.46 | 1.08 | 1.13 | 2.11 | 2.84 | 3.37 | 3.96 |
| Fit results for different jet geometries | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| ζ | |||||||||
| 0.50 | 0.53 | 0.55 | 0.57 | 0.60 | 0.70 | 0.80 | 0.90 | 1.00 | |
| |${\bf \chi ^2_{\nu }}$| | 5.95 | 2.34 | 1.46 | 1.08 | 1.13 | 2.11 | 2.84 | 3.37 | 3.96 |
Results of the VLA-Spitzer/IRAC fits with ish with different ζ values. In all the fits, the number of degrees of freedom is equal to 5.
| Fit results for different jet geometries | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| ζ | |||||||||
| 0.50 | 0.53 | 0.55 | 0.57 | 0.60 | 0.70 | 0.80 | 0.90 | 1.00 | |
| |${\bf \chi ^2_{\nu }}$| | 5.95 | 2.34 | 1.46 | 1.08 | 1.13 | 2.11 | 2.84 | 3.37 | 3.96 |
| Fit results for different jet geometries | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| ζ | |||||||||
| 0.50 | 0.53 | 0.55 | 0.57 | 0.60 | 0.70 | 0.80 | 0.90 | 1.00 | |
| |${\bf \chi ^2_{\nu }}$| | 5.95 | 2.34 | 1.46 | 1.08 | 1.13 | 2.11 | 2.84 | 3.37 | 3.96 |
As mentioned in subsection 3.2, for strongly non-conical geometries, spectra obtained with longer simulation times could produce flatter spectra. We show in Table 4 the results of the fits for three values of ζ, i.e. 0.53, 0.60, and 0.7, with tsimu = 1 × 105 and tsimu = 3 × 105. As expected, the values of |$\chi ^2_\nu$| are generally lowered by increasing the simulation time, with the exception of the fit with ζ = 0.7, which is mostly unaffected by changing tsimu. Therefore, we found out that even with higher simulation times, the best fit is obtained again for ζ around 0.6. While an even higher simulation time would still be physically acceptable (see subsection 3.1), it would very likely only confirm the results presented here with shorter and more feasible tsimu.
Results of the VLA–Spitzer/IRAC fits with ish with different ζ and tsimu. In all the fits, the number of degrees of freedom is equal to 5.
| Fit results for different jet geometries and simulation times | ||||||
|---|---|---|---|---|---|---|
| ζ | ||||||
| 0.50 | 0.60 | 0.70 | ||||
| tsimu (× 105 s) | 1 | 3 | 1 | 3 | 1 | 3 |
| |${\bf \chi ^2_{\nu }}$| | 2.34 | 1.36 | 1.13 | 0.80 | 2.11 | 2.13 |
| Fit results for different jet geometries and simulation times | ||||||
|---|---|---|---|---|---|---|
| ζ | ||||||
| 0.50 | 0.60 | 0.70 | ||||
| tsimu (× 105 s) | 1 | 3 | 1 | 3 | 1 | 3 |
| |${\bf \chi ^2_{\nu }}$| | 2.34 | 1.36 | 1.13 | 0.80 | 2.11 | 2.13 |
Results of the VLA–Spitzer/IRAC fits with ish with different ζ and tsimu. In all the fits, the number of degrees of freedom is equal to 5.
| Fit results for different jet geometries and simulation times | ||||||
|---|---|---|---|---|---|---|
| ζ | ||||||
| 0.50 | 0.60 | 0.70 | ||||
| tsimu (× 105 s) | 1 | 3 | 1 | 3 | 1 | 3 |
| |${\bf \chi ^2_{\nu }}$| | 2.34 | 1.36 | 1.13 | 0.80 | 2.11 | 2.13 |
| Fit results for different jet geometries and simulation times | ||||||
|---|---|---|---|---|---|---|
| ζ | ||||||
| 0.50 | 0.60 | 0.70 | ||||
| tsimu (× 105 s) | 1 | 3 | 1 | 3 | 1 | 3 |
| |${\bf \chi ^2_{\nu }}$| | 2.34 | 1.36 | 1.13 | 0.80 | 2.11 | 2.13 |
We also notice that our conclusion on the best value of ζ should not be considered conclusive, as further investigations in the range between ζ =0.53 and ζ =0.6, with possibly higher tsimu could, in principle, lead to even more accurate estimates of the best ζ. However, a precise estimate of ζ goes beyond the scopes of this work and would likely provide no or very little improvement to the results.
We therefore conclude that a non-conical geometry, with ζ around 0.6 and possibly even below, improves significantly the fit with ishem using the X-ray PDS.
We then tested the other possible scenario, where the observed X-ray variability does not reflect the fluctuations of ejection velocity, using a ‘flicker-noise’ PDS in conical geometry. We assumed an rms fractional amplitude of 30 per cent and included a range of frequencies ranging from f1 = 10−5 Hz to f2 = 103 Hz. Using the same set of parameters shown in Table 2, we obtained a synthetic SED that is in quite good accordance with the data, i.e. |$\chi ^2_\nu$| (d.o.f.) = 0.27(5). The best-fitting model is shown in Fig. 8, in comparison with the data and a pair of best-fitting models obtained with the X-ray PDS and variable values of ζ.
Simulated SEDS with ISHEM for varying geometries and input PDS as compared with the radio-to-IR data set.
4.2 From optical to X-ray: the disc emission
While in the previous section we focused on the part of the SED dominated by the jet, in this section we will focus on the optical-to-X-rays data, which are expected to be dominated by the disc and the hot corona emission.
First of all we used discir (Gierliński et al. 2008), which includes the disc emission, Comptonization from a hot corona, and the X-ray illumination of the disc, relevant in the presence of data coverage in the optical-UV domain (as in our case). The main parameters of the model are the temperature of the disc at its inner radius kTdisc, the Γ index of the power law reproducing the Comptonization spectrum, the electron temperature of the corona kTe, the ratio LC/LD between the luminosity of the Comptonized emission and the disc luminosity, the fraction of the flux of the Compton tail that is thermalized in the inner and outer radius (fin and fout, respectively), the radius of the illuminated disc rirr, the outer disc radius log rout (both in units of the inner disc radius) and the normalization K, which can be used to derive the inner disc radius.
The value LC/LD is generally used as an indicator of the spectral state and it is usually higher than 1 in intermediate/hard states. In the following fits, we fixed fin and rirr to the standard values of, respectively, 0.1 and 1.1 (Gierliński et al. 2008). Furthermore, as 4U 0614 + 091 is known to be an ultracompact binary, we fixed log rout to the value of 3. This value is reasonable considering that, since the system has an orbital period of about 50 min, the orbital separation is expected to be around 3 × 105 km and a factor of 103 guarantees that even with large inner disc radii, say 100 RG, as expected in hard state, it is larger than the extension of the disc. We also included a blackbody model in the spectrum (bbody in Xspec), already present in the fit by Migliari et al. (2010), which accounts for the emission of the NS surface (or the boundary layer). Finally, we included the two Gaussian components used by these authors, i.e. at 0.67 keV for the O VIII line and at 6.6 keV for the Fe–K fluorescence line. We applied to the model the components redden and tbabs to take into account interstellar extinction in both the UV and X-ray band. Finally, we included constant to serve as a cross-calibration constant and we checked that its value was always around 1, i.e. in the range 0.8-1.2.
The fit did not constrain LC/LD, which was therefore tentatively fixed to 10, which is a reasonable value for hard states (see e.g. Del Santo et al. 2008). The fit is unable to constrain the parameters of the gaussian component at 6.6 keV, due to a marginal contribution of this component to the fit, i.e. a 9 per cent probability of improvement by chance (calculated via ftest). This is not surprising, since this is not the first time that the iron line in 4U 0614 + 091 was found very weak (see e.g. Piraino et al. 1999). In the following, we will therefore not include this Gaussian component. On the other hand, we confirmed the presence of a broad line at ∼ 0.67 keV, likely associated with O viii (see Section 1.1 for references). The fit provides E(B − V; from redden) of about 0.5, which corresponds to a value8 of AV slightly lower than the one reported in Migliari et al. (2010), equal to 2. We were able to find only a relatively high lower limit to the corona temperature, i.e. 110 keV, which might be probably due to a the lack of a proper modelling of the high-energy hard tail (see e.g. Di Salvo et al. 2001; Iaria et al. 2001; D’Aí et al. 2007; Del Santo et al. 2013, and references therein) rather than such a high-temperature electron plasma. The results of the fit led as well to a significantly colder disc with respect to Migliari et al. (2010), i.e. kTdisc < 0.1 keV, but correlations with the disc normalization K and/or with the imposed values of LC/LD might be at play. The normalization of the disc K is bound to the inner radius of the disc Rin by the relation: |$K=(R^2_{\rm in}/D^2_{\rm 10 kpc})\times \cos {i}$|, with D10kpc is the distance of the system in units of 10 kpc. We found an apparent inner disc radius of ∼ 160 km (∼ 73 RG), which has to be taken as a lower limit since it does not take into account the proper (unknown) inclination of the system and the correction factor (for a more detailed calculation of the inner disc radius based on K, see e.g. Marino et al. 2019, and references therein). The few differences in our results with respect to the results obtained by Migliari et al. (2010) on the same X-ray data are likely due to either the inclusion, in our data set, of the SMARTS–UVOT data or to the different models used in the two papers.
4.3 Global multiwavelength analysis
The spectral analyses conducted in subsections 4.1 and 4.2 allowed us to characterize separately the emission from the jet and from the disc, under the hypothesis that their spectral domains were independent. In the following, we report on the global data set fitted with both accretion and ejection models, in order to confirm the previous assumption and provide a final, multiwavelength study of the spectral energy distribution of 4U 0614 + 091.
Starting from the best-fitting model for the optical-to-X-ray data, whose main parameters are reported in Table 5, we included the VLA and Spitzer data used in SubSection 4.1 and added ish to the spectral model employed. Since it is not possible to exclude a priori that our results on the jet might have been biased by the lack of higher frequencies data, we performed several fits trying different ISHEM models with ζ spanning from from 1.0 to 0.5. Unfortunately each fit results in approximately the same |$\chi ^2_{\nu }$| (d.o.f.) value of 1.49(451). Similarly, using the ISHEM model based on the ‘flicker noise’ PDS, the fit results in a |$\chi ^2_{\nu }$| (d.o.f.) value of 1.43(451). This situation is not surprising because the fit is indeed dominated by the higher statistics data in the X-ray and only slightly affected by the modelling of the few data points in the radio-IR domain. Leaving the same model, i.e. including discir and bbody, but neglecting all the data but radio and IR, results in fits which are similar to the fits performed in SubSection 4.1: among the different fits with X-ray PDS with varying ζ, the best fit is obtained again with ζ = 0.57, i.e. for which the fit goes from χ2/d.o.f. = 671.9/451 (including optical-to-X-rays data) to χ2/d.o.f. = 5.25/5., while choosing the flicker noise PDS the fit goes from χ2/d.o.f. = 644.9/451 to a χ2/d.o.f. of 3.55/5. However, the extension of the data set does have some effects on the ISHEM best-fitting parameters, i.e. the shift frequency and re-normalization factors, which are different with respect to the previous set of fits. In the next section, we use equations (1) and (2) in order to check whether the best-fitting shifts and normalization allows for ’reasonable’ physical parameters of the jet. Finally, we checked if the results are dependent on the values assumed by the cross-calibration constant for the IR and the radio data, which were assumed to be equal to 1.0 for the radio and IR data; leaving the cross-calibration parameters free to vary between 0.8 and 1.2 does not change significantly the values obtained for the |$\chi ^2_{\nu }$| (d.o.f.). In the following, we will therefore only refer to the fits with the parameters fixed to 1.0. We refer to Fig. 9 for the SED, overimposed to the best-fitting ISHEM model. The best-fitting parameters of the accretion flow model found with the lower frequencies extension of the data set are all perfectly compatible with the results reported in Table 5.
Best-fitting unabsorbed spectral energy distribution compared with the whole multiwavelength data set available for 4U 0614 + 091, in both Flux Density (Left) and ν × Fν (Right) representation. The optical-to-X-ray data set has been analysed with an irradiated disc + blackbody model, while the radio-to-IR data set was modelled with ISHEM. In this plot, we show both best-fitting ISHEM models found in this paper: on Top the model built with the X-ray PDS and corresponding to a non-conical geometry (ζ = 0.6) and on (bottom) the model corresponding to a ‘flicker noise’ PDS in conical geometry.
Fit results of the disc-dominated SED region, with data from SMARTS, UVOT, XRT, PCA, and HEXTE. Quoted errors reflect 90 per cent confidence level. The parameters that were kept frozen during the fits are reported between round parentheses.
| Spectral analysis | |||
|---|---|---|---|
| redden | E(B − V) | 0.48|$^{+0.07}_{-0.06}$| | |
| tbabs | NH | × 1022 cm−2 | 0.21 ± 0.02 |
| discir | kTdisc | keV | 0.077 ± 0.003 |
| kTe | keV | >110 | |
| Γ | 2.24 ± 0.02 | ||
| LC/LD | (10) | ||
| fin | (0.10) | ||
| fout | (× 10−3) | 2.0|$^{+9.0}_{-1.3}$| | |
| |$R_{\rm in}/\sqrt{\cos {i}}$| | (RG) | 73|$^{+5}_{-4}$| | |
| Rirr/Rin | (1.01) | ||
| Rout/Rin | (103) | ||
| Gaussian | Eline | (keV) | 0.671 ± 0.009 |
| σ | (keV) | 0.077 ± 0.006 | |
| bbody | kTbb | (keV) | 1.38 ± 0.02 |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(442) | |||
| Spectral analysis | |||
|---|---|---|---|
| redden | E(B − V) | 0.48|$^{+0.07}_{-0.06}$| | |
| tbabs | NH | × 1022 cm−2 | 0.21 ± 0.02 |
| discir | kTdisc | keV | 0.077 ± 0.003 |
| kTe | keV | >110 | |
| Γ | 2.24 ± 0.02 | ||
| LC/LD | (10) | ||
| fin | (0.10) | ||
| fout | (× 10−3) | 2.0|$^{+9.0}_{-1.3}$| | |
| |$R_{\rm in}/\sqrt{\cos {i}}$| | (RG) | 73|$^{+5}_{-4}$| | |
| Rirr/Rin | (1.01) | ||
| Rout/Rin | (103) | ||
| Gaussian | Eline | (keV) | 0.671 ± 0.009 |
| σ | (keV) | 0.077 ± 0.006 | |
| bbody | kTbb | (keV) | 1.38 ± 0.02 |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(442) | |||
Fit results of the disc-dominated SED region, with data from SMARTS, UVOT, XRT, PCA, and HEXTE. Quoted errors reflect 90 per cent confidence level. The parameters that were kept frozen during the fits are reported between round parentheses.
| Spectral analysis | |||
|---|---|---|---|
| redden | E(B − V) | 0.48|$^{+0.07}_{-0.06}$| | |
| tbabs | NH | × 1022 cm−2 | 0.21 ± 0.02 |
| discir | kTdisc | keV | 0.077 ± 0.003 |
| kTe | keV | >110 | |
| Γ | 2.24 ± 0.02 | ||
| LC/LD | (10) | ||
| fin | (0.10) | ||
| fout | (× 10−3) | 2.0|$^{+9.0}_{-1.3}$| | |
| |$R_{\rm in}/\sqrt{\cos {i}}$| | (RG) | 73|$^{+5}_{-4}$| | |
| Rirr/Rin | (1.01) | ||
| Rout/Rin | (103) | ||
| Gaussian | Eline | (keV) | 0.671 ± 0.009 |
| σ | (keV) | 0.077 ± 0.006 | |
| bbody | kTbb | (keV) | 1.38 ± 0.02 |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(442) | |||
| Spectral analysis | |||
|---|---|---|---|
| redden | E(B − V) | 0.48|$^{+0.07}_{-0.06}$| | |
| tbabs | NH | × 1022 cm−2 | 0.21 ± 0.02 |
| discir | kTdisc | keV | 0.077 ± 0.003 |
| kTe | keV | >110 | |
| Γ | 2.24 ± 0.02 | ||
| LC/LD | (10) | ||
| fin | (0.10) | ||
| fout | (× 10−3) | 2.0|$^{+9.0}_{-1.3}$| | |
| |$R_{\rm in}/\sqrt{\cos {i}}$| | (RG) | 73|$^{+5}_{-4}$| | |
| Rirr/Rin | (1.01) | ||
| Rout/Rin | (103) | ||
| Gaussian | Eline | (keV) | 0.671 ± 0.009 |
| σ | (keV) | 0.077 ± 0.006 | |
| bbody | kTbb | (keV) | 1.38 ± 0.02 |
| |${\bf \chi ^2_{\nu }}(d.o.f.)$| = 1.49(442) | |||
5 DISCUSSION
The analysis limited to the radio-IR domain carried out in subsection 4.1 suggests two possible scenarios to describe the jet emission for 4U 0614 + 091 within the internal shocks scenario: on one hand the variability in the Lorentz factors of the ejecta is related to the X-ray variability, i.e. the variability in the accretion flow, but the jet is non-conical (in the following scenario a), on the other hand it is also possible that, at least in this source, a flicker noise power spectrum, unrelated to the X-ray variability, is a better proxy for the fluctuations of the jet Lorentz factor (scenario b). Including the optical-to-X-rays portion of the data set in subsection 4.3 turned out to be furthermore inconclusive in discerning between the two proposed scenarios. This is mostly due to the lack of data around 1011–1012 erg s−1. The availability of data in this region, e.g. from the Atacama Large Millimeter/submillimeter Array (ALMA), would have been crucial to distinguish between the two scenarios. For example, ALMA data for other NS LMXBs in the hard state (see e.g. the SEDs shown by Díaz Trigo et al. 2017, 2018) indicate a small rise in flux density in this region which, if observed also in our data set, would have likely favored scenario a.
We will instead explore and discuss both scenarios a and b in the following section, on the basis of the global fit reported in subsection 4.3. In both cases, there is evidence that it is not possible to completely disentangle the disc contribution to the IR domain or the jet contribution to the optical domain as well. In particular, it results that there is a jet contribution varying from 30 per cent to 6 per cent in the SMARTS wavelengths range within scenario a, while this contribution is less prominent in scenario b, i.e. from 20 per cent to 3 per cent. These results confirm the study led by Russell et al. (2006), according to which in NS LMXBs the main emission process to be taken into account in the optical domain is the X-ray reprocessing from the disc,9 even though they also point out how the jet contribution might not be negligible at all.
5.1 Scenario a: a non-conical jet?
Due to the high statistics in the X-rays, the fit to the whole data set gives comparable |$\chi ^2_\nu$| values for each of the tested geometries. In order to determine the most likely values of ζ we use equations (1) and (2) to convert the values of the shift and renorm factors found by ish in couples of PJ-ϕ. We explored also how these results were influenced by our choice for θ and Γav, allowing for both to change in some physically reasonable ranges. In particular, since the scaling relations used in this work, i.e. equations (1) and (2), are valid in the approximation of ϕ0 ≲ θ (see B2), with ϕ0 the opening angle at the base of the emitting region, we fixed the lower limit for the inclination θ to 8°. On the other hand, Γav was allowed to vary in the 2–10 range. Running these tests draw areas of possible results in a PJ-ϕ plot, corresponding to specific values of ζ. We show in Fig. 10 the resulting area for ζ = 1.0 (red) and ζ = 0.57 (blue). In each of the resulting skewed areas, the bottom of the areas corresponds to θ = 8°, the top corresponds to θ = 90°, while Γav increases from right (where Γav = 2) to left (Γav = 10). The sub-regions where the condition θ < ϕ0 are coloured in grey and they have to be excluded. It is important to notice that for a non-conical geometry, the opening angle depends on distance z along the jet. The emitting region, located at about 1 au from the base of the jet, will have therefore a different opening angle than the values encompassed by the blue area in Fig. 10. Thereby, this angle can violate the condition θ < ϕ. Indeed, for a fixed geometry parameter ζ = 0.57, using the opening angle at a distance of 1 au ϕ0 results in significantly smaller angles, i.e. below 0.1°, as shown in Fig. 11. Such extreme values are not implausible, as small opening angles have been suggested for jets in XRBs (see e.g. Zdziarski et al. 2016). Both the ranges of jet powers and opening angles individuated by the ζ = 0.57 areas can be accepted and a ζ ≈0.6 value appears still consistent with the best physically motivated scenario.
Jet power – (observed) opening angles curve (from equations 1 to2) for a different combinations of ζ and PDS to a range of possible values for the inclination θ and a range of Γav between 2 and 10. In particular, in each ‘region’ (identified by a specific colour) the inclination increases going upwards, while the Lorentz factor increases going from right to left. Curves for specific fixed values of Γ and θ are identified to help the eye. In each area, the sub-region for which the condition θ > ϕ0 is not satisfied are coloured in grey and have to be discarded, as the used scaling relations are no longer valid. In this plot, the black solid line indicates the expected jet power, while the dashed black line points out the upper limit for the opening angle.
Jet power – opening angle curves for the particular case of ζ = 0.57, using the opening angles at the base of the jet (blue) and the opening angles hypothetically observed at 1 au (magenta). We refer to the caption of Fig. 10 for further details.
For ζ = 1.0, on the contrary, the area found by this procedure does not allow for reasonably small opening angles and requires very high jet powers for the whole range of explored Γav. This represents another point in favour of ruling out the conical geometry scenario.
The concept of a non-conical jet is not groundbreaking: it is known that the conical geometry is an approximation of the real geometry, as valid and efficient as it proved to be. The confinement agent necessary can be internal or external. In the first case, it has been proposed by several authors that collimation might be due to a toroidal component of the magnetic field that increases along the axis and that forces the jet to decrease its opening angle (see e.g. Heyvaerts & Norman 1989; Pudritz, Rogers & Ouyed 2006; Pudritz, Hardcastle & Gabuzda 2012). However, this mechanism has been questioned by Spruit (2010), according to which a magnetic self-confinement of the jet is not physically possible as the toroidal magnetic pressure within the jet would force them to expand. On the other hand, the collimation necessary to ‘break’ the (unstable, Martí, Perucho & Gómez 2016) conical geometry might be exerted for instance by the interstellar medium (e.g. Asada & Nakamura 2012) or an external magnetic field kept in place by the disc (e.g. Spruit, Foglizzo & Stehle 1997).
It is interesting to compare these results with those previously found by applying the same ISHEM model to jets in BH-XRBs, where the conical geometry assumption worked correctly. In a few cases (see e.g. the radio residuals for some spectra in Péault et al. 2019, fig. 3) it is reasonable that flatter models could even improve the already acceptable fits and in this sense a non-conical geometry could be necessary. Any comparison between NS and BH XRBs is therefore premature for at least two reasons: a non-conical geometry was not tested for BH XRBs and, also, we need to test more NS LMXBs to draw any conclusion on a possible difference between jets in these two classes of systems.
5.2 Scenario b: X-ray variability is not a good proxy for the Lorentz factor fluctuations
In the second scenario, the dissipation pattern of the shells internal energy in the jet is not related to the X-ray timing properties, i.e. the timing properties of the accretion flow, but it is mainly due to ‘flicker noise’. The plausibility of this scenario is confirmed by the PJ-ϕ diagram in Fig. 10, since the corresponding area encompasses the expected range of jet powers-opening angle.
This is not groundbreaking either (we refer again to, e.g. Jamil et al. 2010; Malzac 2013) but it would be certainly different to the results obtained on the other sources to which the ISHEM model has been applied in the past. In this case, the fact of having an NS instead of a BH might play a role. Under the hypothesis of a disc–jet coupling, the variability in the emission from the NS/boundary layer may not be transmitted to the ejecta in the jet, breaking subsequently the connection between the ejection pattern of the shells and the X-ray PDS. Alternatively, one might also consider differences in the jet launching mechanism in NSs with respect to BHs. For instance Parfrey, Spitkovsky & Beloborodov (2016) shows how the interaction between a fastly rotating low-magnetized NS and the disc may lead to a state where the magnetic field lines are open and provide the energy necessary for the ejection of particles. In this case, we do not expect that the dissipation pattern in the ejecta and the accretion flow fluctuations in the disc to be exactly matched. Such a mechanism could be at work in Accreting Millisecond X-ray Pulsars and analogous systems, which might possibly include 4U 0614 + 091. Indeed, the system, with a 415 Hz frequency spin (see e.g. van Doesburgh & van der Klis 2017, and references therein), belongs to the family of binaries hosting millisecond NSs. However, the NS magnetic field is likely buried, as witnessed by the lack of observed X-ray pulsations, and this would make the attribution of the aforementioned mechanism to the system unlikely.
We also suggest the possibility that the lack of correlation between X-ray variability and ejecta in the jet proposed here for 4U 0614 + 091 may not necessarily be representative for the whole class of NS LMXBs in hard state. As apparent from the VLA data used in this paper, the jet spectrum looks quite flat, which, we recall from subsection 3.1, requires also an almost flat PDS. The X-ray PDS used has evidently not the required shape and it is therefore not surprising that it may not be the best tracer for the variability in the shells velocity. Similar PDS have been observed frequently in the so-called atoll LMXBs, when in island state (IS; see e.g. van Straaten et al. 2002; van Straaten, van der Klis & Méndez 2003). On the other hand, PDS dominated by broad, flat-topped noise, similar to what observed in BH XRBs in hard state, have been found also in several atoll sources at low luminosity (Belloni, Psaltis & van der Klis 2002; Reig, van Straaten & van der Klis 2004; van Straaten, van der Klis & Wijnands 2005). This variability behaviour has been classified as extreme island state (EIS; Méndez & van der Klis 1997; van Straaten et al. 2003)10 and systems in this state tend to populate a horizontal extension of the IS region in a colour–colour diagram, corresponding to harder spectra than sources in IS (Gierliński & Done 2002; Muno, Remillard & Chakrabarty 2002). Indeed spectra of NS LMXBs in EIS are typically described with power laws of Γ ∼ 1.8 (see e.g Barret et al. 2000; Linares et al. 2008), unlike e.g 4U 0614 + 091 that displays spectra usually steeper, Γ ∼ 2.2−2.4, as in this paper and, e.g. Piraino et al. (1999), Migliari et al. (2010), and Ludlam et al. (2019). In addition, the X-ray variability is significantly stronger (i.e. 30–40 per cent rms amplitude, Wijnands et al. 2017) in EIS than sources in IS. The ensemble of these clues suggests that atoll sources in EIS are likely associated with a physical scenario where the disc is truncated far from the compact object and the NS surface is not very hot (Reig et al. 2004; Bult et al. 2018), while in IS the contribution from the disc and/or from the NS increases, cools down the corona, and reduces the X-ray variability. The distinction between sources in IS and EIS is not strict, and some sources, like 4U 0614 + 091 itself, have been found in both states (see panel 1 in fig. 2 of van Straaten et al. 2002).
Assuming that radio jet spectra in NS LMXBs are usually flat,11 it is plausible that X-ray PDS can be used as proxy of the variability in the ejecta for NS LMXBs in EIS. In this sense, AMXPs and/or low-luminosity bursters, usually found in EIS, could be good candidates to test ISHEM in the future.
6 CONCLUSIONS
In this work, we presented the first ever attempt to describe the broad-band emission of an NS LMXB, i.e. 4U 0614 + 091, with a model taking into account both the jet and the accretion flow emission. We took advantage of the same multiwavelength data set presented by Migliari et al. (2010), with the only exception for the Swift/UVOT data, which were re-analysed. We modelled the radio-to-IR spectrum with the ISHEM code, which calculates the expected spectral energy distribution in the low-energy part of the SED taking into account the X-ray variability (connected in turn to the internal shocks temporal pattern). In particular, we used the quasi-simultaneous Swift/XRT PDS as input for the ‘synthetic’ SED. While the ISHEM model has been applied several times in the past to XRBs hosting BHs as the accreting object, this is the first time that the model is applied to a system hosting an NS. In addition, optical-to-X-ray data were modelled with an irradiated disc model.
We found that the compatibility between the SED built using the X-ray PDS and the data set is critically dependent on the geometry of the jet, enclosed in the geometrical parameter ζ. In particular, a highly non-conical geometry, with ζ ≈ 0.6, results in an acceptable fit. Alternatively, an acceptable fit is found within a conical geometry scenario but using in input a ‘flicker-noise’ PDS instead of the X-ray PDS. This scenario might imply that for NS LMXBs the X-ray PDS are not good tracers for the fluctuations in the Lorentz factors of the ejecta, possibly due to some contribution from the boundary layer/NS emission. The scarce statistics does not allow for the moment to choose one scenario over the other. New observations and/or further studies such as the one presented here are definitely necessary to provide an answer to this issue and in general for a better understanding of the accretion–ejection coupling in NS LXMBs.
ACKNOWLEDGEMENTS
This work received financial support from PNHE ("Programme Nationale Hautes Energies") in France and from the OCEVU (Origins, Constituents and EVolution of the Universe) Laboratory of Excellence (ANR-11-LABX-0060) and the "Aix Marseille Université Initiative d’Excellence" (A*MIDEX) project (ANR-11-IDEX-0001-02) funded by the ‘Investissement d’Avenir’ French government program managed by the ANR ("Agence Nationale de la Recherche"). We acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.0 and INAF main-stream (P.I. Belloni). JLM acknowledges the support of a fellowship from ‘La Caixa’ Foundation (ID 100010434). The fellowship code is LCF/BQ/DR19/11740030. MP acknowledges financial support from the Spanish Ministry of Science and Technology through grants PID2019-105510GB-C31, PID2019-107427GB-C33, and AYA2016-77237-C3-3-P, and from the Generalitat Valenciana through grant PROMETEU/2019/071.
Data availability
The data utilized in this article are publicly available at https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl, while the analysis products will be shared on reasonable request to the corresponding author.
Footnotes
The largest available data base of X-ray/radio observations of XRBs is consultable at https://github.com/bersavosh/XRB-LrLx_pub, while for the most recently published plot see Bassi et al. (2019), fig. 7.
Low-luminosity NS LMXBs that exhibit mainly two spectral states, one hard, dubbed ‘island’ state, and one soft, dubbed ‘banana state’ (Hasinger & van der Klis 1989).
Reported in https://www.swift.ac.uk/analysis/uvot.
Note that in the real world, the extension of the jet depends not only on the time since the ejection started, but also on the interaction of the jets with their ambient medium at large scales that is not modelled here.
Which is close to the peak for recycled NSs in the expected NSs mass distribution (Özel et al. 2012).
However, it was recently suggested that jet opening angles in XRBs could be even smaller, below 1° (Zdziarski et al. 2016).
Contrarily to BHs, where the jet emission is usually extended until the optical wavelengths (see e.g. Péault et al. 2019).
We refer the reader to fig. 2 in Wijnands, Degenaar & Page (2017) for a direct comparison between the two types of PDS.
Aside of 4U 0614 + 091, a couple of other examples can be found in Díaz Trigo et al. (2017).
REFERENCES
APPENDIX A: THE UPDATED VERSION OF ISHEM
The ishem code is extensively detailed in Malzac (2014). However, in this paper we use an updated version of ishem in which radiation transfer was improved in order to account more accurately for the geometrical and relativistic aberration effects on the synchrotron process. The main effects of these modifications is to change the normalization of the predicted SED by at most a factor of a few compared to the previous version. The shape of the predicted SED is not significantly affected (see Fig. A1). Although we expect that the new version is more accurate and provides better estimates of the jet parameters when compared to the data, from a qualitative point of view, the resulting parameters are comparable to those obtained with the previous version. The main changes in the code are described below.
Comparison between the SED simulated with the previous (red) and the updated (blue) version of ISHEM. In both simulations, ζ = 1.0 and tsimu = 105 s.
A1 Emission and absorption coefficients
In the version of ishem presented in Malzac (2014), we used the synchrotron emission and absorption coefficient given in Chaty, Dubus & Raichoor (2011). These estimates are for a uniform magnetic field observed with a specific line of sight that is perpendicular to the magnetic field (Rybicki & Lightman 1986). In the new version of the code, it is instead assumed the field is tangled on scales larger than the Larmor radius and smaller than the emitting region. This constitutes a better approximation of the magnetic field in a shocked region and the result is valid for any viewing angle (Crusius & Schlickeiser 1986). We detail below the formulae that we used for the synchrotron emissivity and absorption coefficients in a tangled magnetic field.
These expressions for emission and absorption coefficients are equivalent to those used in a different form by Zdziarski et al. (2012).
A2 Emission from a homogeneous cylinder in motion
In ishem, the jet is discretized into a large number of homogeneous cylinders whose axes are along the jet axis. These cylindric shells travel along the jet while expanding radially according to the fixed jet geometry. The time-dependent emission from each of these cylinders is calculated by ishem. In order to predict the time-averaged SED of the jet, it is necessary to calculate the time-integrated emission of millions of such cylinders. For reasons of computational efficiency, the radiation transfer has to be simplified.
A2.1 In the rest frame
A2.2 In the observer’s frame
A2.3 Effects of absorption by other shells
APPENDIX B: SCALING LAWS FOR PARABOLIC JETS
In this section, we derive the scaling laws for the synchrotron spectral break frequency and flux that we used to estimate the best fits ishem parameters shown in Figs 10 and 11.
B1 Synchrotron emission of self-similar parabolic jets
First, we are have to estimate the SED of standard self-similar parabolic jets, i.e. jets with radius increasing with height like R = R0(z/z0)ζ, where z0 and R0 are the height and radius at the base of the jet emitting region.
In the following, we assume that ζbd > 1 and zf > >z0.
B2 Local cylinder approximation
This approximation is expected to be accurate at large inclination angles (θ ∼ 90°) where the observed emission is dominated by radiation travelling in the radial direction and does not experience significant gradients in the jet. In fact, it turns out to be remarkably accurate even at smaller viewing angles and up to viewing angles comparable to the jet opening angle (i.e. up to tan θ ∼ R0/Z0). For smaller jet inclinations, a different approach must be adopted (see e.g. Zdziarski et al. 2016).
B3 Scaling laws in the internal shock model
These scaling laws are expected to constitute a good approximation for inclinations θ > ϕ0 where ϕ0 is the jet opening angle at the base of the jet emitting region tan ϕ0 = R0/Z0. Since in the case of 4U 0614 + 091 the first internal shocks occur around z0 ∼ 103Rg, and in all our models we have set Rb =10 RG the approximation is expected to be valid for inclinations such that tan θ > 10−2(1−ζ)tan ϕζ, where ϕ is the jet opening angle at Rb. As they are independent of the magnetic field profile b, the scaling relations are valid for a wide range of dissipation profile along the jet. They extend the formulae used in Péault et al. (2019) to non-conical geometries. The angle dependence of the scaling relations are also improved with respect to the formulae of Péault et al. (2019) in order to reflect the refined treatment of the anisotropy of the jet radiation implemented in the new version of ishem. Note, however, that they do not take into account the possible contribution of the counter jet.
