© ESO 2021

1. Introduction

A black hole (BH) X-ray binary (XRB) is a binary system composed of a stellar mass BH accreting matter from a companion star (see Done et al. 2007, for a general review). Almost all known black hole XRBs are transients (BHTs): they spend most of the time in a quiescent state, at low X-ray luminosity (i.e., below ∼10³¹ erg s⁻¹), but can display sudden episodes of increased X-ray activity called outbursts (with X-ray luminosity up to ∼10³⁸ − 10³⁹ erg s⁻¹ at the peak). When in outburst these systems can be found in a variety of accretion states characterized by different broadband X-ray spectral shape and timing properties (see Remillard & McClintock 2006; Dunn et al. 2010; Belloni et al. 2016, for reviews). In particular, we distinguish between the hard and soft spectral states (HS and SS, respectively). The HS is characterized by a cutoff (around 100 keV when detected) power-law-like spectrum extending up to high energies, and a cool (kT_disc ∼ 0.1 − 0.3 keV) multi-color disk blackbody. In the SS the spectra are dominated by a hotter multi-color disk blackbody component (kT_disc ∼ 0.8 − 1.0 keV), and sometimes an additional steeper hard tail is detected as well. In all these states, a further contribution to the overall spectrum is given by the reflection component, a component thought to be due to the hot photons emitted by the corona reprocessing on the accretion disk (see below). The study of the reflection component, and in particular of the relativistically broadened Fe K line at 6.4 keV, can be used as a powerful diagnostics tool to infer the inner radius and ionization of the accretion disk and the inclination of the system, for example (see, e.g., Reynolds & Nowak 2003, for a review).

The variety of spectral states in BHTs has traditionally been interpreted as being due to different physical properties and geometry of the accretion flow around the BH. In the HS it is thought that the emission is due to thermal Comptonization by a hot optically thin electron plasma (corona). The source of soft seed photons is the (possibly truncated) optically thick accretion disk (Shakura & Sunyaev 1973). Conversely, in the SS such a disk is expected to extend down to the innermost stable circular orbit (ISCO, Bardeen et al. 1972) and only a marginal contribution from the corona is observed. However, the consensus on the truncated disk model is not global, and a number of counterarguments have been given as well (see, e.g., Zdziarski & De Marco 2020, for an extensive discussion on the topic). Among the critical issues, spectral modeling of seemingly broad Fe K lines in BHTs indicated that the optically thick disk was close or even coincident with the ISCO well within the hard state (examples of these can be found in Miller et al. 2006; Tomsick et al. 2008; García et al. 2015).

The controversies on the accretion flow geometry across spectral states is one of the open problems regarding BHTs. When transients rise from quiescence, their spectral evolution usually follows a standard path: the system moves from quiescence to the hard state, increasing in luminosity by several orders of magnitude, then evolves through the hard-intermediate (HIMS) and soft-intermediate (SIMS) up to the soft state. Thereafter, the luminosity decreases and the source exhibits the same transition in reverse but at a much lower luminosity. The common behavior of BHTs in making transitions from hard to soft and from soft to hard states at different luminosities is called hysteresis (Miyamoto et al. 1995; Zdziarski et al. 2004). The origin of this pattern is a matter of strong debate. It is also noteworthy that a relevant fraction of BH binaries (about 40%, Tetarenko et al. 2016a) are “hard-only” or have displayed “failed transition” outbursts (see, e.g., Hynes et al. 2000; Brocksopp et al. 2004; Capitanio et al. 2009; Ferrigno et al. 2012; Del Santo et al. 2016; Bassi et al. 2019; de Haas et al. 2021). BHTs are known to launch mildly relativistic jets, which accounts for the emission from these systems over radio and mid-IR frequencies. Jets are observed only in the HS and HIMS, while they are quenched around the transition to the SS (Fender et al. 1999; Corbel et al. 2000; Fender & Belloni 2004; Russell et al. 2020) and a short-lived transient jet can be launched (e.g., Fender & Belloni 2004; Fender et al. 2010; Rushton et al. 2017; Russell et al. 2019). A third critical ingredient is represented by how the X-ray timing properties (i.e., X-rays periodic and aperiodic variability) change according to the accretion state of BHTs. In this context one of the most intriguing challenges consists in the interpretation of quasi-periodic oscillations (QPOs), the peaks observed in the X-ray Fourier power density spectra of BHTs, which are correlated with the spectral evolution of the systems (see Ingram & Motta 2019; Motta et al. 2021, for recent reviews).

Much effort has been dedicated to interpreting the multi-wavelength spectral behavior and the timing properties of BH binaries in outburst. However, taking into account the correlation between accretion and ejection has proven to be difficult. A first attempt to describe globally such complex behavior was made by Esin et al. (1997), who proposed a multi-flow configuration for the disk. According to these authors, the accretion flow around BHs would consist of an outer Shakura-Sunyaev accretion disk and an inner, less dense, and radiatively inefficient plasma phase, called advection dominated accretion flow (ADAF, Narayan & Yi 1995; Yuan 2001). However, the model suffered several weaknesses; in particular, it was not able to reproduce hard states at luminosities comparable to the Eddington limit (see the Introduction of Marcel et al. 2018b, for an extensive discussion on the topic). In the following decade many updates were made to the model (see, e.g. Yuan 2001; Meier 2005; Xie & Yuan 2012) to make it more suitable to describe the increasing number of X-ray observations of BHTs. However, the correlation between accretion flow and jets has not been addressed in any of the updated versions of the model.

An attempt to connect accretion and ejection in a unified model has been recently reported in a number of papers (Ferreira et al. 2006; Marcel et al. 2018a,b, 2019, 2020). This model connects the spectral evolution BHTs in outbursts to changes in the multi-flow configuration of the accretion disk, similarly to the approach of Esin et al. (1997), among others. The main novelty of this paradigm consists in its ability to also explain the appearance and disappearance of the jet and their correlation with the accretion flow. According to this model, the accretion disk is threaded with a large-scale vertical magnetic field B_Z(R). Recent numerical simulations have shown that such magnetic fields become stronger near the compact object (e.g., Liska et al. 2018). As a consequence, the accretion flow is also expected to become more magnetized as we approach the inner edge of the disk. In the following we define the magnetization as , with P_tot the total pressure (radiation plus gas). In the outer regions of the disk μ ≪ 1, and particles are barely affected by the presence of the vertical magnetic field. As a consequence, jet launching is inactive in this region of the disk. In order to launch magneto-centrifugally driven jets (Blandford & Payne 1982), it has been shown that μ≳0.1 must be achieved (Ferreira & Pelletier 1995; Ferreira 1997; Jacquemin-Ide et al. 2019). For μ < 10⁻³ or less, magnetic winds can be launched, possibly carrying away a significant fraction of the mass. However, these winds exert a negligible torque on the underlying disk (Zhu & Stone 2018; Jacquemin-Ide et al. 2021). At these distances from the BH, particles are only subject to the torque due to the internal turbulent viscosity, so that the disk can be nicely described with the classic Shakura & Sunyaev (1973) model. The inner regions are instead highly magnetized, with μ beyond the 0.1 threshold, and powerful jets can be launched. Jets carry away mass, energy, and angular momentum, and thus exert an additional torque (see, e.g., Ferreira & Pelletier 1993, 1995). Subsequently, accretion proceeds here at supersonic velocity, much higher than in the Shakura-Sunyaev accretion disk, and plasma in the internal regions of the accretion flow result more rarefied. In summary, the accretion flow surrounding the BH has a two-flow configuration. It is composed of a Shakura-Sunyaev disk in the outer regions of the flow, defined as a standard accretion disk (SAD), and an inner, less dense, and optically thin jet emitting disk (JED) (Ferreira et al. 2006). The JED shares the same physical properties of the hot corona, but additionally it also drives jets. Quite intuitively, each accretion state could be obtained by mixing these ingredients with different quantities, specifically with a hybrid JED-SAD configuration where the two realms extend over regions of different scale. In this sense, such a configuration is fundamentally determined by two main control parameters: the transition radius R_J between the JED and the SAD, and the inner mass accretion flow Ṁ_in feeding the BH. A detailed description of the model has been presented in Marcel et al. (2018a,b) and we refer the reader to these papers for more details. This spectral model has been successfully used to reproduce the X-ray spectral evolution of the BHT GX 339−4 during four outbursts between 2001 and 2011 (Marcel et al. 2019, 2020). Moreover, using the R_J − Ṁ_in pairs that best described the X-ray spectra in a HS, these authors were also able to reproduce the jet emission observed simultaneously in the radio band. Finally, it has been recently shown that some timing features could be explained within the JED-SAD paradigm framework. A direct proportionality has been observed between the Keplerian frequency of the transition radius R_J and the type C QPO frequency in four different outbursts of GX 339−4 (Marcel et al. 2020). According to these results, this QPO type could originate at the interface between two regions of different values of μ, being then strictly related to the existence of two different types of accretion flow in the HS and HIMS of BHTs (Ferreira et al. 2021).

The JED-SAD paradigm has successfully explained much of the observed accretion and ejection phenomenology in GX 339−4. However, the variety of different behaviors observed in three decades of BHT studies (see, e.g., Dunn et al. 2010; Tetarenko et al. 2016b, for observational reviews) demands other tests of this model. In this paper we report on the application of the JED-SAD model to the BHT MAXI J1820+070 (ASASSN-18ey), hereafter MAXI J1820, exploiting a broad X-ray data set including data from NuSTAR, the Neil Gehrels Swift Observatory (hereafter Swift) and the Neutron Star Interior Composition Explorer (NICER). This is the first time the JED-SAD model is applied directly to the data through spectral fits and that Compton reflection process is taken into account. It therefore represents a major test for the potentialities of such paradigm.

2. MAXI J1820+070

MAXI J1820 is a BHT that was observed for the first time in the optical band by the All-Sky Automated Search for SuperNovae ASSAS-SN (Shappee et al. 2014) on 2018 March 3 (Tucker et al. 2018), and one week later by MAXI in X-rays (Kawamuro et al. 2018). Detailed studies of the optical counterpart revealed that the system hosts a stellar-mass black hole (∼8.5 M_⊙) accreting from a ∼0.4 M_⊙ companion star (Torres et al. 2019, 2020). Since its discovery, the source has undergone a long (approximately one year) and bright outburst, becoming at its peak the second brightest object in the X-ray sky. Due to its brightness, the system was the object of an impressive multi-wavelength observing campaign (see, e.g., Shidatsu et al. 2018; Paice et al. 2019; Hoang et al. 2019; Trushkin et al. 2018; Bright et al. 2020; Tetarenko et al. 2021) and of a large number of studies. The most recent and precise measure of the distance, determined via radio parallax, is 3.0 ± 0.3 kpc (Atri et al. 2020). Furthermore, the system shows X-ray dips (Kajava et al. 2019) but not eclipses, suggesting an inclination between 60° and 80°. Other evidence of the high inclination of the system are provided by optical spectroscopy (Torres et al. 2019) and by the estimate of the inclination of the jet axis (Atri et al. 2020; Wood et al. 2021), which is about 63°. The orbital period of the system has also been estimated to be around 0.68 days (Patterson et al. 2018; Torres et al. 2020). In X-rays, the outburst was studied in detail in HS (see, e.g., Bharali et al. 2019; Buisson et al. 2019; Zdziarski et al. 2021), SS (see, e.g., Fabian et al. 2020) and in its final phase (see, e.g., Xu et al. 2020). The truncation of the disk during the HS is one of the most controversial aspects of the system. On the one hand, a number of spectral-timing works have pointed out that the disk reaches the ISCO already in HS, while a contracting lamppost corona is responsible for the hard X-ray emission (Kara et al. 2019; Buisson et al. 2019; You et al. 2021; Wang et al. 2021). On the other hand, a truncated disk scenario has also been proposed on the basis of further spectral-timing analyses, even on the same sets of data (Zdziarski et al. 2021; De Marco et al. 2021; Axelsson & Veledina 2021).

3. Observations and data reduction

In this paper we focus on the spectral analysis of the source in the high luminosity hard state, i.e., from ≈ MJD 58190 until ≈ MJD 58300. With the aim of modeling broadband X-ray spectra, we included data collected by Swift/XRT, NuSTAR, NICER and the Burst Alert Telescope (BAT) onboard Swift. In order to obtain broad-band spectra from observations close in time, we selected only XRT, NuSTAR and NICER observations which were quasi-simultaneous, i.e., taken not more than one day apart from each other. Furthermore, in order to avoid any discrepancy between BAT and NuSTAR due to the spectral variability of the source, we extracted BAT spectra exactly over the duration of each NuSTAR observation. These criteria led us to narrow down the list of the available NuSTAR data to eight epochs. We note that this set of NuSTAR observations has been already analyzed and described by Buisson et al. (2019). In the following, we will refer to these observations as epochs and we will label them with numbers from 1 to 8 in chronological order. It is worth noticing that NICER observations cover five out of eight epochs and that the BAT spectrum in epoch 6 has been extracted in a longer time interval, due to the low number of visits performed by Swift during the NuSTAR observation. Details on the selected epochs are reported in Table 1.

3.1. XRT

The 2018 outburst of the source was monitored regularly by the XRT telescope on board Swift from MJD 58191 (2018, March 14) until MJD 58428 (2018, November 6), with a total of 75 observations performed in Window Timing (WT) mode. The XRT data were first reprocessed with the task XRTPIPELINE, included in the software package HEASOFT (v. 6.26.1). The source extraction procedure from the cleaned event files was performed with DS9. Since all of the observations had a high count-rate, i.e., always well above the limit for the pile-up correction in WT (100 ct/s), we used an annulus region centered on the source coordinates to extract spectra not affected by pile-up. The outer radius was always chosen as ∼47″, while the inner radius of the annulus was selected based on the registered count-rate in accordance with the Swift/XRT guidelines¹. In particular, we used a ∼18″ inner radius for epoch 1, a ∼24″ inner radius for epochs 2, 3, 5, 6 and a ∼28″ inner radius for epochs 4, 7, and 8. Each spectrum was re-binned with GRPPHA in order to have 150 counts per bin, allowing the use of the χ² statistics.

3.2. NuSTAR

Data were reduced using the standard NUSTARDAS task, incorporated in HEASOFT (v. 6.26.1). We extracted high scientific products (light curves and spectra) using a circular area of 100″ radius centered at RA = 18:20:22.0702, and Dec = +7:10:58.331, as source region. In order to take into account any background non-uniformity on the detector, we extracted the background spectra using four circles of ∼50″ radii in different positions with negligible contamination from the source. We then used NUPRODUCTS to build spectra and light curves. We used data from both the two hard X-ray imaging telescopes on board NuSTAR: the focal plane mirror (FPM) A and B. The extracted spectra were grouped using the optimal binning recipe by Kaastra & Bleeker (2016) in order to have a grouping that reflects the spectral resolution of the instrument in a given energy range and avoid any oversampling issues. We did not sum the two spectra, but rather fitted them together by leaving a floating cross-normalization constant as suggested by the NuSTAR team for bright sources².

3.3. BAT

Data from the BAT survey were also retrieved from the HEASARC public archive. The downloaded data were processed using BAT-IMAGER software (Segreto et al. 2010). This code, dedicated to the processing of coded mask instrument data, computes all-sky maps and, for each detected source, produces light curves and spectra. Light curves in Crab units were extracted in three energy bands, 15–40 keV, 40–80 keV, 80–150 keV, with one-day binning time (Fig. 1). Spectra were extracted in the range 15–195 keV, with logarithmic binning (for a total of 49 bins) and the official BAT spectral redistribution matrix was used.

Fig. 1. From top to bottom, the Swift/XRT light curves and the related soft X-ray HR, Swift/BAT light curves and the hard X-ray HR are shown. Energy ranges are reported in the see legends. The NuSTAR observations are highlighted with black dashed lines.

3.4. NICER

NICER covered the 2018 outburst of MAXI J1820 with almost daily cadence, resulting in hundreds of individual observations (ObsIDs 1200120101-1200120312). In this work we analyzed a sample of the NICER observations which satisfied the requirement of being quasi-simultaneous (i.e., taken within a day) with respect to the NuSTAR data. The data were reduced using NICER-DAS 2019-05-21 v006. We selected good time intervals using NIMAKETIME and then applied them to the data via NIEXTRACT-EVENTS, selecting events with PI channel between 25 and 1200 (0.25–12.0 keV). In order to fix the distortions due to the NICER calibration uncertainties, we re-normalized the spectra by using the residuals of a power-law fit to the Crab nebula (Ludlam et al. 2018). We used the public files nixtiaveonaxis20170601v002.arf and nixtiref20170601v001.rmf as ancillary response file and redistribution matrix file respectively, retrievable from the NICER website³. As a background spectrum we used the public background file nixtiback20190807.pi, also available in the HEASARC archive.

4. Spectral analysis

The XRT and BAT light curves and the related hardness ratios are shown in Fig. 1, while in Fig. 2 the XRT hardness intensity diagram (HID) of the whole 2018 outburst is shown. We used the (absorbed) 0.5–10 keV flux as an indicator of the intensity of the source. In order to estimate the flux, we fitted each XRT spectrum separately in the energy bands 0.5–2 keV and 2–10 keV with a simple power-law model (POWER in XSPEC). We note that a power-law model, although not an appropriate choice to describe adequately the spectrum in the whole XRT energy range, approximates satisfactorily its shape in the two energy bands considered. We associated an error equal to 10% of the estimated value to each of these fluxes.

Fig. 2. XRT HID of the source during the entire outburst. The epochs used in this work are highlighted by the colors and different symbols (Sect. 4.2), i.e., circle for phase 1 (rise), upper triangles for phase 2 (plateau) and lower triangles for phase 3 (decline).

According to their hardness ratio and their timing properties (e.g., De Marco et al. 2021), during epoch 1 the system was in the canonical HS, while it remained in HIMS from epochs 2 to epoch 8. The latter occurred during a short episode of re-hardening, as better highlighted in the BAT light curve (see Fig. 1). In the following we report on the spectral analysis of these epochs.

The NuSTAR data usually range from 3–79 keV. However, a remarkable (and unexpected) mismatch between FPMA and FPMB data is sometimes observed below 4 keV, likely due to a known instrumental issue (Madsen et al. 2020). This mismatch was visible in our data. Therefore, we decided not to include data below 4 keV, which are already covered by XRT. Similarly, the NICER data below 4 keV were also ignored as they showed a soft excess, which could be related to a further Comptonization or a quasi-thermal component (see Fig. 6, Zdziarski et al. 2021), and several lines, both sharp instrumental features (Wang et al. 2020) and broad physical lines. A proper investigation of these features goes beyond the purposes of this manuscript and will be reported elsewhere. Moreover, the BAT data showed an unexpected mismatch with NuSTAR below 30 keV, possibly due to high systematics in this energy range for extremely bright sources. We therefore ignored BAT data below 30 keV so that the final broadband spectra covered the range 0.8–10 keV for XRT, 4–78 keV for NuSTAR, 4–10 keV for NICER, and 30–190 keV for BAT.

Each spectrum was analyzed using the TBABS model in order to take into account the effect of interstellar absorption, with the photoelectric cross sections from Verner et al. (1996) and the element abundances from Wilms et al. (2000). A CONSTANT component was also used to serve as cross-calibration constant. To account for differences in the calibration of NuSTAR and NICER (see also Zdziarski et al. 2021) we included a cross-calibration function which reads CONST × E^ΔΓ, with ΔΓ the discrepancy in photon index between the data sets (Ingram et al. 2017; Ursini et al. 2020). The value of ΔΓ was left free to vary for NICER, while it was fixed to 0 for all the other instruments included in the broadband spectrum. Finally, as each data set was showing apparent systematics, we applied a 1% systematic error to the analyzed spectra.

4.1. Basic ingredients of the JED-SAD spectral model

The eight X-ray observations were analyzed with the JED-SAD model by applying it to the data in the form of a XSPEC table. The following key parameters are used in the model (but see Marcel et al. 2018a, for a full description):

• transition radius R_J between the JED and the SAD, in units of gravitational radii R_G = GM/c²;

• mass accretion rate at the transition radius R_J or inner mass accretion rate Ṁ_in, here expressed in units of Eddington mass accretion rates Ṁ_Edd = L_Edd/c²;

• sonic Mach number m_S of the accretion flow, defined as the ratio of the mass weighted accretion speed to the sound speed;

• ejection efficiency p⁴, which takes into account how effectively the jet extracts angular momentum from the underlying JED, then modifying its structure. Therefore in the JED, the accretion rate Ṁ(R) scales with R according to the formula , with p the ejection efficiency;

• magnetization μ (see Introduction);

• fraction b of the accretion power P_acc that is released in the JED and ultimately powers the jets (i.e., b = 2P_jet/P_acc);

• geometrical dilution factor ω, corresponding to the fraction of the cold photons emitted by the SAD which cool down (through inverse Compton) the JED (we refer to Marcel et al. 2018a, and references therein for further discussion on this topic). In the current version of the model, ω assumes values in the range 0–0.5;

• radius of the ISCO (R_ISCO), which coincides with the innermost radius of the JED⁵;

• normalization factor K_JEDSAD, which determines how the luminosity of the system scales with the distance D_kpc (in units of 1 kpc), following the formula

  (1)

As explored in detail by Marcel et al. (2018b), the JED emission in the archetypal object GX 339–4 is best described by setting μ between 0.1 and 1.0, b = 0.3, ω = 0.2, and p = 0.01⁶. In the following we fix μ to 0.5 and adopt the same set of values for b, ω, and p in the case of MAXI J1820+070. Contrarily to the works on GX 339–4, we leave m_S free to vary in order to investigate its correlation with the high energy cutoff, which can be better constrained in our work due to the NuSTAR and BAT coverage at hard X-rays. The value of R_ISCO in MAXI J1820 is not precisely known. While it has been reported that the ISCO in MAXI J1820+070 could be even lower than 2 R_G (e.g., Kara et al. 2019), it was recently found that the BH in the system is likely to be slowly spinning (Fabian et al. 2020; Zhao et al. 2021; Guan et al. 2021), with a_* not higher than 0.5. We chose therefore to use an initial value of R_ISCO = 4 R_G (corresponding to a_* = 0.55), approximately consistent with the estimates reported in these works. The distance of the source was fixed at 3 kpc according to Atri et al. (2020). We will allow more freedom in R_ISCO and K_JEDSAD and discuss the dependence of our results on the choice of these parameters in Sects. 4.4 and 4.5. A summary of the parameters in the model and the values adopted in this work is presented in Table 2.

4.2. The two reflection models

A first analysis without any reflection component results in extremely poor fits and strong residuals in the iron line and Compton hump regions. This was expected because reflection is not taken into account in the JED-SAD model. We therefore used a reflection table component, which we label REFL for simplicity, based on the XILLVER reflection model (García et al. 2013). This table was produced by simulating spectra for different combinations of R_J − Ṁ_in pairs, fitting them with a simple cutoff power-law model, and finally injecting the corresponding spectral index and high energy cutoff in the XILLVER table. More details on this model are given in Barnier et al. (2021). The parameters of this component are the iron abundance A_Fe, the disk ionization parameter logξ, and a normalization K. Since REFL does not take into account the relativistic blurring effects, we applied the convolution kernel KDBLUR, which smears the reflection spectrum according to the original calculations by Laor (1991). This model considers four additional parameters: the emissivity ϵ, the inner and outer radius of the reflecting disk R_in and R_out, and the inclination i. In order to reduce the number of degrees of freedom we set ϵ to the reasonable value of 3 (see, e.g., Dauser et al. 2013; Xu et al. 2020) and i to 70° (consistent with the estimated source inclination; see, e.g., Kajava et al. 2019). As the reflection component is produced by the disk (SAD) surface reprocessing the hard X-ray photons emitted by the corona (JED), within the JED-SAD framework the inner radius of the reflection R_in should coincide with the transition radius R_J between JED and SAD. In all spectral fits, R_J and R_in were therefore linked together.

The first model we applied, Model 1, reads:

(2)

It is noteworthy that two reflection components were already suggested by Kara et al. (2019) on the basis of a reverberation time lag study with NICER data. Thereafter a number of authors used a double-reflection model in the spectral analysis of the same NuSTAR observations presented here (Buisson et al. 2019; Chakraborty et al. 2020; Zdziarski et al. 2021). They all found a low ionization reflection component and a highly ionized inner reflection component, which both contribute to the observed iron line profile providing respectively a narrow core and a broad base. Therefore, we also tested a second model, Model 2, where a second REFL component has been added:

(3)

When Model 2 is applied, we use the subscripts “1” and “2” to refer to the parameters of the inner and outer reflection components, respectively. The purpose of this model is to reproduce a physical scenario where reflection comes from two distinct regions, characterized by different ionization parameters and therefore presumably placed at different distances from the illuminating corona (the JED here). Therefore, we linked R_in, 1 to R_J and R_out, 1 to R_in, 2 in order to place the two reflection regions at separate but neighboring regions⁷, as shown in Fig. 3. The outer radius of the remote reflection is instead left free to vary. The normalization and ionization parameters were left free as well in both reflection components, so that we have two different normalization and ionization values: K₁, K₂, ξ₁, and ξ₂. All the other parameters (the inclination, the iron abundance, and the emissivity) were linked between the two reflection components.

Fig. 3. Sketch representing the proposed two-reflection geometry underlying Model 2. From left to right: the black hole (in black), the inner JED (violet) and the jets (lilac), the outer SAD (from red to yellow).

4.3. Spectral fits with R_ISCO = 4 R_G

Based on the fit results described below, we can divide our observations in three groups or phases. Phase 1 corresponds to epoch 1; phase 2 includes observations from epoch 2 to epoch 6; and phase 3, which corresponds only to epochs 7 and 8. The grouping of these epochs is consistent with the XRT HID (Fig. 2), where an analogous clustering of observations in three regions can be observed. A similar subdivision was also proposed by De Marco et al. (2021), with phases labeled “Rise,” “Plateau,” and “Bright decline,” respectively. In the following, we treat these phases separately. In all phases, we did not manage to constrain the values of R_in, 2 and R_out, 2, as the fit was basically insensitive to variations of these parameters. We therefore set R_in, 2 to 300 R_G and R_out, 2 to 10⁴ R_G. We found values of A_Fe always between 2 and 3 times the solar abundance⁸.

The final results are reported on Table 3, while the corresponding best-fit models and residuals are shown in Figs. 4–5.

Fig. 4. Energy spectra with the best-fit models: Model 1 for epoch 1 and Model 2 for epochs 2–4, and residuals. Data: XRT (blue), NICER (green), NuSTAR (yellow–black), and BAT (red). Different line styles were adopted to distinguish between the different components: dot for JED, dash-dot-dot-dot for SAD, dash-dot for inner REFL, dash for the outer REFL.

Fig. 5. Energy spectra, best-fit model (Model 2) and residuals for epochs 5–8. Data: XRT (blue), NICER (green), NuSTAR (yellow–black), and BAT (red). Different line styles were adopted to distinguish between the different components: dot for JED, dash-dot-dot-dot for SAD, dash-dot for inner REFL, dash for the outer REFL.

Phase 1: The system in the hard state with one reflection component. Epoch 1 is satisfactorily fitted with both Model 1 and Model 2, with (d.o.f.) of 1.01(749) and 0.90(748), respectively. In both cases the SAD is truncated far away from the BH (i.e., with R_J of ∼60 R_G and ∼150 R_G for Model 1 and Model 2, respectively). However, with Model 2 the value obtained for logξ₁ goes beyond the threshold of 4 (approximately 4.3). We consider such a high value as unphysical for several reasons. First of all, such a highly ionized medium at the edge of the disk is hard to reconcile with it being truncated at 150 R_G. Second, we note that with the very high value found for logξ₁, the inner reflection component has almost no significant fluorescence line and basically serves as an additional continuum component. Finally, the spectral shape of the reflection component computed with XILLVER for logξ > 4 is unreliable at high energies⁹. According to these arguments, the presence of the second reflection component is most likely spurious for epoch 1. Using Model 1, logξ¹⁰ is found in the range 2.0–2.4, which is instead more plausible.

Phase 2: The outer reflection component arises. Epochs from 2 to 6 are poorly fitted by Model 1, with values all above 1.25 and, more importantly, evident unmodeled structures between 6 and 10 keV in the residuals. Instead, the application of Model 2 results in good fits and quite flat residuals. The picture captured by the fits is similar for the epochs in this group. The best-fit value of R_J decreases significantly with respect to the previous phase, attaining values in the range 35–50 R_G. Simultaneously, Ṁ_in evolves substantially, going from 1.8 Ṁ_Edd (epoch 2) to 2.4 Ṁ_Edd (epoch 3), and then decreases again to 1.3 in epoch 6. Throughout this phase, m_S retains a well-determined value of about 1.25 for all the epochs, with the only exception of epoch 6, where the parameter was basically unconstrained and therefore fixed. In all epochs, ξ₁ is higher than ξ₂ by at least an order of magnitude. Interestingly, while ξ₁ is quite stable, an increasing trend can be observed in logξ₂, which goes from an upper limit of 1.7 in epoch 2, to 3 in epoch 6. The possibility of this trend being related to an evolution of the inner radius of the outer reflection will be explored in the Discussion (see Sect. 5).

Phase 3: The role of m_S in the re-hardening. Model 1 also gives an unacceptable fit for epoch 7, but is instead compatible with epoch 8, with a (d.o.f.) of 1.15 (750). Model 2 gives a good fit to epoch 7 and an even better fit to epoch 8, with a probability of improvement by chance of ∼10⁻⁵ with respect to Model 1. We then conclude that a double reflection scenario is also statistically preferred in phase 3. As in the previous phase, R_J decreases again to about 33–34 R_G, while Ṁ_in is constrained at ∼1 Ṁ_Edd. The most apparent evolution with respect to the previous phase consists in a re-hardening of the spectrum, as witnessed by the high energy cutoff of the BAT data moving toward higher energies (see the BAT light curve, Fig. 1). Since some level of degeneracy between R_J and m_S clearly exists (see Fig. 6), this trend can be explained by an increase in both these parameters. For example, epoch 7 can be fitted keeping m_S fixed to 1.3, in consistency with the previous phase, and with R_J of 38–39 R_G or leaving m_S free, a case where m_S attains a value of 1.5 and R_J decreases to 33–34 R_G. Both fits are acceptable, but keeping m_S fixed results in slightly unmodeled residuals in the iron line region, contrarily to the case where m_S is free. We therefore conclude that the spectral hardening in this phase is not reflected in the values of R_J, which actually suggests that the inner edge of the disk keeps moving inward, as expected during the transition to the SS, but seems to be instead produced by an increase in m_S, producing a more rarefied and hotter accretion flow in the JED. Another interesting evolution can be found in the values of ξ₁ and ξ₂. While epoch 7 has values compatible with epoch 6 (and phase 2 in general), in epoch 8 both values logξ₁ and logξ₂ decrease significantly to values of ∼3.5 and 2.0 (which was fixed since it was left unconstrained by the fit), respectively.

Fig. 6. JED solutions showing the impact of mS and RJ on the high energy cutoff. Left panel: different JED solutions obtained keeping RJ to 40 RG for different values of mS. Right panel: different JED solutions are produced for different values of RJ, here expressed as rJ = RJ/RG, and mS fixed to 1.5. The plot is used to display how increasing both mS and increasing RJ pushes the hard X-rays cutoff to higher energies.

A summary of the evolution of the main parameters at play is presented in Fig. 7. We also show the evolution of the reflection components in Fig. 8.

Fig. 7. Evolution of the best-fit parameters for RJ, Ṁin, and mS over the eight analyzed epochs. In all panels the triangles and bold lines are used for fits performed with RISCO at 4 RG. The reported errors are likely underestimated and an additional systematic error of about 10%–20% has to be considered to account for the uncertainty on ω and KJEDSAD (see Sect. 4.5).

Fig. 8. Best-fit models for inner and outer reflection in epochs 2–7. Inner reflection models have been multiplied by a factor of 10 for clarity.

Using the best-fit values found for R_J, Ṁ_in, and m_S, it is possible to produce diagrams for the main physical properties of the corresponding hybrid JED-SAD configuration. These diagrams are presented in Fig. 9 for epochs 1, 5, and 8, each chosen to represent one of the three phases. In these plots the accretion flow has been divided into 31 radial zones, 30 for the JED, one for the SAD. The scale height H, the Thomson optical depth τ_T, the electron temperature kT_e, and the emitted spectrum were calculated and plotted in Fig. 9 for each of these 31 regions (see Appendix C in Ursini et al. 2020, for further details on how these plots were produced). The plots show that the increase in Ṁ_in in phase 2 corresponds to a geometrically thinner, optically thicker, and colder JED with respect to phase 1. The evolution of the accretion flow is also apparent from the softening of the spectra emitted by the single rings forming the JED. The re-hardening episode in phase 3 is witnessed by the slight increase in the spectral high energy cutoff combined with an increase in both kT_e and τ_T.

Fig. 9. JED-SAD solutions for epochs 1, 5, and 8, used to represent phases 1, 2, and 3, respectively, calculated using the best-fit parameters for RJ, Ṁin, and mS reported in Table 3. The following values were adopted: for phase 1 RJ = 57 RG, Ṁin = 0.8 ṀEdd, and mS = 1.3; for phase 2 RJ = 44 RG, Ṁin = 1.7 ṀEdd, and mS = 1.3; and for phase 3 RJ = 33 RG, Ṁin = 1.1 ṀEdd, and mS = 1.5. In all panels, the radial distribution is divided in 31 portions (annuli), 30 for the JED and one for the SAD. Top and middle panels: how the actual shape z = H/R of the disk evolves with R. The distribution of the Thomson optical depth (τT, top) and the electron temperature (kTe, middle) are color-coded (see color bars at right). In the same plots, the event horizon of the BH is represented by a half black circle. The total emitted spectrum is shown in the bottom panels as a solid black line. The contribution from each of the 31 considered annuli is highlighted with a dotted line for the SAD region and dashed lines for the JED regions. The sections of the plot corresponding to the energy ranges that are not covered by XRT, NICER, NuSTAR, and BAT data are shown with a gray background.

4.4. Spectral fits with R_ISCO = 2 R_G

In Sect. 4.3 we kept the innermost radius R_ISCO frozen to 4 R_G, corresponding to a moderately spinning black hole, as suggested by recent spectral analysis (Fabian et al. 2020; Zhao et al. 2021; Guan et al. 2021). In this section we replicated the previous fitting procedure after setting R_ISCO to 2 R_G, corresponding to a rapidly spinning BH (a_* ∼ 0.95).

In this new ensemble of fits, a subdivision in three phases can again be identified. Furthermore, a double-reflection model is still necessary to fit the epochs, with the exception of epoch 1, for which Model 2 results again in oddly high values of logξ₁. However, a lower R_ISCO significantly affects the values obtained for R_J, Ṁ_in, and m_S. In this case as well, R_J decreases through the three phases, but starts from almost 44 R_G to then decrease to ∼20–30 R_G and finally to ∼15–20 R_G in phase 2 and 3, respectively. The trend is significant, since the discrepancy between these best-fit values is larger than their errors. A different trend is observed in Ṁ_in, which remained very stable throughout the eight epochs analyzed, with only a slight increase in phase 3. The stability in Ṁ_in is compensated by m_S evolving blatantly, going from m_S ∼ 1.1 in phase 1 to m_S ∼ 0.75 and finally rising again in phase 3 up to 1.5. If we impose the same m_S for all the epochs, in conformity with the analysis performed on GX 339–4 (Marcel et al. 2018b), the best fits found by XSPEC are strongly reflection-dominated. This happens because for m_S values beyond 1, the JED model is characterized by a rather high electron temperature. The subsequent high energy cutoff is not compatible with the spectral curvature well traced by the BAT data in our spectra, which indicate a cutoff at about 100 keV. The only way for the JED-SAD model to reproduce such spectral curvature would be to decrease R_J and simultaneously increase the temperature of the inner edge of the SAD. However, a hotter SAD should also be visible in the XRT band, which does not show any sign of such a strong disk component. Since the JED-SAD is not suitable, the best solution to describe the broadband spectrum consists in tuning up the inner (highly ionized) reflection component. Such reflection-dominated spectra are challenging to explain, especially in the framework of the proposed JED-SAD geometry. Similarly high reflection fractions have been explained in the past, for example as a large covering fraction of the disk due to the presence of clouds of colder matter (Malzac & Celotti 2002), as a compact corona close to the BH emitting anisotropically due to light bending (Miniutti et al. 2003), or as being due to a mildly relativistic motion of the corona toward the disk (Beloborodov 1999). Unfortunately all of these scenarios are hard to reconcile with the proposed JED-SAD geometry. We conclude that with R_ISCO = 2 R_G the role played by Ṁ_in in driving the evolution of the accretion flow toward intermediate states is somehow replaced by an evolution in m_S.

4.5. Exploring the effect of the other JED-SAD parameters

In order to reduce the number of degrees of freedom, so far we kept some parameters fixed in both Models 1 and 2. In the following we explore in more detail their impact on the results obtained.

As shown by Marcel et al. (2018b), Fig. 10, p (which was labeled ξ by these authors) has little to no effect on the JED-SAD spectral shape. Similarly, the choice of fixing μ to 0.5 should have a negligible impact, especially since the high energy cutoff is basically insensitive to this parameter (see Fig. 6, Marcel et al. 2018b). On the contrary, the jet power b is expected to influence the JED shape at both hard and soft X-rays. Broadly speaking, by increasing b more power is funneled in the jets, making the underlying JED rather cold, while with a lower b, the JED results hotter. This effect would shift the position of the high energy cutoff and therefore act as m_S. However, b is a function of the aspect ratio ϵ, which for a typical JED should be about 0.2–0.3. As shown by Petrucci et al. (2010), Appendix A, we need b values around 0.3–0.5 to obtain such aspect ratios, as higher values would produce disks that are geometrically too thin to be consistent with the expected JED geometry. We therefore do not explore the effect of varying b further.

Potential biases could be ascribed to the dilution factor ω and to the parameter K_JEDSAD, which is solely related to the distance. We recall ω being a parameter regulating the amount of SAD photons penetrating the JED and thereby cooling it. In the current model, it ranges between 0 (the unlikely case where no photons from the SAD are intercepted by the JED) and 0.3 (when up to 30% of the photons can be intercepted). In order to check the effect of ω, we analyzed epoch 5 again, but keeping ω fixed first to 0 and then to 0.3. As expected, when we put ω to zero, the cooling effect played by the SAD is replaced by lowering R_J, which attains a value of 27–31 R_G, while also Ṁ_in and m_S drop to values of about 1.1 Ṁ_Edd and 1.0, respectively. Increasing ω to 0.3 only moderately impacts the outcomes of the fits (i.e., with R_J and m_S compatible with the results obtained with ω = 0.2 and a 15% increase in Ṁ_in). The distance of MAXI J1820 is well constrained between 2.7 and 3.3 kpc (Atri et al. 2020). In the previous sections we froze ω to 0.2 and K_JEDSAD to 11.11, corresponding to a distance of 3 kpc. We then computed again the best fits for epoch 5 freezing K_JEDSAD first to 13.7 (2.7 kpc), and then to 9.2 (3.3 kpc). Analogously to ω, variations in K_JEDSAD do not affect the overall physical scenario which was identified in this phase phase with K_JEDSAD = 11.11, but they result in a slight fluctuation of the best-fit values of Ṁ_in and R_J. For example, R_J goes from ∼41 R_G at d = 2.7 kpc to ∼48 R_G at d = 3.3 kpc, a range comparable with the estimated errors (about 2–5 R_G). On the other hand, Ṁ_in oscillates between ∼1.58 and ∼1.70 Ṁ_Edd (for the lower and higher edge of the distance range, respectively), which is higher than the range of variability associated with the errors. These results point out that some caution must be taken when considering the best-fit values obtained for R_J and Ṁ_in as face values.

5. Discussion

In the previous section we applied the physical model JED-SAD to fit the XRT+NICER+NuSTAR+BAT observations of the source MAXI J1820+070 in HS and HIMS. We included either one (Model 1) or two (Model 2) reflection components. In the latest model, which is more complex, the two reflection spectra originate from different but neighboring regions of the SAD disk. In all the epochs analyzed, with the exception of epoch 1 (the only one in full hard state), we found that only a two-reflection solution gives acceptable fits. In the following we discuss the main results obtained in this work and their interpretation.

To fit our spectra we considered both a moderately (R_ISCO = 4 R_G) and rapidly (R_ISCO = 2 R_G) spinning BH. We note that in the case of a rapidly spinning BH the evolution of the system is mainly driven by variations in the sonic Mach number m_S. This solution is different with respect to the modeling of GX 339–4, in which m_S was fixed to 1.5 and the interplay between R_J and Ṁ_in was sufficient to describe the data (Marcel et al. 2019). We observed that this trend is statistically significant, since imposing constant values of m_S results in unphysical reflection-dominated scenario. A change in m_S arises from a change in the torque due to the jet (through the toroidal magnetic field). We see no clear reason for the jet torque to change on long timescales (> days), even if it is possible in principle. Although we discard this possibility here, one should nevertheless keep it in mind and look for traces of any correlation between these rapid accretion (increase in m_S) events and any jet signature.

On the contrary, in the fits performed keeping R_ISCO equal to 4 R_G, the m_S values are (almost) constant and the interplay between R_J and Ṁ_in again plays the main role. We therefore consider the fits performed in the moderately spinning BH scenario more physically reliable. This evidence is in line with the estimation of the spin value obtained by Fabian et al. (2020) and Guan et al. (2021), among others. Nevertheless, even in the case R_ISCO = 4 R_G, a (slight) increase in m_S is observed going from phase 2 to phase 3 (i.e., from 1.25 to 1.50). Imposing m_S to 1.5 in epoch 6, in continuity with epoch 7, results in a significantly worst fit where the reflection dominates over the direct JED emission. While this is slightly different with respect to the results obtained by Marcel et al. (2019, 2020), it is noteworthy that these authors did not have data above 40 keV and, more importantly, they did not directly perform spectral fits on the data. In future applications of the JED-SAD model, the inclusion of high energy NuSTAR and especially BAT data might therefore be necessary to highlight the role, if any, of m_S in shaping the spectra of BHTs in hard state and, particularly, their high energy cutoff.

5.1. The geometry of the accretion flow

The geometry of the accretion flow in the HS of MAXI J1820 has been inspiring a vivacious debate since its discovery. Several authors have proposed that, when in HS, the hot corona of MAXI J1820 contracts, while the accretion disk remains stable at the ISCO. This scenario was proposed first by Kara et al. (2019) and recently supported by Wang et al. (2021), based on NICER X-ray spectral and timing analysis. These results were also corroborated by modeling the reflection component in NuSTAR spectra, adopting two lamppost corona models (Buisson et al. 2019; You et al. 2021). However, Zdziarski et al. (2021) proposed that the truncated disk scenario can still explain the HS and HIMS of MAXI J1820. According to their spectral analysis performed with NuSTAR data on our epochs 1, 2, 3, and 4, the disk never extends down to the ISCO. A scenario where the disk is not only truncated, but also the inner edge of the disk approaches the BH throughout the hard state is also suggested by the spectral evolution of the quasi-thermal component responsible for the thermal reverberation lags (De Marco et al. 2021) and the evolution of the characteristic variability frequencies observed in the iron line (Axelsson & Veledina 2021). It is noteworthy that a trend for an approaching truncated inner disk is also qualitatively suggested by Type C QPOs (Buisson et al. 2019).

Regardless of the choice for R_ISCO and of the uncertainty in K_JEDSAD, all the fits performed by us show that the accretion disk in MAXI J1820+070 is truncated in hard state and its inner edge moves inward during the transition to the intermediate state. For R_ISCO = 4 R_G, we found R_J to be about 50–60 R_G in phase 1, 35–45 R_G in phase 2, and 30–35 R_G in phase 3. A bit of caution is required in taking these measurements at face value since the exact estimates of the truncation radius are instead critically dependent on our choices for R_ISCO and K_JEDSAD. Our spectral analysis was performed on a data set which includes the data set of Zdziarski et al. (2021) and reaches a similar conclusion on the geometry. However, our work includes four more observations and provides a coherent picture of the spectral evolution of the system, using a physical model. The calculated transition radii R_J are consistent with the lower limits on the truncation radius estimated by De Marco et al. (2021) by analyzing the quasi-thermal spectral component arising from the disk due to X-ray irradiation from the corona (see in particular their Fig. 9).

Finally, a multi-zone corona was invoked on the basis of both spectral (Zdziarski et al. 2021) and temporal (Dziełak et al. 2021) behavior of the system in the same observations. We highlight that the JED model is also a multi-zone model, as it takes into account the single contribution for each ring of matter that composes the JED, each with variable temperature and optical depth. Therefore, our approach is in line with other works, where a single uniform plasma cloud was proven insufficient to adequately describe the spectrum of black hole X-ray binaries (Nowak et al. 2011; Basak et al. 2017). Moreover a multi-zone model can naturally produce a spectral shape that is not exactly a power law. In the energy range of the iron line this can have a direct impact on the broadness of the line, and could explain the absence of extreme red wings that would require an accretion disk down to the ISCO.

5.2. The nature of the second reflection component

In Kara et al. (2019) the line profile revealed by the NICER data was clearly evolving across the HS and HIMS state, with a shrinking narrow core and a constant broad blurred base. This has been interpreted by invoking either a disk reaching the ISCO in hard state, or a contracting lamppost corona, which produces two reflection components, one stable (the broad inner one) and one variable (the narrow outer one). In this geometry, the top of this corona would mainly illuminate a distant region of the disk, which produces the narrow line, contrarily to the bottom of the corona, which would instead illuminate an inner region of the disk, where a strongly blurred line arises. Using two reflection components was also required in the spectral analyses by Buisson et al. 2019; Chakraborty et al. 2020; Zdziarski et al. 2021. In our modeling, the broad component in the iron line is also stable and changes very slightly from epoch 2 to 7 (see Fig. 8). However, the disk does not reach the ISCO and the broadening of the line profile is possibly due to the combination of two effects. First, a significant role is played by the underlying continuum. Previous papers reporting a highly ionized disk reaching close to the ISCO in the hard state consider a single Comptonization continuum, while a multi-zone Comptonization, as considered in the JED-SAD model, is more physically motivated. Furthermore, as a consequence of the high ionization parameter estimated for the plasma in the inner disk (logξ₁ ∼ 3.7 − 3.8), photons in the iron line tend to suffer multiple Compton scatterings leading to a further broadening of the line profile (Matt et al. 1996). In the JED-SAD configuration the corona does not have a lamppost geometry and an explanation analogous to Kara et al. (2019) cannot hold. An alternative possibility could be the inadequacy of our reflection models in dealing with real-life reflection spectra. We note that the reflection spectrum observed in accretion disks always results from the continuum of reflection spectra emitted at different radii. Using only one reflection component could just be an oversimplified approximation in this case. Nevertheless, given the high inclination of the system, some effects could also create an imbalance in the contributions from the different parts of the disk, lowering the contribution from the intermediate regions or increasing that from the outer regions. We suggest two mechanisms potentially responsible for such an imbalance, sketched in Fig. 10.

Fig. 10. Sketches representing two possible scenarios that might generate the double-reflection component used to model X-rays data of MAXI J1820+070 in hard state. In particular, in panel a we show that a sudden jump in the scale height of the disk might shadow part of the reflection coming from the intermediate regions. In panel b we show that the contribution from the on-the-other-side-of-the-BH JED could particularly illuminate the outer regions of the disk, leading to a boost in the outer reflection component.

One possibility could be that the disk, since viewed at such high inclination, might create a “shadow” zone where part of the reflection incoming between R_in, 1 and R_in, 2 is obscured (see Fig. 10, panel a). As highlighted by the Figure, a jump in the scale height of the disk, located sufficiently far away from its inner edge, could in principle generate such self-shielding effect. A disk flared in the outer region was also put forward by Zdziarski et al. (2021) and Axelsson & Veledina (2021) in order to explain the presence of the outer reflection component and the variability associated to the iron line. However, it is not obvious whether a sudden puffing up of the disk could happen or not in a classical Shakura-Sunyaev disk, like the SAD. A jump could arise from a transition in the mechanism contributing mainly to the opacity in the disk. In the outer regions of a Shakura-Sunyaev disk, the opacity is expected to be mainly due to free-free absorption, while Thomson scattering dominates the opacity in the inner regions of the disk (see Eqs. (2.16) and (2.19) of Shakura & Sunyaev 1973). The geometrical thickness of the disk scales with the radius R following a (slightly) different relation; it is proportional to R^−21/20 in the inner Thomson scattering-dominated, regions and to R^−9/8 in the outer regions. At the boundary between these two zones (i.e., at a radius R_jump) a jump is expected. We calculated R_jump according to our results for Ṁ_in and assuming the same value for α in both equations. We found R_jump between 3 × 10⁴ R_G and 7 × 10⁴ R_G in all the epochs considered. In our fits the location of the outer reflector R_in, 2 was kept frozen to 300 R_G. When the parameter is thawed, it is completely unconstrained, as shown in Fig. 11. Notwithstanding the uncertainty on R_in, 2, a boundary at beyond 10⁴ R_G is likely located too far away to produce an extra reflection component. However, a crucial role might be played by irradiation from the inner regions of the disk to the outer regions, an ingredient neglected in Shakura-Sunyaev disks. Beyond some boundary radius, we expect that the upper layers of the disk should be heated up by the impinging photons coming from the inner regions, leading to evaporation and to a larger effective scale height.

Fig. 11. Contour plots for the parameter Rin, 2 when left free and not tied to RJ.

This effect would be efficient only where thermal winds are launched and unfortunately this region is also located beyond 10⁵ − 10⁶ R_G, where the sound speed of the plasma overcomes the Keplerian speed (see, e.g., Higginbottom et al. 2017). This argument would exclude that the optical winds detected in MAXI J1820 in HS (Muñoz-Darias et al. 2019; Sánchez-Sierras & Muñoz-Darias 2020) could play the role of shadowing reflection from intermediate regions. However, the presence of a large-scale vertical magnetic field make real-life accretion disks more layered and puffed up (Zhu & Stone 2018; Jacquemin-Ide et al. 2021). It is therefore possible that the boundary radius at which the disk is geometrically thickened by irradiation may be located even closer and reasonably impacts the reflection spectrum observed from the disk at high inclinations. In order to confirm this, irradiation must be properly taken into account in models of magnetized disks around BHs¹¹. However, this is far beyond the scope of the present paper.

An alternative scenario is sketched in Fig. 10, panel b. In principle, the radiation from the JEDs on the two different sides of the BH with respect to our line of sight could illuminate two differently located regions of the SAD, originating the two reflection components. Photons coming from the other-side-of-the-BH JED may not be able to impinge on the SAD before R_in, 2, producing an excess in the reflection component coming from this radius. The picture presented here is likely oversimplified, as it may require a finely tuned light bending of the radiation coming from the other side of the BH, and it does not take into account any possible occultation from the jet, expected to be optically thin.

In addition, it is unclear how the evolution, mainly in ionization, of the outer reflection component fits within the proposed scenarios (see Fig. 8). This trend could be explained with an approaching of R_in, 2, which was instead kept frozen to the assumed value of 300 R_G. As shown in Fig. 11, R_in, 2 remains essentially unconstrained by the fit, so that an approaching trend may be hidden here. If we consider the irradiation from the inner disk as being responsible for the puffing up of the outer disk and the subsequent self-shielding effect, the fading in X-ray luminosity of the system during epochs 2 to 8 could be somehow correlated to the observed variability in the outer reflection. However, this cannot be confirmed without a proper magnetohydrodynamic modeling of magnetized disks that takes into account how irradiation impacts the scale height of disks.

6. Conclusions

In this paper we reported on the spectral study of the BH transient MAXI J1820+070 in HS exploiting the JED-SAD accretion-ejection paradigm. This is the second object the model has been applied to, and the first time that Compton reflection is taken into account in detail. We investigated the spectral behavior of the system in eight epochs, spanning roughly 100 days, during which the system was always in HS and HIMS. Due to the uncertainty on the spin of the BH, we considered two possible values for R_ISCO, 4 R_G (a_* = 0.55) and 2 R_G (a_* = 0.95), respectively. In the first case only two parameters, R_J and Ṁ_in, drive the largest changes in the spectral modeling: R_J decreases throughout the period, while Ṁ_in first increases, during the rise, and then decreases, during the plateau and the decline (see De Marco et al. 2021, for the definition of these phases). For a rapidly rotating BH, Ṁ_in remains quite constant, while a decrease in the sonic Mach number m_S is observed and seems to drive the overall evolution of the accretion flow. This scenario is unlikely and suggests that the BH in MAXI J1820 is spinning rather moderately, as also found by Fabian et al. (2020) and Guan et al. (2021). The geometry of the system consists in a truncated disk (i.e., R_J never goes below 15 R_G) with the inner radius decreasing during the monitored period, in agreement with Zdziarski et al. (2021) and De Marco et al. (2021), among others. In order to successfully describe the spectra, two reflection components have to be taken into account for all the epochs considered (except for the first one): one component is highly ionized and originates from the edge of the SAD, while the other is less ionized and presumably originates from an outer region of the disk. While the inner reflection component is stable, the outer reflection component evolves in a way that is compatible with a scenario where the region responsible for such a component is approaching the BH. A self-shielding effect, due to both the high viewing angle and a flared outer disk, might offer a viable explanation for such a phenomenon.

It is worth noticing that the exact Ṁ_in − R_J values estimated here depend on the choice of R_ISCO, the uncertainty on the distance of the system, and some simplifying assumptions on the dilution factor ω and the jet power fraction b. However, we argued that these assumption do not affect the main conclusions provided by this work (the truncated disk geometry and the double reflection component), but only provide additional uncertainty on the best-fit values presented here. Further investigations are necessary to obtain more reliable and precise measurements of the truncation radius and the inner mass-accretion rate. One possibility would be to check if the obtained R_J and Ṁ_in values can be used to fit the radio observations of this system, as successfully done by Marcel et al. (2019) for GX 339–4. Results discussing the jet behavior and brightness using the JED-SAD model for MAXI J1820 will be presented in a forthcoming publication.

Acknowledgments

AM, MDS, AS, SM, ADA, TDS and TDR acknowledge a financial contribution from the agreement ASI-INAF n.2017-14-H.0 and from the INAF mainstream grant (PI: T. Belloni, A. De Rosa). AM and TDS acknowledge financial contribution from the HERMES project financed by the Italian Space Agency (ASI) Agreement n. 2016/13 U.O. POP, SB and JF acknowledges financial support for the CNES space national agency and the CNRS PNHE. We thank R. La Placa for fruitful discussion.

References

All Tables

All Figures

	Fig. 1. From top to bottom, the Swift/XRT light curves and the related soft X-ray HR, Swift/BAT light curves and the hard X-ray HR are shown. Energy ranges are reported in the see legends. The NuSTAR observations are highlighted with black dashed lines.
In the text

	Fig. 2. XRT HID of the source during the entire outburst. The epochs used in this work are highlighted by the colors and different symbols (Sect. 4.2), i.e., circle for phase 1 (rise), upper triangles for phase 2 (plateau) and lower triangles for phase 3 (decline).
In the text

	Fig. 3. Sketch representing the proposed two-reflection geometry underlying Model 2. From left to right: the black hole (in black), the inner JED (violet) and the jets (lilac), the outer SAD (from red to yellow).
In the text

	Fig. 4. Energy spectra with the best-fit models: Model 1 for epoch 1 and Model 2 for epochs 2–4, and residuals. Data: XRT (blue), NICER (green), NuSTAR (yellow–black), and BAT (red). Different line styles were adopted to distinguish between the different components: dot for JED, dash-dot-dot-dot for SAD, dash-dot for inner REFL, dash for the outer REFL.
In the text

	Fig. 5. Energy spectra, best-fit model (Model 2) and residuals for epochs 5–8. Data: XRT (blue), NICER (green), NuSTAR (yellow–black), and BAT (red). Different line styles were adopted to distinguish between the different components: dot for JED, dash-dot-dot-dot for SAD, dash-dot for inner REFL, dash for the outer REFL.
In the text

Fig. 6. JED solutions showing the impact of mS and RJ on the high energy cutoff. Left panel: different JED solutions obtained keeping RJ to 40 RG for different values of mS. Right panel: different JED solutions are produced for different values of RJ, here expressed as rJ = RJ/RG, and mS fixed to 1.5. The plot is used to display how increasing both mS and increasing RJ pushes the hard X-rays cutoff to higher energies.

In the text

	Fig. 7. Evolution of the best-fit parameters for RJ, Ṁin, and mS over the eight analyzed epochs. In all panels the triangles and bold lines are used for fits performed with RISCO at 4 RG. The reported errors are likely underestimated and an additional systematic error of about 10%–20% has to be considered to account for the uncertainty on ω and KJEDSAD (see Sect. 4.5).
In the text

	Fig. 8. Best-fit models for inner and outer reflection in epochs 2–7. Inner reflection models have been multiplied by a factor of 10 for clarity.
In the text

Fig. 9. JED-SAD solutions for epochs 1, 5, and 8, used to represent phases 1, 2, and 3, respectively, calculated using the best-fit parameters for RJ, Ṁin, and mS reported in Table 3. The following values were adopted: for phase 1 RJ = 57 RG, Ṁin = 0.8 ṀEdd, and mS = 1.3; for phase 2 RJ = 44 RG, Ṁin = 1.7 ṀEdd, and mS = 1.3; and for phase 3 RJ = 33 RG, Ṁin = 1.1 ṀEdd, and mS = 1.5. In all panels, the radial distribution is divided in 31 portions (annuli), 30 for the JED and one for the SAD. Top and middle panels: how the actual shape z = H/R of the disk evolves with R. The distribution of the Thomson optical depth (τT, top) and the electron temperature (kTe, middle) are color-coded (see color bars at right). In the same plots, the event horizon of the BH is represented by a half black circle. The total emitted spectrum is shown in the bottom panels as a solid black line. The contribution from each of the 31 considered annuli is highlighted with a dotted line for the SAD region and dashed lines for the JED regions. The sections of the plot corresponding to the energy ranges that are not covered by XRT, NICER, NuSTAR, and BAT data are shown with a gray background.

In the text

Fig. 10. Sketches representing two possible scenarios that might generate the double-reflection component used to model X-rays data of MAXI J1820+070 in hard state. In particular, in panel a we show that a sudden jump in the scale height of the disk might shadow part of the reflection coming from the intermediate regions. In panel b we show that the contribution from the on-the-other-side-of-the-BH JED could particularly illuminate the outer regions of the disk, leading to a boost in the outer reflection component.

In the text

	Fig. 11. Contour plots for the parameter Rin, 2 when left free and not tied to RJ.
In the text