format: text

text:
# Mass-Loss Rate Determination for the Massive Binary V444 Cyg using 3-D Monte-Carlo Simulations of Line and Polarization Variability

## 1 Introduction

The short-period ($`P=4.212`$ days, Khaliullin et al. 1984) eclipsing, massive binary V444 Cyg (WN5+O6 III-V) has been the subject of extensive studies since its discovery (e.g., Wilson 1939; Münch 1950; Shore & Brown 1988; Marchenko et al. 1994, 1997). V444 Cyg exhibits variability in polarization, line strength and X-ray flux as a function of orbital phase. The variability arises from occultation of the photosphere, from perturbations induced in the extended atmosphere of the Wolf-Rayet star (W-R) by the O star and its wind, and from the wind-wind interaction region. Despite the complexities, many authors (Hamann & Schwarz 1992; St-Louis et al. 1993; Marchenko et al. 1994; Cherepashchuk et al. 1995; Moffat & Marchenko 1996; Marchenko et al. 1997; Stevens & Howarth 1999) have used this object to determine fundamental parameters of the W-R star.

One of the most important and uncertain parameters in the stellar evolution calculation of a massive star is the mass loss rate ($`\dot{M}`$) during the W-R star phase. Presently there are at least 6 different methods of determining the mass-loss rate of a W-R star: 1. dynamical method using the change in the orbital period of binary (e.g., Khaliullin et al. 1984), 2. polarization variation method (St-Louis et al. 1988), 3. X-ray spectra synthesis method using a hydrodynamical model (e.g., Stevens et al. 1996; Pittard & Corcoran 2002), 4. radio/IR continuum flux method (e.g. Wright & Barlow 1975), 5. radiative transfer method (e.g., Hillier 1989; Schmutz et al. 1989) and 6. photometric variability method (e.g., Lamontagne et al. 1996). The mass-loss rate estimated from methods 4 and 5 are, in general, about $`23`$ times higher than the values estimated from methods 1, 2 and 3.

In the case of V444 Cyg, the mass-loss rate measured from the orbital period change by Khaliullin et al. (1984), Underhill et al. (1990) and Antokhin et al. (1995) are $`\dot{M}_{\mathrm{WR}}=1.0\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$, $`\dot{M}_{\mathrm{WR}}=0.4\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ and $`\dot{M}_{\mathrm{WR}}=0.7\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ respectively. (See also Table 1.) Prinja et al. (1990) found $`\dot{M}_{\mathrm{WR}}=2.4\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ from the free-free radio emission flux. By modeling the infrared lines, Howarth & Schmutz (1992) obtained $`\dot{M}_{\mathrm{WR}}=\left(2.55\right)\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$. More recently, Nugis et al. (1998) published the “clumping-corrected” mass-loss rates for 37 W-R stars from the radio emission power and the spectral index $`\left(\alpha =d\mathrm{ln}f_\nu /d\mathrm{ln}\nu \right)`$. They obtained $`\dot{M}_{\mathrm{WR}}=0.92\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ for V444 Cyg.

Hamann & Schwarz (1992), by using a non-LTE radiative transfer model with a spherically expanding atmosphere, simultaneously fitted a set of helium emission lines and the light curve of V444 Cyg. The mass-loss rate estimated from their analysis is $`\dot{M}_{\mathrm{WR}}=1.26\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$. Similar to Methods 4 and 5 mentioned above, the model assumes no clumps in the stellar wind of the W-R star; hence, the mass-loss rate value is most likely overestimated. An important conclusion of their work is that the light curve is very sensitive to the inclination angle which must be treated as a fitting parameter.

St-Louis et al. (1988, hereafter STL1) derived an analytical expression for the mass-loss rate of a W-R star in a binary system, based on the model of Brown et al. (1978) which predicts the continuum polarization from a binary system as a function of orbital phase. The original model of Brown et al. (1978) assumes point sources and a single scattering atmosphere (optically thin). In the same paper, STL1 estimated $`\dot{M}_{\mathrm{WR}}`$ of V444 Cyg to be $`0.9\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ assuming $`V_{\mathrm{}}=2500\mathrm{km}\mathrm{s}^1`$. Later St-Louis et al. (1993, hereafter STL2) corrected this value to $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ using the new estimate of $`V_{\mathrm{}}=1785\mathrm{km}\mathrm{s}^1`$ (Prinja et al. 1990). STL2 also derived analytical expressions of Q & U Stokes parameters as a function of orbital phase *especially near the secondary eclipse* where the W-R is eclipsed by the O star. The observed polarization curves were fitted with their model, and the mass-loss rate of the W-R star in V444 Cyg was estimated ($`0.75\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$). Although these models are very simple and potentially powerful, the validity is still in question since the models assume that the W-R star has a smooth wind and its atmosphere is optically thin. In addition, the derived expression of $`\dot{M}_{\mathrm{WR}}`$ in STL1 and $`Q`$ & $`U`$ in STL2 involves an integral which diverges. To avoid the infinity, they fixed the lower radial integration limit to $`2R_{}`$.

Another important, but still uncertain parameter for V444 Cyg is the luminosity ratio, $`q=L_{\mathrm{WR}}/L_\mathrm{O}`$. Beals (1944) spectroscopically estimated $`q`$ for the wavelength range $`40005000`$Å. He obtained $`q=0.21`$ using the emission lines, and $`q=0.19`$ using the absorption lines. For the same system, Cherepashchuk et al. (1995) estimated the value of $`q`$ in a similar manner as Beals (1944), but with some modification to the calculation of the equivalent width of the emission/absorption lines. They found $`q=0.60\pm 0.06`$ for $`\lambda \lambda `$4000–6000; however, the values from different lines show a large scatter, ranging from 0.36 to 0.80. From the two principal emission lines used in these studies, they found $`q=0.51`$ using He ii $`5412`$ Å and $`q=0.59`$ using He i $`5876`$ Å (see their Table 1). Lastly, Hamann & Schwarz (1992) estimated $`q=3.1`$ according to their models of the light curve and helium spectrum, but the value is most likely incorrect since the effect of clumping was not included in their model.

The time dependent spectrum of the helium lines obtained by Shore & Brown (1988), Underhill et al. (1988), Marchenko et al. (1994), Marchenko et al. (1997) and Stevens & Howarth (1999) shows very complex variability associated with the orbital motion of the stars and the bow shock. Although the terminal velocities of the O star ($`2540\mathrm{km}\mathrm{s}^1`$: Underhill et al. 1988) and the W-R star ($`1785\mathrm{km}\mathrm{s}^1`$: Prinja et al. 1990) are comparable, the W-R star has a much denser wind, c.f., $`\dot{M}_\mathrm{O}0.6\times 10^6M_{\mathrm{}}\mathrm{yr}^1`$ (Leitherer 1988; Marchenko et al. 1997) and $`\dot{M}_{\mathrm{WR}}\left(0.41.26\right)\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ (Table 1). The momentum carried by the W-R wind is several times larger than that of the O star wind, and the bow shock produced from the colliding stellar winds is folded toward the O star. Further, the orbital speed of the binary is a significant fraction of the terminal velocity of the W-R stellar wind ($`V_{\mathrm{orb}}/V_{\mathrm{}}1/4`$ ); hence, the shape of the bow shock is not expected to be cylindrically symmetric around the axis joining the center of the W-R star and the O star. In addition, the bow shock is not smooth, i.e., the shape of the shock front is distorted by instabilities (see Stevens & Pollock 1994; Gayley et al. 1997; Pittard 1998; Pittard & Stevens 1999) — the Kelvin-Helmholtz and the non-linear thin-shell instabilities discussed in Chandrasekhar (1961), Vishniac (1983), Vishniac (1994) and Blondin & Marks (1996).

The first X-ray evidence of the colliding winds (CWs) in V444 Cyg was found by Moffat et al. (1982) from the flux variability seen in Einstein Observatory data. In addition, Corcoran et al. (1996) and Maeda et al. (1999) found a similar X-ray variability using ROSAT and ASCA respectively. Theoretical aspects of CW related X-ray variability have been developed by, for example, Luo et al. (1990), Usov (1990), Usov (1992), Stevens et al. (1992), Myasnikov & Zhekov (1993) and Pittard & Stevens (1997).

In a previous paper (Kurosawa & Hillier 2001, hereafter, Paper I), we developed a 3-D Monte-Carlo radiative transfer model which predicts the level of continuum and line polarization produced by scattering of light in a stellar atmosphere with arbitrary geometry. The model can predict the variability features associated with orbital motion of a binary system, including the polarization level, the flux level and line profile shapes. The model can treat a finite sized stellar disk, multiple scattering, line absorption of the continuum photons and emission from multiple light sources. To achieve high precision with a minimum of data points (i.e. to save computer memory), an 8-way tree data structure created via a “cell-splitting” method (e.g, Wolf et al. 1999) is used in the model. This method is essential since the model assumes no symmetry in the atmosphere. A logarithmic grid, which is commonly used in spherical and axi-symmetric codes, is not readily implementable in a model of arbitrary geometry.

The goal of this paper is to apply the 3-D Monte-Carlo model developed in Paper I to V444 Cyg, and to interpret the observed variability features (seen in the continuum polarization and line profile shape) as a function of orbital phase. We estimate the mass-loss rate of the W-R component in the binary, the orbital inclination, and the luminosity ratio of the two stars, by fitting the observed He i $`5876`$ Å and He ii $`5412`$ Å line profiles, and the continuum light curves of three Stokes parameters $`(\text{I},Q,U)`$ in the $`V`$ band consistently. Unlike the radio/IR/spectral methods discussed earlier, the derived mass-loss rate determined by the polarization method is insensitive to the amount of clumps in the stellar wind. The continuum polarization level is proportional to the electron number density while the thermal emissivity (of radio/IR) is proportional to the square of the density. With a complicated density distribution, a Monte Carlo simulation provides the only method to realistically predict variability in polarization, continuum flux and line strength. An original aim was also to obtain improved estimates of the stellar parameters such as the radii of the W-R star and the O star, the monochromatic luminosity ratio of the two stars and the orbital inclination angle. Due to both model and observation uncertainties, this proved not to be feasible.

In § 2, we discuss our models, examine the variability of He i $`5876`$ Å and He ii $`5412`$ Å lines, and determine the mass-loss rate of the Wolf-Rayet star in V444 Cyg. The discussion on the results will be given in § 3, and the conclusion in § 4.

## 2 Analysis

The five most important parameters which we initially hoped to constrain in our analysis were: the mass-loss rate of the W-R star ($`\dot{M}_{\mathrm{WR}}`$), the luminosity of the W-R star ($`L_{\mathrm{WR}}`$), the radius of the W-R star ($`R_{\mathrm{WR}}`$), the volume filling factor ($`f`$) and the monochromatic luminosity ratio of the W-R and O stars, $`q\left(\lambda \right)=L_{\mathrm{WR}}\left(\lambda \right)/L_\mathrm{O}\left(\lambda \right)`$, at $`\lambda =5630`$Å. The volume filling factor is defined as the fractional volume which contains material (the higher density regions or clumps), and it controls the amount of clumps in the atmosphere (see e.g., Abbott et al. 1981). It can also be thought of as the fractional length along any random line of sight which contains clumps.

Since the procedure of the model fitting is rather complicated, we summarize the steps for determining the parameters below:

1. Understand the basic behavior of the variability seen in He i $`5876`$ Å and He ii $`5412`$ Å lines to determine at which binary phase the spectrum is least affected by the bow shock and by the presence of the O star atmospheric absorption. In other words, determine at which phase the spectrum is modeled best by a spherical radiative transfer model (CMFGEN).
2. At the phase chosen to be the best in Step 1, fit the observed helium spectrum for different W-R star radii, and find the corresponding value of $`q`$. In general, $`L_{\mathrm{WR}}=2\times 10^5L_{\mathrm{}}`$ is used111This luminosity value is chosen from the mass-luminosity relation of pure Helium stars given in Meynet et al. (1994), and $`M_{WR}10M_{\mathrm{}}`$ (Marchenko et al. 1994). .
3. Fit the continuum I, Q, U light curves around the secondary eclipse ($`\varphi 0.5`$) to constrain the radius of the W-R star. (*$`q`$ and $`R_{\mathrm{WR}}`$ are constrained*.)
4. With $`R_{\mathrm{WR}}`$ and $`q`$ values determined in Step 3, re-fit the helium spectrum with different values of $`\dot{M}`$.
5. Compute the continuum $`I`$, $`Q`$, $`U`$ light curves for the $`R_{\mathrm{WR}}`$ , $`q`$ and $`\dot{M}_{\mathrm{WR}}`$ values used in Step 4, and fit only the $`I`$ light curve by adjusting the orbital inclination ($`i`$) and $`R_\mathrm{O}`$.
6. Now there are several models which fit the observed $`I`$ light curve and the helium spectrum simultaneously, but only a certain value of $`\dot{M}_{\mathrm{WR}}`$ will fit the $`Q`$ and $`U`$ light curves; therefore, the value of $`\dot{M}_{\mathrm{WR}}`$ is constrained. (*$`i`$, $`R_\mathrm{O}`$ and $`\dot{M}_{\mathrm{WR}}`$ are constrained*.)

### 2.1 Line Variability and Geometric Configuration of the Binary

Since the 3-D Monte Carlo polarization model introduced in Paper I is not a self-consistent model (i.e., the radiation does not modify the opacity of the atmosphere), the opacity and the emissivity must be supplied as the inputs of a model. There are three possible sources of the input field: 1. parametrized functions from a semi-analytical model, 2. output from a spherical non-LTE radiative transfer model with some modifications, and 3. output from a hydrodynamical calculation in conjunction with method 2. Method 1 is useful for testing models, but not suited for modeling a real object. In our analysis, method 2 is mainly used.

The first step, in method 2, is to compute the object’s spectrum using a non-LTE radiative transfer code (CMFGEN) which includes thousands of lines to treat the line-blanketing effect (Hillier & Miller 1998). For example, see Hillier & Miller (1999) and Herald et al. (2001) for a detailed discussion of W-R star spectrum models across a wide wavelength range. CMFGEN provides the opacity and the emissivity of the W-R star and the O star separately, and is designed for a single star with a “spherically” expanding atmosphere; hence, the output fields (opacities etc.) are functions of radius only. The second step is to compute the variability of the Stokes parameters ($`I`$, $`Q`$, $`U`$) and the line profiles by using the 3-D Monte Carlo model as a function of the binary phase.

The velocity field in a CMFGEN model is assumed to have the following form (Hillier & Miller 1999):

$$v\left(r\right)=\frac{V_0+V_1+V_2}{1+\left(V_0/V_{\mathrm{core}}1\right)\mathrm{exp}\left[\left(rR_{}\right)/h_{\mathrm{eff}}\right]}$$
(1)
where $`V_1=\left(V_{\mathrm{}}V_{\mathrm{ext}}V_0\right)\left(1R_{}/r\right)^{\beta _1}`$, $`V_2=V_{\mathrm{ext}}\left(1R_{}/r\right)^{\beta _2}`$, $`h_{\mathrm{eff}}`$ is an isothermal effective scale height in the inner atmosphere, $`V_{\mathrm{}}`$ is the terminal velocity, $`V_{\mathrm{core}}`$ is the expansion velocity of the core, and $`V_{\mathrm{ext}}`$ & $`\beta _2`$ are the parameters for the outer parts of the wind. $`\beta _1`$ is equivalent to the $`\beta `$ in the classic beta-velocity law, $`V\left(r\right)=V_{\mathrm{}}\left(1R_{}/r\right)^\beta `$, for the “inner part” of the stellar wind. See Najarro et al. (1997) and Hillier & Miller (1999) for more explanations on this velocity law. Although the line-driving of O-star winds is well understood, the multiple-scattering of W-R winds complicates theoretical models, and the standard CAK theory (Castor, Abbott, & Klein 1975) does not apply for W-R stars.

Fig. 1 shows the basic configuration of the model geometry. The W-R star with its extended atmosphere is located at the center. Surrounding the O star, the bow shock due to the CWs has a paraboloidal shape, and points away from the W-R star since its wind is much stronger than that of the O star. The head of the shock front is tilted from the purely radial direction because of the rapid orbital motion. Using the ratio of the terminal velocity of the W-R component ($`V_{\mathrm{}}=1785\mathrm{km}\mathrm{s}^1`$, Prinja et al. 1990) and the orbital speed ($`V_{\mathrm{orb}}460\mathrm{km}\mathrm{s}^1`$)222We simply estimated the value by assuming a circular orbit. If the radial velocity measurements of Marchenko et al. (1994) and the orbital inclination $`i=78.8^{}`$ (Robert et al. 1990) were used, we would obtain $`V_{\mathrm{orb}}440\mathrm{km}\mathrm{s}^1`$. Consequently, the tilt angle would be $`\delta 14^{}`$. , the tilt angle of the bow shock region can be estimated as $`\delta \mathrm{tan}^1(V_{\mathrm{orb}}/V_{\mathrm{}})15^{}`$. Extensive discussions on the morphology of the bow shock can be found in Marchenko et al. (1997), including the issue of whether the W-R wind is impacting onto the O star surface or not. (See also the hydrodynamical calculation of Pittard & Stevens 1999).

At $`\varphi `$ (orbital phase) $`=0`$, the W-R star is in front of the O star, and an observer is looking at the system from the top of Fig. 1 since the inclination angle $`i`$ is about $`80^{}`$ (almost edge-on view). At $`\varphi =0.5`$, the O star is in front of the W-R star. The directions of an observer at $`\varphi =0,\mathrm{\hspace{0.17em}0.25},\mathrm{\hspace{0.17em}0.5},\mathrm{\hspace{0.17em}0.75}`$ are indicated in the same figure, and they infer the sense of the orbital motion (counter-clockwise in the figure). The strong stellar wind from the W-R star is interrupted by the shock front, and the density behind the shock is assumed to be negligibly small for simplicity. In other words, we do not include the O star wind explicitly, and we place the bow shock zone rather artificially.

Next, we propose a simple model of the gas flow in the CW zone and the geometry of the CW zone. We assume the shape of the CW zone to be paraboloid333See Huang & Weigert (1982), Girard & Willson (1987) and Canto et al. (1996) for the predicted geometry of the thin-shell shocks created by the interaction of the two stellar winds in a highly radiative binary.  with thickness $`d`$. The center of the paraboloid is located along the y-axis but displaced by $`y_o`$ from the origin as shown in Fig. 2. The surface of the paraboloid can be written as:

$$y=a_yl^2+y_o$$
(2)
where $`l^2=x^2+z^2`$. We define the bow shock region to be the volume between this surface and the same surface displaced by $`d`$ (thickness) in $`+y`$ direction.

The flow of the gas in the CW zone is assumed to be tangential to the surface of the paraboloid; hence, there is no azimuthal component. Then, the velocity ($`𝐯_{\mathrm{cw}}`$) of the gas flow at point $`P`$ in the CW zone can be written as:

$$𝐯_{\mathrm{cw}}=v_y\widehat{𝐲}+v_l\widehat{𝐥}$$
(3)
where $`\widehat{𝐥}=\mathrm{cos}\varphi _c\widehat{𝐳}+\mathrm{sin}\varphi _c\widehat{𝐱}`$. $`v_y`$ and $`v_l`$ are the velocity components which are parallel and perpendicular to the y axis, respectively. $`a_y`$ is a parameter to control the asymptotic open angle of the bow shock. Further, by taking the derivative of Eq. 2 and using the absence of an azimuthal component of $`𝐯_{\mathrm{cw}}`$, we obtain:

$$\frac{dy}{dl}=2a_yl=\frac{v_y}{v_l}.$$
(4)
If the magnitude of the velocity ($`𝐯_{\mathrm{cw}}`$) at point $`P`$ is similar to that of the spherical flow ($`v_r`$) of the pre-shock W-R wind at $`P`$ , i.e., $`v_{\mathrm{cw}}v_r`$, then

$$v_r^2v_y^2+v_l^2.$$
(5)
Since the flow of gas in the bow shock region is usually slower than the surrounding flow according to the hydrodynamics model of Pittard & Stevens (1999), we introduce a free parameter, and rewrite the equation above as:

$$v_{\mathrm{cw}}=sv_r=\left(v_y^2+v_l^2\right)^{1/2}.$$
(6)
In general $`s<1`$, and is a function of the position (e.g., $`s=s\left(y\right)`$). However, in most of the analysis, $`s=1/2`$ is assumed for simplicity.

Using Eqs. 4, 5 and 6, the components of $`𝐯_{\mathrm{cw}}`$ can be expressed, in terms of the position $`(y,l)`$ on the paraboloid, as

$$v_y=\frac{2a_ylsv_r}{\left[1+\left(2a_yl\right)^2\right]^{1/2}},$$
(7)
and

$$v_l=\frac{sv_r}{\left[1+\left(2a_yl\right)^2\right]^{1/2}}.$$
(8)
These expressions are then corrected for the tilt of the bow shock due to the fast orbital motion of the system, as mentioned earlier. A more accurate velocity field in the bow shock region is needed to properly model the line profile and the line polarization, but it is not important for the continuum calculations.

Considering the radial dependence of the density, we define:

$$\rho _{\mathrm{shock}}=\rho _\mathrm{o}\left(\frac{y_\mathrm{o}}{r}\right)^n$$
(9)
where $`r`$ is the distance between the origin ($`0`$) and a point $`p`$ on the surface of the parabola, and $`y_\mathrm{o}`$ is the distance between the origin and the head of the paraboloid as shown in Fig. 2. The index $`n`$ is a scaling parameter, and $`\rho _\mathrm{o}`$ is the density at the head of the parabola (at $`y=y_\mathrm{o}`$). This formulation is similar to one in Stevens & Howarth (1999). $`\rho _\mathrm{o}`$ is estimated from the amount of the excess emission seen in the He i $`5876`$ Å line, and is discussed in the later section (§ 2.2.3). In our analysis, we use the bow shock density expressed in Eq. 9 with $`n=2`$, and the thickness ($`d`$) of the bow shock $`0.1R_\mathrm{O}`$, estimated from Pittard & Stevens (1999). In predicting the continuum polarization level, the global distribution of the electron density is more important than small structure. The column density or the electron scattering opacity through the bow shock is related to the amount of gas in the cavity which will be absent if there is no bow shock.

Some sensitivity of the line strengths and polarization may arise from assuming whether the helium is singly or doubly ionized. He i $`5876`$ Å and He ii $`5412`$ Å lines are often referred to as “diagnostic lines” for modeling a WN star atmosphere (e.g., Hamann et al. 1995). The latter is an optically thick line; hence, the absorption and scattering terms in radiative transfer equations must be treated properly. The gas in the bow shock is rapidly cooled by radiative processes, and it could potentially cause a significant amount of absorption in the He i $`5876`$ Å line. Figure. 3 shows the approximate location where He i $`5876`$ Å and He ii $`5412`$ Å lines originate according to a CMFGEN model of a typical WN5 star with $`R_{\mathrm{WR}}=5R_{\mathrm{}}`$. The figure shows that the He ii emission peaks at around $`3R_{\mathrm{WR}}`$ from the center of the W-R star, and the He i emission peaks around the location where the O star resides. This explains why the emission from the CW zone is more clearly seen in the He i line. In the following, the variability of the He i line due to the presence of the bow shock will be examined.

#### 2.1.1 Effect of the Colliding Winds

As explained earlier, the emission from the bow shock region is prominent for some lines (see Marchenko et al. 1994). In this section, the expected effect of the bow shock emission on the He i $`5876`$ Å line will be demonstrated using the 3-D Monte Carlo model. He i $`5876`$ Å is primarily a recombination line; hence, to first order we assume the line emissivity in the bow shock region is simply proportional to the square of the density. The tilt angle $`\delta =15^{}`$ and the $`\mathrm{open}\mathrm{angle}=90^{}`$ (Marchenko et al. 1994) are used for the geometry of the bow shock. The left graph in Fig. 4 shows a sequence of He i $`5876`$ Å emission lines as a function of phase ($`01`$ with steps of $`0.1`$ from the bottom) arising from a V444 Cyg model which includes an *unrealistically strong* bow shock emission. In fact, the line emission is dominated by the bow shock emission, and the underlying line emission from the W-R atmosphere is barely visible in this figure. This plot demonstrates the basic behavior of the emission from the bow shock region. (See also Lührs 1997 and Bartzakos et al. 2001 for similar models and discussion.) At $`\varphi =0.81.0`$, the bow shock is pointing away from an observer and tilted to the light of sight (Fig. 1); therefore, the emission peak appears on the red side of the flat-top He i emission line. The shape of the CW emission changes very rapidly from $`\varphi =0`$ to $`0.1`$ because a rapid change in the velocity component of the bow shock arms along the line of sight occurs near $`\varphi =0.1`$. The emission from the bow shock flattens out during $`\varphi =0.10.4`$ and $`\varphi =0.70.8`$ since the bow shock arms are almost perpendicular to the direction of the observer. At $`\varphi =0.40.5`$, the peak appears on the blue side since the bow shock arms are pointing towards an observer. A sequence of He i $`5876`$ Å line profiles with a more realistic level of bow shock emission is shown in the right panel of Fig. 4. The figure shows that the line strength increases around $`\varphi =0`$ and $`0.5`$ because the continuum level drops during the primary and secondary eclipses.

#### 2.1.2 O star Absorption Effect

The observed spectrum of V444 Cyg shows absorption lines due to the O star atmosphere on top of the broad emission lines from the W-R atmosphere. The position of absorption shifts as the relative speed of the stars changes. The distortion caused by this effect is prominent for some lines at certain binary phases. Since the model does not take into account the photospheric spectrum of the O star, the O star absorption lines have been added to the CMFGEN model in an approximate manner. In principal they could easily be incorporated into the Monte Carlo method by changing the appropriate boundary conditions.

The O star component of V444 Cyg has spectral type O5-O6.5 III-V (Marchenko et al. 1994). The O star (O6) spectrum is modeled separately by CMFGEN, and the effect of the rotational broadening of lines is added to the spectrum. In this process, the rotational speed $`V_{\mathrm{rot}}\mathrm{sin}i=215(\pm 13)\mathrm{km}\mathrm{s}^1`$ (Marchenko et al. 1994) is used.

Fig. 5 shows the combined spectrum of the W-R and the O star computed by CMFGEN as a function of the orbital phase. The sinusoidal motion of the O star absorption component is clearly seen in both lines. To examine only the effect of the O star, *the emission from the CW zone is not included in the model spectra shown in this figure*. A specific luminosity ratio ($`q=L_{\mathrm{WR}}/L_\mathrm{O}=0.45`$ at $`\lambda =5630`$ Å) was used when the two spectra were combined.

#### 2.1.3 Combined Effects Compared with Observation

In order to use CMFGEN for the model fitting, we must choose a binary phase when the lines are least affected by the bow shock emission and the O star absorption. Fig. 5 shows that the absorption by the O star does not affect the overall shape of He i $`5876`$ Å, but it does affect the shape of He ii $`5412`$ Å. At $`\varphi 0.2`$ and $`\varphi 0.7`$, the He i line shape is least affected by the O star, but at these phases the emission from the bow shock region is spread over the He i line as seen in Fig. 4. As a result, the line flux level is likely slightly higher than that of the spherical model. On the other hand, at $`\varphi 0`$ and $`\varphi 0.5`$, the emission from the bow shock region is well concentrated near the blue or the red edge of the He i line, but at these phases, the stars are eclipsed; hence, the line strength will be underestimated in a simple spherical model (CMFGEN). Meanwhile at $`\varphi 0.3`$ and $`\varphi 0.8`$, the bow shock emission is concentrated only on either the blue or the red side of the He i line, and the O star absorption affects the same half side of the line profile. Unfortunately, at $`\varphi 0.3`$, the O star absorption is on the red side of the helium line, and the bow shock emission is on the blue side of the line. So the line is contaminated on both sides. On the other hand at $`\varphi 0.8`$, both the O star absorption and the emission from CW are on the red side of a line; hence we can reliably use the blue part of He i $`5876`$ Å and He ii $`5412`$ Å for the model fitting. Figure 6 shows the sample CMFGEN model spectrum compared with the observed spectrum of Marchenko at six different phases. The agreement between the model and the observation at $`\varphi 0.8`$ is excellent. Near the eclipsing phases ($`\varphi 0`$ and $`0.5`$), the spherical model CMFGEN does not fit the observation well since the contribution from the bow shock is not included in the model. The two narrow Na i interstellar absorption lines are present in the red side of observed He i $`5876`$ Å lines.

### 2.2 Mass-Loss Rate Estimation of the W-R Component

The idea behind the mass-loss estimation performed here is very similar to the one in STL1 and STL2. We assume that there is a correspondence between the amount of continuum polarization and the mass-loss rate of the W-R component. This is supported by the analytic expression of $`\dot{M}`$ in terms of $`A_p`$ (a semi-major axis of polarization on Q-U plane) derived by STL1, using the basic results from the classic paper on binary polarization by Brown et al. (1978).

#### 2.2.1 Investigation of the W-R Radius

For an assumed mass-loss rate, $`\dot{M}`$, of the W-R star, there are several CMFGEN models which can fit the observational spectrum of He ii $`5412`$ Å and He i $`5876`$ Å simultaneously. As discussed in the previous section, the spectrum at $`\varphi 0.8`$ is used in the model fitting since the contamination from the O star atmospheric absorption and the CW emission is expected to be small at this orbital phase. Fig. 7 shows the fits of the observed spectrum at $`\varphi =0.826`$ for different values of $`R_{\mathrm{WR}}`$ (=$`2.5`$, $`4.0`$, $`5.0R_{\mathrm{}}`$)444The different solutions correspond to different assumed distances for V444 Cyg. The W-R radius ($`R_{\mathrm{WR}}`$) used in this paper corresponds to the inner boundary of the model atmosphere where the Rosseland optical depth $`\tau _\mathrm{R}30`$.  with the mass-loss rate of the W-R star fixed ($`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$). The corresponding volume filling factors of the models are $`0.08`$, $`0.04`$, $`0.05`$. All three models fit the observation very well; therefore, the spectral solutions are degenerate. The He ii profiles with different values of $`R_{WR}`$ are not significantly different from each other when the O star spectrum is added to the W-R spectrum with an appropriate value of $`q`$ ($`=0.30,\mathrm{\hspace{0.17em}0.35},\mathrm{\hspace{0.17em}0.45}`$ for $`R_{\mathrm{WR}}=2.5,\mathrm{\hspace{0.17em}4.0},\mathrm{\hspace{0.17em}5.0}R_{\mathrm{}}`$ models respectively). In the case of a single W-R star spectrum, a small but noticeable difference in the profile shape of the He ii line, usually becoming more triangular for a larger $`R_{\mathrm{WR}}`$, should be present. Unfortunately, it is hard to see this effect in the V444 Cyg model because the continuum flux of the O star component is stronger ($`q0.5`$); hence, the line strength is weakened.

In an attempt to remove this degeneracy, the light curves ($`I`$, $`Q`$, $`U`$) are computed by the 3-D Monte-Carlo model (Kurosawa & Hillier 2001) for the spectral models shown in Fig. 7 for fixed $`R_\mathrm{O}`$ and $`\dot{M}_{\mathrm{WR}}`$ values. First, models without the bow shock (Figs. 812) are considered for clarity. The effect of the bow shock on the light curve will be examined later. A small difference is expected to be seen (c.f., STL2) in the shape of the polarization curve near the secondary eclipse ($`\varphi 0.5`$, i.e., when the W-R star is behind the O star).

The results, using $`R_\mathrm{O}=7.2R_{\mathrm{}}`$ and $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$, are shown in Fig. 8. The orbital inclination angles used for the fits are $`i=78.8^{}`$, $`78.5^{}`$ and $`78.0^{}`$ for the models with $`R_{\mathrm{WR}}=2.5R_{\mathrm{}}`$, $`4.0R_{\mathrm{}}`$ and $`5.5R_{\mathrm{}}`$ respectively. The inclination angles are consistent with $`i=78^{}\pm 1^{}`$ (Cherepashchuk 1975) obtained from the light curve solution. STL2 estimated the inclination angle by fitting the double sine curves of $`Q`$ and $`U`$ polarization components with the formula given by Brown et al. (1978). Their result is $`i=80.8^{}\pm 1.6^{}`$. In Fig. 8, the three models fit the observed $`I`$ light curve moderately well. Although not shown here, all three models also fit well the observed $`I`$ light curve outside of the phase range shown in Fig. 8. The deviations from perfect anti-symmetry around $`\varphi =0.5`$ seen in the observational $`I`$ light curve and the anomaly seen around $`\varphi 0.530.54`$ in the U observational light curve are not well fitted, and may be related to the presence of the bow shock in the system.

The width of the $`I`$ light curve around $`\varphi =0.5`$ is slightly wider for the model with the larger W-R star radius. The $`R_{\mathrm{WR}}=5.0R_{\mathrm{}}`$ model fits the overall shape of the light curve better than the two other models with smaller radii. However, the width of the secondary eclipse depends also on the mass-loss rate of the W-R star, and we discuss this dependency later. Unfortunately, no significant difference in the Q and U light curves are seen for different values of W-R star radii, unlike that found by STL2. In fact, the difference in Q and U light curves seen by STL2 for different $`R_{\mathrm{WR}}`$ (see their Fig. 8) is not purely due to changing $`R_{\mathrm{WR}}`$, being also coupled with $`\dot{M}_{\mathrm{WR}}`$ and the adopted inclination.

In this analysis, we were not able to constrain the value of $`R_{\mathrm{WR}}`$. Models with a wide range of $`R_{\mathrm{WR}}=2.55.0R_{\mathrm{}}`$ can fit the observed $`I`$, $`Q`$ and $`U`$ light curves reasonably well, and there is no significant difference among them as long as appropriate values of the orbital inclination angle and the O star radius are chosen for each model. Despite this uncertainty in the W-R star radius, we continue to derive the mass-loss rate of the W-R star.

#### 2.2.2 Determination of the W-R Mass-Loss Rate

Next, the observed helium spectrum is modeled for different values of $`\dot{M}_{\mathrm{WR}}`$ by separately adjusting the values of $`f`$ and $`q`$ for $`R_{\mathrm{WR}}=5.0`$ and $`2.5R_{\mathrm{}}`$. Fig. 9 shows the results for $`\dot{M}_{\mathrm{WR}}=0.3`$, $`0.6`$ and $`1.2\times 10^5`$$`M_{\mathrm{}}\mathrm{yr}^1`$ with the fixed radius of the W-R star: $`R_{\mathrm{WR}}=5.0R_{\mathrm{}}`$. No significant difference in the profile shapes is seen in these models. To remove this degeneracy, once again, the polarization curves were computed; however, this time the amplitude of the “double-sine” shaped polarization curves were fitted. The polarization curves for each value of $`\dot{M}_{\mathrm{WR}}`$ corresponding to the spectral models in Fig. 9 are shown in Fig. 10 along with the light curves at $`\lambda =5630`$ Å. The model parameters used for the fits in Fig. 10 are summarized in Table 2. The figure clearly shows that the model with $`\dot{M}_{\mathrm{WR}}=1.2\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ overestimates and the one with $`\dot{M}_{\mathrm{WR}}=0.3\times 10^5\dot{M}_{\mathrm{}}\mathrm{yr}^1`$ underestimates the amplitude of the observed $`Q`$ and $`U`$ curves, while $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5\dot{M}_{\mathrm{}}\mathrm{yr}^1`$ model is consistent with the observations. The width of the $`I`$ light curve for $`\lambda =5630`$ Å around the primary eclipse ($`\varphi =0.0`$) with $`\dot{M}_{\mathrm{WR}}=0.3\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ is slightly narrower than the observation, and that with $`\dot{M}_{\mathrm{WR}}=1.2\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ is slightly wider than the observation. The observed $`I`$ light curve at $`\lambda =5630`$ Å is fitted with the $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ model very well.

Since there is a large uncertainty in the $`R_{\mathrm{WR}}`$ value, the results of similar calculations with $`R_{\mathrm{WR}}=2.5R_{\mathrm{}}`$ are shown in Figs. 11 and 12. They are very similar to the results of the $`R_{\mathrm{WR}}=5.0R_{\mathrm{}}`$ models, and again the model with $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5\dot{M}_{\mathrm{}}\mathrm{yr}^1`$ fits the observed $`I`$, $`Q`$, and $`U`$ light curves at the wavelength around $`5630`$ Å consistently. A summary of the model parameters is given in Table 2.

#### 2.2.3 Effects of the Bow Shock on Light Curves

We now consider the possible effects of the presence of the bow shock region on our computed light curves (Figs. 10 and 12). To do so, the observed He i $`5876`$ Å line at the phase where the effect of the CW emission is most prominent is modeled first. The spectrum taken at $`\varphi =0.419`$ (see Fig. 6) is chosen for this purpose. The profile at this phase is fitted by adjusting $`\rho _o`$ in Eq. 9 to produce a similar amount of the excess emission as seen on the top of the He i $`5876`$ Å line. Since the true line opacity and emissivity in the bow shock region are not known, this was done just to produce the bow shock line emission seen in He i $`5876`$ Å. We do not intend to perform a very rigorous fit of this line although that would be a good topic for a future investigation. The result from the 3-D Monte Carlo calculation for Model B is shown in Fig. 13. A tilt of the bow shock, $`\delta =15^o`$, and open angle $`=90^{}`$ (Marchenko et al. 1994) were used in this model. The fit shown in Fig. 13 is reasonably good except for the flux level around $`V=200500\mathrm{km}\mathrm{s}^1`$ where the observed spectrum is affected by the O star absorption. For simplicity, O star atmospheric absorption effects are not included in this model; however, an appropriate $`q`$ was used to produce the correct O star continuum flux. The ratio the bow shock density to the background density ($`\rho _{\mathrm{cw}}/\rho _{\mathrm{bg}}`$), due to the W-R stellar wind, is found to be about 40 according to this model. This ratio is similar to that of Pittard & Stevens (1999).

With this bow shock model, the light curves ($`I`$, $`Q`$, and $`U`$) at $`\lambda =5630`$ Å of Model B ($`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ and $`R_{\mathrm{WR}}=5.0R_{\mathrm{}}`$) in Fig. 10 are recomputed. Fig. 14 shows the result. There is no significant effect of the bow shock region on the continuum $`I`$, $`Q`$ and $`U`$ light curves. The two models are identical to each other within the range of the error. The same conclusion was obtained for the other models shown in Figs. 10 and 12.

#### 2.2.4 Effect of the O-Star Wind using a Hydrodynamical Model

To check the validity of ignoring the O star wind and the geometrical bow shock model, the polarization light curves were computed using the hydrodynamical model of Pittard & Stevens (1999), which includes the orbital motion of the stars and a parameterized form of radiative driving. First, the opacity of the W-R star atmosphere was computed with CMFGEN using the parameters adapted by Pittard & Stevens (1999): $`R_{\mathrm{WR}}=4.0R_{\mathrm{}}`$, $`R_\mathrm{O}=10.0R_{\mathrm{}}`$, $`L_{\mathrm{WR}}=1.40\times 10^5L_{\mathrm{}}`$, $`L_\mathrm{O}=2.4\times 10^5L_{\mathrm{}}`$, $`\dot{M}_{\mathrm{WR}}=9.72\times 10^6M_{\mathrm{}}\mathrm{yr}^1`$, $`\dot{M}_\mathrm{O}=5.78\times 10^7M_{\mathrm{}}\mathrm{yr}^1`$, $`V_{\mathrm{}}^{\mathrm{WR}}=1800\mathrm{km}\mathrm{s}^1`$ and $`V_{\mathrm{}}^\mathrm{O}=2900\mathrm{km}\mathrm{s}^1`$. The CMFGEN model with these parameters does not fit the observed helium spectrum, but the opacity from this model will be used in the polarization calculation to maintain the consistency with the hydrodynamical model.

Using this W-R star atmosphere model, the following two models are computed: 1. a model without the bow shock and the O star wind, 2. a model with the O star wind and the bow shock from the hydrodynamical calculation of Pittard & Stevens (1999). The density distribution for the second case is shown in Fig. 15. To assign the opacity and emissivity in the O star wind and the bow shock regions, we assumed that the thermal opacity and the emissivity are proportional to the square of the density, and the electron scattering opacity is proportional to the density. Then the opacity and the emissivity are scaled with those from the W-R star atmosphere model. The light curves computed for these models are shown in Fig. 16. No significant differences between the two models are seen in this figure. This validates one of our main model assumptions, i.e., the O star wind does not significantly affect the continuum polarization curves. It also assures that the details of the bow shock structure are not important for predicting the *continuum* polarization curves though they are very important for a line calculation.

The model light curves in Fig. 16 do not fit the observed light curves. By increasing the inclination angle, the depth of the eclipses (at $`\varphi =0`$ and $`0.5`$) will increase; however, without adjusting other parameters, it is impossible to decrease the depth of the secondary eclipse. It is clearly seen from the location of minima in the eclipsing $`Q`$ and $`U`$ light curves ($`\varphi 0.46`$) that the O star radius is too large. Surprisingly, the amplitudes of the $`Q`$ and $`U`$ double sine curves fit the observation well with $`\dot{M}_{\mathrm{WR}}=9.72\times 10^6M_{\mathrm{}}\mathrm{yr}^1`$. This seems contradictory to our earlier analysis, but is due to the inconsistent W-R model, which does not fit the observed helium spectrum, being used in the polarization calculations. Since the parameters from the model of Pittard & Stevens (1999) are fixed, $`q=L_{\mathrm{WR}}/L_\mathrm{O}=0.583`$ is also fixed in these models. If $`q=0.38`$ is used by increasing $`L_\mathrm{O}`$ instead, the spectral model will fit the observed helium spectrum, but the higher $`L_\mathrm{O}`$ results in larger amplitudes of the double-sine polarization curves which then no longer fit the observation.

In summary, from this analysis, we conclude that the mass-loss rate of the W-R star in V444 Cyg is $`0.6\left(\pm 0.2\right)\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$. The size of the error in $`\dot{M}`$ is estimated from the light curves in Figs. 10 and 12. The mass-loss rate determined here is consistent with the values of STL2, and it is slightly higher than that of Underhill et al. (1990). (See Table 1 for other published values.) The range of the volume filling factor for the W-R star atmosphere is estimated to be $`0.0500.075`$ with the corresponding range of the W-R star radius, $`5.02.5R_{\mathrm{}}`$. We also found that the presence of the bow shock, the exact details of the bow shock, and the presence of the O star wind are unimportant for the continuum fits, and they do not affect the mass-loss rate estimate given here.

## 3 Discussion

### 3.1 Infrared Light Curves

In order to check the consistency of our best fit models (Models B and E with $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$) in § 2.2 with some other observational aspects, the infrared light curves (at $`\lambda =2.2\mu \mathrm{m}`$) were also computed. In order to calculate the IR light curves, the monochromatic luminosity ratio $`q=L_{\mathrm{WR}}/L_\mathrm{O}`$ at $`\lambda =2.2\mu \mathrm{m}`$ must first be estimated. In Model B, $`q`$ was determined to be $`0.45`$ at $`\lambda =5630`$ Å in the previous section. By normalizing the O star flux with this value of $`q`$ at $`\lambda =5630`$ Å, $`q`$ for the wavelengths between $`3000`$ Å and $`33,000`$ Å is plotted in Fig. 17 for three different O star effective temperatures. The figure shows that the O star is brighter than the W-R star at $`\lambda =5630`$ Å, but the W-R star is brighter than the O star after $`\lambda 20,000`$ Å. In model B, an O star temperature of $`T_\mathrm{O}=35,000\mathrm{K}`$ was used; hence, the corresponding value of $`q`$ at $`\lambda =2.2\mu m`$ is 1.10 according to the graph. For models with higher O star temperatures, the values of $`q`$ are slightly higher.

The resulting model IR light curves are shown in the right hand side of Figs. 10 and 12 along with Q and U polarization light curves for $`\dot{M}_{\mathrm{WR}}=0.3\times 10^5M_{\mathrm{}}`$, $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$ and $`\dot{M}_{\mathrm{WR}}=1.2\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$. The figure also includes the observed IR light curve from the broadband ($`K`$) photometry obtained by Hartmann (1978). The depth of the primary eclipse (at $`\varphi =0`$ or $`1`$ when the W-R star is in front of the O star) is fitted well for both $`R_{\mathrm{WR}}=2.5R_{\mathrm{}}`$ and $`5.0R_{\mathrm{}}`$ models, but the secondary eclipse modeled (at $`\varphi =0.5`$ when the O star is in front of the W-R star) is shallower than the observed IR light curve. We also notice that the maximum $`I`$ flux level ($`I_{\mathrm{max}}`$) in the observed light curve is not well defined. Cherepashchuk et al. (1984) and Hamann & Schwarz (1992) had a similar problem – the primary and secondary eclipses are too shallow in their light curve model. While the light curve model of Hamann & Schwarz (1992) did not fit the depth of both eclipses, our model is consistent with at least the primary eclipse.

Cherepashchuk et al. (1984) discussed two possible causes for this problem: a). the existence of “third light” and b). the result of ignoring the reflection effect. They noted that one possible source of third light is the contribution of emission lines in the broad-band photometry at $`2.2\mu \mathrm{m}`$ as originally mentioned by Hartmann (1978), but we now expect this effect to be generally small (unless the He i $`2.06\mu \mathrm{m}`$ emission line is abnormally large). To examine this possibility, we have run a model which includes the reflection by the disk of the O star. However we find that it does not affect the infrared and optical light curves. Cherepashchuk et al. (1984) also conjectured that the increase in the depth and the width of the secondary eclipse (from optical to IR) is a result of clumps in the WR stellar wind.

If the wrong flux ratio at $`\lambda =2.2\mu \mathrm{m}`$ is used in our model, the secondary eclipse could become shallower at this wavelength. The differences in the depth of the secondary eclipse in Figs. 10 and 12 are approximately $`8`$% and $`5`$%. This could be explained by the wrong O star temperature used in the model. According to Fig. 17, $`q`$ at $`2.2\mu \mathrm{m}`$ also increases about $`8`$% if we use $`T_\mathrm{O}=45,000\mathrm{K}`$ although this would not fit the absorption components of the observed spectrum.

### 3.2 Polarization Curves near the Secondary Eclipse

In both Model B (Fig. 10) and Model E (Fig. 12), the amplitude of the $`U`$ polarization variation near the secondary eclipse ($`\varphi =0.5`$) underestimates the observed level. A larger polarization amplitude can be obtained if the radius of the O star is increased, although this will result in an inconsistent depth of the secondary eclipse. A possible cause of this inconsistency is the inaccurate stellar wind structure very close to the core of the W-R star. For all the models in Table 2, $`\beta _1=1`$ was used in the modified beta-velocity law (Hillier & Miller 1999) (see Eq. 1) where $`\beta _1`$ is equivalent to $`\beta `$ in the classic beta-velocity law for the “inner” part of the stellar wind. To see if changing the structure of the inner W-R stellar wind affects the computed continuum polarization level, a model with $`\beta _1=2`$ (a slower wind) is also calculated. Unfortunately, the amplitude of the polarization curves (not shown here) near $`\varphi =0.5`$ did not change significantly for the new $`\beta _1`$ value. Changing the value of $`\beta _1`$ also did not influence the depth of the secondary eclipse in the IR light curve.

We note that STL2 were able to fit the observed $`Q`$ and $`U`$ variations around $`\varphi =0.5`$. There are two reasons for their success: 1. they did not fit the observed $`I`$ light curve consistently, and 2. they use a fixed orbital inclination angle. In their fitting procedure, the radius of the O star was used as a fitting parameter, but changing the O star radius changes the depth of the eclipse around $`\varphi =0.5`$. Consequently, the inclination angle must also be adjusted when the O star radius is changed, to maintain a consistency with the observed $`I`$ light curve. We should also note that the amplitudes of polarization curves around $`\varphi =0.5`$ are very sensitive to the orbital inclination. STL2 adopted $`i=78.9^{}`$ from Robert et al. (1990). However, if instead $`i=78^{}`$ from Cherepashchuk (1975) were used, they would obtain a fairly different set of parameters from the ones derived in their paper. Fortunately, the amplitudes of the double sine curves are not very sensitive to the adopted inclination angle. Considering the relative sensitivity of different methods on the inclination angle, the mass-loss rate estimation by fitting the double sine curves is more reliable than the one obtained by fitting the polarization variations near $`\varphi =0.5`$.

### 3.3 Magnitudes and Distance Moduli

For each model in Figs. 10 and 12, the absolute magnitude ($`M_v`$) of the system (W-R + O) is calculated, and placed in Table 3. $`M_v=5.60`$ from Model B is very similar to $`M_v=5.7`$ of Lundström & Stenholm (1984), but is far from the value of Hamann & Schwarz (1992), $`M_v=4.2`$. Model E gives $`M_v=5.15`$ which is slightly different from the values estimated by Lundström & Stenholm (1984) and Hamann & Schwarz (1992). Using the observed $`m_v=8.27`$ and $`A_v=2.79`$ from Lundström & Stenholm (1984), we have also estimated the distance modulus ($`DM`$) for each model, and summarized them in the same table. The values of $`DM`$ from Model B and Model E are $`11.1`$ and $`10.6`$ respectively. The former (Model B) agrees with the value given by Lundström & Stenholm (1984) ($`DM=11.01\pm 0.20`$), but it is slightly smaller than the value estimated by Marchenko et al. (1997) ($`DM=11.2711.47`$).

The absolute magnitudes of the O star corresponding to Models B and E are also listed in Table 3. They are $`M_v=5.20`$ and $`M_v=4.85`$ respectively. According to Howarth & Prinja (1989), the mean absolute magnitude of O6 V stars is $`M_v=5.2`$, and that for O6 III is $`M_v=6.6`$. Our magnitude of Model B agrees with the mean of their O6 V stars. Similarly, the average absolute magnitude of observed WN5 stars is $`M_v=4.1\left(\pm 0.8\right)`$ (van der Hucht 2001) which also agrees with our values $`M_v=4.33`$ (Model B) and $`M_v=3.58`$ (Model E).

The derived O star radii are $`R_\mathrm{O}=7.2R_{\mathrm{}}`$(Model B) and $`R_\mathrm{O}=6.9R_{\mathrm{}}`$(Model E). Both values are slightly smaller than the recent value $`R_O=8.5\pm 1.0R_{\mathrm{}}`$ estimated by STL2, and are not consistent with $`R_\mathrm{O}=10.0R_{\mathrm{}}`$ (Cherepashchuk et al. 1984) and $`R_\mathrm{O}=8.49.3R_{\mathrm{}}`$ (Moffat & Marchenko 1996). The rather small O star radii of the models and the model $`M_v`$ values compared to the mean observed $`M_v`$ suggest that the O star is more likely to be a main-sequence star rather than a giant star. A similar conclusion about the O star luminosity class was reached by STL2.

Vanbeveren et al. (1998) quoted Cherepashchuk (1975), ”The observed mass of the WNE component is $`9M_{\mathrm{}}`$; its progenitor must therefore have had a mass of $`30M_{\mathrm{}}`$. Thus, the age of the binary is about 7 million years. The OB companion is an O6 star. Since the age of a normal O6 star is about 1-2 million years, the OB star must have been rejuvenated, i.e. mass transfer must have occurred.” This argument suggests that *the O star is not a normal O star*. According to Vanbeveren et al. (1998), a quasi-conservative Roche lobe overflow model can produce a system which resembles V444 Cyg. At the end of the mass transfer, when the secondary has been entirely mixed, the O6 component may be an *over-luminous* star with $`M_\mathrm{O}=26M_{\mathrm{}}`$. Contrary to the last remark, our best fit models (Models B and E) suggest that the absolute magnitude of the O star is similar to the average of O6 V stars, as discussed earlier.

Unfortunately, we can not strongly conclude that either Model B or E is a more favorable model based on the comparisons with the observed magnitudes, distance moduli and the results from earlier works.

## 4 Conclusions

Using the 3-D Monte Carlo model developed in Paper I, combined with the multi-line non-LTE radiative model of Hillier & Miller (1998), we first investigated the basic behavior of the WN star diagnostic lines (He ii $`5412`$ Å and He i $`5876`$ Å) as a function of the phase of V444 Cyg (Figs. 46). Since the emission lines from the W-R star envelope are contaminated by the O star atmospheric absorption lines and the line emission from the shock-heated region, we first examined at which phase the profile shape of helium lines, especially He i $`5876`$ Å, is least affected by the contaminations. We then fit the He i $`5876`$ Å profile with the multi-line non-LTE radiative transfer model with spherical geometry to obtain a realistic model for the W-R star. The best orbital phase for this purpose is found to be around $`\varphi =0.8`$.

A summary of the derived parameters of V444 Cyg with our best fit models (Models B & E) is given in Table 3. Both models show excellent agreement with the observed light curves of $`I`$, $`Q`$ and $`U`$ Stokes parameters (Figs. 10 and 12). They also fit the spectrum of He ii $`5412`$ Å and He i $`5876`$ Å very well (Figs. 9 and 11). The effect of ignoring the presence of the bow shock in the continuum polarization calculation was also considered. This was done by first approximating the opacity in the bow shock region via fitting He i $`5876`$ Å line at $`\varphi =0.419`$ where the emission from the bow shock is most prominent (Fig 13). According to this model, the bow shock region has little or no effect on the continuum light curves at $`\lambda =5630`$Å for all $`I`$, $`Q`$ and $`U`$ (Fig. 14). Similarly, the presence of the O-star wind and the exact details of the bow shock are found to be unimportant for the continuum fits (Fig. 16). We also found that the spectroscopic solution, when combined with the strong O star spectrum, is very insensitive to the value of the W-R core radius (c.f., Fig. 9 and Fig. 11).

Small discrepancies are seen in the oscillation of the $`Q`$ and $`U`$ light curves (Figs. 10 and 12) near the secondary eclipse ($`\varphi =0.5`$), which are possibly due to the lack of knowledge of the exact opacity and emissivity in the bow shock region. This discrepancy could be a result of the wrong combination of the orbital inclination angle and the O star radius in our model. The shape of the polarization curves near $`\varphi =0.5`$ is very sensitive to the inclination angle, and is difficult to fit consistently with other observational features.

The most noticeable difference between our best model and the observations is in the IR light curves ($`2.2\mu \mathrm{m}`$) in Figs. 10 and 12. The model correctly predicts the depth of the primary eclipse, but that of the secondary eclipse is too shallow. We argued that the difference could be caused by the use of the wrong monochromatic luminosity ratio ($`q`$) at $`\lambda =2.2\mu \mathrm{m}`$.

The mass-loss rate of the W-R component determined in our analysis is $`\dot{M}_{\mathrm{WR}}=0.6\left(\pm 0.2\right)\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$, and this conclusion is insensitive to the model radius of the W-R star. The mass-loss rate of the W-R star obtained by STL2, $`\dot{M}_{\mathrm{WR}}=0.6\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$, is in good agreement with our value derived here. Our derived mass-loss rate lies between the values obtained by the orbital period change method (Table 1). The fits did not allow a unique value for the radius of the W-R star to be derived. The range of the volume filling factor for the W-R star atmosphere is estimated to be $`0.0500.075`$ for the corresponding range of the W-R star radius, $`5.02.5R_{\mathrm{}}`$.

Future work includes: 1. obtaining improved IR light curves for both models and observations, 2. fitting and modeling the UV spectrum and polarization, 3. checking the consistency between the polarization model and the hydrodynamical model more carefully, 4. investigating the cause for the asymmetry of the polarization light curve around the secondary eclipse, and 5. obtaining a tighter constraint on $`q`$. Future interferometer observations of this system will also prove invaluable to our further understanding.

###### Acknowledgements.

We would like to thank Dr. S. Marchenko and Dr. A. F. J. Moffat for providing us the helium spectrum sequence data and for giving us very constructive comments on the manuscript. This research was supported by NASA grant NAGW-3828.


text:
# The Brownian Web

## 1 Brownian Web: Characterization

We begin by defining three metric spaces: $`(\overline{}^2,\rho )`$, $`(\mathrm{\Pi },d)`$ and $`(,d_{})`$. The elements of the three spaces are respectively: points in space-time, paths with specified starting points in space-time and collections of paths with specified starting points. The BW will be an $`(,_{})`$-valued random variable, where $`_{}`$ is the Borel $`\sigma `$-field associated to the metric $`d_{}`$.

$`(\overline{}^2,\rho )`$ is the completion (or compactification) of $`^2`$ under the metric $`\rho `$, where

$$\rho ((x_1,t_1),(x_2,t_2))=\left|\frac{\mathrm{tanh}(x_1)}{1+|t_1|}\frac{\mathrm{tanh}(x_2)}{1+|t_2|}\right||\mathrm{tanh}(t_1)\mathrm{tanh}(t_2)|.$$
(1.1)
$`\overline{}^2`$ may be thought as the set of $`(x,t)`$ in $`[\mathrm{},\mathrm{}]\times [\mathrm{},\mathrm{}]`$ with all points of the form $`(x,\mathrm{})`$ identified (and similarly for $`(x,\mathrm{})`$). More precisely, it is the image of $`[\mathrm{},\mathrm{}]\times [\mathrm{},\mathrm{}]`$ under the mapping

$$(x,t)(\mathrm{\Phi }(x,t),\mathrm{\Psi }(t))(\frac{\mathrm{tanh}(x)}{1+|t|},\mathrm{tanh}(t)).$$
(1.2)
For $`t_0[\mathrm{},\mathrm{}]`$, let $`C[t_0]`$ denote the set of functions $`f`$ from $`[t_0,\mathrm{}]`$ to $`[\mathrm{},\mathrm{}]`$ such that $`\mathrm{\Phi }(f(t),t)`$ is continuous. Then define

$$\mathrm{\Pi }=\underset{t_0[\mathrm{},\mathrm{}]}{}C[t_0]\times \{t_0\},$$
(1.3)
where $`(f,t_0)\mathrm{\Pi }`$ then represents a path in $`\overline{}^2`$ starting at $`(f(t_0),t_0)`$. For$`(f,t_0)\mathrm{\Pi }`$, we denote by $`\widehat{f}`$ the function that extends $`f`$ to all $`[\mathrm{},\mathrm{}]`$ by setting it equal to $`f(t_0)`$ for $`t<t_0`$. Then we take

$$d((f_1,t_1),(f_2,t_2))=(\underset{t}{sup}|\mathrm{\Phi }(\widehat{f_1}(t),t)\mathrm{\Phi }(\widehat{f_2}(t),t)|)|\mathrm{\Psi }(t_1)\mathrm{\Psi }(t_2)|.$$
(1.4)
$`(\mathrm{\Pi },d)`$ is a complete separable metric space.

Let now $``$ denote the set of compact subsets of $`(\mathrm{\Pi },d)`$, with $`d_{}`$ the induced Hausdorff metric, i.e.,

$$d_{}(K_1,K_2)=\underset{g_1K_1}{sup}\underset{g_2K_2}{inf}d(g_1,g_2)\underset{g_2K_2}{sup}\underset{g_1K_1}{inf}d(g_1,g_2).$$
(1.5)
$`(,d_{})`$ is also a complete separable metric space.

Before stating our characterization theorem for the Brownian web, we need some definitions. For an $`(,_{})`$-valued random variable $`\overline{W}`$ (or its distribution $`\mu `$), we define the finite-dimensional distributions of $`\overline{W}`$ as the induced probability measures $`\mu _{(x_1,t_1;\mathrm{};x_n,t_n)}`$ on the subsets of paths starting from any finite deterministic set of points $`(x_1,t_1),\mathrm{},(x_n,t_n)`$ in $`^2`$. There are several ways in which the Brownian web can be characterized; they differ from each other primarily in the type of extra condition required beyond the finite-dimensional distributions. The characterization of the next theorem, or more precisely a variant discussed later in Remark 1.3, is the one most directly suited to the convergence results of Section 2; an alternative characterization in which the extra condition is a type of Doob separability property (see, e.g., Chap. 3 of ) is discussed in Remark 1.2. For the next theorem, we also define, for $`t0`$ and $`ab`$, the ($`\{0,1,\mathrm{},\mathrm{}\}`$-valued) random variable $`\eta (t_0,t;a,b)`$ as the number of distinct points in $`\times \{t_0+t\}`$ that are touched by paths in $`\overline{W}`$ which also touch some point in $`[a,b]\times \{t_0\}`$.

###### Theorem 1.1

There is an $`(,_{})`$-valued random variable $`\overline{W}`$ whose distribution $`\mu `$ is uniquely determined by the following two properties:

* its finite-dimensional distributions are those of coalescing Brownian motions (with unit diffusion constant), and
* for $`\mathrm{}<t_0<\mathrm{}`$, $`0<t<\mathrm{}`$, $`\mathrm{}<ab<\mathrm{}`$,

$$E_\mu (\eta (t_0,t;a,b))=1+\frac{ba}{\sqrt{\pi t}}.$$
(1.6)

###### Remark 1.2

Implicit in condition (i) of the theorem is that starting from any deterministic point, there is almost surely only a single path in $`\overline{W}`$. Condition (ii) can be replaced by the separability property that there is a deterministic dense countable set $`𝒟`$ of space-time starting points, such that almost surely, $`\overline{W}`$ is the closure in $`(\mathrm{\Pi },d)`$ of the set of paths starting from the points of $`𝒟`$. It should be noted that $`(,_{})`$-valued random variables, satisfying condition (i) but not condition (ii) or its separability alternative, can occur naturally. Such a process (closely related to the “Double Brownian Web” of Remark 1.4 below), where the counting variable $`\eta `$ is infinite with strictly positive probability, will be studied elsewhere and shown to arise as the scaling limit of stochastic flows, extending earlier work of Piterbarg .

Sketch of Proof of Theorem 1.1. The construction of the Brownian web (i.e., the existence of such a $`\overline{W}`$) begins as in  with the construction of a set $`W`$ of coalescing Brownian paths starting from a deterministic dense countable set $`𝒟`$ of space-time starting points. This skeleton $`W=\{\stackrel{~}{W}_1,\stackrel{~}{W}_2,\mathrm{},\}`$ is a random subset of $`\mathrm{\Pi }`$ that is constructed by deterministically ordering the points of $`𝒟`$ as $`(x_1,t_1),(x_2,t_2),\mathrm{}`$, then defining $`W_j=(x_j+B_j(tt_j),t_j)\mathrm{\Pi }`$ where the $`B_j`$’s are independent standard Brownian motions, and finally using the ordering to inductively define $`\stackrel{~}{W}_j\mathrm{\Pi }`$ by following $`W_j`$ until it meets some $`\stackrel{~}{W}_k`$ with $`k<j`$ after which point it follows $`\stackrel{~}{W}_k`$.

The next several steps of the construction are to show that the closure $`\overline{W}`$ in $`(\mathrm{\Pi },d)`$ of this BW-skeleton is compact, that the distribution of $`\overline{W}`$ does not depend on the choice of $`𝒟`$ or its ordering, and that $`\overline{W}`$ satisfies (i) and (ii) of Theorem 1.1 above. The compactness can be proved in a number of ways; one of these is to verify a condition, as in (2.1) below, but with $`\mu _\delta `$ replaced by the distribution of $`\{\stackrel{~}{W}_1,\mathrm{},\stackrel{~}{W}_m\}`$ and the sup over $`\delta `$ replaced by a sup over $`m`$ (and then argue as at the beginning of the proof of Theorem 2.1 below, eventually invoking the Arzelà-Ascoli theorem). To verify the said condition, one argues as in the last two paragraphs of the proof of Theorem 2.3 below. The argument actually involves only a single bound like (2.6), which is obtained in the same way as in the proof of Theorem 2.3. We remark that by considering the quantity $`\stackrel{~}{g}(t,u)`$, as in (2.1), but with $`u=t^\xi `$, one can show not just compactness of $`\overline{W}`$, but also Hölder continuity with any exponent $`\xi <1/2`$ for all the paths of $`\overline{W}`$.

The lack of dependence of the distribution on the choice of $`𝒟`$ or its ordering follows fairly directly after verifying property (i) for $`\overline{W}`$. Property (i) itself follows by a trapping argument about a deterministic point $`(\overline{x},\overline{t})`$ (and similarly for finitely many points) and any sequence $`(\overline{x}_i,\overline{t}_i)=(x_{j(i)},t_{j(i)})`$ from $`𝒟`$ converging to $`(\overline{x},\overline{t})`$ as $`i\mathrm{}`$: that (even if $`j(i)`$ is nondeterministic and regardless of whether $`t_{j(i)}>\overline{t}`$ or $`t_{j(i)}\overline{t}`$) for large $`i`$, $`\stackrel{~}{W}_{j(i)}`$ is (with probability very close to one) trapped between $`\stackrel{~}{W}_k`$ and $`\stackrel{~}{W}_k^{}`$ for deterministic $`k,k^{}`$ with $`x_k<\overline{x}<x_k^{}`$, $`t_k<\overline{t}`$, $`t_k^{}<\overline{t}`$, $`x_k`$ and $`x_k^{}`$ close to $`\overline{x}`$ and $`t_k,t_k^{}`$ even closer to $`\overline{t}`$ so that $`\stackrel{~}{W}_{j(i)}`$ (with probability very close to one) quickly coalesces with both $`\stackrel{~}{W}_k`$ and $`\stackrel{~}{W}_k^{}`$; thus $`\stackrel{~}{W}_{j(i)}`$ converges almost surely as $`i\mathrm{}`$ to a path (independent of the specific sequence $`(\overline{x}_i,\overline{t}_i)`$) that is distributed as a Brownian motion starting from $`(\overline{x},\overline{t})`$.

Verifying property (ii) is somewhat indirect. First, one shows that the random variable $`\eta `$ for our constructed $`\overline{W}`$ is almost surely finite with a finite mean which, by the translation invariance in space and time that results from the lack of dependence on $`𝒟`$, must be of the form $`\mathrm{\Lambda }(ba,t)`$. Second, the specific evaluation of $`\mathrm{\Lambda }`$ as given on the righthand side of (1.6) is carried out. As the explicit expression for $`\mathrm{\Lambda }`$ will not actually be used in our convergence results of the next section (see Remark 1.3), we can and will use those convergence results in the evaluation of $`\mathrm{\Lambda }`$.

The first part of the verification of (ii) is a consequence of an inequality,

$$P(\eta (t_0,t;a,b)k)[P(\eta (t_0,t;a,b)2)]^{k1}=[\mathrm{\Theta }(ba,t)]^{k1},$$
(1.7)
where $`\mathrm{\Theta }(ba,t)`$ is the probability that two independent Brownian motions starting at a distance $`ba`$ apart at time zero will have met by time $`t`$ (which itself can be expressed in terms of a single Brownian motion). The inequality in (1.7) is first derived for finite subsets $`\{\stackrel{~}{W}_1,\stackrel{~}{W}_2,\mathrm{},\stackrel{~}{W}_m\}`$ of the skeleton, and thus for the whole skeleton. For the whole skeleton and its closure $`\overline{W}`$, the equality in (1.7) is seen to be valid by choosing the countable set $`𝒟`$ so that its first two points are $`(a,t_0)`$ and $`(b,t_0)`$. Then the inequality is extended to $`\overline{W}`$ from the skeleton by a limit/approximation argument which uses that $`\{K:\stackrel{~}{\eta }(K)k\}`$ is open in $`(,d_{})`$, where $`\stackrel{~}{\eta }`$ is the modification of $`\eta `$ that counts points in $`\times \{t_0+t\}`$ touched by paths in $`K`$ which touch $`(a\stackrel{~}{\epsilon }_1,b+\stackrel{~}{\epsilon }_1)\times \{t_0+\stackrel{~}{\epsilon }_2\}`$ and start earlier than $`t_0+\stackrel{~}{\epsilon }_2`$.

The second part of the verification of (ii), in which $`\mathrm{\Lambda }`$ is explicitly evaluated, is a consequence of all the following: a result of Bramson and Griffeath  on the large-time asymptotics of mean interparticle distance in coalescing random walks, the conversion of that result by standard arguments to asymptotics for the mean of the rescaled random walk version of the counting variable $`\eta `$, convergence of the distribution of $`\eta `$ in the scaling limit (see Remark 2.2), and finally the analogue of (1.7) for coalescing walks (see (2.5)) which implies uniform integrability of $`\eta `$ as $`\delta 0`$ and hence convergence of the mean of $`\eta `$.

It remains to show that conditions (i) and (ii) for a measure $`\mu ^{}`$ on $`(,_{})`$ together imply that $`\mu ^{}`$ equals the distribution $`\mu `$ of the constructed Brownian web $`\overline{W}`$. Let us denote by $`X^{}`$ the $`(,_{})`$-valued random variable distributed by $`\mu ^{}`$ and by $`\eta ^{}`$ the counting random variable appearing in condition (ii) for $`\mu ^{}`$. Choose some deterministic dense countable subset $`𝒟`$ and consider the countable collection $`W^{}`$ of paths of $`X^{}`$ starting from $`𝒟`$. By condition (i), $`W^{}`$ is equidistributed with our constructed Brownian skeleton $`W`$ (based on the same $`𝒟`$) and hence the closure $`\overline{W^{}}`$ of $`W^{}`$ in $`(\mathrm{\Pi },d)`$ is a subset of $`X^{}`$ that is equidistributed with our constructed Brownian web $`\overline{W}`$. To complete the proof, we will use condition (ii) to show that $`X^{}\overline{W^{}}`$ is almost surely empty by using the fact that the counting variable $`\eta ^{}`$ for $`\overline{W^{}}`$ already satisifies condition (ii) since $`\overline{W^{}}`$ is distributed as a Brownian web. If $`X^{}\overline{W^{}}`$ were nonempty (with strictly positive probability), then there would have to be some rational $`t_0,t,a,b`$ for which $`\eta ^{}>\eta ^{}`$. But then

$$P(\eta ^{}(t_0,t;a,b)>\eta ^{}(t_0,t;a,b))>\mathrm{\hspace{0.17em}0}$$
(1.8)
for some rational $`t_0,t,a,b`$ and so condition (ii) for $`\eta ^{}`$ would not be valid for that $`t_0,t,a,b`$.

###### Remark 1.3

The proof of Theorem 1.1 makes clear that the idea behind (i) and (ii) together implying uniqueness of the distribution is that (i) implies sufficiently many paths and (ii) implies no extraneous ones. Thus condition (i) can be weakened to the existence of a subset of paths distributed as the coalescing Brownian motions of the skeleton $`W`$ (for any deterministic dense countable $`𝒟`$) and condition (ii) can also be modified, e.g., by replacing the equality in (1.6) by an inequality ($``$) and by replacing an (in)equality for the mean by one for the distribution. Similarly, in applying our characterization results to obtain convergence criteria as we do in Theorem 2.1, an explicit expression for the mean as given in the righthand side of (1.6) or an explicit expression for the distribution is not needed; i.e., to verify that an $`X^{}`$ is equidistributed with our explicitly constructed Brownian web $`\overline{W}`$, condition (ii) for the $`\eta ^{}`$ of $`X^{}`$ can be replaced by the condition that the distribution of $`\eta ^{}`$ equal (or only is stochastically dominated by) the distribution of the $`\eta `$ of $`\overline{W}`$.

###### Remark 1.4

In the graphical representation of Harris for the one-dimensional voter model , coalescing random walks forward in time and coalescing dual random walks backward in time (with forward and backward walks not crossing each other) are constructed simultaneously (see, e.g., the discussion in ). The simultaneous construction of forward and (dual) backward Brownian motions was emphasized in  and their approach and results can be applied to extend both our characterization and convergence results to the Double Brownian Web (DBW) which includes simultaneously the forward BW and its dual backward BW. We note that in the DBW, the $`\eta `$ of (1.6) equals $`1+\eta ^{\text{dual}}`$, where $`\eta ^{\text{dual}}`$ is the number of distinct points in $`[a,b]\times \{t_0\}`$ touched by backward paths which also touch $`\times \{t_0+t\}`$.

###### Remark 1.5

As in , space-time points $`(x,t)`$ can be characterized by the number of locally disjoint paths $`m_{\text{in}}`$ (resp., $`m_{\text{out}}`$) of the BW entering (resp., leaving) that point from earlier (resp., to later) times. The corresponding dual BW characterization has $`m_{\text{in}}^{\text{dual}}=m_{\text{out}}1`$ and $`m_{\text{out}}^{\text{dual}}=m_{\text{in}}+1`$. Generic (e.g., deterministic) points have $`(m_{\text{in}},m_{\text{out}})=(0,1)`$. Almost surely, there are nongeneric points of type $`(0,2),(0,3),(1,1),(1,2)`$ and $`(2,1)`$ but no others. We note that as in , ruling out points of higher type uses improvements of (1.7) for $`k>2`$. Type $`(2,1)`$ (resp., $`(0,3)`$) points are those where coalescing (resp., dual coalescing) occurs. Type $`(1,2)`$ points are particularly interesting in that the single incident path continues along exactly one of the two outward paths — with the choice determined intrinsically rather than by some convention.

## 2 Convergence to the Brownian Web

Let $`X_\delta `$ be an $`(,_{})`$-valued random variable indexed by $`\delta (0,1]`$, with distribution $`\mu _\delta `$. We present criteria sufficient to insure convergence in distribution as $`\delta 0`$ of $`X_\delta `$ to the Brownian web $`\overline{W}`$, in the setting where the $`X_\delta `$’s have coalescing paths; for simplicity, we will not present here more general criteria that do not require the coalescing property. We next introduce the various conditions on $`\mu _\delta `$ which together will imply convergence.

The first condition will guarantee tightness of the $`\mu _\delta `$’s. Let $`R(x_0,t_0;u,t)`$ denote the rectangle $`[x_0u/2,x_0+u/2]\times [t_0,t_0+t]`$ in $`^2`$. We call $`\{x_0\pm u/2\}\times [t_0,t_0+t]`$ its right and left boundaries. For $`t>0,u>\sqrt{t}`$, define $`A_{t,u}(x_0,t_0)`$ to be the event (in $`_{}`$) that $`K`$ (in $``$) contains a path touching both $`R(x_0,t_0;\sqrt{t},t)`$ and (at a later time) the left or right boundary of the bigger rectangle $`R(x_0,t_0;u,2t)`$; see Figure 2. Our tightness condition is

$$(T_1)\stackrel{~}{g}(t,u)t^{3/2}\underset{\delta >0}{sup}\underset{x_0,t_0}{sup}\mu _\delta (A_{t,u}(x_0,t_0))0\text{ as }t0+\text{ for fixed }u>0.$$
(2.1)
Our second condition will guarantee a weakened version of (i) in Theorem 1.1 (see Remark 1.3) for any limit $`\mu `$ of $`\mu _\delta `$. Let $`𝒟`$ be any deterministic countable dense set of points in $`^2`$. The condition concerns the existence for each $`\delta >0`$ and $`y𝒟`$ of measurable (on the probability space of $`X_\delta `$) single-path valued random variables $`\theta _\delta ^y(\omega )X_\delta (\omega ):`$

$`(I_1)`$ There exist such $`\theta _\delta ^yX_\delta `$ satisfying: for any deterministic $`y_1,\mathrm{},y_m𝒟`$, $`\theta _\delta ^{y_1},\mathrm{},\theta _\delta ^{y_m}`$ converge in distribution as $`\delta 0`$ to coalescing Brownian motions (with unit diffusion constant) starting at $`y_1,\mathrm{},y_m`$.

Our next two conditions will together guarantee (when $`X_\delta `$ is coalescing) a version of (ii) in Theorem 1.1 (see Remark 1.3). For $`\mathrm{}<t_0<\mathrm{}`$ and $`0<t<\mathrm{}`$,

$`(B_1)`$ $`\underset{\delta 0}{lim\; sup}\underset{a}{sup}\mu _\delta (\eta (t_0,t;a,a+ϵ)2)0\text{ as }ϵ0+;`$ (2.2)
$`(B_2)`$ $`ϵ^1`$ $`\underset{\delta 0}{lim\; sup}\underset{a}{sup}\mu _\delta (\eta (t_0,t;a,a+ϵ)3)0\text{ as }ϵ0+.`$ (2.3)
###### Theorem 2.1

Suppose $`X_\delta `$ for $`0<\delta 1`$ are $`(,_{})`$-valued random variables with coalescing paths. If $`T_1,I_1,B_1`$ and $`B_2`$ all hold, then the distributions $`\mu _\delta `$ of $`X_\delta `$ converge weakly as $`\delta 0`$ to the distribution $`\mu _{\overline{W}}`$ of the standard Brownian web.

Sketch of Proof of Theorem 2.1. We first explain why $`T_1`$ implies tightness. Let $`g_\delta (t,u)`$ denote the sup over $`x_0,t_0`$ of $`\mu _\delta (A_{t,u})`$ as in (2.1). This represents an upper bound on the $`\mu _\delta `$-probability that there is some path $`(f,t^{})`$ passing through some point $`(x^{},t^{})=(f(t^{}),t^{})`$ in the deterministic $`\sqrt{t}\times t`$ rectangle $`R(x_0,t_0;\sqrt{t},t)`$ located at any $`(x_0,t_0)`$, such that for some $`t^{\prime \prime }[t^{},t^{}+t]`$ the spatial increment $`|f(t^{\prime \prime })f(t^{})|u`$ even though the time increment is $`t`$. Now taking a large $`L\times T`$ space-time rectangle centered at the origin and covering it with $`O(LT/t^{3/2})`$ $`\sqrt{t}\times t`$ small rectangles, we see that $`LT\stackrel{~}{g}(t,u)`$ represents an upper bound on the $`\mu _\delta `$-probability (for any $`\delta `$) that some path has $`|f(t^{\prime \prime })f(t^{})|u`$ while $`t^{\prime \prime }t^{}t`$ with $`(f(t^{}),t^{})`$ anywhere in the large rectangle. We next choose sequences $`u_n0`$, $`L_n\mathrm{}`$, $`T_n\mathrm{}`$ and then $`t_n0`$ sufficiently rapidly that $`L_nT_n\stackrel{~}{g}(t_n,u_n)`$ is summable. Now moving to the compactified space-time $`\overline{}^2`$ (and using the notation of (1.2)), it follows that there are sequences $`\varphi _n,\psi _n0`$ so that for large enough $`n`$, with $`\mu _\delta `$-probability close to one (for any $`\delta `$), $`|\mathrm{\Psi }(t^{\prime \prime })\mathrm{\Psi }(t^{})|\psi _n`$ implies $`\mathrm{\Phi }(f(t^{\prime \prime }),t^{\prime \prime })\mathrm{\Phi }(f(t^{}),t^{})|\varphi _n`$. This equicontinuity with probability close to one (for any $`\delta `$) combined with a version of the Arzelà-Ascoli theorem leads to the paths, as elements of $`(\mathrm{\Pi },d)`$, belonging to a compact subset $`K_{\stackrel{~}{\epsilon }}`$ of $`(\mathrm{\Pi },d)`$ with $`\mu _\delta `$-probability $`1\stackrel{~}{\epsilon }`$ (for any $`\delta `$), which implies tightness because the collection of compact subsets of $`K_{\stackrel{~}{\epsilon }}`$ is itself a compact set in $`(,d_{})`$.

Tightness implies that every subsequence of $`\mu _\delta `$ has a sub-subsequence converging weakly to some $`\mu `$. To complete the proof, we need to show that any such $`\mu `$ equals $`\mu _{\overline{W}}`$. To do this, we will show that $`\mu `$ satisifies the two characterization properties of Theorem 1.1, as modified in Remark 1.3. Combining condition $`I_1`$ with convergence in distribution (along a subsequence) of $`X_\delta `$ to some $`X`$ distributed by $`\mu `$, we see that we can realize $`X`$ on some probability space so that it contains paths starting from the points of $`𝒟`$ distributed as coalescing Brownian motions. This is just the desired (weakened version of) property(i) in Theorem 1.1. Indeed, this shows that $`X`$ contains an $`X^{}`$ that has the Brownian web distribution.

To complete the proof we use conditions $`B_1`$ and $`B_2`$. Note first that by limit/approximation arguments these two conditions (without the $`lim\; sup`$ over $`\delta `$) are valid with $`\mu `$ replacing $`\mu _\delta `$. For fixed $`t_0,t,a,b`$, we now consider $`M+1`$ equally spaced points, $`z_j=(a+j(ba)/M,t_0)`$ for $`j=0,\mathrm{},M`$. For the random $`X`$ we will denote the counting variable $`\eta (t_0,t;a,b)`$ by $`\eta `$, and the corresponding variable for $`X^{}`$ by $`\eta ^{}`$. We also want to count the number of points on $`\times \{t_0+t\}`$ that are touched by paths of $`X`$ that also touch $`\{z_0,\mathrm{},z_M\}`$ and we will denote these variables for $`X`$ and $`X^{}`$ by $`\eta _M`$ and $`\eta _M^{}`$. Of course, $`\eta \eta _M`$ and $`\eta \eta _M^{}`$. By condition $`B_1`$ (for $`\mu `$) applied to small intervals about each of the $`z_j`$’s, it follows that $`\eta _M=\eta _M^{}`$ almost surely. Applying condition $`B_2`$ (for $`\mu `$) to the $`M`$ spatial intervals $`[z_{j1},z_j]`$ of length $`ϵ=(ba)/M`$, and using the coalescing (or at least non-crossing) property of $`X`$ that it inherits from the $`X_\delta `$’s, it follows that

$$P(\eta >\eta _M^{})=P(\eta >\eta _M)0\text{ as }M\mathrm{}.$$
(2.4)
Thus $`P(\eta >\eta ^{})=0`$ so that the distribution of $`\eta ^{}`$ is stochastically dominated by (and hence equal to) the distribution of $`\eta `$. This gives the desired (modified version of) property(ii) in Theorem 1.1 and completes the proof.

###### Remark 2.2

The arguments used in the proof of Theorem 2.1 also show that under the same four conditions $`T_1,I_1,B_1`$ and $`B_2`$, the counting random variables $`\eta _\delta (t_0,t;a,b)`$ for $`X_\delta `$ converge in distribution to the Brownian web counting variable $`\eta _{\overline{W}}(t_0,t;a,b)`$. If one also has uniform integrability as $`\delta 0`$, then the means converge to the mean of $`\eta _{\overline{W}}`$.

To apply Theorem 2.1 to random walks, we begin by precisely defining $`Y`$ (resp., $`\stackrel{~}{Y}`$), the set of all discrete (resp., continuous) time coalescing random walks on $``$. The sets of rescaled walks, $`Y^{(\delta )}`$ and $`\stackrel{~}{Y}^{(\delta )}`$, are then obtained by the usual rescaling of space by $`\delta `$ and time by $`\delta ^2`$. The (main) paths of $`Y`$ are the discrete-time random walks $`Y_{y_0,s_0}`$, as described in the Introduction and shown in Figure 1, with $`(y_0,s_0)=(i_0,j_0)\times `$ arbitrary except that $`i_0+j_0`$ must be even. Each random walk path goes from $`(i,j)`$ to $`(i\pm 1,j+1)`$ linearly. In addition to these, we add some boundary paths so that $`Y`$ will be a compact subset of $`\mathrm{\Pi }`$. These are all the paths of the form $`(f,s_0)`$ with $`s_0\{\mathrm{},\mathrm{}\}`$ and $`f\mathrm{}`$ or $`f\mathrm{}`$. Note that for $`s_0=\mathrm{}`$ there are two different paths starting from the single point at $`s_0=\mathrm{}`$ in $`\overline{}^2`$.

The continuous time $`\stackrel{~}{Y}`$ can be defined similarly, except that here $`y_0`$ is any $`i_0`$ and $`s_0`$ is arbitrary in $``$. Continuous time walks are normally seen as jumping from $`i`$ to $`i\pm 1`$ at the times $`T_k^{(i)}(\mathrm{},\mathrm{})`$ of a rate one Poisson process. If the jump is, say, to $`i+1`$, then our polygonal path will have a linear segment between $`(i,T_k^{(i)})`$ and $`(i+1,T_k^{}^{(i+1)})`$, where $`T_k^{}^{(i+1)}`$ is the first Poisson event at $`i+1`$ after $`T_k^{(i)}`$. Furthermore, if $`T_k^{(i_0)}<s_0<T_{k+1}^{(i_0)}`$, then there will be a constant segment in the path before the first nonconstant linear segment. If $`s_0=T_k^{(i_0)}`$, then we take two paths: one with an initial constant segment and one without.

###### Theorem 2.3

Each of the collections of rescaled coalescing random walk paths, $`Y^{(\delta )}`$ (in discrete time) and $`\stackrel{~}{Y}^{(\delta )}`$ (in continuous time) converges in distribution to the standard Brownian web as $`\delta 0`$.

Sketch of Proof of Theorem 2.3. By Theorem 2.1, it suffices to verify conditions $`T_1,I_1,B_1`$ and $`B_2`$. We will save the tightness condition $`T_1`$ for last as it is the messiest to verify, at least in the continuous time case of $`\stackrel{~}{Y}^{(\delta )}`$.

Condition $`I_1`$ is basically a consequence of the Donsker invariance principle, as already noted in the Introduction. Conditions $`B_1`$ and $`B_2`$ follow from the coalescing walks version of the inequality of (1.7), which is

$$\mu _\delta (\eta (t_0,t;a,a+ϵ)k)[\mu _\delta (\eta (t_0,t;a,a+ϵ)2)]^{k1}.$$
(2.5)
Taking the sup over $`a`$ and the $`lim\; sup`$ over $`\delta `$ and using standard random walk arguments produces an upper bound of the from $`C_k(ϵ/\sqrt{t})^{k1}`$ which yields $`B_1`$ and $`B_2`$ as desired.

It remains to verify $`T_1`$. We will sketch the arguments for the continuous time $`\stackrel{~}{Y}^{(\delta )}`$; the discrete time $`Y^{(\delta )}`$ is easier and corresponds to a portion of the continuous time arguments. As in the proof of Theorem 2.1, we denote by $`g_\delta (t,u)`$ the sup over $`x_0,t_0`$ of $`\mu _\delta (A_{t,u})`$, where $`\mu _\delta `$ now denotes the distribution of $`\stackrel{~}{Y}^{(\delta )}`$. For the continuous time case (and for $`u1`$ and $`\sqrt{t}`$ much smaller than $`u`$), we will obtain a $`\delta `$-independent bound on $`g_\delta (t,u)`$ that will yield $`T_1`$ by first obtaining, as we explain below, separate upper bounds in three regions of $`\delta `$-values that depend on $`t,u`$:

$`C_1\mathrm{exp}(C_2u/\sqrt{t})`$ $`\text{for }D_0/\sqrt{t}\delta ^1,`$ (2.6)
$`(2t/\delta ^2)^{u/(3\delta )}`$ $`\text{for }6/u\delta ^1D_0/\sqrt{t},`$ (2.7)
$`C_3(t/u^2)^2`$ $`\text{for }\delta ^16/u.`$ (2.8)
Together, these bounds (with $`D_0`$ chosen appropriately) yield

$$g_\delta (t,u)C_4t^2/u^4,$$
(2.9)
which gives $`T_1`$ as desired.

The first region of $`\delta `$-values corresponds to a spatial interval of width $`\sqrt{t}`$ being multiple lattice spacings $`\delta `$ wide and a spatial interval of width $`u`$ being multiple $`\sqrt{t}`$-intervals wide. The bound (2.6) comes about because the event $`A_{t,u}`$ is prevented if between the small rectangle and both the left and right boundaries of the larger rectangle (see Figure 2), there is a random walk path that stays within some spatial $`\sqrt{t}`$-interval between times $`t_0`$ and $`t_0+2t`$. The second region corresponds to two (or more) spatial lattice sites between the small rectangle and the left (or right) boundary of the larger rectangle. The bound here comes from preventing $`A_{t,u}`$ by having a random walk path stay between two adjacent spatial lattice sites between times $`t_0`$ and $`t_0+2t`$. The third bound comes from preventing $`A_{t,u}`$ by not having the Poisson process occurences at adjacent spatial lattice sites, $`T_k^{(i)}`$ and $`T_k^{}^{(i\pm 1)}`$, too close together in time. This completes our sketch of the proof.

Acknowledgments. Research partially supported by FAPESP and CNPq (Brasil), MURST (Italia), and NSF (U. S. A. ). The authors thank S. R. S. Varadhan for useful discussions.