A broadband view on microquasar MAXI J1820+070 during the 2018 outburst
J. Rodi, A. Tramacere, F. Onori, G. Bruni, C. Sánchez-Fernández, M. Fiocchi, L. Natalucci, P. Ubertini
DDraft version January 27, 2021
Typeset using L A TEX twocolumn style in AASTeX63
A broadband view on microquasar MAXI J during the 2018 outburst
J. Rodi, A. Tramacere, F. Onori,
1, 3
G. Bruni, C. S´anchez-Fern´andez, M. Fiocchi, L. Natalucci, andP. Ubertini INAF - Istituto di Astrofisica e Planetologia Spaziali, via Fosso del Cavaliere 100, 00133 Roma, Italy Department of Astronomy, University of Geneva, Ch. d’Ecogia 16, 1290, Versoix, Switzerland INAF - Osservatorio Astronomico d’Abruzzo, via M. Maggini snc, I-64100 Teramo, Italy European Space Astronomy Centre (ESA/ESAC), Science Operations Department, 28691 Villanueva dela Ca˜nada, Madrid, Spain (Received -; Revised -; Accepted -)
Submitted to ApJABSTRACTThe microquasar MAXI J1820 + 070 went into outburst from mid-March until mid-July 2018 withseveral faint rebrightenings afterwards. With a peak flux of approximately 4 Crab in the 20 −
50 keV,energy range the source was monitored across the electromagnetic spectrum with detections from radioto hard X-ray frequencies. Using these multi-wavelength observations, we analyzed quasi-simultaneousobservations from 12 April, near the peak of the outburst ( ∼
23 March). Spectral analysis of thehard X-rays found a kT e ∼
30 keV and τ ∼ CompTT model, indicative of an accreting blackhole binary in the hard state. The flat/inverted radio spectrum and the accretion disk winds seen atoptical wavelengths are also consistent with the hard state. Then we constructed a spectral energydistribution spanning ∼
12 orders of magnitude using modelling in
JetSeT . The model is composedof an irradiated disk with a Compton hump and a leptonic jet with an acceleration region and asynchrotron-dominated cooling region.
JetSeT finds the spectrum is dominated by jet emission upto approximately 10 Hz after which disk and coronal emission dominate. The acceleration regionhas a magnetic field of B ∼ . × G, a cross section of R ∼ . × cm, and a flat radiospectral shape naturally obtained from the synchroton cooling of the accelerated electrons. The jetluminosity of > × erg/s ( > . L Edd ) compared to an accretion luminosity of ∼ × erg/s, assuming a distance of 3 kpc. Because these two values are comparable, it is possible the jet ispowered predominately via accretion with only a small contribution needed from the Blanford-Znajekmechanism from the reportedly slowly spinning black hole. Keywords:
Black Holes: individual (MAXI J1820 + 070) — X-rays: binaries — radiation mechanisms:non-thermal INTRODUCTIONThe term ”microquasar” was first applied to the per-sistent black hole candidate (BHC) 1E 1740 . − Corresponding author: J. [email protected] tions between radio and X-ray luminosities (Gallo et al.2003), indicating a relationship between the emissionmechanisms despite a large physical separation betweenthe two. Additionally, this correlation holds also forsupermassive BHs in AGN when accounting for mass(Merloni et al. 2003), thus linking the mechanisms instellar mass and supermassive BHs. Therefore under-standing the jet and X-ray components in microquasarscan shed light on AGN.The low-mass X-ray binary MAXI J1820 + 070(=ASASSN-18ey) was first detected on 6.59 March a r X i v : . [ a s t r o - ph . H E ] J a n Rodi et al.
Figure 1. (Top) The light curves for
Swift /BAT 15 −
50 keV(black diamonds), 0 . −
10 keV
Swift /XRT (green squares),and 4.7 GHz RATAN (blue triangles) for the initial phase ofthe outburst. The time span of observations analyzed in thiswork are denoted in red. with the All-Sky Automated Survey for Super-Novae (Shappee et al. 2014) and was detected ∼ MAXI /GSC at 11 March 2018 19:48 UTC(Kawamuro et al. 2018). With a peak flux of ∼ −
50 keV energy band (Roques & Jourdain 2019)and a long decay, the source was a good candidate for nu-merous observing campaigns across the electromagnetic(EM) spectrum (e.g. Mu˜noz-Darias et al. 2019; Tuckeret al. 2018; Stiele & Kong 2020; Bright et al. 2020) toexplore various aspects of the source. Combining ob-servations from various campaigns enables studying thevarious emission processes together.Therefore we compiled quasi-simultaneous observa-tions from public archives, Astronomer’s Telegrams, andGamma-ray Coordination Network Circulars to con-struct the widest possible frequency range, we were ableto find detections covering nearly 12 orders of magni-tude from the meter-wavelength frequencies to hard X-rays on 12 April 2018 (MJD 58220). With this spectralenergy distribution (SED), we studied the spectral com-ponents independently before investigating them jointlyby constructing a model consisting of a leptonic jet, an irradiated disk, and a corona, using the JetSeT soft-ware . OBSERVATIONSFigure 1 shows the initial period of MAXI J1820+070outburst in several wavelengths across the spectrumusing data
Swift /BAT (black diamonds),
Swift /XRT(green squares) (Stiele & Kong 2020), and 4.7 GHzRATAN (blue triangles) (Trushkin et al. 2018). TheXRT and RATAN data have been normalized to be onthe same scale of the BAT data. The period of the ob-servations used in this work are bracketed in red.In the following, we give information about the differ-ent simultaneous observations collected from archives,covering different bands from radio to gamma-ray onApril 12, 2018. Further details can be found in Table 1.2.1.
JVLA
We retrieved calibrated Karl G. Jansky Very LargeArray (VLA) data for experiment VLA/18A-470 fromthe National Radio Astronomy (NRAO) online archive.The JVLA antennas, in A configuration, were split in3 subarrays in order to obtain simultaneous observa-tions at 6 different central frequencies (4.7 GHz, 7.5GHz, 8.5 GHz, 11 GHz, 20.7 GHz, 25.5 GHz). Datawere imaged using
CASA ( Common Astronomy SoftwareApplications package) version 5.6.2 following stan-dard procedures. 2.2. ALMA
Atacama Large Millimiter/submillimiter Array(ALMA) data for project 2017.1.01103.T were retrievedfrom the ESO archive, and pipelined at the the Italiannode of the European ALMA Regional Centre (INAF-Istituto di Radioastronomia, Bologna). Imaging wasperformed with
CASA version 5.1.1, separately for eachone of the 4 spectral windows (spw) present in the data(Band 7, spw 5, 7, 9, 11) corresponding to the followingcentral frequencies: 336 GHz, 338 GHz, 348 GHz, 350GHz (Bonato et al. 2018).2.3.
VLT/X-shooter
A number of observations of MAXI J1820 + 070were performed with the X-shooter spectrograph (Ver-net et al. 2011) in the framework of the ESO program0101.D-0356(A). We retrieved the processed spectra ob-tained during the 2018 outburst on 12 April from theEuropean Southern Observatory (ESO) archive science https://jetset.readthedocs.io/en/latest/ https://casa.nrao.edu roadband View of MAXI J (cid:48)(cid:48) ×
11, 1.2 (cid:48)(cid:48) × (cid:48)(cid:48) ×
11 for the UVB, VIS and NIR arm, respec-tively. This configuration yields a spectral resolutionR= λ /∆ λ of 4100, 6500 and 4300 for the UVB, VIS andNIR arm, respectively. The observing conditions wheregood with a seeing of 0.47 (cid:48)(cid:48) and an average airmass of thesource during the acquisition of 1.3 (cid:48)(cid:48) . The total exposuretimes are 1640 s, 1300 s and 1520 s for the UVB, VIS,and NIR arm, respectively. The reduced spectra havebeen corrected for the foreground extinction using theCardelli function (Cardelli et al. 1989) with R(V)=3.1and A V =0.627 mag (Schlafly & Finkbeiner 2011, via theNASA/IPAC Extragalactic Database (NED)).In order to estimate the slit loss effect in the X-shooterspectra, we first applied standard aperture photometryon the i (cid:48) acquisition image using the iraf task phot .The zero point was calibrated using the stars in thePanoramic Survey Telescope and Rapid Response Sys-tem (Pan-STARRS1 Flewelling et al. 2016) catalog.From the aperture photometry we obtain an i (cid:48) bandapparent magnitude of m AB = (12.20 ± λ F i = (2.01 ± × − erg s − cm − which is in agreement with the averageflux measured from the spectrum in the 7300-7600 ˚Awavelength range: λ F i =(1.8 ± × − erg s − cm − .2.4. XMM-Newton/EPIC-pnXMM-Newton
ToO observations were carried outfrom 2018-04-12 07:27:58 to 09:39:28 UTC (ob-sid 0820880501) using burst mode. The Eu-ropean Photon Imaging Camera (EPIC)-pn datawere analyzed using the standard procedures withthe Science Analysis System (SAS) software version xmmsas 20190531 1155-18.0.0 .2.5. INTEGRAL
The
INTErnational Gamma-Ray Astrophysics Labo-ratory ( INTEGRAL ) observed MAXI J1820 + 070 ev-ery 2–3 days between March 16, and May 8, via a se-ries of Target of Opportunity (ToO) observations. Forthis work, we selected the data covering the interval(UTC) 11 April 2018 23:41:01 to 12 April 2018 1114:00:21 (
INTEGRAL revolution 1941). Here we focuson the analysis of data provided by the Integral SoftGamma-Ray Imager (ISGRI; 18 − Table 1.
MAXI J1820 + 070 observations logInstrument Start Time (UTC) Stop Time (UTC)VLITE 07:25:00 12-04-2018 13:07:00 12-04-2018.JVLA 07:15:00 12-04-2018 13:14:50 12-04-2018ALMA 08:13:18 12-04-2018 09:20:16 12-04-2018X-shooter 07:41:08 12-04-2018 08:15:48 12-04-2018EPIC-pn 07:27:58 12-04-2018 09:39:28 12-04-2018OMC 23:41:01 11-04-2018 14:00:21 12-04-2018ISGRI 23:41:01 11-04-2018 14:00:21 12-04-2018 Frequency (Hz) F l u x d e n s i t y ( J y ) VLITEJVLAALMAOMCX-shooter
Figure 2.
Jet emission between radio and optical bands, asreconstructed from VLITE, JVLA, ALMA, and X-shooterobservations. The blue dashed line is a broken power law,used to identify the synchrotron peak frequency and flux den-sity. The red line is a power-law fit of the most expandedregion of the jet, considered as a physically separated com-ponent. The optical flux from OMC is also shown, althoughnot considered for the fit. the upper layer of the detector plane of the Imager onBoard the INTEGRAL Satellite (IBIS) telescope (Uber-tini et al. 2003) and by the Optical Monitoring Camera(OMC) (500 −
600 nm) instruments. The data were ana-lyzed using the Offline Science Analysis software (OSA)v11.0 available at the
INTEGRAL
Science Data Center(ISDC). We followed standard analysis procedures. RESULTS AND DISCUSSION3.1.
The compact jet emission Rodi et al.
Table 2.
Collected flux densities for the jetmodelling.Instrument Frequency Flux density(Hz) (Jy)VLITE 3.39E+08 0.033 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± With the collected flux densities between radio andUV bands, we built the jet SED. In addition to theJVLA and ALMA data mentioned above, we consideredVery Large Array Low-band Ionosphere and TransientExperiment (VLITE, Clarke et al. 2016), also collectedon 12 April 2018 (Polisensky et al. 2018). We showin Figure 2 data from JVLA, ALMA, X-shooter, andOMC. For X-shooter, we considered only the part of thespectrum not affected by absorption/emission features,and averaged values for each of these intervals to obtaincontinuum flux density values. In this way, we calcu-lated 12 photometric measurements (Table 2) coveringthe interval from NIR to UV. The overall shape of theX-shooter spectrum shows a break, resulting in a changeof the spectral index, at about 2 × Hz. This is mostprobably the frequency at which the broadband SED isno longer dominated by the jet synchrotron emission, while the accretion disk thermal emission increases (seeSec. 4). The single value from OMC is in good agree-ment with the X-shooter photometry.We fitted the radio to optical data set with a broken-power law to identify the synchrotron peak frequencyand the spectral slopes of the optically thin and thickregions (Russell et al. 2013). We used the Astropy
BrokenPowerLaw1D function (see Astropy Collabora-tion et al. 2013; Price-Whelan et al. 2018) adoptingthe
LevMarLSQFitter routine to perform a Levenberg-Marquardt least squares statistic. We excluded fromthe fit the X-shooter points above the third one (Op-tical and UV ranges) since they show a turn-up of theflux density liekly due to accretion disk emission. Sim-ilarly, we did not consider data points below 20 GHz,since they show a different slope, probably belonging toa physically distinct (and more expanded) radio jet com-ponent. We obtained a synchrotron peak frequency of1.6 ± × Hz. The estimated slope for the opticallythick part of the spectrum is α thick = 0 . ± .
02 and α thin = − . ± .
01 for the optically thin one (adoptingthe convention S ∝ ν α , where S is the flux density, ν thefrequency, and α the spectral index). The slope of thelower frequency radio SED (0.3 −
10 GHz), fitted withthe
PowerLaw1D
Astropy function, is α = 0 . ± .
02. InSec. 4, a detailed physical model of this source is pre-sented and provides a more precise estimate of the jetparameters. 3.2.
Accretion disk winds
Mu˜noz-Darias et al. (2019) discovered accretion diskwinds in MAXI J1820 + 070 during both the 2018 hardstate rise and decay. Absorption wind signatures weredetected in the blue wings of He I λ λ v t ) of1200 km/s in one of the epochs. For H α , both a blue-wing broadened emission line profile, implying a windcomponent of 1800 km/s, and a superimposed absorp-tion trough with a v t =1200 km/s were found. Thoseauthors collected several epochs from 15 March to 4November 2018, allowing them to follow the evolutionof the winds from the hard state to the disappearanceduring soft state, and back.The VLT/X-shooter data presented in this work addan epoch to monitoring in Mu˜noz-Darias et al. (2019),falling in an uncovered time window of one month be-tween 26 March and 23 April. We normalized thesespectra by dividing them for the continuum emission,fitted with a third order spline3 function by using the IRAF task continuum . These spectra are rich withemission lines from the UVB to the NIR arms. Follow-ing the Mu˜noz-Darias et al. (2019) analysis, we explored roadband View of MAXI J N o r m a li s e d f l u x He I - 5876 Å H - 6563 Å He I - 6678 Å Figure 3.
Accretion disk wind absorption features in the VLT/X-shooter optical spectrum. Normalised flux errors are shownas a cyan shaded area. Troughs deeper than 3 × RMS are highlighted in purple. The RMS value is 0.0027 for the He I λ α and He I λ ± RMS intervals. the presence of wind signatures linked to the mentionedemission lines (He I and H α ), and found a significantabsorption on the left wing of He I λ v t of880 km/s. Further blue-ward absorption features arevisible, but since they are narrower and not connectedto the previous ones, we consider those as not related tothe accretion disk wind. The same is true for the nar-row absorption features detected blue-wards of the H α emission line (Fig. 3, central panel). For the He I λ v t of -825km/s.During this period, we observe strong asymmetries inthe emission lines which are commonly observed in linesemitted from the disc, particularly in the He I λ α . Therefore, we explored line profile propertiesby applying multi-component Gaussian fits using the python packages curvefit and leastsq . In Figure 4 weshow the result of this analysis for H α (left panel) andHe I λ α line analysis has beenperformed in the wavelength region 644 −
670 nm, whichincludes the feature of interest and the local continuum.In this case we have ignored from the fit the He I emis-sion line, which falls at the end of the analyzed wave-length range. The H α profile is well modelled by twonarrow Gaussian components, and only one broad Gaus-sian component is needed to fit the red wing of the emis-sion line. We note the absence of a blue-shifted broadwing, which has been observed in Mu˜noz-Darias et al. (2019), as well as p-cygni profiles signatures. Howevera forest of narrow absorption lines is clearly visible inthe blue region of the H α . The two narrow componentsare characterized by a central wavelength ( λ c ) and a fullwidth at half maximum (FWHM) of λ c = (6574.4 ± ±
20) km/s and λ c = (6559.8 ± ±
18) km/s, respectively. While thebroad red wing is centered at λ c = (6583.0 ± ± α rest frame wave-length we derive an outflow velocity of v = (923 ± ∼
667 km/s.For the He I λ −
672 nm, which includes also thelocal continuum but excludes the wavelength range inwhich the H α falls. The He I λ λ c = (6679.2 ± ± ∼
450 km/s with respect to the first com-ponent, and are characterized by λ c = (6669.4 ± ±
12) km/s and λ c = (6689.7 ± ±
29) km/s, respectively.As a whole, the detected optical disk wind featuresshow properties with in between what was found in thehard and soft state epochs collected by Mu˜noz-Dariaset al. (2019), confirming the decreasing trend of the op-tical wind between the two states of the source.3.3.
Soft X-ray
We fitted the
XMM Newton /EPIC-pn spectrumin the 0 . −
12 keV energy range in
XSPEC usinga(
Tbabs*Powerlaw model. With this model we de-rive the following parameters: nH = 0 . ± . × Rodi et al. . . . . . N o r m a li s e d F l u x H α
645 650 655 660 665 670Wavelength [nm] − . . . . . . . . N o r m a li s e d F l u x He I λ
664 666 668 670 672Wavelength [nm] − . . . Figure 4.
Fits of the emission line profiles for H α (left panel) and He I λ cm − , kT bb = 0 . ± .
03 keV, and Γ = 1 . ± . χ /ν = 1 .
00. 3.4.
Hard X-ray
We fitted the
INTEGRAL /IBIS/ISGRI spectrum inthe 30 −
400 keV energy range. A systematic 1.5 % er-ror was added to the data, following OSA 11 standardrecommendations . A power-law fit to the data in XSPEC found a photon index of Γ = 2 . ± .
01 and a χ /ν =71 .
80. The spectrum deviates from a simple power-lawmodel, especially at high energies, with residuals sug-gesting a Comptonized spectrum. Fitting the data witha
CompTT model using a photon temperature of 0.24 keVfixed to the kT bb value from XMM finds a better fit, with k T = 36 . ± . τ = 1 . ± . , and χ /ν = 6 . Reflect ) withthe reflection fraction fixed to 1, the fit improves to χ /ν = 3 .
45 with k T = 38 ± τ = 1 . ± . Reflect*(CompTT)+cutoff with Γ = 1 . . k T = 27 ± τ = 2 . . −
400 keV using the model
Tbabs*Reflect*(diskbb+CompTT)+Tbabs*cutoff found best-fit parameters of kt BB = 0 . ± .
01 keV, kT = 27 ± τ = 2 . ± . χ /ν = 0 . Table 3.
Irradiated disk fit parameters. diskir diskir+po kT disk (keV) 0 . ± .
007 0 . ± . . ± .
02 1 . ± . kT e (keV) 58 ± ± L C /L D . ± . . ± . f out (1 ± × − (4 ± × − log( r out ) 3 . ± .
04 3 ± po − . ± . po − . ± . χ /ν .
30 0 . Using this joint spectrum, we calculated the accretionluminosity in the 1 −
200 kev energy range and found avalue of ∼ × erg/s for a distance 3 kpc.3.4.1. IR − Hard X-ray Spectrum
Subsequently, we fit our data from the near-IR to hardX-rays using an irradiated disk model to compare withShidatsu et al. (2018), which analyzed a similar energyrange using observations from 24 March. The irradiateddisk model accounts for the effects of the Comptonizedemission on the accretion disk and the soft-excess thatis seen in the hard state (Gierli´nski et al. 2009). Fig-ure 5 shows the spectrum from 0 . −
400 keV with the diskir model shown as a solid red line and the power-law component of the diskir+po model as a dashedblack line. The power law is used to model the high-energy cutoff power-law component in the previous sec-tion. Table 3 contains the fit parameters using a diskir roadband View of MAXI J
Figure 5.
MAXI J1820 + 070 spectrum from 0 . − diskir model is shown with a solid red line andthe po component is shown as a black dashed line. model with and without an additional power-law compo-nent. We found that including a power-law componentimproved the fit at high energies and reduced the χ /ν from 1.30 to χ /ν = 0 .
97 with an f-test probability of3 . × − .The origin of the power-law component is unclear. Asshown below in Figure 7, the expected jet flux is toolow for the component to be jet emission at those en-ergies. However, the emission could possibly be fromComptonization of non-thermal electrons as in the caseof GRS 1716 −
249 ( ? ). BROADBAND SED MODELLINGWe modelled the broadband SED of MAXI J1820+070using a combination of jet leptonic models and irradi-ated disk and corona model implemented in Jets SEDmodeler and fitting Tool (
JetSeT ) (Tramacere 2020;Tramacere et al. 2011, 2009). A more accurate descrip-tion of the model is discussed in Tramacere (in prep.).We assume that the optical/UV up to keV energies isdominated by disc irradiation and coronal emission. Theemission in the mm to optical region is dominated by thenon-thermal emission of leptons accelerated in the jetby shock and/or stochastic acceleration, and we assumethat the break at ≈ . × Hz is due to the transitionfrom the optically thin to the optically thick synchrotronemission. The radio emission is dominated by the ter-minal part of the jet that starts beyond the accelerationregion and extends up to a distance of ≈ × cmaccording to Bright et al. (2020). A schematic view ofthe model is provided in Fig. 64.1. Individual Model Components Description
In the following we describe the implementation ofeach model component. https://jetset.readthedocs.io/en/latest/ Irradiated Disk and Hot Corona
To model the UV to hard-X-ray emission we haveused the disk Comptonization plus disk irradiationmodel,
DiskIrrComp implemented in
JetSeT . The
DiskIrrComp is based on the diskir model (Gierli´nskiet al. 2009) and the Comptonization model of Zdziarskiet al. (2009). In detail, we assume that a classicalmulti-temperature disk with an inner temperature T Disk and an extension R in = 3 R S to R out , expressed bythe dimensionless parameters r in = R/R in and r out = R/R out . The disk spectrum is modified due to the re-processing of irradiated emission from the disk itself andfrom the corona Compotonization tail. The corona emis-sion is described by a power law with an exponentialcutoff with a photon index Γ
Comp and a cut-off energy E Comp = kT e where kT e is the corresponding electrontemperature. The Compton hump component is de-scribed by a power-law with exponential cut-off witha photon index Γ hump and a cut-off energy E hump . Werefer to this model as Comp. hump . The normaliza-tion of the Compton tail component is parameterizedas a fraction of the disk luminosity L Disk according to L ratioComp = L C /L Disk . The total bolometric flux will be L bol = L Disk + L rep + L C , where L rep represents the ther-malized fraction f in of L C thermalized within r in and r irr = R irr /R in , where R irr is the radius of the innerdisk irradiated by the Compton tail. A fraction f out ofthe bolometric luminosity will irradiate the outer disk.The irradiation creates a shoulder with a spectral trend f out ∝ L bol ν − that extends between ν = 3 kT ( r out )and ν = 3 kT ( r t ), where r t is the transitional radius be-tween gravitational and irradiation energy release. Thiseffect depends strongly on r out and f out , and it is presenteven without corona Comptonization, because it repre-sents the disk self-irradiation. The presence of a Comp-tonization component will provide a further heating ofthe disk in the inner part modifying the pure gravita-tional temperature profile.4.1.2. Pre-acceleration and Acceleration Region
We assume that electrons in the pre-acceleration re-gion close to the base of the jet are described by athermal plasma with cooling dominated by adiabaticlosses. Once the particles approach the acceleration re-gion they are accelerated under the effect of diffusiveshock acceleration and/or stochastic acceleration andthe corresponding energy distribution can be modeledby a power-law with a high-energy cutoff N e,acc ( γ ) = N γ − s exp( − γ/γ cut ) (1)where the value of γ cut takes into account the balancebetween cooling and acceleration terms. The index s Rodi et al. is dictated by the competition of the acceleration timescales and escape time scales (Tramacere et al. 2011).We assume that the acceleration region extends from z startacc to z endacc , with cross section R acc equal to the aver-age cross section of the jet at z = ( z startacc + z endacc ) /
2, with z endacc − z startacc = 2 R acc . The emission from the acceler-ation region is reproduced using the jet leptonic model Jet implemented in
JetSeT , and we refer to it as
JetAcc (Tramacere in prep.).4.1.3.
Radio Jet
To model the radio jet emission we have used the
JetSeT multi-zone radio jet model
RadioJet . Thismodel implements a continuous jet as a sum of N c singlezones, following the approach of Kaiser (2006), where foreach zone the values of R and B are ruled by Eq. 3, andthe particle density scales as N s,i = N s, ( z s, /z i ) m N (2)where N s, is the initial density of emitters at the start-ing point of the radio jet z s, , z i is the average positionof the i th component, and m N is the index of the par-ticle density law fixed to 2. The initial particle densityis a fraction N frac of that present in the accelerationregion and we fix it to 1. The radio jet extends from z startradio = ( z acc + R acc ) K startR to z endradio ( z acc + R acc ) K endR ,where K startR and K endR are free parameters. In thepresent analysis we fix K startR = 1 , and K endR is fixed inorder to match the value of 1 × cm according to theanalysis presented in Bright et al. (2020). The particledistribution in each region has the same spectral law asin the acceleration region, but we decrease the value of γ cut to take into account the effect of the cooling whenthe particles leave the acceleration region. In our analy-sis we take into account only synchrotron cooling and weevolve γ cut according to Eq. 27 in Kaiser (2006). Moredetails about the connection between the accelerationand radio are discussed in Sec. 4.34.2. Phenomenological Model Setup
As a first step we set the geometrical properties of thejet, i.e. we define the extent of the pre-acceleration, ac-celeration, and radio emission sites, and the values ofthe magnetic field. We assume that the jet is launchedat distance z from the BH, with an initial cross sec-tion R , and that the bulk Lorentz factor (Γ jet ) of thejet is constant over the full jet extent. The accelera-tion region starts at a distance z startacc with a width equalto jet cross section diameter R acc = 2 R ( z acc ), and wetreat it as a spherical region. The radio region startsat z startradio = K startR ( z acc + R acc ) and ends at a distance z endradio = K endR ( z acc + R acc ) (a scheme of the model is pre-sented in Fig. 6). According to Bright et al. (2020) we fix the distance of jet from the observer to the value of d = 3 kpc, the termination of the radio jet to the valueof z end ≈ × cm, and the value of the beamingfactor to δ = [Γ jet (1 − β jet cos( θ obs )] − ≈ .
2, using thevalues of β jet = 0 .
89 and θ obs = 63 ◦ reported in Brightet al. (2020). We assume a ballistic jet model (Kaiser2006; Begelman et al. 1984) characterized by B ( z ) ∝ B ( z /z ) m B R ( z ) ∝ R ( z/z ) m R (3) N ( z ) ∝ R ( z /z ) m N with m B ≈ m R = 1, and m N = 2 .
0. Thischoice assumes that the jet is very close to the ballisticregime, with a magnetic field dominated by the toroidalcomponent and justifies the assumption that the bulkLorentz factor is constant along the jet.We use a black hole mass of M BH = 8 M (cid:12) . The jetluminosity L jet is linked to the Eddington luminosity( L Edd ) according to L jet = 12 q jet L Edd (4)where L Edd ≈ . × ( M BH /M (cid:12) ) erg s − (Rybicki& Lightman 1986). It is worth noting that our q jet parameter is not linked directly to the accretion effi-ciency process, because the jet powering could, in prin-ciple, be supported also by other mechanisms such asthe Blandford–Znajek mechanism (Blandford & Znajek1977), that predicts electromagnetic extraction of en-ergy and angular momentum from magnetized accretiondisc surround a black hole. Hence, our q jet parametershould not be used to infer or constrain the accretion effi-ciency, and will be discussed in a more accurate physicalcontext in Tramacere (in prep.).We assume that the jet is launched at a distance z = z from the BH with z = 50 R S ≈ . × cm,where R S = (2 GM BH ) /c . The launching jet position z = z , in the current analysis is assumed constant, toreduce the model complexity, and it is chosen accord-ing to reference values published in previous analysis(Vila & Romero 2010). The initial radius of the jet isset to R ( z ) = 0 . z , resulting in an opening angle of θ open ≈ . ◦ . We impose that in the launching regionthe entire jet power is in the form of magnetic energy L jet = L B ( z ) = πU B ( z ) R ( z ) Γ jet β jet c (5)where U B = B / (8 π ), and setting q jet = 0 . B ≈ . × G.The value of m B can be constrained from the spectralindex of radio jet emission, α R ≈ .
15, according tothe Eq. 39 in Kaiser (2006) that refers to the case of roadband View of MAXI J
Figure 6.
A schematic view of the jet model setup. The purple region identifies the pre-acceleration region, the cyan regionidentifies the acceleration region, and the green region identifies the radio jet. The z axis on the bottom shows the starting andend end point of each region. The acceleration region is assumed to be spherical with a radius equal to the jet cross section.The vertical black lines in the radio jet region marks qualitatively the division of region in slices. strong radiative cooling and almost constant value ofthe electron distribution high-energy cutoff. Accordingto this scenario, that is very similar to what we expectin our case, we can rearrange Eq. 39 in Kaiser (2006)as: m B = 1 + m R − α R (6)that is similar to the trend of the thick radio spectrumdiscussed in Pe’er & Casella (2009a), and we obtain avalue of m B ≈ . m B and m R considered in the RadioJet emission.To constrain the value of z acc we impose that R acc = R ( z acc ), B acc = B ( z acc ) and N e,acc correspond to a syn-chrotron self-absorption frequency of ν t ≈ . × Hz.This value of ν t is obtained from the phenomenologicalfit of the optically-thin to optically-thick synchrotronemission between mm and optical data shown in Fig. 2.In order to solve this problem we combine the analyticalexpression of the synchrotron self-absorption frequency( ν t ) (Rybicki & Lightman 1986), evaluated at the peaki.e. α ν = 0 ν t = ν L (cid:104) π √ π qR acc N e,acc B acc f k ( s ) (cid:105) s +4 , (7)and of that the synchrotron emissivity (Rybicki & Light-man 1986) (cid:15) s ( ν ) : (cid:15) s ( ν ) = F ν d L V = 3 σ T cN e,acc U accB π (cid:112) ( π ) ν L f (cid:15) ( s ) , (8) where q is the electron charge, U accB is the value of themagnetic field, V is the volume of a spherical geome-try of volume V of radius R acc , s is the slope of theelectron distribution power-law and ν L = qB πm e c is theLarmor frequency, and where the functions f k ( s ) and f (cid:15) ( s ) are approximated to percent accuracy as reportedin Ghisellini (2013). The value of s is obtained usingthe optically thin spectral index ≈ . s = 2 α + 1 ≈ . N acce and then substitute in Eq.7, and we insert the functional form of B = B ( z acc ) and R = R ( z acc ) according to Eq. 3. The final equationsolved with respect to z acc reads: z acc = (cid:104)(cid:16) ν t πm e cqB z m B (cid:17) s +42 σ T ( B R ) f (cid:15) ( s ) z m r e π f k ( s ) F ν d L (cid:105) ψ (9)where r e = q / ( m e c ) is the classical electron radius,∆ m = m B − m R , and ψ = m − m B ( s +4) .Consequently, the starting position of the radio jet isset to z startradio = z endacc = z acc + R acc ≈ . × cm, withan extent derived from Bright et al. (2020) of z end (cid:38) z startradio The value of the cut-off of the electron distributionis set to γ cut = 60, in order to produce the peak ofthe synchrotron emission above the IR frequencies fora magnetic field B acc ≈ . × G, with a power-lawslope s ≈ . z acc can be used to derive thehadroninc content of the jet energetic in form of coldprotons. Following Vila & Romero (2010) we imposethat in the acceleration region of the jet the magnetic0 Rodi et al.
Table 4.
Phenomenological Setup ParametersInput Parameterspar. name units input value z cm 1 . × r cm 1 . × M BH M (cid:12) q jet F tν Jy 0 . ν t Hz 1 . × s ρ accp,B > m B m R B G 6 . × B acc G 1 . × − L accp erg s − > . × z startacc cm 2 . × z endacc cm 2 . × z acc cm 2 . × R acc cm 2 . × z startradio cm 2 . × z endradio cm (cid:38) × energy of the jet is in subequipartition with the bulkkinetic energy of the cold protons, a condition that ismandatory to allow the mechanical compressibility ofthe plasma (Komissarov et al. 2007). We define theparameter ρ accp,B = U p ( z acc ) /U B ( z acc ), where U p ( z ) = n p ( z ) m p c , and we require that ρ p,B >
1. This choicesets a value of cold proton luminosity in the accelerationregion L p ( z acc ) > . × erg s − .4.3. Model Fit and Results
Initial model setup
To optimize the model we use the composite model in-terface
FitModel provided by
JetSeT , that allows com-bining different models in a global model. This modelcan be optimized by inserting it to the
ModelMinimizerJetSeT plugin. In the current analysis we use a frequen-tist approach and we use the
Minuit ModelMinimizer option. We have used the
Data and
ObsData JetSeT tools to import the observed data, and we have added a5% systematic error in the range [1 × , × ] Hz,to avoid that the large inhomogeneity on the fractional error between radio and optical/UV data, could bias thefit convergence. For the error estimate we provide onlyerrors derived form the MIGRAD module of
Minuit , amore reliable estimate based on a Markov chain MonteCarlo (MCMC) will be presented in Tramacere (in prep.)The
DiskIrrComp model, the
Comp. hump model,and the
JetAcc are independent, on the contrary,
JetAcc and radio
RadioJet are bound.The initial values of the parameters for the
DiskIrrComp model are chosen according to the analysispresented in Sec 3.3 and Sec. 3.4. In detail, we set theinitial values of L disk = 1 × erg s − , of r out = 5000,of f out = 0 .
01, and L ratioComp = 4 . T Disk = 1 . × K, and the pa-rameters r irr = 1 . f in = 0 .
1, the choice adoptedin Gierli´nski et al. (2009), when the Comptonization ofthe outer disk is included in the irradiated disk.For the
JetAcc model, we fix θ jet = 63 ◦ , Γ jet = 2 . R acc =2 . × cm, and B acc ≈ . × G, we freeze theinitial value of z acc = 2 . × cm, and we leave freethe parameters for the electron distribution.The initial setup of the parameters of the RadioJet is more complex and we need to take into account thephysical connection with the acceleration region and thecooling process. This effect plays a crucial role, in-deed, as already discussed in Kaiser (2006) and Pe’er& Casella (2009a), the combination of synchrotron cool-ing and jet expansion (assuming a negligible contribu-tion from adiabatic cooling) will result in an asymp-totic value of γ cut ( t ), that can naturally explain the flatradio spectrum without the need to introduce signifi-cant particle re-acceleration in the radio jet. We fol-low the approach reported in Pe’er & Casella (2009a)(in the case of negligible adiabatic cooling) and we set m radioB = m radioR = m jet . The particle cut-off evolutionin the radio jet will evolve according to (Kaiser 2006): γ cut ( t ) = γ cut σ T B m e cπ ( f ) γ cut t − f ( t f − t finj ) (10)where f = 1 − m jet , and t = z /β jet c Γ jet , t inj = z inj /β jet c Γ jet and t = z / β jet c Γ jet , are the comovingtime scales. We freeze a starting value of z inj = z startacc ≈ . × cm.Another effect to take into account is the fact thatfor z > z acc the structure of the jet could change, forthis reason we leave free the parameters m jet with a fitboundary of [0.5,1.5], with an initial value of 1.18, thatis slightly larger than the value used for the phenomeno-logical constraining, but gives a better agreement withradio-to-optical data. roadband View of MAXI J Table 5.
JetSeT best fit model parametersmodel name par. name units best fit value error starting value fit boundaries frozen
CompHump E hump keV 26 14 20 [ 15 ; 35] False” Γ hump -0.5 2 -1.2 [ -2 ; 2] False DiskIrrComp T Disk
K 1 . × True” L Disk erg s − . × . × × [ 1 × ; 1 × ] False” r out . × . × × [ 1 ; – ] False” r irr Comp E Comp keV 150 100 140 [ 20 ; 200 ] False” L ratioComp f in f out × − × − DiskIrrComp r out . × . × . × [ 1 ; – ] False” f out . × − . × − × − [ 0 ; – ] False” L ratioComp JetAcc N e,acc cm − . × . × . × [0 ; – ] False” s γ cut . . R acc cm 2 . × . × . × [ 1 . × ; 3 . × ] False” z acc cm 2 . × True” B acc G 17986 1 . × − θ jet deg 63 True” Γ jet RadioJet z inj . × True” N frac K startR K endR m jet The density of emitters at the base of the
RadioJet ,is bound to be equal to the density of emitters in theacceleration region N e calculated according to Eq. 3, at z = z startradio , by fixing N frac = 1 .
0. We fix the values of K endR = 3000 and of K startR = 1.A list of the free and frozen and of the bounds is re-ported in Table 5 in the columns ‘staring values’, ‘fitboundaries’, respectively.4.3.2. Model fit results for the Disk and Corona emission
We fit first the
DiskIrrComp and
Comp. hump com-ponents restricting the fit range to ν = [5 × , ] Hzand we get χ = 152 for 98 degree of freedom ( N dof ),corresponding to reduced χ red = 1 .
55. The parametersvalues are reported in the upper part of Table 5.The best-fit parameters resulting from
JetSeT aresimilar to those obtained form the
XSPEC analysis forthe diskir model. In particular the r out value (3 . × R in and 3 . × R in ) and L C /L D (4.1 and 4.6), for JetSeT and
XSPEC respectively. The f out parame-ter, is unconstrained both for XSPEC and the
JetSeT .However, a well constrained value is obtained when thejet component is added as shown in section 4.3.3. Be-cause the
JetSeT model for the Comptonized emissionis phenomenological, the high-energy range of the irra-diated disk is fit as a cutoff power-law and thus is notdirectly comparable to the diskir parameters. For thatportion of the spectrum,
JetSeT found Γ = 1 .
64 (com-pared to 1.78 from diskir ) and E C = 150 keV (com-pared to kT e = 58 keV from diskir ). We do note that E C /kT e ≈ .
6, which falls within the predicted range ofratios between cutoff energy and electron temperature(Petrucci et al. 2000, 2001), suggesting that the valuesare in agreement, even though the incertitude on the
JetSeT value is quite large.4.3.3.
Model fit results for the jet emission Rodi et al. l o g ( F ) ( m J y ) l o g ( F ) ( ( e r g c m s H z )) Disk. Irr. Comp.radio jetbest fit modeljet acc.deabs jetset6 4 2 0 2 4 6log(E) (eV) 1012141618 l o g ( nu F ) ( m J y H z ) l o g ( F ) ( e r g c m s ) Disk. Irr. Comp.radio jetbest fit modeljet acc.deabs jetset8 10 12 14 16 18 20log( ) (Hz)505 r e s Figure 7.
The best-fit
JetSeT model of the broadband SED. Top panel: the F ν representation of the global model fit. Bottompanel: the νF ν representation. The red line represents the global model, the dashed lines correspond to the single components,the color is reported in the legend. The best fit parameters are reported in Tab. 5. The residuals plot is evaluated with respectto the νF ν representation. To fit the full band SED we freeze all the parametersin the
CompHump and
DiskIrrComp components, exceptfor r out , f out , and L C /L D , and we fit the global modelover the full SED band in the range ν = [5 × , ]Hz .The model fit converged with a final χ = 181for 122 degree of freedom ( N dof ), corresponding to a χ red = 1 .
48. The parameters values are reported inthe bottom part of Table 5, and the parameters derivedfrom the best-fit model are reported in Table 6. Re-garding the
DiskIrrComp , we note that adding the jetcomponent results in a better constraint on the valueof f out =(7 . ± . × − , that is in the expectedrange of other black hole binaries in the hard state roadband View of MAXI J Table 6.
Model parameters evaluated from the best-fitmodelpar. name units value setup value q jet > .
15 0.20 U e /U B N acce /N coldp <
94 – L accjet erg s − > . × – L accrad erg s − . × – L accB erg s − . × . × L acce erg s − . × – L accp erg s − > . × > . × (Gierli´nski et al. 2009). Moreover, restricting the fitstatistics to the same interval used in the XSPEC anal-ysis, we get a χ = 157 with N dof = 107, correspond-ing to a χ red = 1 .
6. Regarding the jet component, wenote that final best-fit model parameters did not changesignificantly from the input values, suggesting that thephenomenological setup was able to find a configurationvery close to the optimal one, even though the fit mightbe biased by the degeneracy among some parameters.We will investigate this problem in a forthcoming workTramacere (in prep.), by means of Bayesian approachbased on the MCMC technique. In general, we find thatour assumption based on the connection between a com-pact acceleration region feeding the extended radio jet isable to model self-consistently the UV-to-optical emis-sion, reproducing the observed flat radio spectrum. Inparticular, we find that, according to our best fit model,the particles in the radio region reach the asymptoticvalue of γ cut ≈
8, and keep it almost constant as resultof the decrease in the cooling synchrotron cooling ratedue to the jet expansion. This behaviour is in agree-ment with the results of Pe’er & Casella (2009b) andKaiser (2006) for the case of synchrotron cooling dom-inating over the adiabatic one. It is worth discussingsome specific parameters in detail: • q jet > .
15. This value is compatible with the in-put value q jet = 0 .
2. As already discussed in theprevious section, the q jet parameter is not linkeddirectly to the accretion efficiency because the jetpowering could, in principle, be supported also byother mechanisms such as the Blandford–Znajekmechanism Blandford & Znajek (1977), that takesinto account advection of magnetic flux from anaccretion disk surrounding the Black Hole. Hence,our q jet parameter can not be used to infer or con-strain the accretion efficiency. • U e /U B = 0 .
18. The U e /U B is not far fromequipartition, and it is obtained without provid-ing any constraint. This proves that the combina-tion of the phenomenological model setup and theminimization of the global model converged natu-rally toward a configuration close to the physicalequipartition of U e and U B , giving further supportto the choice of a compact acceleration region thatis connecting the pre-acceleration region to the ra-dio jet. • N acce /N coldp <
94. Since our model is leptonic, thecontent of cold protons can be derived from an-cillary conditions, as the condition that the mag-netic energy of the jet has to be in subequipartitionwith the bulk kinetic energy of the cold protons,in order to allow the mechanical compressibility ofthe plasma (Komissarov et al. 2007), and forma-tion of shocks/turbulent acceleration sites in theacceleration region. From the best fit model weget that to respect the condition ρ accp,B > . N e /N coldp < N e /N coldp = 10 (Celotti & Ghisellini 2008) usedin the case of relativistic jets with a leptonic ra-diative domination. Moreover, we note that thevalue of B acc obtained from the best fit did notrequire a significant change in the value of L B asderived from the phenomenological model setup,and demonstrating that constraining z acc based onthe value of ν t is naturally in agreement with for-mation of mechanical compression in the jet when U p > U B . • m jet = 1 .
2. The value of m jet is one the most crit-ical, indeed it dictates the topology and intensityof the magnetic field beyond the acceleration re-gion, and it is interesting to compare to the valueof m B that is used to model the jet below the ac-celeration region. The initial guess based on thevalue of α R has required a small modification inorder to reproduce the observed radio spectrum,and the final model naturally explains the almostflat radio spectrum as emission of the cooled elec-tron leaving the acceleration region. DISCUSSION AND CONCLUSIONSAs MAXI J1820 + 070 was observed numerous timesacross the EM spectrum during its outburst, thereare multiple works relevant to portions of our multi-wavelength analysis, though to date none study the4
Rodi et al. source in such a complete picture as is presented withour model from
JetSeT .Shidatsu et al. (2018) provides the most direct com-parison to the analysis in this work, though the sourcebehavior is different before and after 26 March (MJD58206). They found that the optical and near-IR emis-sion is not entirely from disk emission and thus in-cluded a power-law to their diskir model with a spec-trum described by fixed parameters kT disk = 0 .
35 keV, L C /L D = 70, f out = 5 × − , and R out = 10 R in . Ourfit found a considerably lower kT disk (0.12 keV), L C /L D (4.7), and R out (10 R in ), but a higher f out (4 × − ).We note that the source behavior between the two obser-vations is different with changes in the spectral hardnessin hard X-rays (Roques & Jourdain 2019), the develop-ment of type-C QPOs (Stiele & Kong 2020), and a re-duction in the size of the corona (Kara et al. 2019) thatcan possibly explain the differences.Following the work of Mu˜noz-Darias et al. (2019), weexplored the presence of disk wind signatures in ourVLT/X-shooter optical spectrum, as this data-set fallsbetween epochs 11 and 12 of Mu˜noz-Darias et al. (2019)campaign and adds an epoch in their uncovered timewindow (between 26 March and 23 April). We focus ourspectral analysis on the He I λ λ α wavelength regions. We found shallow p-cygni profilesand strong line asymmetries in all the three mentionedlines, while a broad outflow component is detected onlyin the red wing of the H α . Among the observed absorp-tion troughs, the one detected in the blue wing of He I λ v t =880 km/s, which is consistentwith the outflow velocity of v ∼
900 km/s, derived fromthe H α redshifted broad component. These propertiesindicate that at this epoch optical disk winds are stillpresent, although with slower velocities with respect towhat found in Mu˜noz-Darias et al. (2019). The authoralso report an evolution of the line profiles during theirmonitoring campaign and our observation confirms thistrend. In particular the observed H α profile can be in-terpreted as a continuation in the evolving pattern ofthe line between the epochs 9 and 12 shown in Figure2 of Mu˜noz-Darias et al. (2019). Similar spectral vari-ations were previously reported by Tucker et al. (2018)and ascribed to the orbital motion of the system. In-terestingly, some of the most conspicuous optical winddetections in Mu˜noz-Darias et al. (2019) occur in epochscorresponding to the hard state of the source, when ra-dio emission and strong jet activity are present (Brightet al. 2020) and the peak of the optical outburst of thesource is reported. This led the authors to the conclu-sion that the optical wind detected in MAXI J1820+070 is simultaneous with the jet. Our wind signatures de-tection, together with the results from our broad bandspectral analysis are consistent with this scenario.Our phenomenological analysis of the compact jetfound the data could be modelled by a broken power-lawwith α thick = 0 . ± .
02 and α thin = − . ± .
01. Com-bining observations from late March and early April,Russell et al. (2018) performed a similar analysis andfound spectral indices of α thick ∼ . α thin ∼ − . ∼ × Hz and acorresponding flux density of ∼ . B ∼ × G and R ∼ × cmusing equations from Shidatsu et al. (2011). Our modelpeaks at 1.6 ± × Hz with a flux density of ∼ . JetSeT .In particular the
JetSeT best-fit model gives a magneticfield in the acceleration region of ≈ . × G, and aregion radius of ≈ . × cm.The corresponding energy density of the magneticfield is ≈ . × erg/cm compared to the values of8 × erg/cm from Shidatsu et al. (2018).Additionally, we identify a separate radio spectralcomponents at frequencies below ∼
10 GHz, showing aninverted power-law spectrum with slope α = 0 . ± . JetSeT broad-band model by the
RadioJet component, and stems nat-urally from the cooling of the accelerated particle leavingthe acceleration region. Interestingly we find that thebest-fit index m jet ≈ . α = 1 − /m jet ≈ .
166 that is close to the value foundin the power-law fit. We note, that the small differencebetween the two values, is due to the fact the
JetSeTRadioJet model takes into account the data range fromradio-to-mm frequencies, differently from the power-lawfit, whose range extends up ≈ Hz. roadband View of MAXI J
JetSeT broadband model was able toreproduce the full SED taking into account both thedisk/corona emission, and the leptonic radiativelly dom-inated relativistic jet contribution. We found that therelativistic jet required a total energy of L jet ≥ . × erg/s, corresponding to 0.15 L Edd . This value representsa lower limit, since we assume that the hadronic contentof the jets is only in terms of cold protons, without a sig-nificant radiative contribution. The flat radio spectralshape stems naturally from the synchroton cooling of theelectrons in the acceleration regions, in agreement withprevious analyses (Kaiser 2006; Pe’er & Casella 2009a).In comparison, the accretion luminosity (6 × erg/s)is comparable to the lower limit of the jet luminosity.Thus in MAXI J1820 + 070, it is possible for the jetto be powered predominately via accretion with onlya small contribution from the Blanford-Znajek mech-anism, which in this case cannot provide much powersince the black hole spin is reported to be low (Bassi etal. 2020; Zhao et al. 2020). ACKNOWLEDGMENTSWe thank the Italian node of the European ALMARegional Centre (ARC) for the support. JR and GBacknowledge financial support under the INTEGRALASI-INAF agreement 2019-35-HH.0 and ASI/INAF n.2017-14-H.0. FO acknowledge the support of theH2020 European Hemera program, grant agreementNo 730970. The research leading to these resultshas received funding from the European Union’s Hori-zon 2020 Programme under the AHEAD2020 project(grant agreement n. 871158) F.O. acknowledges thesupport of the H2020 European Hemera program,grant agreement No 730970, and the support of theGRAWITA/PRIN-MIUR project: ” The new frontier ofthe Multi-Messenger Astrophysics: follow-up of electro-magnetic transient counterparts of gravitational wavesources ”. Based on observations with INTEGRAL, anESA project with instruments and science data centrefunded by ESA member states (especially the PI coun-tries: Denmark, France, Germany, Italy, Switzerland,Spain) and with the participation of Russia and theUSA. This research has made use of the services of theESO Science Archive Facility. Based on observationscollected at the European Southern Observatory underESO programmes 2017.1.01103.T (ALMA) and 0101.D-0356(A) (VLT). ALMA is a partnership of ESO (rep-resenting its member states), NSF (USA), and NINS(Japan), together with NRC (Canada) and NSC andASIAA (Taiwan), in co- operation with the Republicof Chile. The Joint ALMA Observatory is operated byESO, AUI/NRAO, and NAOJ. The National Radio As-tronomy Observatory is a facility of the National Sci-ence Foundation operated under cooperative agreementby Associated Universities, Inc.REFERENCES