© The Authors 2025

1. Introduction

Nonthermal emission from gamma-ray binaries – systems composed of a massive star and a compact object, typically a pulsar or black hole – has been robustly detected at high energies (HEs; giga-electronvolts) and very high energies (VHEs; tera-electronvolts) (see e.g., Dubus 2006, 2013; Maraschi & Treves 1981; Huber et al. 2021; Bogovalov et al. 2008; Bednarek 2009; Bosch-Ramon et al. 2012; Hinton et al. 2009; De Sarkar et al. 2022). In addition, the recurrent nova binary system RS Ophiuchi, comprising of a white dwarf and a red giant companion, has recently been observed in both giga-electronvolt and tera-electronvolt gamma rays (see e.g., Tatischeff & Hernanz 2007; Acciari et al. 2022; H.E.S.S. Collaboration 2022; Cheung et al. 2022; Diesing et al. 2023; De Sarkar et al. 2023), further indicating the presence of efficient particle acceleration. A distinct subclass of binary systems, colliding wind binaries (CWBs), consists of two massive stars whose powerful stellar winds interact, forming a wind collision region (WCR) capable of accelerating particles to relativistic energies (see e.g., Eichler & Usov 1993; Dougherty & Williams 2000; De Becker & Raucq 2013; Benaglia et al. 2001; Benaglia & Romero 2003; Reimer et al. 2006; Reitberger et al. 2014; Grimaldo et al. 2019). These systems are promising candidates for gamma-ray emission via hadronic and/or leptonic processes.

The only firmly established giga-electronvolt-emitting CWB to date is η Carinae (see e.g., Farnier et al. 2011; Bednarek & Pabich 2011; Ohm et al. 2015; Gupta & Razzaque 2017; Reitberger et al. 2015; Balbo & Walter 2017), which has also been detected at tera-electronvolt energies (H.E.S.S. Collaboration 2020). Another notable system is γ² Velorum (WR 11), the closest known CWB at a distance of ${336}_{-7}^{+8}$ pc (North et al. 2007), consisting of a WC8 star and an O7.5 companion (Henley et al. 2005) in a 78.53 ± 0.01 day orbit (Schmutz et al. 1997) with a moderate eccentricity of 0.334 ± 0.003 (North et al. 2007). This system exhibits both thermal (Purton et al. 1982; Benaglia et al. 2019) and nonthermal (Chapman et al. 1999; Benaglia 2016) radio emission, explained by synchrotron emission of accelerated electrons at the WCR. De Becker & Raucq (2013) explains the nonthermal radio emission, including its absorption in the Wolf–Rayet (WR) star wind. The WCR has been observed in ultraviolet (St Louis et al. 1993) and X-rays by ROSAT (Willis et al. 1995), ASCA (Rauw et al. 2000), Chandra (Skinner et al. 2001), and XMM-Newton (Schild et al. 2004), including a wide shock-cone opening half angle of approximately ≈85° revealed by Chandra observations (Henley et al. 2005).

While γ² Velorum has previously been associated with a gamma-ray source (Pshirkov 2016), only recent Fermi-LAT data analysis has revealed a significant high-energy source, 4FGL J0809.5−4714, spatially consistent with γ² Velorum, with hints of orbital variability (Martí-Devesa et al. 2020), where flux increases around apastron and decreases near periastron. This behavior differs from η Carinae, where maximum γ-ray emission occurs near periastron and vice versa. This phenomenon thus motivates a detailed modeling of the phase-resolved emission.

Although a leptonic scenario has been proposed to explain the gamma-ray emission from γ² Velorum (Tatischeff et al. 2004), comparatively more recent literature favors a hadronic origin, with the emission primarily arising from neutral pion (π⁰) decay following proton acceleration in the WCR (Reitberger et al. 2017). However, the earlier studies did not particularly focus on the orbital modulation of the gamma-ray emission, which was only identified relatively recently. In this paper, we present a simple semianalytic treatment of the phase-dependent hadronic emission for γ² Velorum, incorporating both the orbital geometry and γ − γ absorption. We demonstrate that reproducing the observed modulation indeed necessitates an intrinsic geometric variation of the acceleration region, while attenuation effects along the line of sight play a comparatively minor role.

In Sect. 2, we outline the theoretical framework adopted in this study. Sect. 2.1 describes the physical environment of the γ² Velorum system. The proton injection model is detailed in Sect. 2.2, while the treatment of γ − γ absorption is presented in Sect. 2.3. The results of our phase-resolved modeling are discussed in Sect. 3, followed by our main conclusions in Sect. 4.

2. The model

2.1. CWB environment

In CWBs such as γ² Velorum, each massive star drives a powerful, accelerated stellar wind characterized by high mass-loss rates and supersonic terminal velocities. As these winds expand outward, they eventually collide, forming a wind interaction zone bounded by two strong shocks and a contact discontinuity that separates the shocked plasma of each star. This structure, i.e., the WCR, is shaped by the wind momentum balance and orbital geometry, and typically adopts a curved morphology that wraps around the star with the weaker wind momentum flux; namely, the O-star companion in the case of γ² Velorum. Across these shocks, the initially supersonic winds are abruptly decelerated and compressed, resulting in the formation of high-temperature plasma and an enhancement of the local magnetic field. These shocked regions are ideal sites for diffusive shock acceleration (DSA), wherein charged particles can gain energy through repeated crossings of the shock front. Both electrons and protons can be accelerated to relativistic energies in this environment, producing observable nonthermal emission via leptonic and hadronic processes, respectively. The efficiency of acceleration and the resulting emission depend sensitively on the local physical conditions – such as the magnetic field strength, ambient density, and orbital geometry – which vary with orbital phase due to the system’s eccentricity.

As was discussed earlier, a wide shock-cone opening half angle of approximately 85° was observed in γ² Velorum (Henley et al. 2005), which can be attributed to the effects of radiative braking (Reitberger et al. 2017). The sudden radiative braking effect in CWBs, especially in WR + O type binaries, has previously been extensively studied both analytically and numerically in Gayley et al. (1997). The authors showed that radiative braking may prove to be a necessary effect, as in some systems, not considering this effect may render the hydrodynamical wind-wind ram balance impossible. In this scenario, the intense radiation field of the O star decelerates the approaching WR wind before the collision, effectively reducing its velocity significantly from the terminal value near the WCR. This radiative braking alters the local momentum balance and leads to a wider shock cone than what would be expected from unimpeded wind interactions alone. Motivated by this, we adopted an effective WR wind velocity that is a reduced fraction of its nominal terminal value. This choice not only reproduces the observed shock-cone opening half angle, but also provides a more realistic estimate of the wind momentum flux in the vicinity of the collision region.

The key physical parameters for the two companion stars of the γ² Velorum system used in this study are summarized in Table 1. The orbital parameters of the source are given in Table 2. Based on these parameters, we first estimated the wind momentum balance ratio,

$$\begin{array}{c}\hfill {\eta }_{\varphi }=\frac{{\dot{M}}_{\mathrm{O}}{V}_{\mathrm{O},\infty }}{{\dot{M}}_{\mathrm{W}}{V}_{\mathrm{W},\mathrm{brak}}},\end{array}$$(1)

where Ṁ_W and Ṁ_O are the mass-loss rates of the WR and O stars, respectively. The mass-loss rate of the WR star in γ² Velorum remains uncertain, as different observational methods yield significantly different values. Polarimetric analyses suggest a rate of approximately 8.0 × 10⁻⁶ M_⊙ yr⁻¹, while estimates based on radio emission indicate a higher value of ∼3.0 × 10⁻⁵ M_⊙ yr⁻¹. This discrepancy is likely a consequence of wind clumping, which can significantly enhance the free-free radio emission and thereby inflate mass-loss estimates derived from radio observations. In contrast, polarimetric measurements are less sensitive to small-scale density inhomogeneities and are generally considered more robust in systems where clumping is expected. For these reasons, and to avoid overestimating the wind momentum flux of the WR star, we adopted the polarimetric value in this study.

Next, V_O, ∞ is the O-star wind terminal velocity, whereas V_W, brak is the radiatively braked velocity of the WR wind, obtained by considering a fraction of the WR wind terminal velocity, V_W, ∞. To obtain this fraction, we solved the shock-cone opening half angle equation from Eichler & Usov (1993), i.e.,

$$\begin{array}{c}\hfill \theta (\mathrm{in}\mathrm{deg})\approx 2.1(1-\frac{{\eta }_{\varphi }^{2/5}}{4}){\eta }_{\varphi }^{1/3}\left(\frac{{180}^{°}}{\pi }\right),\mathrm{for}{10}^{-4}\le {\eta }_{\varphi }\le 1,\end{array}$$(2)

for θ ≈ 85°. We found that the wide shock-cone opening half-angle can be reconciled for the braked WR wind velocity, being around 5% of its nominal terminal velocity at the collision region. The braking radius, where wind braking becomes significant, can be calculated from the analytical formulation provided in Gayley et al. (1997). We found that for the reduction in WR wind velocity due to radiative braking, the braking radius falls very close to the stagnation point radius, which is expected in this scenario. We further found that without the inclusion of radiative braking, the ram balance could not be achieved during the periastron, as was found in Gayley et al. (1997). We found that the reflection factor, S (see Gayley et al. 1997), representing the radiative momentum from the O star, proves to be significant in this case, further confirming that sudden braking of WR wind due to the O-star radiative momentum is highly expected in γ² Velorum. We considered the braked WR wind velocity, V_W, brak, for further calculations.

The shock apex, or stagnation point, defined as the location along the line connecting the two stars where the ram pressures of the winds equilibrate, lies closer to the O star in the case of γ² Velorum. The distance of the shock from each star – and hence the overall shape of the WCR – can be derived from the aforementioned ratio. In an elliptical orbit, the separation between two stars can be given by

$$\begin{array}{c}\hfill R_{\mathrm{sep}}=a(1-e^2)/(1+ecos(\varphi )).\end{array}$$(3)

Using this equation, and assuming the collision of two spherical winds, one can calculate the distances, R_sh, O and R_sh, W, respectively, from the O and WR stars to the wind interaction region along the line joining the O and WR stars following

$$\begin{array}{c}\hfill R_{\mathrm{sh},\mathrm{O}}=\frac{\sqrt{{\eta }_{\varphi }}}{1+\sqrt{{\eta }_{\varphi }}}{R}_{\mathrm{sep}}\mathrm{and}{R}_{\mathrm{sh},\mathrm{W}}=\frac{1}{1+\sqrt{{\eta }_{\varphi }}}{R}_{\mathrm{sep}}.\end{array}$$(4)

Furthermore, due to the eccentric orbit of γ² Velorum, the stellar separation – and thus the location of the WCR – varies with orbital phase. This variation leads to a changing shock geometry, impacting the local physical conditions such as the density, and the apex magnetic field strength. To properly account for these phase-dependent effects, we calculated the local wind velocity of the O star at the position of the WCR using a β-velocity law. We adopted an index of β_O = 0.8, as is appropriate for O stars (Gayley et al. 1997). Given the O-star radius, R_O, the local wind velocity at the standoff distance, R_sh, O, was computed as

$$\begin{array}{c}\hfill V_{\mathrm{O},\mathrm{local}}=V_{\mathrm{O},\infty }{(1-\frac{R_{\mathrm{O}}}{R_{\mathrm{sh},\mathrm{O}}})}^{{\beta }_{\mathrm{O}}}.\end{array}$$(5)

For the WR star, we adopted a constant local wind velocity, V_W, local = V_W, brak, corresponding to the effective braked velocity. This choice was motivated by the assumption that, due to a much smaller radius of the WR star, its wind acceleration takes place in a region largely unaffected by the influence of the companion O star (Gayley et al. 1997). It is important to note that radiative braking acts predominantly along the line connecting the two stars, where the winds directly interact with each other.

These local wind velocities were used to derive the magnetic field strength from each star at the WCR apex and the WR wind number density. The kinetic luminosity of the WR wind entering the shock, as well as the advection timescale, was also estimated using the local velocities.

To compute the magnetic field strength at the location of the WCR, we adopted a structured magnetic field prescription following the approach of Eichler & Usov (1993), in which the field transitions from a dipole-dominated regime near the stellar surface to a toroidal configuration at larger radii. We first calculated the Alfvén radius by defining the wind magnetic confinement parameter, η_*, as (Ud-Doula & Owocki 2002)

$$\begin{array}{c}\hfill {\eta }_{\ast }=\frac{B_{\ast }^2{R}_{\ast }^2}{{\dot{M}}_{\ast }{V}_{\ast ,\mathrm{local}}},\end{array}$$(6)

where B_* is the surface magnetic field strength, R_* is the stellar radius, Ṁ_∗ is the mass-loss rate, and V_*, local is the local wind velocity (or the braked velocity in the case of the WR). The Alfvén radius, R_A, was estimated using the strong confinement scaling (Ud-Doula et al. 2008):

$$\begin{array}{c}\hfill R_{\mathrm{A}}\approx R_{\ast }(0.3+{\eta }_{\ast }^{1/4}).\end{array}$$(7)

Given the stellar surface rotational velocity, V_rot (∼0.2V_*, ∞), the magnetic field strength, B(r), at a radial distance, r, from the center of the star is then defined in three regimes:

• A dipole regime (R_* ≤ r <  R_A),

  $$\begin{array}{c}\hfill B(r)=B_{\ast }{\left(\frac{R_{\ast }}{r}\right)}^3,\end{array}$$(8)

• A radial regime ($R_{\mathrm{A}}\le r<R_{\ast }\left(\frac{V_{\ast ,\mathrm{local}}}{V_{\mathrm{rot}}}\right)$),

  $$\begin{array}{c}\hfill B(r)=B_{\ast }\frac{R_{\ast }^3}{R_{\mathrm{A}}{r}^2},\end{array}$$(9)

• A toroidal regime ($r\ge R_{\ast }\left(\frac{V_{\ast ,\mathrm{local}}}{V_{\mathrm{rot}}}\right)$),

  $$\begin{array}{c}\hfill B(r)=B_{\ast }\left(\frac{V_{\mathrm{rot}}}{V_{\ast ,\mathrm{local}}}\right)\frac{R_{\ast }^2}{R_{\mathrm{A}}r}·\end{array}$$(10)

This three-zone prescription was applied separately to each star in the binary. The magnetic field at the WCR location was calculated considering the contributions from both stars.

Note that the shock can undergo a strong magnetohydrodynamic (MHD) compression in the region. The value of the corresponding compression ratio depends on whether the shock is in a strong adiabatic regime or has transitioned into the radiative regime. For a strong, adiabatic shock with adiabatic index γ = 5/3, the compression ratio is given by χ ≈ 4. However, if the shock has transitioned into the radiative regime, then the compression ratio can increase way above the value of χ ≈ 4. To determine the thermal regime of the shock, we adopted the diagnostic criterion outlined in Stevens et al. (1992) and Pittard & Parkin (2010), which introduces the parameter ξ, which can be calculated with

$$\begin{array}{c}\hfill \xi \approx \frac{{(V_{\ast ,\mathrm{local}}/{10}^8\mathrm{cm}{\mathrm{s}}^{-1})}^4(R_{\mathrm{sh},\ast }/{10}^{12}\mathrm{cm})}{({\dot{M}}_{\ast }/{10}^{-7}{M}_{\odot }{\mathrm{yr}}^{-1})}·\end{array}$$(11)

The parameter ξ acts as a characteristic measure of the importance of cooling. If ξ ≳ 1, then the wind can be assumed to be adiabatic, whereas if ξ ≪ 1, the wind can be assumed to be in the radiative regime. Using the parameters of γ² Velorum, we find that while the O-star shock remains in the adiabatic regime (ξ >  1) across all orbital phases, the WR shock consistently transitions into the radiative regime, with ξ ≪ 1 at all orbital phases. A radiative WR shock indicates a stronger post-shock compression and amplified magnetic field. In contrast, the adiabatic nature of the O-star shock leads to weaker compression and consequently less efficient magnetic field amplification. To incorporate this into our modeling, we adopted a higher compression ratio for the WR wind, χ_W ≈ 10, consistent with its radiative nature. For the O-star wind, we retained the standard adiabatic compression ratio, χ_O ≈ 4.

Assuming both shocks compress the upstream magnetic field, the total post-shock magnetic field at the WCR apex is given by

$$\begin{array}{c}\hfill B_{\mathrm{apex}}={[{(B_{\mathrm{sh},\mathrm{O}}·{\chi }_{\mathrm{O}}^{2/3})}^2+{(B_{\mathrm{sh},\mathrm{W}}·{\chi }_{\mathrm{W}}^{2/3})}^2]}^{1/2}.\end{array}$$(12)

This combined magnetic field strength, B_apex, was then used in calculating the maximum particle energies in the WCR.

In our model, we consider relativistic protons interacting primarily with the shocked WR wind material at the WCR. Crucially, given the presence of a WC8 subtype star in γ² Velorum, WR wind is hydrogen-poor, but primarily composed of heavier elements such as helium (He), carbon (C), and oxygen (O), as is observed from their strong line emissions (see Crowther 2007 and references therein). The pre-shock target material number density at a radial distance, R_sh, W, from the WR star, assuming a steady, spherically symmetric wind, is given by the standard continuity equation: $n_{\mathrm{w}}={\dot{M}}_{\mathrm{W}}/(4\pi R_{\mathrm{sh},\mathrm{W}}^2{V}_{\mathrm{W},\mathrm{local}}\mu m_{\mathrm{p}})$, where Ṁ_W and V_W, local are the same parameters discussed above, μ is a factor signifying the mean molecular weight, and m_p is the proton mass. However, since the gamma-ray-producing protons are assumed to interact in the post-shock region on the side of the WR wind, we adopted a shocked material density enhanced by a compression factor, χ_W(≈10), as is expected for the radiative nature of the WR wind shock. This yields a post-shock density of

$$\begin{array}{c}\hfill n_{\mathrm{w},\mathrm{shocked}}={\chi }_W·n_{\mathrm{w}}=\frac{{\chi }_W{\dot{M}}_{\mathrm{W}}}{4\pi R_{\mathrm{sh},\mathrm{W}}^2{V}_{\mathrm{W},\mathrm{local}}\mu m_{\mathrm{p}}},\end{array}$$(13)

which we used as the effective target density for the hadronic interactions. For a representative composition of 55% He, 40% C, and 5% O by mass (see e.g., Crowther 2007), we derived the value of μ ≈ 5.75, which we adopted when computing the number density of targets for hadronic interactions. This quantity plays a key role in setting the normalization of the hadronic gamma-ray emission, as it directly governs the pion production rate at the WCR apex.

2.2. Proton injection

To model the hadronic gamma-ray emission from γ² Velorum, we computed the phase-dependent injection of relativistic protons at the WCR. This injection process is governed by a combination of physical quantities: the maximum energy attainable by shock-accelerated protons; the WR wind kinetic energy available for particle acceleration over the advection timescale; the fraction of the WR wind intercepted by the WCR; the fraction of accelerated protons retained in the emission region before they can escape; and the effective volume of the acceleration region, which controls the number of particles that can be accelerated. Most of these quantities are modulated by the binary’s orbital geometry and wind interaction dynamics, and together they determine both the normalization and shape of the proton injection spectrum at each phase. Despite the hydrogen-poor nature of the WR wind, the presence of even a trace population of protons can suffice for DSA. Additionally, protons can be supplied at the WCR by the hydrogen-rich O-star wind through mixing or shock interface instabilities. For simplicity, we assumed that the accelerated population is predominantly composed of protons, consistent with standard cosmic-ray acceleration theory. The target material for hadronic interactions was modeled using the WR wind composition. To account for the heavy-element-enriched nature of the WR wind in γ² Velorum, we adopted a nuclear enhancement factor of ϵ_enh ≈ 2 in our hadronic gamma-ray calculations. This accounts for the increased pion production efficiency in proton interactions with heavier nuclei (e.g., Kafexhiu et al. 2014), and is appropriate for the energy range considered (≳1 GeV).

The maximum energy cutoff of the proton distribution, E_p, max, was set by comparing the characteristic acceleration time with energy loss and escape timescales. Following the formalism in Bednarek & Pabich (2011), we considered two distinct regimes that limit acceleration: the loss-limited regime, in which losses due to hadronic interaction dominate, and the escape-limited regime, in which particles advect out of the WCR before significant energy losses occur.

In the loss-limited scenario, the maximum energy is given by

$$\begin{array}{c}\hfill E_{p,\max }\approx 6.3{\xi }_{-5}{B}_{\mathrm{apex}}{R}_{13}^2{v}_3/{\dot{M}}_{-4}\mathrm{TeV},\end{array}$$(14)

where ξ₋₅ = (V_W, local/c)²/10⁻⁵ is the acceleration parameter, B_apex is the post-shock magnetic field strength at the apex of the WCR, R₁₃ = R_sh, W/10¹³ cm is the WR shock standoff distance scaled to 10¹³ cm, v₃ = V_W, local/10⁸ cm s⁻¹ is the wind speed in units of 1000 km s⁻¹, and Ṁ₋₄ = Ṁ_W/(10⁻⁴ M_⊙yr⁻¹) is the normalized mass-loss rate.

In the escape-limited regime, in which the dominant limitation is the spatial confinement inside the WCR, the maximum energy is

$$\begin{array}{c}\hfill E_{p,\max }\approx 30{\xi }_{-5}{B}_{\mathrm{apex}}{R}_{13}/v_3\mathrm{TeV}.\end{array}$$(15)

The choice between these regimes was made dynamically at each orbital phase by evaluating the critical condition

$$\begin{array}{c}\hfill {\dot{M}}_{-4}>0.2R_{13}{v}_3^2.\end{array}$$(16)

If this inequality was satisfied, the loss-limited formula was used; otherwise, the escape-limited one applied. This approach ensured that the maximum energy of protons responded sensitively to phase-dependent variations in magnetic field strength, wind velocity, and shock standoff distance.

Even after successful acceleration, relativistic protons may escape the WCR via diffusion before undergoing hadronic interactions. This escape reduces the hadronic gamma-ray yield. We accounted for this effect by using the retention fraction, f_ret(E), defined as the ratio of the diffusion timescale to the total timescale including both diffusion and hadronic interactions:

$$\begin{array}{c}\hfill f_{\mathrm{ret}}(E)=\frac{{\tau }_{\mathrm{diff}}}{{\tau }_{\mathrm{diff}}+{\tau }_{\mathrm{loss}}}·\end{array}$$(17)

The isotropic diffusion timescale was estimated as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{diff}}=\frac{R_{\mathrm{sh},\mathrm{W}}^2}{6D(E)}·\end{array}$$(18)

The diffusion coefficient, D(E), was modeled using a power-law energy scaling:

$$\begin{array}{c}\hfill D(E)=D_0{\left(\frac{E}{1\mathrm{MeV}}\right)}^{\delta },\end{array}$$(19)

with δ = 0.33, corresponding to Kolmogorov turbulence, and D₀ ∼ 10¹⁹ cm² s⁻¹, representative of turbulent conditions in shocked stellar winds. The choice of diffusion coefficient is consistent with those explored in Reitberger et al. (2017), and orders of magnitude less than that found in the interstellar medium (see e.g., Ptuskin et al. 2006; Putze et al. 2010; Génolini et al. 2019; De Sarkar et al. 2021). The hadronic interaction timescale is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{loss}}=\frac{1}{n_{\mathrm{w},\mathrm{shocked}}\sigma c},\end{array}$$(20)

where n_w, shocked is the post-shock WR wind material number density, σ ≈ 3 × 10⁻²⁶ cm² is the cross section for hadronic interactions, and c is the speed of light. This retention factor modulates the injected proton spectrum, particularly at higher energies where escape becomes more probable. In denser or more extended shocks, f_ret → 1; in smaller or more diffusive regions, it drops significantly.

The fraction of the WR wind that impinges upon the WCR was determined by the wind momentum ratio, η_ϕ, given in Eq. (1), which controls the opening half angle of the conical shock surface. The shock-opening half angle, θ, was approximated in Eq. (2). The fraction of the WR wind intercepted by the WCR was then computed via the solid angle:

$$\begin{array}{c}\hfill f_{\mathrm{shock}}=\frac{2\pi (1-cos\theta )}{4\pi }·\end{array}$$(21)

This factor was used to compute the fraction of wind energy diverted into the shock.

To account for the phase-dependent size of the acceleration region, we introduced a volume scaling factor:

$$\begin{array}{c}\hfill f_{\mathrm{vol}}={\left(\frac{R_{\mathrm{sh},\mathrm{W}}}{R_{\mathrm{sep},\mathrm{peri}}}\right)}^{{\gamma }_{\mathrm{vol}}},\end{array}$$(22)

where R_sh, W reflects the extent of the WCR on the WR side, and R_sep, peri = a(1 − e) is the periastron separation used as a reference scale. A cubic exponent, γ_vol = 3, reflects the volumetric scaling of a conical or quasi-spherical acceleration region. As the stellar separation increases toward apastron, the WCR expands, enlarging the acceleration volume and enhancing proton conversion. Near periastron, the WCR contracts, reducing both the acceleration volume and the particle residence time, thereby suppressing acceleration efficiency. By incorporating f_vol alongside the conical interception fraction, f_shock, and the retention efficiency, f_ret(E), we consistently capture the phase-dependent geometry and energetics of the system in the resulting hadronic gamma-ray emission from our simplified model.

The total kinetic power of the WR wind is given by

$$\begin{array}{c}\hfill L_{\mathrm{wind},\mathrm{W}}=\frac{1}{2}{\dot{M}}_{\mathrm{W}}{V}_{\mathrm{W},\mathrm{local}}^2.\end{array}$$(23)

However, not all of this energy is available instantaneously for particle acceleration. The relevant timescale is the advection time:

$$\begin{array}{c}\hfill t_{\mathrm{adv}}=\frac{3R_{\mathrm{sh},\mathrm{W}}}{V_{\mathrm{W},\mathrm{local}}},\end{array}$$(24)

accounting for the finite residence time of the particles in the acceleration zone. The total wind energy budget at each phase is then

$$\begin{array}{c}\hfill E_{\mathrm{wind}}=L_{\mathrm{wind},\mathrm{W}}·t_{\mathrm{adv}}.\end{array}$$(25)

The injection spectrum of relativistic protons was assumed to follow a power law in energy with an exponential cutoff at the maximum energy, E_p, max,

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=\mathcal{A}{E}^{-\alpha }exp(-\frac{E}{E_{p,\max }}),\end{array}$$

with spectral index α ≈ 2.2, as is expected from standard DSA in shocks. To ensure that the total kinetic energy injected into protons was properly normalized, the normalization constant, 𝒜, was computed as

$$\begin{array}{c}\hfill \mathcal{A}=\frac{{\eta }_p·f_{\mathrm{shock}}·f_{\mathrm{vol}}·f_{\mathrm{ret}}(E)·E_{\mathrm{wind}}}{{\displaystyle \int E^{-\alpha +1}}exp(-\frac{E}{E_{p,\max }})\mathrm{d}E},\end{array}$$(26)

where η_p is the proton acceleration efficiency (typically 10⁻⁴ − 10⁻¹). Finally, the total energy injected into protons at a given orbital phase was obtained by integrating

$$\begin{array}{c}\hfill E_{p,\mathrm{inj}}={\displaystyle \int E}\frac{\mathrm{d}N}{\mathrm{d}E}\mathrm{d}E.\end{array}$$(27)

This formalism yields a proton injection energy that is naturally modulated by orbital phase, with higher values near apastron (due to increased WCR volume and residence time) and suppressed values near periastron (due to compression and escape). Notably, the effective kinetic power available for proton acceleration in our model is governed by the product f_shock ⋅ f_vol ⋅ f_ret(E), which encapsulates the geometric, volumetric, and confinement-related modulation of the injection process. This formulation closely resembles the effective injection factor, f_eff, introduced in De Becker & Raucq (2013). While f_eff is defined in a system-specific context, in our simplified model it is effectively realized through this phase-dependent product, providing a physically motivated analog that captures the evolving energetics of the wind interaction zone.

Finally, the phase-dependent hadronic γ-ray emission from γ² Velorum was computed by convolving the spectrum of injected relativistic protons with the local post-shock WR wind density, which serves as the target for hadronic interactions. At each orbital phase, the proton injection spectrum was folded with the shocked WR wind target density, n_w, shocked, to calculate the π⁰-decay γ-ray emissivity, while taking into account the nuclear enhancement factor, ϵ_enh. The resulting emission reflects both the energy distribution of the accelerated protons and the evolving ambient density conditions across the orbit. These calculations were performed using the GAMERA library (Hahn 2015; Hahn et al. 2022), which implements the hadronic interaction formalism and pion production cross sections parameterized in Kafexhiu et al. (2014).

As was discussed earlier, the O-star shock remains in the adiabatic regime throughout all orbital phases. Nevertheless, the higher unbraked velocity of the O-star wind, combined with its comparatively lower mass-loss rate, leads to a shorter advection timescale, reduced kinetic energy available for particle acceleration, and a lower post-shock density for hadronic interactions. Consequently, the O-star wind contributes negligibly to the overall gamma-ray emission when compared to the denser and more energetically dominant WR wind.

2.3. γ − γ absorption

High-energy gamma rays traversing the intense radiation fields of massive stars in the CWB system can be absorbed via pair production process (γ + γ → e⁺ + e⁻), resulting in orbital-phase-dependent attenuation. In γ² Velorum, we computed the total γ − γ optical depth, τ_γγ, by considering the finite sizes of both stars, their anisotropic blackbody photon fields, orbital inclination, and the location of the emission region at the WCR.

Both stars were treated as finite-sized blackbody emitters. The spectral number density of target photons at a given energy, ϵ, is described by the Planck distribution:

$$\begin{array}{c}\hfill n(ϵ)=\frac{2}{h^3{c}^3}·\frac{{ϵ}^2}{exp(ϵ/k_{\mathrm{B}}{T}_{\ast })-1},\end{array}$$(28)

where T_* is the effective temperature of the star, and h, c, and k_B are the usual physical constants.

The γ − γ pair production cross section was calculated using the Breit–Wheeler formula, which depends on the interaction angle, ψ, between the gamma-ray and target photon (Breit & Wheeler 1934; Gould & Schréder 1967):

$$\begin{array}{cc}& s=\frac{E_{\gamma }ϵ(1-cos\psi )}{2{(m_e{c}^2)}^2},\hfill \end{array}$$(29)

$$\begin{array}{cc}& {\sigma }_{\gamma \gamma }(s)=\frac{3{\sigma }_{\mathrm{T}}}{16}(1-{\beta }^2)[(3-{\beta }^4)ln\left(\frac{1+\beta }{1-\beta }\right)-2\beta (2-{\beta }^2)],\hfill \end{array}$$(30)

where $\beta =\sqrt{1-1/s}$, and σ_T is the Thomson cross section. The cross section is nonzero only above the pair production threshold (s ≥ 1).

The optical depth was computed by integrating along the gamma-ray path starting from the shock apex and extending outward along the line of sight. At each integration step, the angle, ψ, and distance, d, from the point to the star were used to compute the solid angle subtended by the stellar disk:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\ast }=2\pi (1-\sqrt{1-{\left(\frac{R_{\ast }}{d}\right)}^2}),\end{array}$$(31)

and the minimum target photon energy for pair production was

$$\begin{array}{c}\hfill {ϵ}_{min}=\frac{2{(m_e{c}^2)}^2}{E_{\gamma }(1-cos\psi )}·\end{array}$$(32)

The differential optical depth along the line of sight over the path length, dz, is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\tau }_{\gamma \gamma }}{\mathrm{d}z}=\frac{1-cos\psi }{4\pi }{\mathrm{\Omega }}_{\ast }{\displaystyle {\int }_{{ϵ}_{min}}^{\infty }}n(ϵ){\sigma }_{\gamma \gamma }(E_{\gamma },ϵ,\psi )\mathrm{d}ϵ,\end{array}$$(33)

which was then integrated over the path to obtain the total optical depth, τ_γγ. If the gamma-ray path intersects a stellar surface (d <  R_*), we set τ_γγ ≫ 1 to simulate full absorption. Note that the geometry of the stars and the interaction region were projected into the observer’s frame by applying an inclination angle rotation about the orbital plane.

The total optical depth is the sum of contributions from both stars:

$$\begin{array}{c}\hfill {\tau }_{\gamma \gamma }(E_{\gamma })={\tau }_{\gamma \gamma ,\mathrm{O}}(E_{\gamma })+{\tau }_{\gamma \gamma ,\mathrm{W}}(E_{\gamma }),\end{array}$$(34)

where τ_γγ, O and τ_γγ, W are the individual optical depths from the O- and WR-star photon fields, respectively. The final attenuation factor is given by exp(−τ_γγ), and varies with orbital phase and gamma-ray energy for a constant inclination angle. This treatment captures the anisotropic and phase-dependent absorption relevant for modeling the observed modulation of high-energy emission from γ² Velorum.

3. Results

In this section, we present the results of our phase-dependent hadronic emission model for γ² Velorum, which captures the evolving interplay between stellar wind geometry, proton acceleration, and γ − γ absorption across the binary orbit. Our simulations sample 20 evenly spaced orbital phases and account for phase-dependent variations in the WCR magnetic field strength, the shocked WR wind density, and the total energy injected into relativistic protons.

The key physical conditions in the acceleration region – the magnetic field strength, post-shock WR wind density, and maximum attainable proton energy – exhibit systematic variations across orbital phase due to the evolving geometry of the WCR. The magnetic field at the WCR apex increases toward periastron due to enhanced compression of the WR wind. Similarly, the post-shock WR wind density reaches its peak near periastron as the standoff distance decreases, leading to a denser target for hadronic interactions. However, the efficiency of gamma-ray production is suppressed in this phase because the compact WCR limits both the spatial volume available for particle acceleration and the retention of high-energy protons. The maximum energy of accelerated protons also varies with phase, with higher values near apastron, where a larger WCR volume and longer escape timescales facilitate more efficient acceleration. Note that the maximum proton energy calculated from our model comes out just above ∼1 TeV, which is consistent with that calculated in Reitberger et al. (2017). At each orbital phase, the total energy injected into relativistic protons was computed by combining the WR wind kinetic energy with the conical shock interception fraction, volume scaling factor, and retention efficiency. Along with the available WR wind energy budget, E_wind, the injection energy in protons, E_p, inj, exhibits a clear enhancement near apastron (ϕ ≈ 0.5). This increase in E_p, inj arises primarily from the geometric expansion of the WCR, which enhances proton conversion by increasing the acceleration volume, extending the diffusion timescale, and increasing the probability of proton interaction before escape. In contrast, the periastron phase (ϕ ≈ 0) is characterized by a compact, dense region with reduced residence time and acceleration volume, leading to a suppression of the injected proton energy budget and hadronic gamma-ray emission. The orbital variation of these parameters is presented in Fig. 1 of the paper.

Fig. 1. Variation in the magnetic field (top left), post-shock WR wind number density (top right), proton maximum energy (bottom left), and injected proton energy budget (bottom right) with orbital phases.

Using the computed proton injection spectra at each orbital phase, we evaluated the resulting hadronic spectral energy distribution (SEDs), both with and without γ − γ absorption. The required proton acceleration efficiency to match the observed Fermi-LAT SED is found to be η_p ∼ 10⁻⁴, which is somewhat lower than the value adopted in Reitberger et al. (2017), yet remains within a physically plausible range. We note that η_p is sensitive to the choice of diffusion normalization, D₀ (∼10¹⁹ cm² s⁻¹ in this paper); larger D₀ values lead to lower retention efficiency and thus require higher η_p to compensate. However, in the absence of tighter observational constraints, a broader parameter exploration is not currently feasible. The γ − γ absorption significantly suppresses the high-energy flux, particularly at periastron, where the photon field density is highest due to the smaller binary separation. This attenuation is illustrated in the left panel of Figure 2, which shows the energy-dependent transmission factor exp(−τ_γγ) across the orbital phase. The right panel of Figure 2 compares the resulting phase-dependent SEDs to Fermi-LAT observations from Martí-Devesa et al. (2020) and Pshirkov (2016). Our model, with and without γ − γ absorption, reproduces both the spectral shape and flux normalization reasonably well. While Martí-Devesa et al. (2020) interpret the Fermi-LAT spectrum as arising entirely from either apastron or phase-averaged emission, we assume that the observed γ-ray flux originates predominantly near apastron. In this scenario, the suppressed emission at periastron falls below the LAT sensitivity threshold, consistent with its non-detection during those phases. Finally, given the limited maximum energy of accelerated protons (E_p, max ∼ a few tera-electronvolts), our model shows a steep cutoff in the gamma-ray SED above ∼100 GeV, suggesting that γ² Velorum is unlikely to be detected by next-generation tera-electronvolt instruments such as the Cherenkov Telescope Array (CTA) or the Southern Wide-field Gamma-ray Observatory (SWGO).

Fig. 2. The left panel shows the attenuation factor exp(−τγγ) as a function of energy and orbital phase used in this model. The right panel shows the phase-dependent model SEDs varying across the orbital phases, along with the Fermi-LAT data points (Pshirkov 2016; Martí-Devesa et al. 2020). The CTA, SWGO, and Fermi-LAT sensitivities are also shown in the figure for comparison purposes. The Fermi-LAT sensitivity corresponds to the 10-year benchmark (P8R3_SOURCE_V3), assuming a point source with a power-law spectrum in a uniform diffuse background. It serves as a reference only and does not represent a strict detection threshold. Variations in background intensity and exposure at the location of γ2 Velorum, along with low test statistics (TS) or large uncertainties, can cause some data points to appear below this curve. The solid lines represent the unabsorbed SEDs, whereas the dashed lines signify the SEDs affected by γ − γ absorption.

Next, we computed the phase-dependent gamma-ray flux integrated in energy between 100 MeV and 500 GeV across the orbit of γ² Velorum. The full phase-energy-integrated flux, accounting for γ − γ absorption, is found to be ∼8.6 × 10⁻⁷ MeV cm⁻² s⁻¹, corresponding to a photon flux of ∼1.5 × 10⁻⁹ ph cm⁻² s⁻¹. The values are a few factors lower than those reported in Martí-Devesa et al. (2020), with the discrepancy likely arising from differences in the assumed proton spectral index, as our diffusion-limited model yields a harder slope compared to their softer assumption. To assess flux modulation across the orbit, we adopted the phase ranges ϕ = 0.25 − 0.75 (apastron) and ϕ = 0.0 − 0.25 and 0.75 − 1.0 (periastron), following Martí-Devesa et al. (2020). The resulting apastron-to-periastron flux ratio of ∼2.4 aligns well with the non-gated (∼2.3) and gated (∼2.8) values from that study, indicating that our model successfully reproduces the observed orbital modulation, which is dominated by WCR geometry, with absorption playing a secondary role.

To visually illustrate the orbital modulation, we constructed phase-binned light curves using the photon flux integrated between 100 MeV and 500 GeV. The orbital phase interval [0, 1] was divided into N = 4, 6, 8, and 20 equal bins. For each bin, we integrated the energy-integrated photon flux over the corresponding phase interval and normalized the result by the bin width, yielding the average photon flux per bin. These values were then used to generate the model gamma-ray light curves shown in Figure 3. The resulting profiles clearly exhibit a flux minimum near periastron (ϕ ≈ 0) and a maximum near apastron (ϕ ≈ 0.5), in agreement with the observed modulation reported by Martí-Devesa et al. (2020). While we did not attempt to fit the light curve data due to limited observational statistics, this exercise demonstrates that incorporating orbital-phase-dependent variations in WCR geometry is essential to reproduce the flux contrast between apastron and periastron observed in Fermi-LAT data.

Fig. 3. Variation in average photon flux per bin with orbital phase for a single orbit is shown in this figure. The orbit is divided into 4 bins (green), 6 bins (blue), 8 bins (yellow), and 20 bins (red). As usual, the periastron occurs at phase ϕ = 0, and the apastron occurs at ϕ = 0.5.

The phase-dependent behavior of gamma-ray emission contrasts sharply between γ² Velorum and η Carinae, and can be understood in terms of both wind properties and shock cooling regimes. In the moderately eccentric system γ² Velorum, the relatively low mass-loss rate leads to longer hadronic interaction timescales, making diffusion a significant factor compared to the losses. The WR shock, in this case, is radiative, resulting in moderately enhanced post-shock compression and target densities. However, near periastron, the compact WCR limits the available acceleration volume and reduces the residence time of relativistic protons, lowering the efficiency of wind energy conversion into protons, leading to the suppression of hadronic gamma-ray emission. At apastron, the expanded WCR increases both the acceleration volume and the retention efficiency, enhancing the probability of hadronic interactions and boosting the gamma-ray output. In contrast, the highly eccentric system η Carinae exhibits peak gamma-ray emission near periastron due to its extreme mass-loss rate, which generates a highly dense post-shock region. This condition significantly shortens the hadronic interaction timescale, making losses dominant over diffusion and enabling efficient emission despite the compactness of the WCR. Furthermore, the shock is expected to be much more strongly radiative in η Carinae compared to γ² Velorum, producing a higher compression ratio and an elevated post-shock density that further amplify the gamma-ray yield, though such a strongly radiative condition will also make the shock structure more complicated, which in turn could impact the acceleration process. At apastron, the increased orbital separation reduces wind densities and interaction probability, leading to suppressed emission. Thus, the key physical difference between the two systems lies not only in orbital geometry and eccentricity, but also in the prevailing shock cooling regimes and their influence on particle acceleration and target densities throughout the orbit.

Leptonic processes are not included in our model, as the strong stellar radiation field in γ² Velorum imposes severe inverse-Compton losses that limit electron energies to ∼10 MeV at the WCR apex and ∼100 MeV in the outer regions (Reitberger et al. 2017), making their contribution to giga-electronvolt emission negligible. Thus, the high-energy output is expected to be hadron-dominated. Although the absorbed nonthermal radio component suggests the presence of accelerated electrons (Purton et al. 1982; Chapman et al. 1999; De Becker & Raucq 2013), we focus on the orbital modulation of the giga-electronvolt emission, where leptonic contributions are likely subdominant.

4. Concluding remarks

We have developed a physically motivated, phase-dependent model of hadronic gamma-ray emission in the colliding-wind binary γ² Velorum, incorporating variations in WCR geometry, proton acceleration and escape, and γ − γ absorption. This framework successfully reproduces the observed SED and explains the hints of orbital variability seen in the Fermi-LAT data.

Our results indicate that the primary driver of the flux variability is the evolving geometry of the WCR with orbital phase. The widening of the WCR at apastron enhances the acceleration volume and proton retention, increasing the efficiency of hadronic interactions and boosting gamma-ray output, and vice versa during periastron. While γ − γ attenuation contributes to emission suppression near periastron, it alone cannot explain the apastron flux enhancement, underscoring the importance of intrinsic modulation.

Although based on a simple semianalytic treatment, our model demonstrates the key role of geometric effects in shaping high-energy emission from CWBs. Future efforts involving 3D MHD simulations will be essential to capture additional complexities – such as shock curvature, instabilities, magnetic field evolution, and dynamic WCR disruption – that can further influence particle acceleration and gamma-ray production. The present study provides a foundation for such advanced modeling and emphasizes the relevance of orbital geometry in interpreting the high-energy behavior of eccentric binary systems such as γ² Velorum.

Acknowledgments

A.D.S. thanks the anonymous reviewer for helpful suggestions and constructive criticism, which vastly improved the quality of this work. A.D.S. is supported by the Juan de la Cierva JDC2023-052168-I grant, funded by MICIU/AEI/10.13039/501100011033 and the FSE+. This work is also supported by the grant PID2021-124581OB-I00 of MCIU/AEI/10.13039/501100011033 and 2021SGR00426 and the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and by MCIU with funding from European Union NextGeneration EU (PRTR-C17.I1).

References

All Tables

All Figures

	Fig. 1. Variation in the magnetic field (top left), post-shock WR wind number density (top right), proton maximum energy (bottom left), and injected proton energy budget (bottom right) with orbital phases.
In the text

Fig. 2. The left panel shows the attenuation factor exp(−τγγ) as a function of energy and orbital phase used in this model. The right panel shows the phase-dependent model SEDs varying across the orbital phases, along with the Fermi-LAT data points (Pshirkov 2016; Martí-Devesa et al. 2020). The CTA, SWGO, and Fermi-LAT sensitivities are also shown in the figure for comparison purposes. The Fermi-LAT sensitivity corresponds to the 10-year benchmark (P8R3_SOURCE_V3), assuming a point source with a power-law spectrum in a uniform diffuse background. It serves as a reference only and does not represent a strict detection threshold. Variations in background intensity and exposure at the location of γ2 Velorum, along with low test statistics (TS) or large uncertainties, can cause some data points to appear below this curve. The solid lines represent the unabsorbed SEDs, whereas the dashed lines signify the SEDs affected by γ − γ absorption.

In the text

	Fig. 3. Variation in average photon flux per bin with orbital phase for a single orbit is shown in this figure. The orbit is divided into 4 bins (green), 6 bins (blue), 8 bins (yellow), and 20 bins (red). As usual, the periastron occurs at phase ϕ = 0, and the apastron occurs at ϕ = 0.5.
In the text