DYNAMICAL MASS CONSTRAINTS ON THE ULTRALUMINOUS X-RAY SOURCE NGC 1313 X-2

Based on the observed sinusoidal light curves and the black hole's radial velocity amplitude for NGC 1313, the dynamical parameters, including the black hole mass and the secondary properties, can be derived with the X-ray irradiated binary model as described in Section 2. For this purpose, we define $\chi ^2 = \sum _i{(m_i - m_d)^2/e_i^2} + (P_d - P)^2/\sigma _P^2 + (K^\prime _d - K^\prime)^2/\sigma ^2_{K^\prime }$. Here, m_i and e_i are the observed HST magnitudes and errors, P = 6.12 days and σ_P = 0.16 days are the detected period and error (Liu et al. 2009), 2K' = 200 km s⁻¹ and $\sigma _{K^\prime }=40$ km s⁻¹ are the radial velocity semi-amplitude of the black hole and error based on the Gemini observations (Roberts et al. 2011; J. C. Gladstone et al. 2011, in preparation), while m_d is the predicted magnitude at the observed phase, and P_d and K'_d are the corresponding quantities predicted by the model. Five F555W observations are excluded from the fit because they exhibited abnormally red colors and abnormally large deviations from the sinusoidal light curve (Liu et al. 2009). The Padova stellar model grid (Girardi et al. 2002), with Z = 0.1 Z_☉ as appropriate for NGC 1313 (Ryder 1993), is used to determine the secondary properties with two parameters, the age t and the initial mass m_ini = m_f · m_{ini, max} with m_f defined as a fraction of the maximum initial mass allowed for that age. In addition to the 10 parameters to describe the X-ray irradiated binary model, we assume an extinction E(B − V) with the standard galactic extinction law to convert model magnitudes to observed magnitudes.

Acceptable models are searched over the 11-parameter phase space by minimizing the above-defined χ² with 26 degrees of freedom (37 constraints and 11 free parameters). To not miss any acceptable models, we start with very loose constraints on the parameters, with 3 < M_•/M_☉ < 3000, $6.6< \lg (t) < 8.4$, 0.3 < m_f < 1, 0.05 < f_R < 1, 0 < Θ < 90, 30 < L_X/10³⁸ < 300, 1 < X_h < 10, 0.05 < X_d < 0.9, 0.01 < β < 0.99, 0.01 < γ < 0.99, and 0.11 < E(B − V) < 0.44. Here the extinction value is bounded within 0.11 mag (Galactic extinction) and 0.44 mag, i.e., n_H = 2.7 × 10²¹ cm⁻² as inferred from the fit to the X-ray spectrum (Miller et al. 2003). The X-ray luminosity is bounded within 3 × and 30 × 10³⁹ erg s⁻¹ with an average of ∼4 × 10³⁹ erg s⁻¹ based on previous X-ray observations (Mucciarelli et al. 2005). The average X-ray luminosity corresponds to an accretion rate of $\dot{M} = 2\times 10^{-7} (L/10^{39}) (0.1/\eta) \sim 8\times 10^{-7} \, M_\odot \, {\rm yr^{-1}}$ from $L_{{\rm acc}} = \eta \dot{M} c^2$ regardless the black hole mass. Based on the accretion rate consideration, models are discarded unless the secondary is filling (>98% of) its Roche lobe, or the secondary has already evolved off the main sequence into its red/blue giant or asymptotic giant phases and fills a significant fraction of its Roche lobe, so that the heavy stellar winds can be focused onto the primary to provide the required accretion rate.

Given the vast parameter space, the χ² minimization is an extremely difficult problem that demands carefully designed intensive computations, which are carried out using the Odyssey cluster at Harvard University. As the first step, 10⁸ models are randomly drawn over the 11-dimensional phase space. This is equivalent to an average of 10^8/11 ∼ 5.3 sampled values over the parameter range for each dimension, but with the advantage of continuous coverage due to randomization. Most of the 10⁸ models have extremely large χ² values, including about 20 models with χ² < 260 (i.e., χ²_ν < 10), and about 25,000 models with χ² < 10, 000. As shown in Figures 2(a)–(f), M_•, Θ, $\lg (t)$, m_f, and E(B − V) are clearly constrained even at the χ² < 10, 000 level, while L_X, X_h, X_d, f_R, β, and γ are not constrained. At the χ² < 10, 000 level, models with black hole masses above 1000 M_☉ or inclination angles below 20° are not acceptable. Only certain combinations of the secondary age and initial mass in the two thin stripes are possibly acceptable. Secondaries in the lower stripe are starting the helium core and hydrogen shell burning phase, and secondaries in the upper stripe are starting the carbon core and helium shell burning phase, while all have similar average densities appropriate for the detected orbital period ($P = \sqrt{110/\rho }$ with ρ as the mean matter density within the secondary's Roche lobe). While we start with a uniform distribution of E(B − V) within 0.11 mag and 0.44 mag, it is found that models with larger E(B − V) are more likely to be unacceptable.

Zoom In Zoom Out Reset image size

Local χ² minima are then searched using the AMOEBA method as described in Numerical Recipe (Press et al. 2003) adapted to computing with the Odyssey cluster. Selected as starting points are about 6000 models as shown in Figures 2(a)–(f), including all models with χ² < 2000 and those with lowest χ² within fine intervals of each parameter. This large number of initial guesses provides a fine coverage of the parameter space, and yet is still doable with the available computing power. For each AMOEBA search, we construct the initial simplex with a length scale of 20% of the whole range for each parameter, in hope to have a complete coverage of the whole parameter space. The simulated annealing method is applied in the first 300 model evaluations of an AMOEBA search in order to increase the chances to obtain a "global" χ² minima in the large parameter space to be searched. The AMOEBA search claims a χ² minima if the fractional tolerance ftol <10⁻⁵, and stops if the iteration number exceeds 1000. To find the "true" χ² minima, we restart another AMOEBA search from a claimed χ² minima, but with a 4 × smaller length scale. In the end, two rounds of AMOEBA searches with about eight million model evaluations lead to local χ² minima from χ² ∼ 52 to ⩾1000. As shown in Figure 3, models with black hole masses above 100/300 M_☉ are excluded at the χ²_ν = 3/4 level. Due to the large photometric fluctuations present in the observations, even the "best" model has χ²_ν > 2. Therefore, we do not claim one "true" model for NGC 1313 X-2, but instead consider models with χ² < 78 (i.e., χ²_ν ⩽ 3) as "acceptable" models.

Zoom In Zoom Out Reset image size

The acceptable models are grouped into three categories based on the secondary properties as shown in Figures 4(a)–(e). The secondaries in the first group of models are very young (⩽5 Myr) and massive in the carbon core and helium shell burning phase filling about 20% of their Roche lobes. The black hole mass and the inclination angle cluster tightly in small ranges of (30–40 M_☉, 80°–90°). The X-ray luminosities are within 1–2 × 10⁴⁰ erg s⁻¹ (about 2–6 times the Eddington luminosity, or ∼2–6 × L_E), but the accretion disks are quite small (<15% of the tidal-disruption radius) and as thin as the standard thin disk. The accretion disk albedo is in the range of 0–0.5, and the X-ray-to-thermal conversion factor is above 0.8. The best model in this group (model A1, χ² = 70.88) contains a 45.3 M_☉ black hole and a young (5 Myr) and massive (initially 44 M_☉, currently 17 M_☉, 5.1 R_☉) secondary that fills 20% of its Roche lobe. As shown in Figure 5(a), the optical spectrum is dominated by the hot (∼70, 000 K) secondary, which is only slightly affected by the X-ray heating. Although the X-ray luminosity is high (L_X ∼ 1.6 × 10⁴⁰ erg s⁻¹ ∼3 × L_E), the optical light from the irradiated accretion disk is only 1% of that from the secondary in the B/V bands because the accretion disk is small (X_d ∼ 0.05) and thin (X_h ∼ 1). The binary system is viewed almost edge-on (θ ∼ 82°), and the predicted sinusoidal light curves have amplitudes of 0.04/0.04 mag in the B/V bands (Figure 5(d)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The secondaries in the second group are 40–50 Myr old in the helium core and hydrogen shell burning phase filling about 40%–80% of their Roche lobes. As compared to the first group, the secondaries are 10 times older and less massive. The black hole mass and the inclination angle occupy much larger parameter ranges of (3–100 M_☉, 20°–80°). The X-ray luminosities are 3–30 × 10³⁹ erg s⁻¹ (3–20 × L_E), and the accretion disks can be rather large (>20% of the tidal-disruption radius) and about 2–6 times thicker than the standard thin disk. The accretion disk albedo and the X-ray-to-thermal conversion factor on the secondary surface range both from 0.01 to 0.99. The best model in this group (model B1, χ² = 54.94) contains a 15.3 M_☉ black hole and an old (40 Myr) intermediate-mass (initially 8.08 M_☉, currently 8.07 M_☉, 9.0 R_☉) secondary that fills 54% of its Roche lobe. The X-ray luminosity is L_X ∼ 1.2 × 10⁴⁰ erg s⁻¹, or ∼6 × L_E. The secondary has a surface temperature of ∼18, 500 K and is significantly affected by the X-ray heating, leading to sinusoidal light curves with amplitudes of 0.15/0.15 mag in the B/V bands for the inclination angle of θ ∼ 38° (Figure 5(d)). The accretion disk is 55% of the tidal-disruption radius and 3.5 times thicker than the standard thin disk, and the X-ray heated disk contributes roughly equally as the secondary in the B/V bands (Figure 5(b)).

The secondaries in the third group are 80–160 Myr old in the helium core and hydrogen shell burning phase filling (>98% of) their Roche lobes. As compared to the second group, the secondaries are 2–3 times older and less massive. The black hole mass and the inclination angle cluster tightly in small ranges of (3–10 M_☉, 50°–70°). The X-ray luminosities are within 1–2 × 10⁴⁰ erg s⁻¹ (or 7–20 × L_E), and the accretion disks can be rather large (40%–70% of the tidal-disruption radius) and about 6–9 times thicker than the standard thin disk. The accretion disk albedo is constrained within 0.7–0.8, while the X-ray-to-thermal conversion factor on the secondary surface ranges from 0.01 to 0.8. The best model in this group (model F1, χ² = 61.26) contains a 6.7 M_☉ black hole and a very old (150 Myr) and rather low-mass (initially 4.28 M_☉, currently 4.27 M_☉, 10.6 R_☉) secondary that fills its Roche lobe. The X-ray luminosity is L_X ∼ 1.9 × 10⁴⁰ erg s⁻¹ (∼20 × L_E), the accretion disk is 40% of the tidal-disruption radius and seven times thicker than the standard thin disk. The X-ray heated disk contributes about four times more light than the secondary in the B/V bands (Figure 5(c)). The secondary has a surface temperature of ∼9500 K and is significantly affected by the X-ray heating, leading to sinusoidal light curves with amplitudes of 0.17/0.17 mag in the B/V bands for the inclination angle of θ ∼ 56° (Figure 5(d)).

The model χ² for assumed 2K' = 100 ± 40 km s⁻¹ is defined in the same way as described in Section 3.1, and χ² minima are sought in a large phase space with two rounds of AMOEBA searches following the same procedures. The resulted models are grouped into three categories based on the secondary properties as in Section 3.1, and the "acceptable" models with χ² < 78 are plotted in Figures 6(a)–(f). As compared to the case of 2K' = 200 ± 40 km s⁻¹, the first group of models with secondaries in their carbon core and helium shell burning phase are no longer "acceptable." Even the best model (model A2) for the first group has χ² = 84.0, an apparently bad fit to the observed light curve as shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Models in the second group occupy similar phase space as compared to the case of 2K' = 200 ± 40 km s⁻¹, while the best model is quite different. As shown in Figures 6 and 7, the best model (model B2, χ² = 51.9) contains a 39 M_☉ black hole and a 50 Myr old intermediate-mass secondary (initially 7.18 M_☉, currently 7.16 M_☉, 8.2 R_☉) in its helium core and hydrogen shell burning phase filling 51% of its Roche lobe. The X-ray luminosity is L_X ∼ 1.2 × 10⁴⁰ erg s⁻¹, or ∼2.4 × L_E. The accretion disk is 15% of the tidal-disruption radius and five times thicker than the standard thin disk, and the X-ray heated disk contributes two times that of the secondary in the B/V bands. The secondary has a surface temperature of ∼17, 500 K and is significantly affected by the X-ray heating, leading to sinusoidal light curves with amplitudes of 0.155/0.155 mag in the B/V bands for the inclination angle of θ ∼ 45° (Figure 7).

Models in the third group occupy roughly similar phase space as compared to the case of 2K' = 200 ± 40 km s⁻¹, while the secondaries are older and less massive. The best model (model F2, χ² = 55.8) contains an 8.2 M_☉ black hole and a 200 Myr old secondary (initially 3.75 M_☉, currently 3.74 M_☉, 9.9 R_☉) in its helium core and hydrogen shell burning phase filling its Roche lobe. The X-ray luminosity is L_X ∼ 9 × 10³⁹ erg s⁻¹ (∼8.5 × L_E), the accretion disk is 30% of the tidal-disruption radius and six times thicker than the standard thin disk. The X-ray heated disk contributes about four times that of the secondary in the B/V bands. The secondary has a surface temperature of ∼8700 K and is significantly affected by the X-ray heating, leading to sinusoidal light curves with amplitudes of 0.156/0.156 mag in the B/V bands for the inclination angle of θ ∼ 38° (Figure 7).

The model χ² for assumed 2K' = 25 ± 15 km s⁻¹ is defined in the same way as described in Section 3.1, and χ² minima are sought in a large phase space with two rounds of AMOEBA searches following the same procedures. The resulted models are grouped into three categories based on the secondary properties as in Section 3.1, and the "acceptable" models with χ² < 78 are plotted in Figures 8(a)–(f). As compared to the case of 2K' = 200 ± 40 km s⁻¹, the first group of models with secondaries in their carbon core and helium shell burning phase are no longer "acceptable" like the case of 2K' = 100  ±  40 km s⁻¹. Even the best model (model A3) for the first group has χ² = 96.5, an apparently bad fit to the observed light curve as shown in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Models in the second group occupy similar phase space as compared to the case of 2K' = 200 ± 40 km s⁻¹, except that the black hole mass can be above 1000 M_☉. The best model (model B3, χ² = 52.43) contains a 275 M_☉ black hole and a 40 Myr old intermediate-mass secondary (initially 8.09 M_☉, currently 8.08 M_☉, 7.9 R_☉) in its helium core and hydrogen shell burning phase filling 44% of its Roche lobe. The X-ray luminosity is L_X ∼ 1.3 × 10⁴⁰ erg s⁻¹, or ∼0.4 × L_E. The accretion disk is 12% of the tidal-disruption radius and 4 times thicker than the standard thin disk, and the X-ray heated disk contributes 1.7 times that of the secondary in the B/V bands. The secondary has a surface temperature of ∼20, 000 K and is significantly affected by the X-ray heating, leading to sinusoidal light curves with amplitudes of 0.157/0.157 mag in the B/V bands for the inclination angle of θ ∼ 73° (Figure 9).

Models in the third group occupy roughly similar phase space as compared to the case of 2K' = 200 ± 40 km s⁻¹, while the secondaries are older and less massive. The best model (model F3, χ² = 53.89) contains a 6.1 M_☉ black hole and a 220 Myr old secondary (initially 3.58 M_☉, currently 3.57 M_☉, 9.8 R_☉) in its helium core and hydrogen shell burning phase filling its Roche lobe. The X-ray luminosity is L_X ∼ 1.6 × 10⁴⁰ erg s⁻¹ (∼20 × L_E). The secondary has a surface temperature of ∼8200 K and is significantly affected by the X-ray heating, leading to sinusoidal light curves with amplitudes of 0.156/0.156 mag in the B/V bands for the inclination angle of θ ∼ 15° (Figure 9). The accretion disk is 28% of the tidal-disruption radius and six times thicker than the standard thin disk, and the X-ray heated disk contributes about two times that of the secondary in the B/V bands.

The best models for the three cases are listed in Table 1. To summarize, there are three groups of acceptable models based on the secondary properties depending on whether the secondary is in the carbon core and helium shell burning phase or in the helium core and hydrogen shell burning phase, and whether the secondary is filling the Roche lobe or just part of the Roche lobe. The first group of models with secondaries in the carbon core and helium shell burning phase, while acceptable for the highest 2K' = 200  ±  40 km s⁻¹, become unacceptable when 2K' drops. Models of the other two groups, with secondaries in the helium core and hydrogen shell burning phase, are acceptable for all cases of assumed 2K', but the black hole mass can become as massive as 1000 M_☉ as the assumed 2K' becomes lower as shown in Figure 10. We note that, while the extinction value is started to range from 0.11 mag to 0.44 mag, it is always close to E(B − V) = 0.11 mag for all three groups of "acceptable" models regardless the assumed 2K' as shown in Figures 4(f), 6(f), and 8(f).

Zoom In Zoom Out Reset image size

We follow the formulation of Kopal (1959) to compute the Roche lobe for given M_•, M_*, and a (i.e., primary and secondary masses and separation). In our computation, we define a Cartesian coordinate system with a as the unit length to place the secondary at (0, 0, 0) and the primary at (1, 0, 0). The Roche potential at each point (x, y, z) is

$$\begin{eqnarray*} &&\Phi (x,y,z) = -GM_*/r_1 - GM_\bullet /r_2 - G(M_\bullet +M_*)r^2/2, \end{eqnarray*}$$

where $r_1 = \sqrt{x^2+y^2+z^2}$, $r_2 = \sqrt{(x-1)^2+y^2+z^2}$, and r² = (x − M_•/(M_• + M_*))² + y². The L1 point is located at (x₁, 0, 0), with x₁ solved numerically from x⁻¹₁ − x₁ = M_•/M_*[(1 − x₁)⁻² − (1 − x₁)]. To compute the Roche lobe, we set up a polar coordinate system (r, θ, ϕ), and solve numerically for each of the θ, ϕ grid the radius r at which the Roche potential is the same as at the L1 point. The effective radius R_{L, *} (in unit of a) for the Roche lobe around the secondary is computed from the volume (V = 4πR³_{L, *}/3) enclosed by the equipotential surface. The effective radius R_{L, •} is also computed for the Roche lobe around the primary black hole.

A secondary with mass M_*, radius R_*, surface temperature T_*, and gravity g_* (=GM_*/R²_*) is filled into its Roche lobe by equating a · R_{L, *}f_R = R_*, with f_R as the fraction the secondary fills its Roche lobe. This immediately sets the physical dimension for a and the orbital period with

$$\begin{eqnarray*} &&P^2 = 2\pi a^3/G(M_\bullet +M_*). \end{eqnarray*}$$

The gravity at each surface element of the Roche lobe can be computed as

$$\begin{eqnarray*} &&\vec{g} \equiv g\hat{g} = \nabla \Phi (x,y,z) = ({\partial /\partial x},{\partial /\partial y},{\partial /\partial z})\Phi (x,y,z) \end{eqnarray*}$$

and $-\hat{g}$ is the surface normal vector. This gravity varies over the stellar surface. It is higher near the poles and lower near the equators, and the lowest at the parts of the star nearest to the L1 point. In hydrostatic equilibrium, the gravity variation modifies the local effective temperature over the stellar surface. To describe this gravitational darkening (brightening) effect, we adopt the scaling relation T_g = T_*(g/g_*)^b with b = 1/4 for radiative stellar atmospheres (von Zeipel 1924) and b = 0.08 for convective atmospheres (Lucy 1967).

The accretion disk around the black hole is described as a standard accretion disk with the γ prescription (Shakura & Sunyaev 1973; Frank et al. 2002, hereafter FKR2002). For such an accretion disk with radius r, the inner radius r_i = 6GM_•/c² assuming a non-spinning black hole, the luminosity $L_d = GM_\bullet \dot{M_\bullet }/2r_i$, and the temperature profile T_a(r) = T₀(r_i/r)^3/4f with f⁴ = 1 − (r_i/r)^1/2 and T⁴₀ = 3L_d/4πσr²_i. To avoid tidal disruption, the disk must be within the tidal radius R_Tides≈0.9R_{L, •}. We assume the outer radius r_o = X_r · R_Tides with X_r ⩽ 1. The disk can be divided into the inner disk, the middle disk, and the outer disk, with the transitions at $r_{{\rm im}} = 3.6\times 10^2 \gamma ^{2/21} \dot{m}^{16/21} m^{2/21} f^{64/21}$ and at $r_{{\rm mo}} = 2.6\times 10^4 \dot{m}^{2/3}f^{8/3}$. Here, m = M_•/M_☉, $\dot{m} = \dot{M}_\bullet /\dot{M}_{{\rm Edd}}$, $\dot{M}_{{\rm Edd}} = L_{{\rm Edd}}/0.1c^2$, and L_Edd = 4πGM_•c/κ = 1.3 × 10³⁸m erg s⁻¹. The height of the disk $h(r) = 7.5\dot{m} f^4$ for the inner disk, $h(r) = 1.6\times 10^{-2}\gamma ^{-1/10}\dot{m}^{1/5}m^{-1/10}r^{21/20}f^{4/5}$ for the middle disk, and $h(r) = 7.2\times 10^{-3}\gamma ^{-7/10}\dot{m}^{3/20}m^{-1/10}r^{9/8}f^{3/5}$ for the outer disk. In our model, we assume the "real" disk height H(r) = X_h · h(r) to account for the fattened disk at high X-ray luminosities. With the above formulae, we compute the surface temperature T(r), height H(r), and surface normal vector $\hat{s}(\vec{r})$ for the disk under the reasonable assumption that L_X = L_d, i.e., X-ray dominates the disk radiation in ULXs.

X-ray irradiation of the accretion disk is treated following FKR2002 (Section 5.10). The X-ray flux crossing the disk surface at $\vec{r}$ is

$$\begin{eqnarray*} &&F = {L_X\over 4\pi d^2} (1-\beta) \cos \psi, \end{eqnarray*}$$

with $\vec{d} = d\hat{d} = \hat{d}\sqrt{r^2+H^2(r)}$ as the vector from the black hole to the surface element, β as the albedo, and $\cos \psi = -\hat{s}\cdot \hat{d}$ representing the angle between the direction of the incident radiation and the disk surface normal vector. The effective temperature resulting from irradiation can be defined as

$$\begin{eqnarray*} &&T_{{\rm Irr}}^4 = F/\sigma = {1\over 3} T_0^4 (r_i/d)^2 (1-\beta) \cos \psi \end{eqnarray*}$$

and the total temperature for the disk surface element

$$\begin{eqnarray*} &&T_d^4 = T_a^4(r) + T_{{\rm Irr}}^4. \end{eqnarray*}$$

The X-ray irradiation of the secondary is treated in a similar fashion. With respect to the black hole, a stellar surface element at (x, y, z) raises above the orbital plane with an angle of $\delta = \arcsin {z \over d}$ to the black hole with $\vec{d} = d\hat{d} = (x-1,y,z)$, while the accretion disk edge at the outer radius r_o raises above the orbital plane with an angle of $\Delta = \arctan {H(r_o)\over r_o}$. The effective temperature from irradiation is defined as

$$\begin{eqnarray*} &&T_{{\rm Irr},*} = {1\over 3} Q T_0^4 (r_i/d)^2 \gamma \cos \psi, \end{eqnarray*}$$

with $\cos \psi = -\hat{d}\cdot -\hat{g}$, $-\hat{g}$ as the surface normal vector, Q as the transmission function, and γ representing the fraction of intercepted X-ray flux converted to thermal energy. The stellar surface element with cos ψ < 0 is blocked by the star itself and thus not irradiated. The transmission function describes how the X-ray is blocked by the accretion disk, with Q = 0 if δ < Δ (i.e., completely blocked), and $Q=1 - e^{\delta -\Delta \over \Delta }$ if δ > Δ to account for the photoelectric absorption by the extended disk atmosphere. The total temperature at each element of the secondary surface is then computed as

$$\begin{eqnarray*} &&T_s^4 = T_g^4 + T_{{\rm Irr},*}^4. \end{eqnarray*}$$

Whether a surface element of the accretion disk or the secondary is visible to the observer at an inclination of Θ is determined as follows. For each binary phase Φ with the observing vector $\hat{v}(\Theta,\Phi)$ from the binary to the observer, treated first is the component closer to the observer, which is the accretion disk at $-{\pi \over 2} < \Phi < {\pi \over 2}$ and the secondary at ${\pi \over 2} < \Phi < {3\pi \over 2}$ under the convention that Φ = 0 if the accretion disk is right in between the secondary and the observer. For the accretion disk, the disk surface may for some Θ be blocked by the disk edge, which is treated as a cylinder with a radius of r_o, a half-height of H(r_o), a surface temperature of T_a(r_o), and a surface normal vector $\hat{m}$. We thus project the visible parts of the disk edge (with $\hat{m}\cdot \hat{v} \ge 0$) onto a plane normal to $\hat{v}$, and compute its silhouette as a polygon. For the disk surface elements, visible are those with $\hat{s}\cdot \hat{v} \ge 0$ and outside the edge's silhouette. A disk silhouette is computed for the disk edge and surface as a whole. For the secondary, we examine the surface elements from pole to pole, and compute a silhouette polygon of all previously examined elements for each latitude. A secondary surface element is visible if it is outside the previous silhouette polygon and with $-\hat{g} \cdot \hat{v} \ge 0$. We then treat the component farther away from the observer in the same way, except that a surface element is invisible if it is within the silhouette polygon for the closer component. The ray-tracing method developed by Haines (1994) is used to check whether a surface element, after projecting to the plane normal to $\hat{v}$, is within a silhouette polygon.

Both the spectrum and the wide-band photometry are computed for the emergent emission from the X-ray irradiated binary model. For the ith visible surface element, the flux is computed as

$$\begin{eqnarray*} &&\Delta F_i(\lambda) = I_i(\lambda) {A_i\cos \Psi _i\over D^2}, \end{eqnarray*}$$

with I_i(λ) as the specific intensity, A_icos Ψ_i as the surface element area projected along the line of sight, and D as the distance to the observer. Here, $\cos \Psi = \hat{v}\cdot -\hat{g}$ for the secondary, $\cos \Psi = \hat{v}\cdot \hat{s}$ for the disk surface, and $\cos \Psi = \hat{v}\cdot \hat{m}$ for the disk edge. The disk edge and surface are treated simply as a blackbody, and I_i(λ) = B_λ(T_d). The secondary is treated as Kurucz's stellar atmosphere model (Kurucz 1993), which gives the emergent intensity spectrum I(λ; μ, T, g) for a grid of the incident angle μ, surface temperature T, and gravity g (ref). In our model, the intensity spectrum I_i(λ) = I(λ; μ_i, T_{s, i}, g_i) is interpolated from the model grid with $\mu = \hat{v}\cdot -\hat{g}$. Summing up all visible surface elements of the secondary and the accretion disk, we obtain the total flux spectrum as

$$\begin{eqnarray*} &&F(\lambda) = \sum _i \Delta F_i(\lambda). \end{eqnarray*}$$

The flux observed through a filter is computed as

$$\begin{eqnarray*} &&F = \int F(\lambda) T(\lambda) d\lambda \end{eqnarray*}$$

and the magnitude is computed as

$$\begin{eqnarray*} &&m = -2.5\lg (F/F_0). \end{eqnarray*}$$

Here, T(λ) is the filter transmission, and F₀ is the zeroth magnitude flux, which are taken from the stsdas.synphot package. In our model illustration, we calculate the magnitudes for Johnson UBV, Cousins RI, and Bessel JHK, but any number of filters can be easily added into the code.

Given the model parameters, the model first sets up the binary geometry, computes the secondary and disk temperature profiles with and without X-ray irradiation, then computes the emergent spectrum and the magnitudes for a binary phase Φ. For each binary phase, we put the model output in one figure as shown in Figure 1 for a base model with M_• = 30 M_☉, a B0V secondary (M_* = 17.5 M_☉, R_* = 7.4 R_☉, and T_* = 30,000 K with solar abundance) filling its Roche lobe (f_R = 1), L_X = 10³⁹ erg s⁻¹, X_h = 1, X_d = 0.5, β = 0.5, γ = 0.5, and Θ = 60°. In such a figure, the upper left panel shows the binary components projected on the plane normal to the viewing vector $\hat{v}$, and the upper right panel shows the UBVRIJHK light curves for a whole period with the magnitudes marked for this phase. The lower panel shows the emergent spectra from the secondary with and without X-ray irradiation, from the disk itself and the X-ray irradiation contribution, and from the binary as a whole. As an illustrative product for a given model, we create an animated movie by stacking up a sequence of such figures for binary phases over a period.

To illustrate how the model depends on the key parameters, we start from the base model shown in Figure 1. In such a model, the star contributes 70%–60% of the light in the U − I bands, the un-X-ray-irradiated disk contributes 15%–20%, and the X-ray irradiation of the disk contributes the rest 15%–20%. X-ray heating of the L1 side of the secondary boosts it to be the hottest side, yielding a sinusoidal light curve instead of an ellipsoidal one. We vary one parameter in its reasonable range with all other parameters fixed at the base model, and investigate the according changes in the spectrum, the magnitudes, the light curve shapes, and variation amplitudes with respect to the base model.

We first investigate the five parameters on X-ray irradiation. The X-ray luminosity is varied from 10³⁸ to 10⁴¹ erg s⁻¹. With increasing X-ray luminosity, the X-ray heated disk and the L1 side of the secondary become hotter and brighter. Accordingly, the light curve shape changes from ellipsoidal to sinusoidal, with brighter mean magnitudes and larger variation amplitudes. The effects are illustrated in Figure 12.

Zoom In Zoom Out Reset image size

The disk fattening factor is varied from 0.1 to 10. With increasing fattening, the disk becomes thicker and intercepts more X-ray radiation, thus becomes slightly brighter. The disk casts a larger shadow on the secondary, and the secondary becomes less X-ray heated. Accordingly, the light curve shape changes from sinusoidal to ellipsoidal, with slightly brighter mean magnitudes but smaller variation amplitudes. The effects are illustrated in Figure 13.

Zoom In Zoom Out Reset image size

The disk size factor is varied from 0.01 to 1. With a larger radius and emitting area, the disk contributes slightly more in the optical, leading to slightly brighter mean magnitudes for the light curves. The slightly larger disk shadow on the secondary reduces the X-ray heating only slightly, and the light curve shape remains sinusoidal yet with smaller variation amplitudes. The effects are illustrated in Figure 14.

Zoom In Zoom Out Reset image size

The albedo for the disk surface is varied from 0.01 to 0.99. The increasing albedo reduces the heating of the disk, leading to much smaller contribution from the X-ray irradiation of the disk. X-ray heating of the secondary remains the same, and the light curve shape remains sinusoidal. The light curve mean magnitudes decrease slightly, and the variation amplitudes increasing slightly, both due to the diminishing of the X-ray irradiation of the disk. The effects are illustrated in Figure 15.

Zoom In Zoom Out Reset image size

The X-ray-to-thermal conversion factor is varied from 0.01 to 0.99. We do not link this factor to the disk albedo because the secondary albedo may be different due to temperature and ionization differences, and some of the X-ray energy may be converted to kinetic energy of the materials blown away from the secondary surface. With increasing conversion efficiency, the L1 side of the secondary becomes hotter. The light curve shape changes from ellipsoidal to sinusoidal with slightly brighter mean magnitudes and much larger variation amplitudes. The effects are illustrated in Figure 16.

Zoom In Zoom Out Reset image size

Now we investigate the parameters that affect the binary configuration. The black hole mass is varied from 3 M_☉ to 3000 M_☉. With the increasing black hole mass, the Roche lobe of the black hole gets larger, the separation gets larger, and the secondary star is less distorted and less gravitationally darkened. The heating of the secondary becomes less and less significant so that the light curves change from sinusoidal to ellipsoidal as illustrated in Figure 17(a). Because the secondary star is less gravitationally darkened at higher black hole mass, the variation amplitudes for the light curves become smaller and smaller (Figure 17(b)). With a fixed disk size factor, the physical size of the disk gets larger at higher black hole mass, and its contribution in the optical gets higher as illustrated in Figure 17(a). This also leads to higher total optical light, and brighter mean magnitudes of the light curves (Figure 17(b)). At black hole mass ⩽10 M_☉, the separation becomes so small that the secondary star and the disk begin to block each other when they line up, leading to extra shallow drops in the light curves around phases 0° and 180° (Figure 17(a)).

Zoom In Zoom Out Reset image size

The fraction the secondary fills its Roche lobe is varied from 0.01 to 1. With a smaller filling fraction and a fixed secondary, the Roche lobe of the secondary becomes larger. Subsequently, the separation, the Roche lobe of the black hole, and the physical size of the accretion disk all become larger. The larger disk contributes more in the optical, leading to brighter mean magnitudes of the light curves, more so at longer wavelengths. X-ray heating of the secondary becomes less significant, and the light curve shape changes from sinusoidal to ellipsoidal with reduced variation amplitudes. The tidal distortion and gravitational darkening of the secondary also become less significant, adding to the decrease of the variation amplitudes. The effects are illustrated in Figure 18.

Zoom In Zoom Out Reset image size

The inclination angle is varied from 0° (face-on) to 90° (edge-on). With the inclination angle increasing, the projected area of the disk becomes smaller, leading to smaller disk contribution in the optical, and dimmer mean magnitudes of the light curves. The projected area of the L1 side of the secondary also becomes larger, leading to larger variation amplitudes for the sinusoidal light curves. When the inclination angle increases to 70°, the disk begins to block parts of the secondary, which enhances the contrast between the heated side and its antipodal side, and leads to much larger variation amplitudes. When the inclination angle increases to above 80°, the disk begins to block the hottest part of the L1 side of the secondary, which lowers the contrast between the heated side and its antipodal side, and leads to reduced variation amplitudes. The effects are illustrated in Figure 19.

Zoom In Zoom Out Reset image size

We experiment with the secondary properties by varying spectral types and luminosity classes as listed in Table 2. The quoted mass, radius, and temperature are taken from the MK spectral type calibration (Cox 2001) for stars in the solar neighborhood. For main-sequence stars, later spectral types are less massive and smaller in size, with higher densities, lower temperatures, and lower luminosities. For an O star secondary, the separation is the largest, and X-ray heating of the secondary is a minor effect, leading to an ellipsoidal light curve. X-ray heating becomes significant for a B star secondary, and the light curve becomes sinusoidal with much larger variation amplitude. For still later stellar types, the variation amplitudes become smaller despite the smaller separations and more significant X-ray heating effects, because the accretion disk contribution to the optical light exceeds the star contribution by larger and larger factors. The effects are illustrated in Figure 20.

Zoom In Zoom Out Reset image size

For supergiant stars, later spectral types are less massive yet larger in size, with lower densities and lower temperatures. For a B5I supergiant with a high temperature, X-ray heating of the L1 side of the secondary is not enough to raise it to be the hottest side, and the light curve is an ellipsoidal one. For later spectral types with much lower temperatures, X-ray heating is enough to raise the L1 side of the secondary to be the hottest side, and the light curves become sinusoidal. The physical disk size, along with the Roche lobe of the black hole, becomes larger for later spectral types. Because of its low temperature, the expanded outer disk does not contribute much to the optical at shorter wavelengths (U/B bands), but contributes significantly at longer wavelengths. In the U/B bands, the mean magnitudes of the light curves become dimmer for later spectral types because the secondary becomes dimmer while the disk contribution remains largely the same. In the V − K bands, the mean magnitudes become brighter for later spectral types, more so in the redder bands, because both the disk and the secondary become brighter. The effects are illustrated in Figure 21.

Zoom In Zoom Out Reset image size

To summarize, five parameters are required to set up the binary model without X-ray irradiation: (1) the black hole mass M_•, (2, 3) the secondary mass (M_*), radius (R_*), and temperature (T_*) as determined by its initial mass and age, (4) the fraction f_R the secondary fills its Roche lobe, and (5) the inclination angle Θ. Five extra parameters are required to incorporate the X-ray irradiation: (6) the X-ray luminosity L_X, (7) the disk fattening factor X_h, (8) the disk size factor X_d, (9) the disk albedo β, and (10) the X-ray-to-thermal conversion factor γ. In total, 10 key parameters are required to fully describe an X-ray irradiated binary model for ULXs.