Supercritical colliding wind binaries

Context. Particle-accelerating colliding-wind binaries (PACWBs) are systems that are formed by two massive and hot stars and produce nonthermal (NT) radiation. The key elements of these systems are fast winds and the shocks that they create when they collide. Binaries with nonaccreting young pulsars have also been detected as NT emitters, again as a consequence of the wind-wind interaction. Black holes (BHs) might produce NT radiation by this mechanism if they accrete at super-Eddington rates. In such cases, the disk is expected to launch a radiation-driven wind, and if this wind has an equatorial component, it can collide with the companion star yielding a PACWB. These systems are supercritical colliding wind binaries (SCWBs). Aims. We aim to characterize the particle acceleration and NT radiation produced by the collision of winds in binary systems composed of a superaccreting BH and an early-type star. Methods. We estimated the terminal velocity of the disk-driven wind by calculating the spatial distribution of the radiation fields and their effect on disk particles. We then found the location of the wind collision region and calculated the timescales of energy gain and losses of relativistic particles undergoing diffusive acceleration. With this information, we were able to compute the associated spectral energy distribution of the radiation. Results. We find that the interaction of winds can produce NT emission from radio up to tens of GeV, with luminosities in the range of $\sim 10^{33}-10^{35} \, {\rm erg \, s^{-1}}$, which for the most part are contributed by electron synchrotron and inverse Compton radiation. Conclusions. We conclude that SCWBs, such as some ultraluminous X-ray sources and some Galactic X-ray binaries, are capable of accelerating cosmic rays and producing NT electromagnetic emission from radio to $\gamma$-rays, in addition to the thermal components.


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 3 km s −1 (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 π 0 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 34 to more than 10 37 erg s −1 .The most extreme cases are WR89, WR98, and WR140, with powers of between 6 and 8 times 10 37 erg s −1 .Less than 10 −7 of this power is finally radiated through synchrotron radio emission.The most luminous nonthermal radioemitting CWB is WR140, with a total radio luminosity of ∼2.6 × 10 30 erg s −1 .
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 31 erg s −1 .This luminosity amounts to ∼6 × 10 −6 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 6 L .The mass-loss rate of the primary is extremely high, reaching up to 10 −3 M yr −1 .The binary was detected in hard X-rays by INTE-GRAL (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 34 erg s −1 .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 π 0 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) andV404 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 37 erg s −1 .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.

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): R * lob = a 0.49 q2/3 0.6 q 2/3 + ln (1 + q 1/3 ) , where M * is the mass of the star and M BH the mass of the BH.Hence, the star overflows its Roche lobe R * 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.

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 with L Edd the Eddington luminosity 1 , η ≈ 0.1 the accretion efficiency, and c the speed of light.
The critical or spherization radius, given by 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 2 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 = f Q 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 2 .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 where √ c 3 = H/r d = tan δ, with H the scale height of the disk, δ the disk opening angle, and Here, c 3 (along with c 1 and c 2 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 radiationdominated gas in the inner disk.These values lead to a diskopening angle of δ = 30 • .

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 Q = (r d , φ d , H) of the disk surface, at any point P = (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).
The relativistic Doppler factor D provides the transformation between the reference frames (McKinley 1980): where z red is the redshift factor given by (Watarai & Fukue 1999) Here, D is the distance between P and Q, v φ = c 2 v K is the azimuthal velocity and v r = −c 1 αv K is the radial velocity, with v K = √ GM BH /r 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 −3 d .
The radiation-field tensor is given by (Rybicki & Lightman 1986) 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 Q: where dS = √ 1 + c 3 r d dr d dφ d .

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) 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 2 /r 3 , 0, 0), with l = r 2 d ω 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: Azimuthal coordinate: Height coordinate: A9, page 3 of 12 where u µ denotes the four-velocity of the particles and γ the Lorentz factor, which is given by The free parameter of these equations of motion is the launching radius of the particles, r 0 , and we assume as initial condition that the particles co-rotate with the disk at this radius, u α 0 = (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 .

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 /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 massloss 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 diskdriven 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.
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.

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, 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 whereas the density of the disk-driven wind reads 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.

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): where R * is the stellar radius, r A the Alfvén radius, and v rot * ∼ 0.1v * w the surface rotation velocity.

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): 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 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): 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 A9, page 4 of 12 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 magneticinhomogeneities 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 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 where T sh is given by 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 B T 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.

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: where the total cooling rate is 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 max e,p < E max Hillas .
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.

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.
-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.

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 −3 or less).In the shell, there should be sufficient material for hadronic interactions with the protons diffused or transported from the central source 3 .
The relativistic particles have a distribution given by dN = n(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: where t esc ∼ ∆x ac /v sh is the advection time, and the particle injection function, 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 firstorder acceleration mechanism (see e.g., Drury 1983).The normalization constant Q 0 is obtained from where ∆V is the volume of the acceleration region, and E max (e,p) 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.

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).

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): 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 = √ 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): 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 ).

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 0 γ (E γ ), the attenuated luminosity reads: where e −τ is the attenuation factor.The targets of the ambient radiation fields are photons from the star and from the diskdriven wind photosphere.The process of annihilation is possible only above a kinematic energy threshold given by 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): 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): where r e is the classical radius of the electron, and The blackbody density radiation of the star and the photosphere of the disk-driven wind is given by 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

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 2 ṀEdd and 10 3 ṀEdd , respectively; in scenarios S3 and S4 we consider a BH with mass M BH = 20M and again accretion rates of 10 2 ṀEdd and 10 3 ṀEdd , respectively.The complete set of parameters is summarized in Table 2.

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.
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 ).
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.
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 −1 , 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.

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 −2 (see Sect. 3.3).The shock velocity is ≈4.4×10 9 cm s −1 and the temperature of the shocked gas reaches ≈4.8 × 10 10 K.
We calculate the energy gain and losses of the shockaccelerated 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)).

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 37 erg s −1 , whereas for S3 and S4 the peak is ≈10 38 erg s −1 .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.
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 semiaxis 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 33 erg s −1 for S1 and S2, and ∼10 35 erg s −1 for S3 and S4.At high energies, the emission peak is ∼10 32 erg s −1 (S1 and S2) and ∼10 34 erg s −1 (S3 and S4).The gamma-ray absorption due to γγ annihilation is total for energies >10 GeV in all scenarios4 .

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  We plot the nonattenuated inverse Compton contributions in gray.The emission peak at high energies is ∼10 33 erg s −1 for S1 and S2, and ∼10 34 erg s −1 for S3 and S4.The gamma-ray absorption due to γγ annihilation is total for energies >10 GeV. 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.

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.0keV energy range: L X ∼ 3−8 × 10 39 erg s −1 .
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, 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.

Collision of winds
The terminal velocity of the disk-driven wind is v dw = 4.95 × 10 9 cm s −1 , and therefore L dw K = 1.5 × 10 39 erg s −1 , while L * K = 2.17 × 10 34 erg s −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 13 cm and the size of the acceleration region is ∆x ac = 6.68 × 10 10 cm; therefore, the shock is adiabatic and the acceleration efficiency of the process is η ac = 10 −2 , 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.
show the sensitivity of the instruments ALMA, VLA (sub-mm waves), Fermi, and CTA (gamma rays), and observational data from XMM-Newton.
The luminosity in the IR band is ∼10 34 erg s −1 , which is relatively strong, though still undetectable at megaparsec distances.The luminosity in gamma-rays also reaches ∼10 34 erg s −1 .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.

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 −1 .These protons might interact with ambient clouds at some distance from the system, producing gamma rays through pp → π 0 + pp interactions and the subsequent pion decays π 0 → γγ.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  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.(Aharonian & Atoyan 1996).Recent modeling by Pittard et al. (2021) of particle acceleration in colliding wind binaries with wind velocities of a few 10 3 km s −1 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 40 erg s −1 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 40 erg s −1 .Begelman et al. (2006) suggest rates of ∼5 × 10 3 ṀEdd ∼ 5 × 10 −4 M yr −1 , 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.

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 A9, page 11 of 12 can be significant nonthermal sources at mm wavelengths and GeV energies.

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

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

Fig. 3 .Fig. 4 .
Fig. 3. Contour maps of the spatial distribution of the normalized radiation energy density in the rz plane above the accretion disk.Both axes 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).

Fig. 5 .
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 r 0 = 20r s .

Fig. 6 .
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 ∼10 33 erg s −1 for S1 and S2, and ∼10 34 erg s −1 for S3 and S4.The gamma-ray absorption due to γγ annihilation is total for energies >10 GeV.

Fig. 8 .
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).

Table 1 .
from Reynoso et al.  2011).Parameters adopted in the model for the star of type O.5V.

Table 2 .
Parameters of the different scenarios calculated for the model.
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 .