© The Authors 2023

1. Introduction

Early-type stars are very hot and their radiation fields can launch powerful particle winds (Lamers & Cassinelli 1999). Such winds quickly reach supersonic velocities and accelerate to terminal velocities in the range (2 − 4)×10³ km s⁻¹ (Abbott 1978; Muijres et al. 2012). When two massive stars with powerful winds form a binary system, the winds collide producing shocks separated by a contact discontinuity from where matter is evacuated (e.g., Stevens et al. 1992). A reverse shock moves in the wind of each star. When such shocks are adiabatic, they can accelerate suprathermal particles up to relativistic energies (Eichler & Usov 1993; Pittard et al. 2020. These particles, in turn, cool mainly by synchrotron radiation and inverse Compton upscattering of stellar photons, emitting nonthermal radiation (Eichler & Usov 1993; Benaglia & Romero 2003; Reimer et al. 2006; De Becker 2007; Reitberger et al. 2014; del Palacio et al. 2016; Pittard et al. 2021). Proton acceleration can also lead to gamma-ray emission through pp collisions and the subsequent π⁰ decays (e.g., Balbo & Walter 2017; Grimaldo et al. 2019).

The actual fraction of particle-accelerating colliding-wind binaries (PACWBs) among massive colliding wind binaries (CWBs) is not well known. De Becker & Raucq (2013) list 43 confirmed cases, mostly detected at radio wavelengths. These authors mention several other candidates, and new sources have been found since the publication of this latter work (e.g., Benaglia et al. 2015; del Palacio et al. 2016). The total kinetic power of these systems ranges from ∼10³⁴ to more than 10³⁷ erg s⁻¹. The most extreme cases are WR89, WR98, and WR140, with powers of between 6 and 8 times 10³⁷ erg s⁻¹. Less than 10⁻⁷ of this power is finally radiated through synchrotron radio emission. The most luminous nonthermal radio-emitting CWB is WR140, with a total radio luminosity of ∼2.6 × 10³⁰ erg s⁻¹.

Contrary to the radio emission, high-energy radiation has been more difficult to detect in CWBs. At X-rays, the thermal component usually dominates and hinders the detection of nonthermal components. In the gamma-ray domain, only two systems have been detected so far: η Carinae and WR11. The latter is the nearest known CWB. At d ∼ 340 pc, it shows a gamma-ray luminosity in the Fermi-LAT energy range of L_γ = (3.7 ± 0.7)×10³¹ erg s⁻¹. This luminosity amounts to ∼6 × 10⁻⁶ of the total wind kinetic power (Pshirkov 2016). Similar fractions for other, more distant PACWBs yield fluxes that are undetectable with the currently available instrumentation. The notable exception is the mentioned η Carinae.

η Carinae is a heavily obscured and peculiar object. The system includes a luminous blue variable (LBV) star of about 90 solar masses and a secondary Wolf-Rayet (WR) star of ∼30 solar masses. η Carinae is the most luminous binary in the Galaxy, with a bolometric luminosity of about 5 × 10⁶ L_⊙. The mass-loss rate of the primary is extremely high, reaching up to 10⁻³ M_⊙ yr⁻¹. The binary was detected in hard X-rays by INTEGRAL (Leyder et al. 2008) and Suzaku (Okazaki et al. 2008), suggesting the presence of relativistic electrons in the system. AGILE detected gamma rays from η Carinae for the first time (Tavani et al. 2009). The system was subsequently detected by Fermi (Abdo et al. 2010) with a luminosity of ∼10³⁴ erg s⁻¹. The observations reveal the presence of a hard component in the spectrum around periastron, which disappears near apastron. Such a component has been explained through the decay of π⁰ produced by relativistic protons interacting with the dense stellar wind (Farnier et al. 2011). There is a clear variability with the orbital phase. Different behaviors are observed at low (0.3 − 10 GeV) and high (> 10 GeV) gamma-ray energies. The low-energy component is likely produced by inverse Compton scattering of stellar photons (Balbo & Walter 2017).

The case of η Carinae suggests that super-Eddington systems might be particularly powerful PACWBs. When a compact object such as a black hole accretes with rates that exceed the Eddington rate, the radiation pressure on the surface of the disk will overcome the gravitational attraction and matter will be expelled from the surface of the disk in the form of a strong wind. Such winds can rival and even surpass those of the most luminous CWBs in terms of kinetic power. When the donor star is a hot early-type star also endowed with a wind, a supercritical colliding wind binary (SCWB) can be formed. Such systems should have strong shocks and are potential particle accelerators and nonthermal emitters.

In our Galaxy, there are some examples of black hole X-ray binaries with disks that launch strong outflows. Two examples are GRS 1915+105 (Mirabel & Rodríguez 1994; Neilsen & Lee 2009) and V404 Cygni (Muñoz-Darias et al. 2016; Tetarenko et al. 2017). However, the donor star in both of these systems is a low-mass star. Another well-known supercritical source is the Galactic microquasar SS433, which is a confirmed nonthermal emitter and might be a possible example of a SCWB in our Galaxy (see Fabrika 2004 for an extensive review). Many ultraluminous X-ray sources (ULXs) detected in nearby galaxies might also belong to this category of sources.

In this paper, we explore the CWB scenario where one of the winds is launched by a supercritical disk around a black hole. We start by characterizing the disk model and the radiation fields it produces (Sects. 2.1 and 2.2). We then investigate the motion of particles under the radiation pressure in such fields (Section 2.3). This allows us to get reasonable estimates of the terminal velocities expected for the matter ejected in the direction of the companion star. We then proceed to study the wind interactions, shock adiabaticity, and other relevant issues for particle acceleration in Sect. 3. This is followed by estimates of energy losses for accelerated particles, particle distributions, and calculations of the nonthermal output (Sect. 4). In Sect. 5 we present results for some specific models, with different choices of the accretor mass and the accretion power. The donor star is supposed to be a hot O.5V with a temperature of 41 500 K and a kinetic power of a few times 10³⁷ erg s⁻¹. We finally apply our model to the extragalactic binary system NGC 4190 ULX 1. After a discussion (Sect. 7), we close with a summary and our conclusions.

2. The accretion disk and its wind

We assume that the X-ray binary is composed of a Population I star and a nonrotating stellar mass black hole (BH) in a close orbit.

The orbital semi-axis a, the stellar radius, and the mass ratio of the system, q = M_*/M_BH, satisfy (Eggleton 1983):

$$\begin{array}{c}\hfill R_{\mathrm{lob}}^{\ast }=\frac{a0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})},\end{array}$$(1)

where M_* is the mass of the star and M_BH the mass of the BH. Hence, the star overflows its Roche lobe $R_{\text{lob}}^{*}$, transfers mass to the BH through the Lagrange point, and an accretion disk is formed due to the angular momentum of the system.

In this section, we describe the semi-analytical models we use to study the accretion disk, the spatial distribution of the radiation fields produced by the disk, and the wind ejected from its surface. We assume a Newtonian potential for the gravity field, because we are interested in weak-field processes.

2.1. Accretion disk

We adopt cylindrical coordinates with axial symmetry along the z-axis, neglect the self-gravity of the disk gas, and consider a nonmagnetized disk with a super-Eddington accretion rate at the outer part of the disk, ṁ_input = Ṁ_input/Ṁ_Edd ≫ 1, where Ṁ_input is the input of mass per time unit in the accretion disk. The Eddington rate is given by

$$\begin{array}{c}\hfill M_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{\eta c^2}\approx 2.2\times {10}^{-8}{M}_{\mathrm{BH}}{\mathrm{yr}}^{-1}=1.4\times {10}^{18}\frac{M_{\mathrm{BH}}}{M_{\odot }}{\mathrm{g}\mathrm{s}}^{-1},\end{array}$$(2)

with L_Edd the Eddington luminosity¹, η ≈ 0.1 the accretion efficiency, and c the speed of light.

The critical or spherization radius, given by

$$\begin{array}{c}\hfill r_{\mathrm{crit}}\sim 40{\dot{m}}_{\mathrm{input}}{r}_{\mathrm{g}},\end{array}$$(3)

separates the disk in two regions: a standard outer disk (Shakura & Sunyaev 1973) and a radiation-dominated inner disk with advection (Fukue 2004). In relation (3), r_g = GM_BH/c² is the gravitational radius of the BH, with G the gravitational constant. In the disk model, the advection is parameterized as a fraction f of the viscous heating, Q_adv = fQ_vis, and the disk becomes geometrically thick in the inner region, where the ejection of winds by the radiation force helps to regulate the mass-accretion rate onto the BH (Ṁ_acc) at the Eddington rate².

As the disk is optically thick, we assume that it radiates locally as a blackbody. The radiation intensity of a plasma element in the comoving frame of the outer and inner disk, at a radius r_d measured on the equatorial plane, is

$$\begin{array}{c}\hfill I_0=\frac{1}{\pi }\sigma T_{\mathrm{eff}}^4=\{\begin{array}{c}\frac{1}{\pi }\frac{3G{M}_{\mathrm{BH}}{\dot{M}}_{\mathrm{input}}}{8\pi r_{\mathrm{d}}^3}{f}_{\mathrm{in}},r_{\mathrm{d}}>r_{\mathrm{crit}}\hfill \\ \hfill \\ \frac{1}{\pi }\frac{3}{4}\sqrt{c_3}\frac{L_{\mathrm{Edd}}}{4\pi r_{\mathrm{d}}^2},r_{\mathrm{d}}\le r_{\mathrm{crit}},\hfill \end{array}\end{array}$$(4)

where $\sqrt{c_3}=H/r_{\mathrm{d}}=tan\delta$, with H the scale height of the disk, δ the disk opening angle, and f_in = 1 − r_in/r_d ≈ 1 (as r_d > r_crit, then r_d ≫ r_in). Here, c₃ (along with c₁ and c₂ used in the following section) is a coefficient that depends on the advection parameter, the adiabatic index of the gas γ, and the viscosity α (see Appendix in Fukue 2004). We adopt a disk with f = 0.5 and α = 0.5; that is, we assume equipartition between advection and viscous heating. The index γ = 4/3 corresponds to a radiation-dominated gas in the inner disk. These values lead to a disk-opening angle of δ = 30°.

2.2. Radiation fields

The wind launched from the radiation-dominated region of the disk will be determined by the radiation forces acting upon the particles on the disk surface and along their subsequent trajectories. These forces will have contributions from different parts of the disk in relative motion with respect to the particles. Some radiation will be blueshifted and some will be redshifted, resulting in differential azimuthal forces onto the particles and then transferring angular momentum from the disk to the wind.

In order to obtain the radiative contribution of each plasma element 𝒬 = (r_d, ϕ_d, H) of the disk surface, at any point 𝒫 = (r, ϕ, z) above or below the disk, we make a transformation of the intensity between the inertial and comoving reference frames (see Fig. 1). Azimuthal symmetry allows us to perform the calculations for any constant value of ϕ; therefore, we do it in the rz plane (ϕ = 0).

Fig. 1. Geometry of the present disk model. The radiation fields are calculated in the rz plane, where ϕ = 0. Here, 𝒬 is the position of the plasma element of the disk and 𝒫 the point of calculation on the rz plane. The scale height of the disk is H, and D is the distance between 𝒬 and 𝒫. The short arrow is the direction cosine jμ. This figure is adapted from Watarai & Fukue (1999).

The relativistic Doppler factor 𝒟 provides the transformation between the reference frames (McKinley 1980):

$$\begin{array}{c}\hfill I={\mathcal{D}}^4{I}_0=\frac{I_0}{{(1+z_{\mathrm{red}})}^4},\end{array}$$(5)

where z_red is the redshift factor given by (Watarai & Fukue 1999)

$$\begin{array}{c}\hfill z_{\mathrm{red}}=-\frac{(rcos{\varphi }_{\mathrm{d}}-r_{\mathrm{d}})v_r-(rsin{\varphi }_{\mathrm{d}})v_{\varphi }+(z-H)v_r{c}_3}{\mathit{cD}}.\end{array}$$(6)

Here, D is the distance between 𝒫 and 𝒬, v_ϕ = c₂v_K is the azimuthal velocity and v_r = −c₁αv_K is the radial velocity, with $v_{\mathrm{K}}=\sqrt{G{M}_{\mathrm{BH}}/r_{\mathrm{d}}}$ the Keplerian velocity. We note that we only consider the inner part of the disk for these calculations, because the intensity decays with $r_{\mathrm{d}}^{-3}$.

The radiation-field tensor is given by (Rybicki & Lightman 1986)

$$\begin{array}{c}\hfill R^{\mu \nu }=\left(\begin{array}{cc}E& \frac{1}{c}{F}^{\alpha }\\ \frac{1}{c}{F}^{\alpha }& P^{\alpha \beta }\end{array}\right)=\frac{1}{c}{\displaystyle \int }I{j}^{\mu }{j}^{\nu }\mathrm{d}\mathrm{\Omega }.\end{array}$$(7)

This is a symmetric tensor of rank 2 and therefore we calculate ten elements in total: one for the energy density E, three for the flux vector F^α, and six for the stress tensor P^αβ. In Eq. (7), j^μ and j^ν are the direction cosines in Cartesian coordinates, and Ω is the solid angle subtended by 𝒬:

$$\begin{array}{c}\hfill j^{\mu }=(\frac{r-r_{\mathrm{d}}cos{\varphi }_{\mathrm{d}}}{D},\frac{-r_{\mathrm{d}}sin{\varphi }_{\mathrm{d}}}{D},\frac{z-H}{D}),\end{array}$$(8)

$$\begin{array}{c}\hfill \mathrm{d}\mathrm{\Omega }=\frac{-(rcos{\varphi }_{\mathrm{d}}-r_{\mathrm{d}})sin\delta +(z-H)cos\delta }{D^3}\mathrm{d}S,\end{array}$$(9)

where $\mathrm{d}S=\sqrt{1+c_3}{r}_{\mathrm{d}}\mathrm{d}{r}_{\mathrm{d}}\mathrm{d}{\varphi }_{\mathrm{d}}$.

2.3. Particles in the photon field

We now calculate the trajectory and velocity of the particles ejected from the disk when they interact with photons of the ambient radiation field.

The equation of motion under a relativistic, radiation treatment, is given by (Kato & Fukue 2020)

$$\begin{array}{c}\hfill f_{\mu }=-\frac{\partial {\mathrm{\Phi }}_{\mathrm{e}}}{\partial x^{\nu }}+R_{\mu ;\nu }^{\nu },\end{array}$$(10)

where f_μ is the four-force per unit volume. The effective potential Φ_e is the sum of gravitational (Φ_g) and centrifugal (Φ_c) potentials. The semicolon (; ) in the second term refers to the covariant differentiation of the energy-momentum tensor.

As we consider a disk with axial symmetry, the gravitational potential cancels out in the azimuthal coordinate: ∂Φ_g/∂x^α = (∂Φ_g/∂r, 0, ∂Φ_g/∂z). Furthermore, the centrifugal potential acts only in the radial direction: ∂Φ_c/∂x^α = (l²/r³, 0, 0), with $l=r_{\text{d}}^2{\omega }_{\text{K}}$ being the specific angular momentum of the disk, and ω_K the angular velocity.

The equations of motion of the ejected particles can be found working with Eq. (10). In terms of the nondimensional form of the radiation-field tensor elements ϵ, f^α, and p^αβ, the system of differential, tensorial, and coupled equations is as follows (equations originally derived by Watarai & Fukue 1999 Eq. (42)–(44), but now extended to second order in velocity):

Radial coordinate:

$$\begin{array}{cc}& \frac{\mathrm{d}{u}^r}{\mathrm{d}\tau }=-\frac{\partial {\mathrm{\Phi }}_{\mathrm{g}}}{\partial r}+\frac{l^2}{r^3}\hfill \\ \hfill & +\frac{1}{2}[\gamma f^r-p^{r\beta }{u}_{\beta }-{\gamma }^2ϵu^r+u^r(2\gamma f^{\beta }{u}_{\beta }-p^{\beta \delta }{u}_{\beta }{u}_{\delta })].\hfill \end{array}$$(11)

Azimuthal coordinate:

$$\begin{array}{cc}& \frac{1}{r}\frac{\mathrm{d}l}{\mathrm{d}\tau }=\frac{1}{2}[\gamma f^{\varphi }-p^{\varphi \beta }{u}_{\beta }-{\gamma }^2ϵ(l/r)\hfill \\ \hfill & +(l/r)(2\gamma f^{\beta }{u}_{\beta }-p^{\beta \delta }{u}_{\beta }{u}_{\delta })].\hfill \end{array}$$(12)

Height coordinate:

$$\begin{array}{cc}& \frac{\mathrm{d}{u}^z}{\mathrm{d}\tau }=-\frac{\partial {\mathrm{\Phi }}_{\mathrm{g}}}{\partial z}\hfill \\ \hfill & +\frac{1}{2}[\gamma f^z-p^{z\beta }{u}_{\beta }-{\gamma }^2ϵu^z+u^z(2\gamma f^{\beta }{u}_{\beta }-p^{\beta \delta }{u}_{\beta }{u}_{\delta })],\hfill \end{array}$$(13)

where u^μ denotes the four-velocity of the particles and γ the Lorentz factor, which is given by

$$\begin{array}{c}\hfill \gamma =\sqrt{1+u^r{u}^r+l^2/r^2+u^z{u}^z}.\end{array}$$(14)

The free parameter of these equations of motion is the launching radius of the particles, r₀, and we assume as initial condition that the particles co-rotate with the disk at this radius, $u_0^{\alpha }=(0,l_0/r_0,0)$.

We solve this system of equations numerically and assume that the kinematics of the disk-driven wind is roughly described by the trajectory and terminal velocities obtained for the test particles. As the accretion rate in the inner region of the disk is regulated at the Eddington rate, the mass loss in the wind is of the order of the super-Eddington accretion rate, Ṁ_dw ∼ Ṁ_input.

3. Collision of winds

The wind ejected from the disk collides with the stellar wind at the interaction region, where shocks are generated giving rise to particle acceleration (Fig. 2). An important quantity that characterizes the wind is the kinetic luminosity, L_K = Ṁv²/2, where Ṁ is the mass-loss rate and v the velocity of the fluid. A small fraction of the total kinetic power of the wind is transferred to relativistic particles, L_rel ∼ 0.1L_K, where we assume equipartition between relativistic protons and electrons (L_e = L_p). The mass-loss rate and velocity of the stellar wind are set according to the parameters found in the literature for the type of star we have chosen (e.g., Kobulnicky et al. 2019). In the case of the disk-driven wind, the velocity is obtained following the procedures described in the previous section. Given the orbital separation, the disk inclination, and the stellar size, we estimate that ∼10% of the original kinetic power reaches the acceleration region. We assume a circular orbit, that is, the geometry associated with the collision of winds does not depend on the orbital phase.

Fig. 2. Scheme of the wind collision seen in the rz plane (not to scale), adapted from Abaroa et al. (2021).

In this section, we describe the models for the collision region, the magnetic ambient field, and the shocks. We adopt a one-zone approximation for these calculations.

3.1. Contact discontinuity

The winds collide at a surface called the contact discontinuity (CD). The stagnation point (SP) is the closest position of the CD to the star, and is located where the ram pressures of the winds are in equilibrium,

$$\begin{array}{c}\hfill P_{\mathrm{ram}}(r_{\mathrm{BH}})={\rho }_{\mathrm{dw}}{v}_{\mathrm{dw}}^2={\rho }_{\ast \mathrm{w}}{v}_{\ast \mathrm{w}}^2=P_{\mathrm{ram}}(r_{\ast }).\end{array}$$(15)

Here, r_BH and r_* are the distances to the SP from the BH and from the center of the star, respectively. The density of the spherical stellar wind at this location is given by

$$\begin{array}{c}\hfill {\rho }_{\ast \mathrm{w}}=\frac{{\dot{M}}_{\ast }}{4\pi r_{\ast }^2{v}_{\ast \mathrm{w}}},\end{array}$$(16)

whereas the density of the disk-driven wind reads

$$\begin{array}{c}\hfill {\rho }_{\mathrm{dw}}=\frac{{\dot{M}}_{\mathrm{dw}}}{\mathrm{\Omega }{r}_{\mathrm{BH}}^2{v}_{\mathrm{dw}}},\end{array}$$(17)

where Ω = 2π(1 − cos θ) is the solid angle of the wind and θ the semi-opening angle of the wind. Solving these equations we obtain the position of the SP.

3.2. Magnetic field

The strength of the magnetic field at the CD is essentially determined by the stellar surface magnetic field B_*. The intensity of B_CD and its topology –dipole (i), radial (ii), or toroidal (iii)–, is given by (Eichler & Usov 1993):

$$\begin{array}{c}\hfill B_{\mathrm{CD}}\approx B_{\ast }\times \{\begin{array}{c}{R}_{\ast }^3/r_{\ast }^3,R_{\ast }<r_{\ast }<r_{\mathrm{A}},(\mathrm{i})\hfill \\ \hfill \\ R_{\ast }^3/r_{\mathrm{A}}{r}_{\ast }^2,r_{\mathrm{A}}<r_{\ast }<R_{\ast }(v_{\ast \mathrm{w}}/v_{\ast }^{\mathrm{rot}}),(\mathrm{ii})\hfill \\ \hfill \\ R_{\ast }^2{v}_{\ast }^{\mathrm{rot}}/r_{\mathrm{A}}{r}_{\ast }{v}_{\ast \mathrm{w}},R_{\ast }(v_{\ast \mathrm{w}}/v_{\ast }^{\mathrm{rot}})<r_{\ast },(\mathrm{iii}),\hfill \end{array}\end{array}$$(18)

where R_* is the stellar radius, r_A the Alfvén radius, and $v_{*}^{\text{rot}}\sim 0.1v_{\text{*w}}$ the surface rotation velocity.

3.3. Particle acceleration and shock

Particles are accelerated up to relativistic energies in the collision region through a first-order diffusive shock mechanism. Two shock fronts are generated: a forward shock (FS) that propagates through the stellar wind, and a reverse shock (RS) that propagates through the wind of the disk. The diffusive acceleration rate of the particles is given by (e.g., Protheroe et al. 1999):

$$\begin{array}{c}\hfill t_{\mathrm{ac}}^{-1}={\eta }_{\mathrm{ac}}\frac{eZc{B}_{\mathrm{CD}}}{E},\end{array}$$(19)

where e is the electric charge, Z the atomic number, and E is the energy of the particle. The acceleration efficiency, η_ac, depends on the diffusion coefficient of the particles, the shock velocity, and the angle between the magnetic field and the normal to the shock plane. We assume that the shock propagates perpendicular to the magnetic field and that diffusion occurs in the Bohm regime. Thus, the acceleration efficiency is

$$\begin{array}{c}\hfill {\eta }_{\mathrm{ac}}\approx \frac{3}{8}{\left(\frac{v_{\mathrm{sh}}}{c}\right)}^2,\end{array}$$(20)

where the shock velocities in the reference frame where one of the fluids is at rest, v_*w = 0, and the other one moves with a velocity v_dw, are given by (Lee et al. 1996):

$$\begin{array}{c}\hfill v_{\mathrm{RS}}=-\frac{4}{3}\frac{1}{1+\sqrt{n_{\ast \mathrm{w}}/n_{\mathrm{dw}}}}{v}_{\mathrm{dw}},\end{array}$$(21)

$$\begin{array}{c}\hfill v_{\mathrm{FS}}=\frac{4}{3}\frac{1}{1+\sqrt{n_{\mathrm{dw}}/n_{\ast \mathrm{w}}}}{v}_{\mathrm{dw}}.\end{array}$$(22)

Here, n_*w and n_dw are the numerical densities of the winds (n_w = ρ_w/m_p, with m_p the mass of the proton). The pressure and density of the shocked medium are calculated following the Rankine-Hugoniot relations (e.g., Lamers & Cassinelli 1999).

As we are interested in the nonthermal particle distribution, we investigate only adiabatic shocks; that is, where radiative losses are negligible. This is because in radiative shocks the gas in the shocked region emits large amounts of thermal radiation; the system therefore loses energy, the entropy increases, and the medium becomes increasingly homogeneous. If magnetic-inhomogeneities disappear, the acceleration efficiency decays abruptly, aborting the formation of nonthermal distributions.

The shock is adiabatic if the thermal cooling length R_Λ is larger than the size of the acceleration region Δx_ac (McCray & Snow 1979). The cooling length reads

$$\begin{array}{c}\hfill R_{\mathrm{\Lambda }}=\frac{5.9\times {10}^{11}\mu {(v_{\mathrm{sh}}/\mathrm{km}{\mathrm{s}}^{-1})}^3}{(n_{\mathrm{w}}/{\mathrm{cm}}^{-3})[\mathrm{\Lambda }(T_{\mathrm{sh}})/\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-3}]}\mathrm{cm}.\end{array}$$(23)

Here, n_w is the number density of the undisturbed medium, μ is the average molecular weight (μ = 0.6 for a fully ionized plasma), and Λ(T_sh) is the cooling function, which depends on the shock temperature (Raymond et al. 1976; Myasnikov et al. 1998; Wolfire et al. 2003). This latter function can be written as

$$\begin{array}{c}\hfill \mathrm{\Lambda }(T_{\mathrm{sh}})=\{\begin{array}{cc}4\times {10}^{-29}{T}_{\mathrm{sh}}^{0.8},\hfill & 55\mathrm{K}\le T_{\mathrm{sh}}<{10}^4\mathrm{K}\hfill \\ 7\times {10}^{-27}{T}_{\mathrm{sh}},\hfill & {10}^4\mathrm{K}\le T_{\mathrm{sh}}<{10}^5\mathrm{K}\hfill \\ 7\times {10}^{-19}{T}_{\mathrm{sh}}^{-0.6},\hfill & {10}^5\mathrm{K}\le T_{\mathrm{sh}}<4\times {10}^7\mathrm{K}\hfill \\ 3\times {10}^{-27}{T}_{\mathrm{sh}}^{0.5},\hfill & T_{\mathrm{sh}}\ge 4\times {10}^7\mathrm{K},\hfill \end{array}\end{array}$$(24)

where T_sh is given by

$$\begin{array}{c}\hfill T_{\mathrm{sh}}=18.21\mu {\left(\frac{v_{\mathrm{sh}}}{\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2\mathrm{K}.\end{array}$$(25)

We note that this temperature has a maximum value in a collisional plasma: it is self-regulated by the pair-creation, satisfying in any case k_BT_sh < 1 MeV (k_B is the Boltzmann constant).

We assume that the size of the acceleration region is a fraction of the distance from the BH to the SP, Δx_ac ∼ 0.1r_BH. As we consider a one-zone model, the acceleration region must be narrow enough to generate near-homogeneous conditions.

4. Radiative processes

Particles accelerated at the shock can cool through different processes and produce nonthermal radiation. The timescales associated to this cooling are related to the total energy-loss of the particles:

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t}\approx \frac{-E}{t_{\mathrm{cool}}},\end{array}$$(26)

where the total cooling rate is

$$\begin{array}{c}\hfill t_{\mathrm{cool}}^{-1}={\displaystyle \sum _i{t}_i^{-1},}\end{array}$$(27)

where t_i corresponds to each timescale of the involved cooling processes.

We assume advective escape; that is, particles are removed from the acceleration region by the bulk motion of the fluid. If the timescales of cooling are shorter than those of escape, particles radiate before they escape from the acceleration region. The maximum energy for each kind of particle can be inferred by looking at the point where the acceleration rate is equal to the total cooling or escape rate. This energy cannot exceed the maximum energy imposed by the Hillas criterion, $E_{\mathrm{e},\mathrm{p}}^{\max }<E_{\mathrm{Hillas}}^{\max }$.

As we are interested in nonthermal processes, we work at scales smaller than the size of the binary system and assume that rotation effects are negligible there. Effects caused by the orbital motion, such as Coriolis or centrifugal forces, could be relevant on larger scales and lead to strong disturbances in the flow and thermal processes. The analysis of such effects usually requires numerical simulations and is beyond the scope of this work.

4.1. Energy losses

We consider adiabatic and radiative losses. Adiabatic cooling is related to the work done by the particles of the wind to expand the shocked gas. Radiative cooling is caused by nonthermal processes as a consequence of the interaction of the wind particles with ambient fields and matter.

Our model is lepto-hadronic, and so we calculate the following radiative processes numerically:

– Synchrotron: interaction of protons and electrons with the ambient magnetic field, which will be amplified by a factor of 4 in the shocked region due to Rankine-Hugoniot relations.

– Inverse Compton (IC): collision of relativistic electrons with photons of the ambient radiation field.

– Bremmstrahlung: Coulombian interactions between relativistic electrons and cold matter.

– Photo-hadronic interactions: interaction of highly relativistic protons with photons of the ambient radiation field.

– Proton-proton: collision of relativistic protons with cold matter.

In addition, we take into account inelastic collision of particles with atoms of the dense medium; that is, ionization losses, which can be relevant in the 1–100 MeV range. We note that in this energy range, ionization losses largely dominate over Coulomb scatterings (see e.g., Fig. 7 from O’C Drury et al. 1996), and so the latter are not included in our analysis. The reader is referred to Romero & Paredes (2011), Romero & Vila (2014), and Müller & Romero (2020) plus references therein for additional details on radiative processes.

4.2. Particle distribution

We investigate the evolution of particles that are accelerated at the shock and injected into the surrounding medium. The medium around the shock is the shocked gas of the winds. In this paper, we restrict our analysis to this region. Beyond the binary, the surrounding medium has been influenced by the effects of the stellar winds, and so the system is expected to be located inside a bubble inflated by the winds and surrounded by a shell formed with the swept-up material at distances of a few to several parsecs, depending on the mass of the black hole progenitor. Inside the bubble, where the advected protons will be injected, the density is expected to be lower than that of the standard interstellar medium (e.g., around 0.01 cm⁻³ or less). In the shell, there should be sufficient material for hadronic interactions with the protons diffused or transported from the central source³.

The relativistic particles have a distribution given by dN = (r,E,t)dEdV, where n is the number density of particles, t the time, r the position, V the volume, and E the energy. The evolution of this distribution is determined by the transport equation (see e.g., Ginzburg & Syrovatskii 1964; Romero & Paredes 2011). We solve this equation numerically in steady state and in the one-zone approximation:

$$\begin{array}{c}\hfill \frac{\partial }{\partial E}\left[\frac{\mathrm{d}E}{\mathrm{d}t}N(E)\right]+\frac{N(E)}{t_{\mathrm{esc}}}=Q(E),\end{array}$$(28)

where t_esc ∼ Δx_ac/v_sh is the advection time, and the particle injection function,

$$\begin{array}{c}\hfill Q(E)=Q_0{E}^{-p}exp(-E/E_{\max }),\end{array}$$(29)

is a power-law in the energy with an exponential cutoff and a spectral index p = 2.2, which is characteristic of the Fermi first-order acceleration mechanism (see e.g., Drury 1983). The normalization constant Q₀ is obtained from

$$\begin{array}{c}\hfill L_{(\mathrm{e},p)}=\mathrm{\Delta }V{\displaystyle {\int }_{E_{(\mathrm{e},\mathrm{p})}^{\min }}^{E_{(\mathrm{e},\mathrm{p})}^{\max }}\mathrm{d}{E}_{(\mathrm{e},\mathrm{p})}{E}_{(\mathrm{e},\mathrm{p})}{Q}_{(\mathrm{e},\mathrm{p})}(E_{(\mathrm{e},\mathrm{p})}),}\end{array}$$(30)

where ΔV is the volume of the acceleration region, and $E_{(\mathrm{e},\mathrm{p})}^{\max }$ the maximum energy reached by protons and electrons, which is found by looking at the point where the acceleration rate is equal to the total cooling or escape rate.

4.3. Nonthermal emission

Once we have the particle distributions, we calculate the spectral energy distribution (SED) for each of the relevant processes involved in cooling. We find that in SCWBs, electrons typically cool by synchrotron and IC mechanisms, and protons escape from the acceleration region without significant cooling. The resultant nonthermal SED usually yields a broadband spectrum from radio waves (due to synchrotron emission) to gamma-rays (due to IC emission).

4.4. Wind emission

We calculate the thermal emission of the photosphere of the disk-driven wind assuming a spherically symmetric wind that expands with constant velocity equal to its terminal velocity. Since the mass-loss rate of the disk is much higher than the critical rate, the wind is optically thick and therefore we assume that it radiates locally as a blackbody. The temperature measured by an observer at infinity is given by (Fukue 2009):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{T}}{T}_{\mathrm{dw}}^4=\frac{\dot{e}{L}_{\mathrm{Edd}}}{{(1-\beta cos\mathrm{\Theta })}^{4}4\pi R^2},\end{array}$$(31)

where ė = Ė/L_Edd is the normalized comoving luminosity, β = v_dw/c the normalized velocity, Θ the angle of the flow with respect to the line of sight, and $R=\sqrt{r^2+z^2}$, with r and z the being cylindrical coordinates. We assume that the comoving luminosity is equal to the Eddington luminosity (ė = 1), as is commonly done in supercritical wind-models (e.g., Fukue 2009).

The apparent photosphere of this wind is defined as the surface where the optical depth τ_photo is unity for an observer at infinity. If the velocity of the wind is relativistic, the optical depth in the observer frame depends in general on the magnitude of the velocity and the viewing angle. The location of the apparent photosphere from the equatorial plane z_photo is (Fukue 2009):

$$\begin{array}{c}\hfill {\tau }_{\mathrm{photo}}={\displaystyle {\int }_{z_{\mathrm{photo}}}^{\infty }{\gamma }_{\mathrm{dw}}(1-\beta cos\mathrm{\Theta }){\kappa }_{\mathrm{co}}{\rho }_{\mathrm{co}}\mathrm{d}z=1,}\end{array}$$(32)

where γ_dw is the wind Lorentz factor, κ_co the opacity in the comoving frame, and ρ_co the wind density in the comoving frame. As we assume a fully ionized wind, the opacity is dominated by free electron scattering (κ_co = σ_T/m_p).

4.5. Absorption

Finally, we calculate the gamma absorption by pair creation from photon–photon annihilation, γ + γ → e⁺ + e⁻. The nonthermal photons in their way out of the acceleration region can find photons of the ambient radiation fields and annihilate. The absorption is quantified by the optical depth of the medium, τ_γγ. If the original luminosity of gamma rays is $L_{\gamma }^0(E_{\gamma })$, the attenuated luminosity reads:

$$\begin{array}{c}\hfill L_{\gamma }(E_{\gamma })=L_{\gamma }^0(E_{\gamma })·{\mathrm{e}}^{-\tau },\end{array}$$(33)

where e^−τ is the attenuation factor. The targets of the ambient radiation fields are photons from the star and from the disk-driven wind photosphere.

The process of annihilation is possible only above a kinematic energy threshold given by

$$\begin{array}{c}\hfill E_{\gamma }{E}_{\mathrm{ph}}>{(m_{\mathrm{e}}{c}^2)}^2,\end{array}$$(34)

in a frontal collision, where E_ph is the energy of the targets. The opacity caused by a photon–photon pair production for a photon created at a distance r from the center of the thermal source can be obtained from (Romero & Vila 2008):

$$\begin{array}{c}\hfill {\tau }_{\gamma \gamma }(E_{\gamma },r)={\displaystyle {\int }_{E_{\min }}^{\infty }{\int }_r^{\infty }{n}_{\mathrm{ph}}(E_{\mathrm{ph}},r^{\prime }){\sigma }_{\gamma \gamma }(E_{\mathrm{ph}},E_{\gamma })\mathrm{d}{r}^{\prime }\mathrm{d}{E}_{\mathrm{ph}},}\end{array}$$(35)

where n_ph is the density of the ambient radiation field. The total cross-section is given by (see e.g., Aharonian et al. 1985):

$$\begin{array}{c}\hfill {\sigma }_{\gamma \gamma }=\frac{\pi r_{\mathrm{e}}^2}{2}(1-{\xi }^2)[(3-{\xi }^4)ln\left(\frac{1+\xi }{1-\xi }\right)+2\xi ({\xi }^2-2)],\end{array}$$(36)

where r_e is the classical radius of the electron, and

$$\begin{array}{c}\hfill \xi ={(1-\frac{{(m_{\mathrm{e}}{c}^2)}^2}{E_{\gamma }{E}_{\mathrm{ph}}})}^{1/2}.\end{array}$$(37)

The blackbody density radiation of the star and the photosphere of the disk-driven wind is given by

$$\begin{array}{c}\hfill n_{\mathrm{ph}}=\frac{2E_{\mathrm{ph}}^2}{h^3{c}^3}\frac{1}{exp(E_{\mathrm{ph}}/k_{\mathrm{B}}T)-1},\end{array}$$(38)

where T is the temperature of the thermal source considered for each case; that is, T_dw or T_eff.

On the other side, free-free absorption (FFA) must also be taken into account. The collision of low-energy photons with particles of the dense medium leads to a cutoff in the SED at radio frequencies. The denser the medium, the higher the energy at which the cutoff occurs. Therefore, FFA will determine the turnover of the synchrotron spectrum in SCWBs, which is expected to be at ∼GHz frequencies (see e.g., Rybicki & Lightman 1986; del Palacio et al. 2016).

Other absorption processes, such as the photoelectric effect, direct Compton, or γ-nucleon pair creation, are not taken into account in this paper. Their cross-sections are not high enough to become relevant in the calculation of opacity given the ambient densities that we consider here (see Fig. 1 from Reynoso et al. 2011).

5. Results

In this section, we apply our model to a generic super-Eddington X-ray binary. We consider a star of spectral type O.5V (Table 1) and investigate four scenarios: in scenarios S1 and S2 we regard a BH with mass M_BH = 5M_⊙ and mass-accretion rates of 10²Ṁ_Edd and 10³Ṁ_Edd, respectively; in scenarios S3 and S4 we consider a BH with mass M_BH = 20M_⊙ and again accretion rates of 10²Ṁ_Edd and 10³Ṁ_Edd, respectively. The complete set of parameters is summarized in Table 2.

5.1. Wind

We calculate the radiation-field tensor (Eq. (7)) and in Fig. 3 we show the distribution of the energy density (ϵ) on the rz plane, where the black zone is the inflated inner disk. We obtain a strong azimuthal flux component of the radiation-field tensor. This distribution is the same in all four scenarios, because in the critical disk the radiation-field tensor depends on advection, viscosity, and adiabatic parameters, which remain the same in all cases.

Fig. 3. Contour maps of the spatial distribution of the normalized radiation energy density ϵ in the rz plane above the accretion disk. Both axes are in units of Schwarzschild radius. The color bar is the intensity of ϵ and the black zone is the inflated disk (f = 0.5, α = 0.5, γ = 4/3).

We solve Eqs. (11)–(13) to find the trajectory and velocity of the particles. Both quantities are determined by R^μν and therefore we obtain the same trajectories and terminal velocities in S1–S4. As an example, in Fig. 4 we show the normalized velocity of a test particle, with a launching radius of 40r_g (≡20r_s), which reaches a terminal velocity of ≈0.16c. This result does not vary much if we vary the launching radius ( ± 0.02c for ±20r_g).

Fig. 4. Normalized velocity of a wind test particle as a function of the Schwarzschild radius. The particle reaches a terminal velocity of ∼0.16c for a launching radius of r0 = 20rs (coincident with the vertical axis).

The particles describe a helical trajectory in the vicinity of the BH for two main reasons (Fig. 5). The first is the presence of the strong azimuthal components of the radiation field, which help to maintain the spiral geometry of the particles in the inner disk. The second reason is the condition imposed for the particle ejection, namely that the particles initially have only azimuthal velocity. The intensity of the radiation field decays rapidly with distance from the BH, and therefore the ejected particles follow a spiral trajectory near the BH, but beyond a certain radius (∼r_crit) they follow a free path with a strong component of the radial velocity.

Fig. 5. Trajectory of a test particle in the Cartesian 3D-space in units of Schwarzschild radius. The particles describe a helical trajectory above the inner disk because of the strong azimuthal radiation fields. The launching radius of this test particle is r0 = 20rs.

The overall result is an equatorial wind with terminal velocities of the order of 0.15c. The kinetic power of these winds is in the range 10^{39 − 41} erg s⁻¹, which is well above the power of the winds of typical WR or OB stars. Therefore, in general, the disk wind is expected to overwhelm the stellar wind.

5.2. Energy gain and losses

We follow the calculations in Sect. 3.1 and find that, in all four scenarios, the SP is located near the stellar surface and the wind of the disk completely sweeps up the stellar wind, as expected. Hence, the forward shock is in the stellar atmosphere, fully radiative, and completely unable to accelerate relativistic particles. Only the reverse shock (RS) is suitable for the task. As r_* ≈ R_*, the magnetic field at the CD is B_CD ≈ B_*.

The cooling length of the RS is greater than the size of the acceleration region in all cases (see Table 2); this is why the shock is adiabatic and the acceleration efficiency of the process is relatively high: η_ac ∼ 10⁻² (see Sect. 3.3). The shock velocity is ≈4.4 × 10⁹ cm s⁻¹ and the temperature of the shocked gas reaches ≈4.8 × 10¹⁰ K.

We calculate the energy gain and losses of the shock-accelerated particles following Sect. 4. Highly relativistic protons escape from the acceleration region without cooling in all scenarios considered here (with energies up to E_p ≈ 1 PeV) and are injected into the interstellar medium (ISM). Protons are advected, that is, they are removed from the collision region by the bulk motion of the fluid. They therefore do not interact with ambient material at scales similar to that of the system. Electrons cool mainly through IC and synchrotron mechanisms, and reach a maximum energy of E_e ≈ 100 GeV. To obtain the electron distribution, we solve the transport equation considering only the dominant IC and synchrotron losses, and a power-law injection function with a spectral index of 2.2 and an exponential cutoff (see Eq. (29)).

5.3. Spectral energy distribution

Figure 6 shows the SEDs of the four scenarios. The only thermal component of the spectrum is the photosphere of the optically thick disk-driven wind. The emission peak of the wind for S1 and S2 is ≈10³⁷ erg s⁻¹, whereas for S3 and S4 the peak is ≈10³⁸ erg s⁻¹. This occurs at energies of ∼100 eV for S1 and S3, and ∼30 eV for S2 and S4. Therefore, if M_BH increases, the luminosity is higher and, if the mass-accretion rate increases, the luminosity peak occurs at lower energies.

Fig. 6. Thermal and nonthermal SEDs of the four scenarios considered, S1–S4, in logarithmic scale, where a face-on inclination is assumed. S1 and S3 are shown in the left plot, whereas S2 and S4 are shown in the right plot. Dashed lines correspond to S1 (left) and S2 (right), solid lines correspond to S3 (left) and S4 (right). We plot the nonattenuated inverse Compton contributions in gray. The emission peak at high energies is ∼1033 erg s−1 for S1 and S2, and ∼1034 erg s−1 for S3 and S4. The gamma-ray absorption due to γγ annihilation is total for energies > 10 GeV.

In the case of the nonthermal spectrum, we calculate the emission due to synchrotron and IC losses. In the latter case, we consider the photon fields of the star and of the wind photosphere as targets. In all cases, the dominant IC contribution is that of the star. The luminosity in S3 and S4 is an order of magnitude greater than that in S1 and S2. This is because of the modification of the orbital parameters when the BH mass varies: to guarantee the overflow of the Roche lobe, the orbital semi-axis varies with M_BH, which results in variation in the size of the acceleration region and the photon density at SP, among other parameters. The emission peak at low energies is ∼10³³ erg s⁻¹ for S1 and S2, and ∼10³⁵ erg s⁻¹ for S3 and S4. At high energies, the emission peak is ∼10³² erg s⁻¹ (S1 and S2) and ∼10³⁴ erg s⁻¹ (S3 and S4). The gamma-ray absorption due to γγ annihilation is total for energies > 10 GeV in all scenarios⁴.

Attenuation due to material between the source and the observer, that is, absorption by external cold gas, is mainly in the optical-to-UV range and at soft X-rays. At radio wavelengths, refractive scintillation on free electrons of the ISM occurs at lower frequencies than predicted here. For high-energy gamma rays, the main absorbers are infrared (IR) fields and the cosmic microwave background (CMB), but their effects are only relevant for cosmological distances.

6. Application to NGC 4190 ULX 1

Ultraluminous X-ray sources (ULXs) are extragalactic point-like objects where the luminosity in the X-ray band appears to be higher than the Eddington luminosity (Bachetti 2016). ULXs are thought to be X-ray binaries with a stellar-mass compact object accreting at super-Eddington rates, where a beaming effect could be responsible for the luminosity observed in the X-ray band: the radiation emitted from the inner part of the accretion disk is geometrically collimated by the ejected wind, which is optically thick except in a narrow region around the black-hole axis and forms a cone-shaped funnel (King et al. 2001; King 2009; Kaaret et al. 2017; Fabrika et al. 2021).

We apply our model to estimate the radiation emitted by the ultraluminous X-ray source NGC 4190 ULX 1 (also known as CXO J121345.2+363754). Although many characteristics of this ULX remain poorly understood, several authors have explored the system and have provided constraints on some of its parameters (see e.g., Liu & Bregman 2005; Gladstone et al. 2013; Koliopanos et al. 2017; Kosec et al. 2018; Ghosh & Rana 2021).

In what follows, we describe the parameterization of the system and its components, and investigate the expected collision of winds. The complete set of parameters used in this section is detailed in Table 3.

6.1. System parameterization

The source is located in the nearby Galaxy NGC 4190 at a distance of d ≈ 3 Mpc (Tully et al. 2013). Observations made in 2010 using the XMM-Newton telescope reveal a long-term spectral variability in the 0.3–10.0 keV energy range: L_X ∼ 3 − 8 × 10³⁹ erg s⁻¹.

The angle i between the line of sight and the z-axis at which the disk of a ULX is observed determines the components of its spectrum: blackbody disk (BB) or Comptonization. If i is small, the observer is able to look into the funnel and see the innermost part of the disk: the spectrum shows only the BB component, which corresponds to thermal emission of the disk. This type of spectrum is called broadened disk (BD). If i is sufficiently large, another effect is observed: the interaction between photons and wind particles near the disk surface induces a Comptonization that produces a hardening in the spectrum. Most ULXs exhibit a combination of both phenomena in their X-ray spectrum.

Ghosh & Rana (2021) investigated the spectral properties of NGC 4190 ULX 1 and suggested that the ULX is in a BD state, and that the compact object is a BH with mass ∼10 − 30 M_⊙ accreting at super-Eddington rates. We fit the XMM-Newton observations (Epoch 3) with the supercritical advection-dominated disk model detailed in Sect. 2.1, assuming a mass-accretion rate of Ṁ_input = 10Ṁ_Edd. We also assume a face-on inclination i ≈ 0°, a BH mass 10 M_⊙ and a geometrical beaming factor b = 0.07. This factor is given by,

$$\begin{array}{c}\hfill b=\mathrm{\Omega }/4\pi =0.5(1-cos\vartheta ),\end{array}$$(39)

where Ω is the solid angle of the emission. The angle ϑ is related to the opening angles of the disk (δ) and its wind (θ): ϑ + δ + 2θ = 90°. Both angles, i and ϑ, can change over time, causing the spectral variability of the object (Fabrika et al. 2021).

On the other hand, Gladstone et al. (2013) provided constraints on the characteristics of the optical counterpart of the system. They suggested that, if M_BH = 10 M_⊙, the mass of the star could be < 50 M_⊙ and its radius < 86 R_⊙. We choose a star of type B2V for our model in light of one of the fittings these latter authors made from Hubble Space Telescope observations. If we apply Eq. (1) and consider the mass ratio M_*/M_BH, and the stellar radius involved (see Table 3), the transfer of mass in the binary system occurs for an orbital semi-axis a ≤ 15.2 R_⊙, which results in a period ≤38 h.

6.2. Collision of winds

The terminal velocity of the disk-driven wind is v_dw = 4.95 × 10⁹ cm s⁻¹, and therefore $L_{\mathrm{K}}^{\mathrm{dw}}=1.5\times {10}^{39}\mathrm{erg}{\mathrm{s}}^{-1}$, while $L_{\text{K}}^{*}=2.17\times {10}^{34}\text{erg}{\text{s}}^{\text{-1}}$. The SP is located near the stellar surface and the wind of the disk completely suppresses the stellar wind. We therefore only take into account the reverse shock (RS). As r_* ≈ R_*, the magnetic field at the CD is B_CD ≈ B_*.

The cooling length of the RS is R_Λ = 2.2 × 10¹³ cm and the size of the acceleration region is Δx_ac = 6.68 × 10¹⁰ cm; therefore, the shock is adiabatic and the acceleration efficiency of the process is η_ac = 10⁻², as in our general models. We calculate the energy gain and losses of the shock particles following Sect. 4. Highly relativistic protons escape from the acceleration region without cooling, as in our previous scenarios (with energies up to E_p ≈ 1 PeV), and are injected into the ISM. Electrons cool mainly through IC and synchrotron mechanisms. Figure 7 shows the timescales of electrons, which reach a maximum energy of E_e ≈ 0.32 TeV. To obtain the electron distribution, we solve the transport equation taking into account only IC and synchrotron losses, and a power-law injection function with a spectral index of 2.2 and an exponential cutoff.

Fig. 7. Timescales in logarithmic scale of the electron acceleration, escape, and cooling at the reverse shock in NGC 4190 ULX 1. Electrons reach a maximum energy of ≈0.32 TeV. The acceleration efficiency is 10−2.

6.3. Total SED

The SED of the ULX spans a broadband energy range. Figure 8 shows the thermal (wind and accretion disk) and nonthermal (colliding-winds shock) contributions of the system. We also show the sensitivity of the instruments ALMA, VLA (sub-mm waves), Fermi, and CTA (gamma rays), and observational data from XMM-Newton.

Fig. 8. Thermal and nonthermal SEDs of NGC 4190 ULX 1 in logarithmic scale (dashed lines). The nonthermal SED is partially attenuated for energies > 1 GeV and totally attenuated for energies > 50 GeV due to annihilation of γ-rays with the photon fields of the star and the photosphere of the disk-driven wind. The gray dashed lines are the nonattenuated IC contributions. The total SED is plotted with a solid black line. Data from XMM-Newton (Epoch 3), and the sensitivity of ALMA, Fermi, VLA, and CTA are also shown (instrument sensitivities were taken from Sotomayor & Romero 2022).

The luminosity in the IR band is ∼10³⁴ erg s⁻¹, which is relatively strong, though still undetectable at megaparsec distances. The luminosity in gamma-rays also reaches ∼10³⁴ erg s⁻¹. The attenuation factor (Fig. 9) has an effect on photons with energies ≳1 GeV. Most of the radiation above 1 GeV and all above 50 GeV is suppressed by the annihilation of the γ rays with the photon fields of the disk-driven wind and the star.

Fig. 9. Attenuation factors due to γγ-annihilation between high-energy nonthermal radiation and photon fields from the star and from the photosphere of the disk-driven wind in NGC 4190 ULX 1. The total attenuation is plotted with a black line.

7. Discussion

Our analysis of supercritical colliding wind binaries shows that these systems should exhibit broadband emission from radio to gamma rays. In this sense, they are similar to CWBs formed by two hot stars, such as O+WR binaries. However, there are important differences as well. If we compare our models with recent models of O+WR CWBs (Pittard et al. 2021), we find that (i) in SCWBs, the wind of the disk is far more powerful than the wind of the star. This results in stagnation points that are very close to the surface of the star. Efficient particle acceleration then can only occur in reverse shocks. (ii) We also see that the disk wind advects protons from the acceleration region before they have time to cool. Only electrons can cool locally. The resulting SED is consequently dominated by synchrotron and IC radiation. (iii) As the acceleration region is close to the star, the local magnetic field is relatively strong. Synchrotron emission reaches energies of hundreds of keV. As the medium is far more dense than in stellar CWBs, free-free absorption causes this radiation to turnover below ∼24 GHz. The total power at millimeter (mm) and submm wavelengths can be between three and five orders of magnitude higher in SCWBs than in stellar CWBs. (iv) IC is the dominant radiation mechanism at high energies. The stronger thermal fields of SCWBs (wind photosphere and star) provide the seed photons, but also impose a high-energy cutoff at ∼1 GeV through γ − γ attenuation. Instead, stellar CWBs can reach energies close to 1 TeV. (v) The strong magnetic fields in the acceleration region cut electromagnetic cascades in SCWBs. (vi) The SED is always dominated by the X-ray component associated with the disk or its wind in SCWBs. Finally, (vii) stellar CWBs have wider orbits and a variable separation between the components of the system. This produces variability related to the orbital period. On the contrary, the orbits of SCWBs should be mostly circularized. In general, CWBs are weaker than SCWBs, although span a broader energy range.

An interesting feature of SCWBs is their potential as cosmic ray sources. As mentioned, the strong wind of the disk drags away the relativistic protons before they cool. These protons, with maximum energies of the order of 1 PeV, are then injected into the ISM where they diffuse. Even if a fraction of just ∼1% of the wind kinetic power goes to relativistic protons, the cosmic ray output of a SCWB would be in the range 10^{37 − 39} erg s⁻¹. These protons might interact with ambient clouds at some distance from the system, producing gamma rays through pp → π⁰ + pp interactions and the subsequent pion decays π⁰ → γγ. The gamma-ray emission from the illuminated clouds can be even stronger than the emission from the binary itself. However, the spectrum should be softer because of propagation effects (Aharonian & Atoyan 1996). Recent modeling by Pittard et al. (2021) of particle acceleration in colliding wind binaries with wind velocities of a few 10³ km s⁻¹ and mG magnetic fields in the acceleration region demonstrate that up to ∼30% of the wind power can be transferred to nonthermal particles. This means that, in some extreme cases, a SCWB might inject up to ∼10⁴⁰ erg s⁻¹ in cosmic rays.

Another type of CWB is the so-called gamma-ray binary (GRB; e.g., LS 5039, PSR B1259-63, LSI +61° 303, PSR J2032+4127, and others; see, e.g., Dubus 2013; Chernyakova & Malyshev 2020). These sources are formed by a massive star (usually a Be star with a dense equatorial decretion disk and a fast wind) and a young pulsar in an eccentric orbit. The pulsar ejects a relativistic pair wind. The wind collision produces a broadband spectrum from electrons accelerated at the shock that cool by synchrotron and IC radiation. The two-peak SEDs are similar to those we estimate for SCWBs, but some differences are also clearly seen: (i) GRBs are less energetic because the spin-down luminosity of the pulsar is much smaller than the power of a supercritical wind. (ii) GRBs are highly variable. This variability is modulated with the orbital period. The orbital modulation of the different components of the broadband spectrum is a consequence of the orbital variability of geometrical parameters, such as the geometry of the contact surface of the stellar and pulsar winds. Absorption effects are also strongly variable. (iii) Hadronic interactions are likely when the pulsar crosses the equatorial disk of the star (e.g., Bykov et al. 2021). (iv) GeV flares have been observed after the periastron passage in sources such as PSR B1259-63 (Abdo et al. 2011; Chernyakova et al. 2014). These flares are attributed to the effects of the unshocked pulsar wind interaction with photons from the stellar disk (e.g., Khangulyan et al. 2012).

We finally mention that some black holes accreting at supercritical rates seem to be capable of launching mildly relativistic jets. A remarkable case in our Galaxy is the notorious microquasar SS433 (Fabrika 2004). This object resembles a ULX source seen edge on (Begelman et al. 2006). The accretion rate should be extremely high in order to explain the large jet power L_K ∼ 10⁴⁰ erg s⁻¹. Begelman et al. (2006) suggest rates of ∼ 5 × 10³ Ṁ_Edd ∼ 5 × 10⁻⁴ M_⊙ yr⁻¹, which are consistent with estimates of equatorial mass outflows inferred from radio observations (Blundell et al. 2001). These outflows, ejected toward either side of the jets, present a thermal spectrum and might well correspond to the radiation-driven wind of the hypercritical disk. The contamination from the jet base makes it impossible to disentangle contributions from colliding winds from those coming from the jet. However, the equatorial outflow might propagate well beyond the system and reveal itself if it collides with any clouds. The shock generated in the collision would convert the kinetic energy of the plasmoids into internal energy and relativistic particles, which might then cool by pp interactions with the cloud material. Such a scenario might explain the detection of a GeV source by the Fermi satellite on the side of SS433 (Bordas 2020; Li et al. 2020). We will explore the details of this hypothesis elsewhere.

8. Summary and conclusions

We explored the consequences of supercritical accretion in binary systems consisting of a hot star and a black hole. We find that a fraction of the kinetic power of the radiation-driven wind released by the accretion disk is transformed into relativistic particles in the region of the wind that collides with the star. Electrons are cooled locally, mainly through synchrotron and inverse Compton radiation. The radiation fields of the star and wind photosphere provide abundant thermal photons for the latter process; they also absorb high-energy radiation above a few GeV. Free-free absorption imposes a high-frequency turnover in the radio regime, suppressing centimeter radio waves, unlike the case of colliding wind binaries. The relativistic protons are blown away by the wind before they can cool down significantly. Once trapped by the outflow, these protons are transported to outer regions where they can interact with ambient gas away from the binary system, producing hadronic gamma-rays. Our most important finding is that, in addition to being strong thermal UV and X-ray sources, supercritical colliding wind binaries can be significant nonthermal sources at mm wavelengths and GeV energies.

Acknowledgments

The authors thank the anonymous referee for a careful and constructive review, and for his/her comments that improved this work. We thank also Daniela Pérez and Jiři Horák for fruitful discussions. This work was supported by grant PIP 0554 (CONICET). LA ackowledges the Universidad Nacional de La Plata for the education received. GER acknowledges the support from the Spanish Ministerio de Ciencia e Innovación (MICINN) under grant PID2019-105510GBC31 and through the Center of Excellence Mara de Maeztu 2020-2023 award to the ICCUB (CEX2019-000918-M).

References

All Tables

All Figures

	Fig. 1. Geometry of the present disk model. The radiation fields are calculated in the rz plane, where ϕ = 0. Here, 𝒬 is the position of the plasma element of the disk and 𝒫 the point of calculation on the rz plane. The scale height of the disk is H, and D is the distance between 𝒬 and 𝒫. The short arrow is the direction cosine jμ. This figure is adapted from Watarai & Fukue (1999).
In the text

	Fig. 2. Scheme of the wind collision seen in the rz plane (not to scale), adapted from Abaroa et al. (2021).
In the text

	Fig. 3. Contour maps of the spatial distribution of the normalized radiation energy density ϵ in the rz plane above the accretion disk. Both axes are in units of Schwarzschild radius. The color bar is the intensity of ϵ and the black zone is the inflated disk (f = 0.5, α = 0.5, γ = 4/3).
In the text

	Fig. 4. Normalized velocity of a wind test particle as a function of the Schwarzschild radius. The particle reaches a terminal velocity of ∼0.16c for a launching radius of r0 = 20rs (coincident with the vertical axis).
In the text

	Fig. 5. Trajectory of a test particle in the Cartesian 3D-space in units of Schwarzschild radius. The particles describe a helical trajectory above the inner disk because of the strong azimuthal radiation fields. The launching radius of this test particle is r0 = 20rs.
In the text

Fig. 6. Thermal and nonthermal SEDs of the four scenarios considered, S1–S4, in logarithmic scale, where a face-on inclination is assumed. S1 and S3 are shown in the left plot, whereas S2 and S4 are shown in the right plot. Dashed lines correspond to S1 (left) and S2 (right), solid lines correspond to S3 (left) and S4 (right). We plot the nonattenuated inverse Compton contributions in gray. The emission peak at high energies is ∼1033 erg s−1 for S1 and S2, and ∼1034 erg s−1 for S3 and S4. The gamma-ray absorption due to γγ annihilation is total for energies > 10 GeV.

In the text

	Fig. 7. Timescales in logarithmic scale of the electron acceleration, escape, and cooling at the reverse shock in NGC 4190 ULX 1. Electrons reach a maximum energy of ≈0.32 TeV. The acceleration efficiency is 10−2.
In the text

Fig. 8. Thermal and nonthermal SEDs of NGC 4190 ULX 1 in logarithmic scale (dashed lines). The nonthermal SED is partially attenuated for energies > 1 GeV and totally attenuated for energies > 50 GeV due to annihilation of γ-rays with the photon fields of the star and the photosphere of the disk-driven wind. The gray dashed lines are the nonattenuated IC contributions. The total SED is plotted with a solid black line. Data from XMM-Newton (Epoch 3), and the sensitivity of ALMA, Fermi, VLA, and CTA are also shown (instrument sensitivities were taken from Sotomayor & Romero 2022).

In the text

	Fig. 9. Attenuation factors due to γγ-annihilation between high-energy nonthermal radiation and photon fields from the star and from the photosphere of the disk-driven wind in NGC 4190 ULX 1. The total attenuation is plotted with a black line.
In the text