The Collimated Radiation in SS 433. Constraints from Spatially Resolved Optical Jets and $\texttt{Cloudy}$ Modeling of the Optical Bullets

The microquasar SS 433 is well-known for its precessing, relativistic baryonic jets. Depending on their heating mechanism, the optical jet bullets may serve as a probe of the collimated radiation coming from the inner region close to the compact object. The optical interferometer VLTI/GRAVITY has allowed to spatially resolved the optical jets in SS 433 for the first time. We present here the second such observation taken over three nights in July 2017. In addition, we use the multi-wavelength XSHOOTER spectrograph at VLT to study the optical bullets in SS 433 in detail. GRAVITY reveals elongated exponential-like spatial profiles for the optical jets, suggestive of a heating mechanism acting throughout a long portion of the jet and naturally explained by photoionization by the collimated radiation. We also spatially resolve the movement of the optical bullets for the first time, detecting extended jet components corresponding to previous ejections. \texttt{Cloudy} photoionization models explain both the spatial intensity profiles measured with GRAVITY and the line ratios from XSHOOTER, and constrain the properties of the optical bullets and the ionizing radiation. We find that the latter must peak in the UV with an isotropic luminosity (as inferred by a face on observer) $\approx 10^{41}$ erg/s. Provided that the X-ray SED is sufficiently hard, the collimated X-ray luminosity could still be high enough so that the face on observer would see SS 433 as an ULX ($L_X \lesssim 10^{40}$ erg/s) and it would still be compatible with the H/He/He+ ionization balance of the optical bullets. The kinetic power in the optical jets is constrained to $3-20 \times 10^{38}$ erg/s, and the extinction in the optical jets to $A_V = 6.7 \pm 0.1$. We suggest there may be substantial $A_V \gtrsim 1$ and structured circumstellar extinction in SS 433, likely arising from dust formed in equatorial outflows.


INTRODUCTION
After more than forty years since its discovery, SS 433 (Stephenson & Sanduleak 1977;Clark & Murdin 1978) remains a unique object primarily due to its relativistic, precessing baryonic jets. They were first discovered through broad emission lines of hydrogen and helium moving across its optical spectrum with extreme Doppler shifts (Margon et al. 1979). Their blue/redshifts follow a kinematic precession model (Fabian & Rees 1979;Margon 1984), according to which the jets precess over a period ≈ 163 days, following a cone of half opening angle ≈ 21 • and a precessional axis that is inclined ≈ 78 • relative to the line of sight, with jet material moving at a very stable speed ≈ 0.26c along radial ballistic trajectories (Eikenberry et al. 2001). From the width of the emission lines and assuming a conical geometry, the optical jets are very collimated, with half opening angle θ 1 • . The jet precession is thought to be driven by gravitational torques by the compact object on the donor star (slaved disk model ;Roberts 1974; Based on observations collected at the European Southern Observatory, Chile, Program ID 099. D-0666(A,B) den Heuvel et al. 1980;Whitmire & Matese 1980), whose spin axis is misaligned with the orbital plane presumably due to the supernova explosion that originated the compact object 10 5 years ago (Zealey et al. 1980;Lockman et al. 2007). SS 433 is also an eclipsing binary, and the behavior of the X-ray and optical eclipses together with radial velocities of the emission lines reveal that the supercritical disk and its outflows dominate the continuum and line radiation at all wavelengths (for a review of the fantastic properties of SS 433, see Fabrika 2004).
The jets in SS 433 were soon also detected in radio (Hjellming & Johnston 1981), confirming the kinematic precession model derived from the optical lines and establishing the position angle of the precessional axis on sky 98.2 • (Stirling et al. 2002). VLBI observations have regularly imaged the movement of individual jet knots (e.g. Vermeulen et al. 1993;Paragi et al. 1999), and larger scale images of the corkscrew structure created by the precessing jets have established a precise distance d = 5.5 ± 0.2 kpc from the aberration induced by the light travel time effect between the two jets (Blundell & Bowler 2004). Later the jets were detected in X-rays (Watson et al. 1986), where they are seen as emission lines from highly ionized metals (e.g. Fe, Article  Ni, S, Si) which follow the same kinematic precession model as the optical jets and have similarly small opening angle (e.g. Kotani et al. 1996;Marshall et al. 2002). They have been modeled with multi-temperature, optically thin collisional ionization models to estimate properties such as temperature, density and kinetic power. Because the X-ray jets are continuous, emissionline diagnostics depend on the assumptions on the geometry of the outflow, usually taken to be a radially outflowing cone, and may also be affected by photoionization from putative collimated radiation (Brinkmann & Kawai 2000). The behavior of the X-ray jets during eclipse constrain their length 10 12 cm (Marshall et al. 2013).
Because of its edge-on orientation, little is known about the radiation in the beam containing the optical jets, nor about the radiation from the inner parts of the accretion flow in general. The supercritical, geometrically and optically thick disk is thought to reprocess the latter to a large radius (∼ R sp , the spherization radius within which radiation pressure leads to a thick disk, Shakura & Sunyaev 1973), thermally downgrading it to the observed blackbody temperature T ∼ 30000 K (Fabrika 2004). As a result, SS 433 is a relatively faint X-ray source L X ∼ 10 35 −10 36 erg/s to observers on Earth, with most of the received X-ray flux originating from isotropic thermal Bremststrahlung from the Xray jets, without any apparent X-ray accretion disk (Watson et al. 1986). It has been proposed that, if viewed face on, so that the inner portions of the accretion disk/jet funnel were directly visible, SS 433 would appear as an extremely bright X-ray source such as an ULX (Fabrika 2004;Begelman et al. 2006). Recent optical spectroscopy of ULX counterparts have shown strong emission lines akin to those seen in SS 433 and likely associated with supercritical disks (Fabrika et al. 2015), as well as the discovery of an Ultraluminous Supersoft Source (ULS) containing a baryonic relativistic jet seen in a moving Hα line (Liu et al. 2015), so far the only other known object to show such a feature. It remains an open question whether a subpopulation of ULXs could be SS433-like objects, or whether SS 433 is an even rarer phenomenon.
The optical bullets that make up the optical jets are thought to form from the collapse of gas from the continuous X-ray jet through a thermal instability as the jet expands and cools (Davidson & McCray 1980;Brinkmann et al. 1988). If the optical bullets are heated mainly by photoionization (Bodo et al. 1985;Panferov & Fabrika 1993), they can serve as a probe of the collimated radiation. On the other hand, their heating has generally been ascribed to external processes related to interaction with the ambient gas, either by direct collisions (Davidson & McCray 1980;Brown et al. 1991) or through shocks with subsequent heating and ionization (Begelman et al. 1980). Spatially resolving the optical jets could reveal the dominant heating process; just the same, the optical jets have been associated with scales ∼ 10 14 − 10 15 cm ↔ 1 − 10 mas (e.g. Marshall et al. 2013), beyond the reach of current diffraction limited single telescopes.
The only way to spatially resolve the optical jets of SS 433 is through interferometry. In Gravity Collaboration et al. (2017b) (Paper I) we presented the first such observations taken during commissioning of the GRAVITY instrument (Gravity Collaboration et al. 2017a) in July 2016 at the Very Large Telescope Interferometer (VLTI), which works in the near-infrared K band. These observations revealed that the optical jets peak very close to the binary and follow an extended exponential radial emission profile with decay constant ≈ 2 mas, suggestive of a continuous heating process throughout the jet. Here, we present a second set of three GRAVITY observations of SS 433 taken over four nights in July 2017 in which we could observe the change in spatial emission profiles of the jets as the emission lines bright and fade. We also present the first XSHOOTER observations of SS 433, where we use up to twenty pairs of jet lines to constrain the properties of the bullets and the ionizing collimated radiation under the assumption of heating by photoionization suggested by the GRAVITY data. This paper is organized as follows. In Section 2, we summarize the observations and data reduction. The GRAVITY and XSHOOTER data analysis are presented in Sections 3 and 4, respectively. Finally, Section 5 contains the conclusions.

OBSERVATIONS AND DATA REDUCTION
GRAVITY SS 433 (K ≈ 8) was observed with GRAVITY (Gravity Collaboration et al. 2017a) with the Unit Telescopes (UT) on VLTI on three nights over a period of four days in July 2017. Half of the K band light of SS 433 itself was directed to the fringe tracker (FT), which operates at > 1000 Hz to stabilize the fringes in the science channel (SC), allowing coherent integration over detector integration times of 10s in high spectral resolution (R ≈ 4000). The FT operates in low resolution (R ≈ 20) with five channels over the K band. The data were obtained in split polarization mode. The adaptive optics (AO) was performed at visual wavelength using SS 433 itself as the AO guide star (V ≈ 14). Table 1 shows the precessional φ prec and orbital φ orb phases of each observation based on the ephemerides in Eikenberry et al. (2001) and Goranskii et al. (1998), respectively. For more details on the observations and data reduction, including the uv coverage of the observations, we refer to the companion paper on the equatorial outflows (Waisberg et al., sub.). In light of a slightly improved jet model, we also reanalyze the 2016 observation (Paper I), which is also shown in Table 1.
The data were reduced with the standard ESO X-Shooter pipeline (version 2.9.0), which includes de-biasing, flat-fielding, wavelength calibration, sky subtraction, order merging and flux calibration, the latter based on nightly response curves from flux standard stars. The line-modeling software Molecfit Kausch et al. 2015) was used to correct for telluric absorption in the VIS and NIR arms. The performance of Molecfit was found to be at least as good as using a telluric calibrator star, with the additional benefit that manual removal of the many H I, He I and additional lines of a telluric calibrator is not needed. Notes.
(a) Based on the kinematic parameters in Eikenberry et al. (2001). Phase zero is when the eastern/western jet is maximally blueshifted/redshifted. (b) Based on the orbital parameters in Goranskii et al. (1998). Phase zero corresponds to the eclipse center of the accretion disk.
Although our observations were not designed for precise flux calibration, the latter has to be taken into account when comparing emission line strengths across a large wavelength range since SS 433 has a complex and variable continuum. Because the slit widths used are smaller than the seeing, it is necessary to correct for the wavelength-dependent slit losses. We have done this by assuming the typical wavelength dependence for seeing s ∝ λ −0.2 , and using the overlapping regions between the spectral arms to fit for the average seeing, which is in agreement with the estimated value from the acquisition image at the start of each observation.

GRAVITY DATA ANALYSIS
The K band Spectrum Figure 1 shows the GRAVITY spectra of SS 433 for the 2017 observations. There are stationary emission lines (Brγ, He I 2.06 µm, He I 2.11 µm and high order (upper levels 19-24) Pfund lines) as well as emission lines from the baryonic jets. For the analysis of the stationary Brγ line, we refer to the companion paper on the equatorial outflows (Waisberg et al., sub.). In this paper, we focus on the jets. In the 2016 observation (Paper I), the precessional phase was such (φ prec ≈ 0.7) that jet lines from Brγ, Brδ and He I 2.06 µm fell in the K band spectrum (see Figure  A.1). In the 2017 observations presented here, the precessional phase was significantly different (φ prec ≈ 0.9) so that other jet lines were observed: Brβ from the approaching jet and Paα from the receding jet (as well as very weak Brδ and Br lines from the receding jet). Another difference in the latter observations is that there are often two components (knots) to the jet lines. The set of three observations over four nights allows to follow the spatial evolution of the jets; while the knots in Epoch 1 had faded by Epoch 2, the knots in Epochs 2 and 3 partially overlap. The jet redshifts were such that the Paα jets are partially blended with the Brγ stationary line and the Brβ jets with stationary Pfund lines. Figure 2 shows the measured jet redshifts along with the kinematic precession model using the parameters from Eikenberry et al. (2001). They agree within the perturbations caused by nutational motion and random jitter (Fabrika 2004). Figure 3 shows the differential visibility amplitudes and phases for Epoch 3 of the 2017 observations. The only strong jet lines present in the GRAVITY spectrum are Paα for the receding jet and Brβ for the approaching jet. As in Paper I, the jets line emitting regions are more extended than the region emitting nearinfrared continuum as shown by the decrease of the visibility amplitude across the lines, and the receding and approaching jets have opposite differential visibility phases (suggesting that they are on opposite sides) relative to the continuum.

Optical Jet Interferometric Model
We fit the jet lines and differential visibilities simultaneously. For the spectrum, we fit the jet knots as Gaussians. The model parameters are: 1. The jet inclination i, which determines the redshift via where β = 0.26 is the jet speed in units of c and γ = 1 . Because we can measure the redshift very precisely for each jet line, we allow for different inclinations of the receding and approaching jets, as well as the different components. We assume a constant velocity since it is very stable (Eikenberry et al. 2001), and the possible small variations in velocity are absorbed in the inclination. 2. The FWHM in km/s, which is assumed to be the same for all the jet lines for all components; 3. The strength of each jet line relative to the normalized continuum.
To fit the visibilities, we model the jets as a 1D structure with a radial emission profile, since they are very collimated (opening angle 1 • ) from the width of the optical emission lines . The model parameters are: 1. The position angle (PA) of jet in the sky plane. It is assumed to be the same for all the jet components in a given epoch because it cannot be measured nearly as precisely as the inclination (the typical error is a few degrees); 2. The radial emission profile is controlled by three parameters θ, α and r 0 : where r is the distance from the center. This model is similar to the one used in Paper I, except for the additional parameter α which allows for more general shapes besides an exponential (α = 1), from Gaussian-like (α > 1) to steeper profiles Epoch 1 (2017-07-07) Epoch 2 (2017-07-09) Epoch 3 (2017-07-10) Jet redshifts (mostly approaching (eastern) jet in blue; mostly receding (western) jet in red) and position angle on sky (green) of the optical jets for each observation as a function of precessional phase. The curves were made using the kinematic model parameters in Eikenberry et al. (2001) and the precessional axis position angle in Stirling et al. (2002). The PA of the optical jets as derived from the GRAVITY data agrees well with the prediction from the radio jets. We plot both the 2016 (φ prec ≈ 0.7) as well as the 2017 (φ prec ≈ 0.9) GRAVITY observations.
(α < 1). The parameters θ (which together with α controls the shape of the profile) and r 0 (the inner edge where the emission starts) are also debiased from the projection effect i.e. they are already divided by sin(i). The errors are estimated from the scatter in line-free regions. We fit for the spectrum and the differential visibilities simultaneously; however, because the former is sensitive to telluric correction and has very small statistical errorbars, we increase the flux errorbars by a factor of two. We find that this scaling led to a comparable reduced χ 2 between flux and visibilities in all observations. Moreover, because of the blending with the Brγ stationary emission line, which also produces differential visibility signatures, it is necessary to perform simultaneous fits for the Brγ line and the jets. For the model and results for the Brγ stationary line we refer to the companion paper on the equatorial outflows (Waisberg et al., sub.). The results presented here correspond to the "outflow" model in the companion paper (which we favor over the "disk" model). The differential visibilities are computed with respect to the best fit continuum model (see companion paper Waisberg et al., sub.). The 2D Fourier transform of the jet model detailed above is computed semi-analytically. The total differential visibility at a given spectral channel is then where V c is the continuum visibility, and f i and V i are the flux ratios relative to the continuum and visibilities for each component i (different jet knots, Brγ stationary line). The fits are done through non-linear least squares minimization with the Levenberg-Marquardt method through the python package LMFIT 1 . The quoted errors correspond to the 1-σ errors from the least squares fit i.e. the estimated derivatives around the optimal solution (scaled by χ 2 red ). We caution, however, that true uncertainties are dominated by (i) degeneracies between the many parameters, which create a complicated multi-dimensional χ 2 map ; (ii) systematic errors from the continuum model; (iii) the assumption of our simple "geometric" models, which cannot capture all the complexities involved. We note, in addition, that in Epoch 1 there is very severe spectral blending of the Paα jet lines with the stationary Brγ line, so that its results are less robust.

Results
Table 2 shows the model fit results for the jet emission lines in all epochs. The data and best fits are shown in Figure 3 for Epoch  There is substantial blending of the Paα lines from the receding jet with the stationary Brγ line, which also has differential visibility signatures across it. The fits are done to all lines simultaneously, but we plot only the model for jet lines for clarity. The baselines projected lengths and position angles on sky are indicated in the labels.
As noted before, there are two components ("knots") in the jet lines in the 2017 observations. The need for two components in 2017 is clear from the much broader lines, as well as from the interferometric data (differential visibility amplitudes and phases which are not aligned). Two components are clearly distinguishable in all epochs for the Paα line from the receding jet but only on Epoch 3 for the Brβ line from the approaching jet. One of the components is compact (emission centroid 2 mas), with an emission profile similar to the one in the 2016 observation, and associated with the most recent or current jet ejection. The Article number, page 5 of 19 A&A proofs: manuscript no. SS433_Jets_2cols other component is significantly more extended (emission centroid 6 mas), associated with a previous jet ejection. This interpretation agrees with the corresponding redshifts from the precession model in Epochs 1 and 2. Figure 4 shows the collection of radial emission profiles of the optical jets for all observations. The compact knots have an exponential-like profile, whereas the extended knots have more rotund shapes. We find that for the more compact jet components, the emission peaks substantially close to the binary system, and is more compact than previous estimates from optical spectroscopy monitoring, e.g.  derived an exponential decay of the jet emission with falloff distance 6.7 × 10 14 cm ≈ 8.2 mas for the Hα jet line. The elongated emission profiles that we measure are strongly suggestive of a continuous heating mechanism along the entire jet, which is naturally accomplished by photoionization by the collimated radiation from the inner regions close to the compact object. As in Paper I, we attempted to fit the jets with more localized emission profiles (such as a point source or Gaussian), but they are inconsistent with the data: an elongated structure for the jets is strongly preferred.
Furthermore, the 2017 observations allow to probe the spatial evolution of the jet profiles from night to night. The jet lines in Epoch 1 have clearly disappeared by Epoch 2 two nights later ( Figure 1). However, the jets lines in Epochs 2 and 3 (which are separated by one night only) partly overlap. The extended component in Epoch 3 could correspond to the compact component in Epoch 2 after it has travelled ≈ 8 mas/day ×1 day ×sin(60 • ) ≈ 7 mas on the sky plane, whereas a new compact component in Epoch 3 clearly appeared at the same redshift where the extended component in Epoch 2 was located. This is the first time the movement of the optical bullets has been spatially resolved, although longer observations during a single night would be needed to trace the motion of an individual component unambiguously.

XSHOOTER DATA ANALYSIS
We fit the jet lines with Gaussian or Lorentzian profiles (the jet lines are usually better fit by a Gaussian, but sometimes a Lorentzian profile is clearly preferred) to estimate their central wavelength, FWHM and total intensity. All line fits were performed with the task splot in IRAF 2 (Tody 1986(Tody , 1993. Whenever jet lines were blended with other jet or stationary lines, deblending was used. For all line identifications and adopted wavelengths, the "Atomic Line List v2.05b21" (van Hoof, P.) was used. The noise for the spectral fits was estimated from the scatter in continuum regions near each line. A particularly important constraint for photoionization models of the optical bullets is the absence of ionized helium in the jets due to the lack of the He II 4686Å line. We estimate upper limits on its flux from the known location of where the line would appear based on the measured redshifts and the FWHM from the other jet lines. Figure 5 shows the spectrum for the first XSHOOTER observation of SS 433 (Epoch X1), which contains the largest number of jet lines of all epochs. Emission lines from the jet are  There are two classes of profiles: compact (solid lines), associated with recent or current knots and which peak close to the central binary, and extended (dashed lines), associated with older knots. Bottom: Emission profiles for selected emission lines calculated from a Cloudy photoionization model using the best fit parameters to the emission line fluxes from Epoch X1 of the XSHOOTER observations. The elongated profiles resolved by GRAVITY suggest photoionization by the collimated radiation as the heating mechanism. The steeper profiles relative to the model are very likely due to screening effects (large area filling factor by the bullets in the jets).
shown in blue and red for the eastern and western jets, respectively (at this epoch, the eastern jet, which is approaching most of the time, is receding). Throughout all observations we identify up to twenty pairs of lines: Hα through H ; Paα through Pa9; Brβ through Br12; and five lines of He I (5875Å, 6678Å, 7065Å, 1.083 µm, 2.056µm). The strongest hydrogen lines (Hα, Hβ, Paα, Brβ) and He I 1.083µm often show one or two additional components (shown in cyan and orange in Figure 5), with redshift corresponding to previous ejections according to the precession model (and which we associate with older jet knots), whereas most emission lines show only one component from the most recent or current ejection (this behavior matches what is seen in the GRAVITY observations). The average redshift for the more recent jet knots in each epoch is shown in Figure 2 and agrees with the kinematic precession model within the random jitter. The strength and number of jet emission lines varies significantly between the five XSHOOTER epochs.

Cloudy Models for Photoionized Optical Bullets
The spatial emission profiles for the optical jets resolved by GRAVITY are strongly suggestive of photoionization as the main heating mechanism of the optical jets. Therefore, we proceeded to fit the jet emission line fluxes measured from XSHOOTER with photoionization models.
All line emission in SS 433 originates from dense gas; a clear lower limit to the density n e 10 7 cm −3 follows from the absence of any forbidden lines in the spectrum (Osterbrock & Ferland 2006, p.60). The true density in the optical bullets is estimated to be ∼ 10 13 cm −3 , and the optical jets take the form of dense bullets distributed throughout the jets with a low volume filling factor ∼ 10 −6 (Fabrika 2004). At such densities, collisional ionization from excited levels and radiative transfer effects become important, causing strong deviations from Case B recombination; still, such densities are not high enough for LTE to apply to all levels, and full collisional radiative models must be used (Ferland et al. 2013(Ferland et al. , 2017. A clear sign that collisional and radiative transfer effects are important in the optical bullets are the rather flat or even inverted Hα/Hβ jet line ratios (after extinction correction).
We model the jet emission lines with the photoionization code Cloudy (Ferland et al. 2013(Ferland et al. , 2017 version c17.01. The ionizing radiation is taken to be a single blackbody. We construct a grid in bullet density, bullet size, blackbody temperature and blackbody intensity. In addition, because SS 433 is known to have appreciable extinction, we expand the grid in A V , assuming an extinction law A λ A(V) following Cardelli et al. (1989) and R V = 3.1. We attempted to fit for the He abundance, since it is likely that the donor star in SS 433 is at an advanced evolutionary stage, but it was unconstrained in our fits; therefore, we fixed a standard He abundance (0.098 by number). The grid parameters and intervals are shown in Table 3. The solutions are iterated until convergence of the optical depths. Each point in the grid gives a model intensity for each line (erg/s/cm 2 ). The total number of bullets is then computed to minimize the χ 2 with respect to the measured line intensities. The bullets are assumed to be identical in both jets, but the number of bullets (and therefore the kinetic power) can be different for the receding and approaching jets. We corrected the line intensities in each jet for Doppler boosting (1 + z) 3 using the average measured redshift.
We use the lines from the most recent jet knots (which contain a much larger number of lines) for the fits. The GRAVITY observations revealed that recent ejections have a typical exponential intensity profile as shown in Figure 4. We integrate our Cloudy models over five radial points covering a factor of 100 in intensity. The reported best fit intensity is the one at the base of the optical jet.
The statistical line flux errors estimated from the spectral fits are very small. The actual errors are dominated by systematic effects in the flux calibration, telluric correction, contamination by weak lines, line blending and model uncertainties (such as collisional cross sections). We estimate a 25% error in flux for all the lines, which results in χ 2 red ∼ 1 for the best fit grid points. The parameter uncertainties are estimated by considering all the models which fall within ∆χ 2 = p of the best χ 2 (Lampton et al. 1976), where p is the number of model parameters.
Article number, page 7 of 19 A&A proofs: manuscript no. SS433_Jets_2cols 1. The volume filling factor of the optical bullets: where V bullets is the total volume in the bullets (calculated from their size and number) and V jet = πψ 2 jet l 3 jet 3 , where ψ jet is the half-opening angle and l jet is the length of the jet; 2. The kinetic power of the optical jets: where v jet = 0.26c, M is the total mass in the bullets (calculated from their density, size and number). 3. The total line luminosity in the optical bullets L bullets . This includes not only the recombination lines within XSHOOTER, but all of the strongest hydrogen and helium lines in the full spectrum output from Cloudy. 4. The total ionizing luminosity within the beam containing the optical bullets: where I(r 0 ) is the intensity at the base of the beam. The total luminosity in the collimated radiation could be higher if it is broader than the beam containing optical bullets, which is likely the case from the presence of older jet knots that can keep radiating for a few days (Panferov & Fabrika 1993). 5. The luminosity inferred by an observer (assuming isotropy) whose line of sight is within the collimated beam: L face−on = L beam 4π π(2ψ jet ) 2 = 2I(r 0 )r 2 0 (7) For these estimates, we assumed ψ jet = 1 • and r 0 = 0.4 mas and l jet = 2 mas from the GRAVITY observations.
Our model does not take into account screening of radiation by the bullets, which plays a role because the jets are very compact (area filling factors are high). To ensure self-consistency in the energetics, we only accept solutions with enough total luminosity in the beam to power the bullets (i.e. L beam > L bullets ). We also note that we assume the bullet radiation is isotropic and identical between the two jets (the only difference we allow between the two jets is in the total number of bullets), whereas differences between the two jets with precessional phase have been interpreted in terms of anisotropic radiation Fabrika 2004).

Results
We note that several related calculations to estimate the properties of the optical jets can be found in previous papers (e.g. Davidson & McCray 1980;Begelman et al. 1980;Bodo et al. 1985;Brown et al. 1991;Panferov & Fabrika 1993), primarily based on the Hα luminosity. These calculations have several uncertainties: unknown heating mechanism/line emissivity, degeneracy between density and volume filling factor, unknown extinction, optical depth effects, unknown spatial emission profile of the optical jets (e.g. jet size). The calculations presented here are based on the first optical interferometric measurements that have spatially resolved the optical jets, and on a large number of jet line species from XSHOOTER spectra. Table 4 shows the model fit results for the five XSHOOTER epochs. Figure 6 shows the measured and model line fluxes for the first XSHOOTER epoch (corresponding plots for the other epochs are shown in Appendix B).

The Bullet Properties
The results confirm that the optical bullets in SS 433 are very dense n H ∼ 10 13 cm −3 and have a size R ≈ 10 6 − 10 7 cm. Even though they are optically thin to electron scattering (τ T ∼ 0.01), many emission lines are optically thick. For the best fit model in Epoch X1, τ Hα ≈ τ Pα ≈ 20, τ Hβ ≈ 4, τ Brγ ≈ 2, τ Br10 ≈ 0.3, τ HeI6578 ≈ 5 and τ HeI1.083 ≈ 0.6, so that optical depth effects in the line emission indeed have to be taken into account. The cooling in the bullets is dominated by hydrogen and helium recombination lines (∼ 80%), with total line luminosity ≈ 2 × 10 37 erg/s. Only around 15% of such luminosity is in lines within the XSHOOTER spectrum, with the strongest lines being the Lyman series in the UV.
A V -6.4 − 6.8 6.2 − 6.8 6.6 − 7.6 6.6 − 7.4 6.6 − 7.6 We constrain the kinetic power in the optical bullets to ∼ 2 − 20 × 10 38 erg/s. This agrees with estimates from emission line modeling in the X-rays, which vary from from ∼ 3 × 10 38 to 5 × 10 39 erg/s (e.g. Marshall et al. 2002;Brinkmann et al. 2005;Medvedev & Fabrika 2010). The collapse of the continuous Xray jets into optical bullets must therefore be an efficient process, both in terms of mass and kinetic energy. Previous estimates from optical spectroscopy give L kin ∼ 10 39 erg/s (e.g. .

The Collimated Radiation
The results constrain the collimated radiation to be relatively soft T ∼ 3 − 4 × 10 4 K and the total luminosity in the 1 • beam containing the optical bullets to L beam ≈ 2 − 6 × 10 37 erg/s. Lower luminosities cannot power the jet emission lines, whereas higher luminosities and temperatures cause too intense heating and ionize helium too much. For an observer looking face on at the col-limated radiation and assuming isotropy, the inferred luminosity would be ≈ 7 − 20 × 10 40 erg/s i.e. SS 433 would appear as an extremely bright UV source. This luminosity is higher than the ∼ few × 10 39 erg/s inferred from the SED (Wagner 1986, although this is a rather uncertain number -see below), suggesting that indeed there is collimated radiation in SS 433 (not only thermal downgrading at low latitudes). The total luminosity in collimated radiation depends on its opening angle, and would be ∼ 2 × 10 39 − 10 41 erg/s for angles 10 • − 50 • . From modeling of the optical filaments in the W50 nebula and assuming photoionization by collimated radiation, Fabrika & Sholukhova (2008) estimated an ionizing luminosity ∼ 10 40 erg/s in an opening angle ∼ 50 • .
Figure 4 (bottom) shows the spatial emission profile for selected emission lines for the best fit Cloudy photoionization model to Epoch X1. It assumes the bullets are distributed homogeneously over the jet with constant properties (density, size), and shows the normalized line emissivity as a function of distance r in the jet. It resembles the spatial profiles directly resolved by GRAVITY, and confirms that the different H I and He I lines have similar emission regions. The elongated profile results from the interplay between the radiation intensity and gas density, and is a good consistency check for our photoionization models. A caveat is that we do not consider screening of the intensity along the jet by the bullets (large area filling factors by the bullets in the jet), which is likely the reason for the steeper measured profiles relative to the model.
Because of the possibility that SS 433 could be an ULX for an observer looking at it face on, we attempted to constrain the X-ray luminosity in the beam containing the optical bullets. We did this through a perturbative approach: we selected the best fit model in Epoch X1 from the analysis above, and added a second X-ray component with varying intensity, and checked at which X-ray intensity the optical emission lines are substantially affected and in clear violation of the data (e.g. by producing a too strong He II 4868 Å line). Because the main effect comes from the soft X-ray photons capable of ionizing H, He and He+, the constraint depends on the relative contribution of soft X-ray photons to the total X-ray luminosity. To this end, we repeated the procedure for two different SEDs: a very soft one corresponding to a Supersoft Ultraluminous (SSUL) source, in the form of a blackbody kT = 0.14 keV, and a harder one corresponding to a Hard Ultraluminous (HL) source, in the form of a blackbody kT = 0.27 keV plus a strong hard component (Kaaret et al. 2017). The corresponding SEDs and limits are shown in Figure  7.
In the case of the SSUL SED, we constrain the X-ray luminosity to be 10 −3 of the UV component i.e. 5 × 10 34 erg/s in the 1 • beam containing the optical bullets or 10 38 erg/s for a face on observer. In this case, SS 433 would not be an ULX. For the hard UL SED, the X-ray luminosity could be much larger, up to 10 −1 of the UV component i.e. 5 × 10 36 erg/s in the 1 • beam containing the optical bullets or 10 40 erg/s for a face on observer. We conclude that SS 433 could be an ULX, as long as its X-ray spectrum is dominated by hard X-rays. Just the same, face on ULXs (with a clear view through the funnel) are generally expected to have dominantly X-ray spectra (Kaaret et al. 2017), whereas our results suggest that SS 433 is UVdominated even in the collimated beam. This might mean that thermal downgrading happens already in the funnel of SS 433 (Begelman et al. 2006). The most promising way to find face on SS 433-like objects in other galaxies might be to look for very bright and variable (due to jet precession) UV sources.
The X-ray luminosity of SS 433 has also been constrained from a putative reflection component ∼ 10 35 erg/s in the hard Xray spectrum 10 keV of SS 433 (Medvedev & Fabrika 2010;Middleton et al. 2018). Middleton et al. (2018) estimates that the intrinsic X-ray luminosity is L X 10 38 erg/s. This is consistent with our upper limits as long as the X-ray radiation is slightly less collimated ( 5 • ) than the beam containing the optical bullets, which is almost certainly the case from the presence of older jet knots that keep radiating for a few days (Panferov & Fabrika 1993).

The Extinction Towards SS 433
The extinction towards SS 433 is known to be large from the very red continuum spectrum but its exact value is difficult to assess. Galactic dust extinction maps give A V = 7.8 towards the direction of SS 433 (Schlegel et al. 1998;Perez M. & Blundell 2010), but that is an upper limit to the total integrated line of sight extinction. More recent 3D dust maps rather give A V = 5.7 ± 0.1 at d = 5.5 kpc towards SS 433 (Green et al. 2018). Strong Diffuse Interstellar Bands (DIBs) are also present in the spectrum of SS 433 (as noted in previous work e.g. Murdin et al. 1980;Margon 1984). The strength of DIBs has been shown to be correlated with interstellar extinction (e.g. Herbig 1975;Friedman et al. 2011;Kos & Zwitter 2013), although with substantial scatter. We measured the EWs of several DIBs (5780, 5797, 5850, 6196, 6202, 6270, 6283, 6379, 6613, 6660 Å) in our XSHOOTER spectra and used the correlations in Friedman et al. (2011) and Kos & Zwitter (2013) to estimate A V = 5.1 ± 1.0 and A V = 5.6 ± 2.1, respectively, where the uncertainties are the 1σ scatter between the different DIBs. Therefore, there is evidence that the interstellar extinction towards SS 433 may be A V 6.0.
On the other hand, A V 7.8 − 8.0 has also been estimated from fitting the spectral energy distribution with a single reddened blackbody (e.g. Murdin et al. 1980;Wagner 1986;Dolan et al. 1997). This approach, however, suffers from several problems (i) in the Rayleigh-Jeans range, the temperature of the blackbody and the extinction are very strongly correlated; (ii) there are numerous and very strong emission lines in SS 433; (iii) it is not clear whether the supercritical disk should look like a single blackbody e.g. the temperature seems to change with precession phase (Wagner 1986). Our XSHOOTER observations are not optimized for SED continuum fitting, but we confirm the strong degeneracy between blackbody temperature and extinction. T 20000 K (a reasonable expectation from the presence of He II 4868 Å stationary line) requires A V 7.5.
Our modeling of the jet lines yields A V ≈ 6.7 ± 0.1, which is intermediate between the lower values inferred from 3D dust maps and DIBs and those estimated from the very reddened SED. We suggest that there may be substantial A V 1 and structured circumstellar extinction in SS 433, affecting the equatorial part of the system more than the optical jets. It may be caused by dust forming from the equatorial outflows seen in radio (Paragi et al. 1999;Blundell et al. 2001) and near-infrared stationary emission lines (see companion paper on equatorial outflows, Waisberg et al., sub.). Mid-infrared observations of SS 433 show evidence of dust from excess emission at λ 20µm (Fuchs et al. 2006). We speculate that mid-infrared interferometric observations with VLTI+MATISSE (Lopez et al. 2014) might resolve an extended dust torus in SS 433.

CONCLUSIONS
We presented a second set of GRAVITY observations of SS 433 after Paper I, as well as the first XSHOOTER observations of this object, focusing on the optical jets. We summarize our main conclusions from the GRAVITY observations as follows: 1. The optical jets have elongated, exponential-like spatial emission profiles, suggestive of a continuous heating process throughout the entire jet; we argue for photoionization by collimated radiation; 2. We have spatially resolved the movement of the optical bullets for the first time, finding more extended jet knots corresponding to previous jet ejections.
Using the up to twenty simultaneous pairs of measured jet line fluxes in the XSHOOTER observations, we have constrained properties of the optical bullets and the putative ionizing radiation with Cloudy photoionization models: 1. The optical bullets are dense ∼ 10 13 cm −3 and have a size ∼ 10 6 − 10 7 cm, from which optical depth effects in the jet emission lines are important; 2. The kinetic power of the optical jets is ∼ 2 − 20 × 10 38 erg/s; 3. The beamed radiation is dominantly UV with a luminosity ≈ 2 − 6 × 10 37 erg/s in the 1 • beam containing the optical bullets. An observer looking directly at the beam would infer an isotropic luminosity ≈ 7−20×10 40 erg/s i.e. SS 433 would appear as an extremely bright UV source; 4. In the photoionization picture, SS 433 could still be an ULX with a face on observer inferring L X 10 40 erg/s, as long as the X-ray SED is dominantly hard, since soft X-ray photons destroy the H/He/He+ ionization balance in the optical bullets; 5. We constrain the extinction in the optical jets A V = 6.7 ± 0.1 and suggest there is substantial and structured circumstellar extinction in this object.   Fig. 7. Potential SS 433 SED for a face on observer who sees the beam radiation as the optical bullets see it. The solid blue curve shows the UV component from the best fit Cloudy photoionization model to the line intensities in Epoch X1 of XSHOOTER. The orange and red lines show the upper limit to a possible X-ray component for a Supersoft Ultraluminous (SSUL) and Hard Ultraluminous (Hard UL) type SEDs, respectively. The dashed lines show the ground state photoionization cross sections for H, He and He+. For SS 433 to look like an ULX for the face on observer, its X-ray spectrum must be significantly hard, since soft X-ray photons break the H/He/He+ ionization balance in the optical bullets.

Appendix A: GRAVITY: Full Data and Model Fits
Here we show as in Figure 3 the full K band spectrum and interferometric data, as well as the best fit jet models, for the remaining GRAVITY observations. Compact and extended jet knots are shown in colored full and dashed lines, respectively. Stationary emission lines are marked in green. In the 2017 observations, in which there is strong blending between the jets and the Brγ stationary line, a combined fit is done, but we show only the jet model for clarity. For the Brγ stationary line model and results, we refer to the companion paper on the equatorial outflows (Waisberg et al., sub.).

Appendix B: XSHOOTER: Full Data and Model Fits
Here we show the data and best fits for the remaining XSHOOTER epochs as in Figure 6.