text:
# The 1997 periastron passage of the binary pulsar PSR B1259-63

## 1 Introduction

The pulsar PSR B1259$``$63 is a member of a unique binary system. Discovered using the Parkes telescope in a survey of the Galactic plane at 1.5 GHz \[Johnston et al. 1992a\], it was shown by Johnston et al. (1992b) to be in a highly eccentric 3.4-yr orbit with a 10th-magnitude Be star, SS 2883.

Extensive observations of this system have been made, particularly near the previous periastron passage on 1994 January 9. These include observations at radio wavelengths \[Johnston et al. 1996, McClure-Griffiths et al. 1998, Wex et al. 1998, Johnston et al. 1999\], optical \[Johnston et al. 1994\], in the X–ray \[Cominsky, Roberts & Johnston 1994, Kaspi et al. 1995, Hirayama et al. 1996\] and in $`\gamma `$-rays \[Grove et al. 1995, Hirayama et al. 1999\]. These observations have led to a better understanding of the physical conditions in the wind and disc of the Be star \[Tavani, Arons & Kaspi 1994, King & Cominsky 1994, Melatos, Johnston & Melrose 1995, Tavani & Arons 1997, Johnston et al. 1999, Ball et al. 1999\].

The understanding of the system is thus as follows: SS 2883 is a (main-sequence) B2e star with mass $``$10 M and radius $``$6 R. Its characteristic emission disc extends to at least 20 R in the optical, similar to the distance between the pulsar and Be star at periastron. It is likely that the disc is highly tilted with respect to the orbital plane; the orbital plane is itself at 36$`\mathrm{°}`$ with respect to the line of sight. Due to the misalignment of the emission disc with the orbital plane, the pulsar goes behind the disc (with respect to the observer) near $`𝒯`$ –16 and emerges from the front of the disc near $`𝒯`$ +14. (Here, $`𝒯`$ is the epoch of periastron). The disc material is extremely dense near the stellar surface but falls off rapidly as a function of distance. The polar wind, by contrast is $``$104 times less dense than the disc material at the stellar surface and has an $`r^2`$ density profile.

The radio observations of the pulsed emission from the 1994 periastron were rather crude in the sense that frequency agility was limited, with frequency switching requiring manual receiver changes and hence only possible on a $``$1 day timescale. Furthermore, only relatively few observations were obtained before and after the eclipse of the pulsar. In 1994 there was no real-time pulse folding software which limited our best observing strategy in that the off-line data reduction significantly lagged behind the observing. However, these observations showed dispersion measure (DM) and scintillation parameter changes on timescales of at least as short as one day and demonstrated that these changes were highly asymmetric with respect to periastron \[Johnston et al. 1996\].

In this paper we report on observations of the pulsed radio emission made around the latest periastron passage of the pulsar on 1997 May 29. Section 2 describes the observations. In Section 3 we present the results and in Section 4 discuss their implications.

## 2 Observations

In 1996, a major upgrade was made to the focus cabin of the Parkes radio-telescope which allowed several receivers to be placed on a movable platform. This enables the choice of receiver to be made from the control room and allows for frequency changes in only a few minutes.

We observed PSR B1259–63 regularly with the 64-m Parkes radio telescope at frequencies of 1.4, 4.8 and 8.4 GHz \[Wex et al. 1998\]. The observations were intensified around the periastron passage of the pulsar on 1997 May 29. At each frequency the receiver consists of a dual-channel, cryogenically cooled system sensitive to two orthogonal linear polarizations.

Data were recorded simultaneously using two recording systems. The Caltech Fast Pulsar Timing Machine (FPTM) is a correlator capable of fast sampling \[Navarro 1994\]. The signals from the two polarizations are down-converted to intermediate frequency and 2-bit digitised by the FPTM which computes correlation functions (CFs) with 128 lags over a 128 MHz bandwidth. The system can be used either to record full Stokes parameters for one 128 MHz bandwidth or total intensity for two independent bandwidths. For the higher frequencies we always recorded the data in full Stokes mode, with the band centered at 4800.25 and 8400.25 MHz. In the 20cm band, we either recorded the full Stokes parameters at 1351.0 MHz or two independent bandwidths centered at 1274.75 and 1396.0 MHz. The CFs are then hardware integrated at the apparent pulsar period ($``$47 ms) and the data recorded every 60 s. A graphical output displays the pulse profile every integration cycle. Thus, changes in observing strategy could be made in real-time.

Simultaneously, we also recorded data using a filterbank detection system in an identical fashion to that described elsewhere (e.g. Johnston et al. 1996). The signals from the receivers were down-converted to an intermediate frequency and then passed into a filterbank which consisted of 64 channels, each 5 MHz wide for a total of 320 MHz. Full Stokes parameters were recorded every 0.6 ms and written to magnetic tape for off-line analysis. The centre frequencies were 1374, 4748 and 8348 MHz.

For flux density calibration, observations were made of Hydra A at 1.4, 4.8 and 8.4 GHz once every few days at the telescope. These measurements yielded system flux densities of $``$28, $``$100 and $``$100 Jy at the three frequencies. Before every observation of the pulsar, a pulsed noise source, injected at 45 degrees to the signal probes, was recorded. This allows for calibration of the system gains and phases.

### 2.1 Derivation of the Dispersion Measure

Over a period of a few hours we obtained observations at 5 different frequencies: 1.2, 1.4, 1.5, 4.8 and 8.4 GHz. The three frequencies in the 20-cm band were obtained simultaneously. The observation would then be followed by one at either 4.8 and 8.4 GHz before again observing the 20-cm band. Typical observations were 1-hr in duration. Pulse time-of-arrivals were obtained in the standard manner (e.g. Wex et al. 1998). The DM was obtained by assuming the DM contribution from the 1.5 kpc path through the interstellar medium is 147.6 cm-3pc and does not change with time. Then, the extra DM contribution from the wind of the Be star could be computed by the relative differences in the residuals between the different frequencies. Having closely spaced frequencies in the 20-cm band allows any ambiguities to be resolved (i.e. where the $`\mathrm{\Delta }`$DM could impose more than 1 phase turn on the residual). Typical errors using this method are small ($``$0.1 cm-3pc), much smaller than those obtained simply from fitting the delays in frequency-phase space over a single observing band.

### 2.2 Derivation of the Scattering Time

Scatter-broadening of a pulse can be approximated by convolving the intrinsic pulse shape with a truncated exponential function. The situation for PSR B1259–63 is somewhat complicated by the presence of two strong pulses, only 20 ms apart in time \[Manchester & Johnston 1995\]. For large scattering times, the scatter tail of the first pulse overlaps the second pulse and this must be taken into account in the derivation of the scattering time. We used the Levenberg-Marquart method (LMM) to perform a least-squares fit to the observed profile given a ‘standard’ (unscattered) profile. The fit returns the 4 free parameters : amplitude, baseline and decay time of the exponential and the offset in time between the profiles. Generally, for scattering times less than 7 ms, the fits are good. Above this, however, the pulses tend to be weak, and the extremely scattered and weak pulses make the fitting less robust. This is reflected in much larger errors bars on the fitted parameters.

## 3 Results

The broad overview of the observational data is as follows. A change in the DM over and above that from the ISM was first observed at $`𝒯`$ –48. Over a period of $``$20 days, both the DM and the pulse scattering showed an overall increase, and strong hourly variations. The pulsed emission was last detected on $`𝒯`$ –17. Observations made on $`𝒯`$ –15 showed no pulsed emission at any frequency, a result confirmed with independent observations on the Australian Telescope Compact Array \[Johnston et al. 1999\]. The pulses were re-detected at Parkes on $`𝒯`$ +16, following a non-detection two days earlier. Post-periastron, no change in DM or pulse-scattering could be measured, although significant changes in rotation measure and polarization properties were recorded. Significant flux variability was seen both prior to and following the eclipse. The general trend though was towards lower fluxes in the few days before and after the eclipse.

Figure 1 shows the DM changes and the scattering time, $`\tau `$, as a function of epoch prior to periastron. Figure 2 shows a shorter time span of only 8 days prior to the eclipse.

The figures show a wealth of detail. Points to note are that the DM changes on short timescales ($``$1 hr or less) by $``$1 cm-3pc. There is a slow increase in the DM from $`𝒯`$ –50 until –27. The DM then slowly drops, before rising steeply again at $`𝒯`$ –17. The pulsar was undetected following these observations until $`𝒯`$ +16. Both the daily variations and the drop in DM towards $`𝒯`$ –20 are evidence for the clumpy nature of the wind.

Flux density variations in PSR B1259–63 occur on a number of timescales. The diffractive scintillation timescale at GHz frequencies close to periastron is expected to be very short ($`1`$ second) and hence unobservable. The refractive timescale is probably of the order of minutes. One also expects changes in the flux density due to free-free absorption as the pulsar moves behind clumps in the Be star wind. In Figure 3, we show flux variations as a function of time starting on 1997 May 8.3 ($`𝒯`$ –21). Each point is an integration over 2 minutes of data. These variations are broad band (at least across the 300 MHz observing bandwidth) at 1.4 GHz. At higher frequencies, the pulsar is too weak to allow reliable flux densities to be obtained in a few minutes.

This extreme flickering is only observed on $`𝒯`$ –21. Observations made the previous day and on the following three days show only slow changes in flux density over hours and a much lower modulation index, consistent with the refractive scintillation timescale. $`𝒯`$ –21 corresponds exactly to the first enhancement in the flux density seen in the unpulsed source \[Johnston et al. 1999\]. This early increase in the flux density subsequently disappeared and does not fit well with the model of the unpulsed source \[Ball et al. 1999\], although Johnston et al. (1999) have argued that it is likely to be the result of the first interaction between the pulsar wind and the stellar disc coupled with optical depth effects.

## 4 Discussion

At $`𝒯`$ –20, the pulsar is $``$50 stellar radii from the Be star. Assuming a wind density at the stellar surface of $`7\times 10^8`$cm-3 which varies with radius as $`r^2`$, the line-of-sight density to the pulsar should be $``$3$`\times 10^5`$cm-3. However, one cannot expect this wind to be homogeneous and there are likely to be significant density variations because of turbulence.

We can use the equation derived by McClure-Griffiths et al. (1998) to relate the electron density in the scattering material, $`\delta n_e^2`$, to the scattering parameters.

$$\delta n_e^2=4.0\mathrm{\Delta }\nu _{\mathrm{d},\mathrm{MHz}}^{5/6}D_{\mathrm{kpc}}^{5/6}\nu _{\mathrm{GHz}}^{11/3}\frac{l_{o,\mathrm{pc}}^{2/3}}{\mathrm{\Delta }L_{\mathrm{pc}}}\left[\frac{(1+x)^2}{x}\right]^{5/6}\mathrm{cm}^6$$
(1)
Here, $`\mathrm{\Delta }\nu `$ is the scintillation bandwidth, $`D`$ the distance to the pulsar and $`\nu `$ the observing frequency. $`l_o`$ is the outer scale length, $`\mathrm{\Delta }L`$ is the thickness of the scattering screen and $`x`$ is the ratio of the observer-screen distance and the pulsar-screen separation. We simplify this by using $`D=1.5`$ kpc, $`\nu =1.4`$ GHz and assuming $`l_o\mathrm{\Delta }L`$, $`x1`$ to obtain

$$\delta n_e^2=9.8\mathrm{\Delta }\nu _{\mathrm{d},\mathrm{MHz}}^{5/6}x^{5/6}\mathrm{cm}^6$$
(2)
Converting the scintillation bandwidth, $`\mathrm{\Delta }\nu `$ to the scattering time, $`\tau `$, in milliseconds via $`2\pi \tau \mathrm{\Delta }\nu =1`$ yields

$$\delta n_e^2=1.43\times 10^4\tau ^{5/6}x^{5/6}\mathrm{cm}^6$$
(3)
Assuming a distance to the screen, $`r_s`$, of $`8.4\times 10^{12}`$cm, (i.e $`x=5.5\times 10^8`$), a scattering time of 3 ms gives $`\delta n_e^2`$ of $`7\times 10^{11}`$cm-6 and hence $`n_e`$ of $`8.3\times 10^5`$cm-3 assuming a fully modulated wind. This is remarkably close to the expected value obtained from a simple scaling of the electron density in the wind and shows that the scattering is caused by the global properties of the fully modulated wind. The size of the scattering disc, given by $`(2\tau cr_s)^{1/2}`$ is rather small, only $`3.9\times 10^{10}`$ cm.

These values are consistent with the DM measurements. A mean wind density of $`8.3\times 10^5`$cm-3 along a pathlength of roughly $`2r_s`$ yields a $`\mathrm{\Delta }`$DM of 4.5 cm-3pc, close to the observed values. We also note that to double the observed scattering tail requires only a change in $`n_e`$ of $``$30 per cent and that this change affects the DM at the $``$1 cm-3pc level. Such changes are seen on roughly a 2-hr timescale.

From equation 3 above it is evident that $`\tau n_e^{12/5}`$ whereas DM $`n_e`$. Providing the same electrons are responsible for the scattering and the dispersion we should observe that $`\tau \mathrm{DM}^{12/5}`$. Figure 2 shows the expected correlation between dispersion measure and scattering time in that they both rise and fall together. In Figure 4 we show a plot of scattering time versus $`\mathrm{\Delta }`$DM. The scatter in the diagram leads to a large uncertainty in any formal straight line fit. For illustrative purposes, we show a line with slope 12/5 as derived above, which is in satisfactory agreement with the data.

Consider the free-free optical depth along the line of sight, given by

$$\tau _{\mathrm{ff}}=8.2\times 10^2T^{1.35}\nu ^{2.1}_0^dn_e^2\mathrm{dl}$$
(4)
where $`T`$ is the temperature in Kelvin, $`\nu `$ the observing frequency in GHz, $`d`$ the distance through the absorbing material in pc and $`n_e`$ the free electron density within the clump. For an electron temperature in the wind of $`10^4`$ K \[Waters 1986\], the wind parameters described above yield $`\tau _{\mathrm{ff}}0.6`$ at 1.4 GHz, fully consistent with the lower flux densities recorded at these epochs. Flux density variations are seen both on timescales of a few minutes and on a roughly similar timescale to the DM and scattering variations. The longer timescale is probably due to variations in the free-free optical depth, either because of temperature variations or small variations in the electron density. The short timescale variations generally have a low modulation index and are likely to be refractive scintillation. Given a timescale of $``$200 s for the variations and knowledge of the size of the scattering disc ($`3.9\times 10^{10}`$cm), the dominant velocity in the system is then $``$2000 kms-1. Such a velocity is consistent with UV observations of Be star winds in general.

In summary, we have shown that a plausible model for the wind of the Be star can account for the dispersion measure and scattering data before periastron and is also consistent with the variability timescale and the flux density variations. We stress that we are discussing the spherical, high-velocity, low-density wind from the star. The pulsar encounters the circumstellar disc material at $`𝒯`$ –18, the density of the disc is so large that the pulsed flux becomes eclipsed shortly thereafter through a combination of free-free absorption and extreme pulse broadening and remains so until the pulsar re-emerges from the disc at $`𝒯`$ +16. Post-periastron, the pulsar lies between the observer and the Be star. The path length through the stellar wind is much less than pre-periastron and the DM and scattering variations are hence much smaller.

## Acknowledgments

We thank the Parkes staff, especially H. Fagg and M. McColl, for their considerable help throughout the observations. We thank M. Britton for useful discussions. The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by the CSIRO.


format: text

text:
# Hedged Monte-Carlo: low variance derivative pricing with objective probabilities

## 1 Introduction

The Black and Scholes model for options has two truly remarkable properties: (i) one can find a ‘perfect’ hedging strategy that eliminates risk entirely, and (ii) the option price does not depend at all on the average return of the underlying asset . This last property shows that the option price is not simply the (discounted) average of the future pay-off over the objective (or historical) probability distribution, as one would have naively expected. This is even more striking in the case of the Cox-Ross-Rubinstein binomial model  where the pricing measure is completely unrelated to the actual distribution of returns. These two models have deeply influenced financial mathematics, and have lead to the development of a rather abstract and general framework for derivative pricing. Within this framework, the fundamental result is that the requirement of absence of arbitrage opportunities is equivalent to the existence of a ‘risk-neutral probability measure’ (a priori distinct from the objective one), such that the price of a derivative is indeed its (discounted) average pay-off, but where the average is performed over the risk-neutral distribution rather than over the objective distribution . It is thus a common belief that the knowledge of the ‘true’ probability distribution of returns is a useless information to price options. The credence is rather that the relevant risk neutral distribution is somehow ‘guessed’ by the market. A possibility is to reconstruct this risk neutral distribution from quoted option prices, and to use this information to price other instruments consistently (for example exotic options) . In many cases, however, exact analytical expressions for these exotic options are not available and numerical methods must be used. A very versatile method, which allows to price complicated path dependent options, is the Monte-Carlo method, where paths are generated with a weight consistent with the risk-neutral distribution. Efficient numerical techniques have been developed, in particular in , to implement this ‘risk neutral Monte-Carlo’ (rnmc).

However, in most models of stock fluctuations, except for very special cases, risk in option trading cannot be eliminated, and strict arbitrage opportunities do not exist, whatever the price of the option. That risk cannot be eliminated is furthermore the fundamental reason for the very existence of option markets. It would thus be more satisfactory to have a theory of options where the objective stochastic process followed by the underlying asset was used to compute the option price, the hedging strategy, and the residual risk. The latter is clearly important to estimate for risk control purposes. A natural framework for this is the risk minimization approach developed by several authors , where the optimal trading strategy is determined such that the chosen measure of risk is minimized (for example the variance of the wealth balance, although other choices are possible ) . The ‘theoretical’ price is then obtained using a fair game argument. Interestingly, this framework allows one to recover exactly the Black-Scholes results when the objective probabilities are log-normal, and when the continuous time limit is taken (this is shown in details in ). In particular, the average trend indeed completely disappears from the price and hedge.

The aim of this paper is to present a very general Monte-Carlo scheme based on this approach, that we call ‘hedged Monte-Carlo’ (hmc). The method, which has been inspired in part by the least square method (lsm) of Longstaff and Schwartz , shares with it the property that it can price a wide variety of exotic options including those with path-dependent or early exercise features. On top of that, the hmc has at least four major advantages over rnmc:

* The most important one is considerable variance reduction. This is related to the fact that the financial risk arising from the imperfect replication of the option by the hedging strategy is directly related to the variance of the Monte-Carlo simulation. When minimizing the former by choosing the optimal strategy, we automatically reduce the latter. The standard deviation of our results are typically five to ten times smaller than with rnmc, which means that for the same level of precision, the number of trajectories needed in the Monte-Carlo is up to a hundred times smaller.
* The hmc method provides not only a numerical estimate of the price of the derivative, but also of the optimal hedge (which may be different from the Black-Scholes $`\mathrm{\Delta }`$-hedge for non Gaussian statistics) and of the residual risk.
* The method does not rely on the notion of risk neutral measure and can be used with any model of the true dynamics of the underlying (even very complex), in particular those for which the risk neutral measure is unknown and/or not uniquely defined.
* Last, but not least, the hmc method allows one to use purely historical data to price derivatives, short-circuiting the modeling of the underlying asset fluctuations. These fluctuations are known to be of a rather complex statistical nature, with fat-tailed distributions, long-range volatility correlations, negative return-volatility correlations, etc. . Within the hmc method, one can directly use the historical time series of the asset to generate the paths. The fact that a rather small number of paths is needed to reach good accuracy means that the length of the historical time series need not be very large.

## 2 Basic principles of the method

Option pricing always requires to work backwards in time. This is because the option price is exactly known at maturity, where it is equal to the pay-off. As with other schemes, we determine the option price by working step by step for maturity $`t=N\tau `$ to the present time $`t=0`$, the unit of time $`\tau `$ being, for example, one day. The price of the underlying asset at time $`k\tau `$ is denoted as $`x_k`$ and the price of the derivative is $`𝒞_k`$. We assume for simplicity that $`𝒞_k`$ only depends on $`x_k`$ (and of course on $`k`$). However, the method can be generalized to account for a dependence of $`𝒞_k`$ on the volatility, interest rate, etc., or to price multi-dimensional options (such as interest rate derivatives). We therefore also introduce the hedge $`\varphi _k(x_k)`$, which is the number of underlying asset in the portfolio at time $`k`$ when the price is equal to $`x_k`$. Within a quadratic measure of risk, the price and the hedging strategy at time $`k`$ is such that the variance of the wealth change $`\delta W_k`$ between $`k`$ and $`k+1`$ is minimized. More precisely, we define the local ‘risk’ $`_k`$ as:

$$_k=\left(𝒞_{k+1}(x_{k+1})𝒞_k(x_k)+\varphi _k(x_k)[x_kx_{k+1}]\right)^2_o,$$
(1)
where $`\mathrm{}_o`$ means that we average over the objective probability measure (and not the risk-neutral one!). As shown in , the functional minimization of $`_k`$ with respect to both $`𝒞_k(x_k)`$ and $`\varphi _k(x_k)`$ gives equations that allow one determine the price and hedge, provided $`𝒞_{k+1}`$ is known. Note that we have not included interest rate effects in Eq. (1). When the interest rate $`r`$ is non zero, one should consider the following modified equation:

$$_k=\left(e^\rho 𝒞_{k+1}(x_{k+1})𝒞_k(x_k)+\varphi _k(x_k)[x_ke^\rho x_{k+1}]\right)^2_o,$$
(2)
where $`\rho =r\tau `$ is the interest rate over an elementary time step $`\tau `$.

In order to implement this numerically, we parallel the lsm method of Longstaff and Schwartz, developed within a risk-neutral approach . We generate a set of $`N_{\mathrm{MC}}`$ Monte-Carlo trajectories $`x_k^{\mathrm{}}`$, where $`k`$ is the time index and $`\mathrm{}`$ the trajectory index. We decompose the functions $`𝒞_k`$ and $`\varphi _k`$ over a set of $`M`$ appropriate basis functions $`C_a(x)`$ and $`F_a(x)`$:111For numerical purposes, these basis functions have been chosen to be piecewise linear for $`F_a`$ and piecewise quadratic for $`C_a`$, with adaptive breakpoints.

$$𝒞_k(x)=\underset{a=1}{\overset{M}{}}\gamma _a^kC_a(x)\varphi _k(x)=\underset{a=1}{\overset{M}{}}\phi _a^kF_a(x).$$
(3)
In other words, we solve the minimization problem with the variational space spanned by the functions $`C_a(x)`$ and $`F_a(x)`$. This leads to a major simplification since now we have a linear optimization problem in terms of the coefficients $`\gamma _a^k,\phi _a^k`$, for which efficient numerical algorithms are available . These coefficients must be such that:

$$\underset{\mathrm{}=1}{\overset{N_{\mathrm{MC}}}{}}\left(e^\rho 𝒞_{k+1}(x_{k+1}^{\mathrm{}})\underset{a=1}{\overset{M}{}}\gamma _a^kC_a(x_k^{\mathrm{}})+\underset{a=1}{\overset{M}{}}\phi _a^kF_a(x_k^{\mathrm{}})[x_k^{\mathrm{}}e^\rho x_{k+1}^{\mathrm{}}]\right)^2$$
(4)
is minimized. Those $`N`$ minimization problems (one for each $`k=0,\mathrm{},N1`$) are solved working backwards in time with $`𝒞_N(x)`$ the known final pay-off function.

Although in general the optimal strategy is not equal to the Black-Scholes $`\mathrm{\Delta }`$-hedge, the difference between the two is often small, and only leads to a second order increase of the risk . Therefore, one can choose to work within a smaller variational space and impose that:

$$\phi _a^k\gamma _a^kF_a(x)\frac{dC_a(x)}{dx}.$$
(5)
This will lead to exact results only for Gaussian processes, but reduces the computation cost by a factor two.

## 3 Numerical results for the Black-Scholes model

We have first checked our hmc scheme when the paths are realizations of a (discretized) geometric random walk. We have priced an at the money three month European option, on an asset with $`30\%`$ annualized volatility and a drift equal to the risk-free rate which we set to 5 % per annum. The number of time intervals $`N`$ is chosen to be $`20`$. The initial stock and strike price are $`x_0=x_s=100`$, and the corresponding Black-Scholes price is $`𝒞_0^{BS}=6.58`$. The number of basis functions is $`M=8`$. We run $`500`$ simulations containing $`500`$ paths each, for which we extract the average price and standard deviation on the price. An example of the result of linear regression is plotted in Fig. 1. Each data point corresponds to one trajectory of the Monte-Carlo at one instant of time $`k`$, and represents the quantity:

$$e^\rho 𝒞_{k+1}(x_{k+1})+\varphi _k(x_k)[x_ke^\rho x_{k+1}],$$
(6)
as a function of $`x_k`$. The full line represents the result of the least-squared fit, form which we obtain $`𝒞_k(x_k)`$. We show in the inset the corresponding hedge $`\varphi _k`$, that was constrained in this case to be the $`\mathrm{\Delta }`$ hedge.

We obtain the following numerical results. For the rnmc (un-hedged) scheme, we obtain $`𝒞_0^{RN}=6.68`$ with a standard deviation of $`0.44`$. For the hmc, we obtain $`𝒞_0^H=6.55`$ with a standard deviation of $`0.06`$, seven times smaller than with the rnmc. This variance reduction is illustrated in Fig. 2, where we show the histogram of the MC results both for the un-hedged case (full bars) and for the hedged case (dotted bars).

Now we set the drift to $`30\%`$ annual. The Black-Scholes price, obviously, is unchanged. A naive un-hedged Monte-Carlo scheme with the objective probabilities would give a completely wrong price of $`10.72`$, $`60\%`$ higher than the correct price, with a standard deviation of $`0.56`$. On the other hand, the hmc indeed produces the correct price ($`6.52`$) with a standard deviation of $`0.06`$. The rnmc scheme in this case simply amounts to setting ‘by hand’ the drift to risk-free rate, and therefore obviously gives back the above figures.

Therefore, we have checked that in the case of a geometric random walk, the hmc indeed gets rid of the drift and reproduces the usual Black-Scholes results, as it should. This allows us to confidently extend the method to other type of options and other random processes.

## 4 American and other exotic options

The hmc method can be used to reduce the Monte Carlo error for all types of exotic options. We illustrate this point by showing how the method can be extended to price an American put option. In order to implement the early exercise condition, one can simply replace $`𝒞_{k+1}(x_{k+1})`$ in Eq. (4) by $`\mathrm{max}(𝒞_{k+1}(x_{k+1}),x_sx_{k+1})`$, where $`x_s`$ is the strike price. We have chosen a slightly different implementation, where we first find the early exercise point $`x_{k+1}^{}`$ and exercise all options for which $`x_{k+1}<x_{k+1}^{}`$.

We have tested the method on a one year American put option on a stock following a standard log-Brownian process. We follow the choice of parameters made in  in order to compare our results with theirs. The initial price and the strike are set to $`x_0=x_s=40`$, the volatility to $`20\%`$ annual and the risk-free rate and the drift to $`6\%`$. As a benchmark price we use the value $`2.314`$ (quoted in ) computed using a very accurate finite difference method. We computed the price within the hmc using $`N_{\mathrm{MC}}=500`$ paths and $`M=8`$ basis functions. To measure the accuracy of the method, we ran the Monte Carlo 500 times with different random seeds. The average price found was $`2.302`$ with a standard deviation (around the true value $`2.314`$) of $`0.032`$.

We also used the least-square method (lsm) of  with the same parameters ($`N_{\mathrm{MC}}=500`$ paths and $`M=8`$ basis functions). Note that the lsm used in  differs from the hmc in the following way: (i) the hedge is not used in the least-square fit Eq. (4), (ii) on subsequent iterations the option price is kept as the discounted pay-off on that particular path (final pay-off or early exercise value) and (iii) the final option price is given as the average pay-off of all paths. In this framework the least-square fit only serves to find the early exercise point. Therefore, for European options, the lsm is identical to the standard rnmc.

The average price within the lsm is found to be $`2.423`$ with a standard deviation around the true value of $`0.170`$, five times larger than for the above quoted $`0.032`$ for the hmc. These numbers are compatible with those found in  where the error quoted is $`0.01`$, i.e. $`17`$ times smaller but with $`200`$ times more paths and $`2.5`$ times more intermediate points.

Obviously, the same variance reduction would hold for other exotic paths, as those discussed in . We have not investigated in depth the optimal values to be given to the parameters $`M`$ and $`N_{\mathrm{MC}}`$, or the choice of the basis functions that minimize the computation cost for a given accuracy. These are implementation issues that are beyond the scope of the present paper.

## 5 Purely historical option pricing

We now turn to the idea of a purely historical hmc pricing scheme. We price a one month (21 business days) option on Microsoft Corp., hedged daily, with zero interest rates. We used 2000 paths of length 21 days, obtained form the time series of Microsoft during the period May 1992 to May 2000. The initial price is always renormalized to $`100`$. We use a set of $`M=10`$ basis functions, and keep with the simple $`\mathrm{\Delta }`$-hedge. From our numerically determined option prices, we extract an implied Black-Scholes volatility by inverting the Black-Scholes formula and plot it as a function of the strike, in order to construct an implied volatility smile. The result is shown in Fig. 3, together with the error bars, obtained from the residual risk of the hedged options222The error on the price is given by the residual risk divided by the square-root of the number of independant paths. Since we over-sample the paths, we use $`2000\times 2/21190`$ as the effective number of paths.. The residual risk itself, divided by the call or the put option price (respectively for out of the money and in the money call options), is given in the inset. We find that the residual risk is around $`42\%`$ of the option premium at the money, and rapidly reaches $`100\%`$ when one goes out of the money. These risk numbers are comparable to those obtained on other options of similar maturity (see ), and are much larger than the residual risk that one would get from discrete time hedging effects in a Black-Scholes world.

The smile that we obtain has a shape quite typical of those observed on option markets. However, it should be emphasized that we have neglected the possible dependence of the option price on the local value of the volatility. This could be taken into account by letting the function $`𝒞_k`$ depend not only on $`x_k`$ but also on the value of some filtered past volatility $`\sigma _k`$.

## 6 Conclusion and prospects

We have presented what we believe to be a very useful Monte-Carlo scheme, which closely follows the actual history of a trader hedged portfolio. The inclusion of the optimal hedging strategy allows one to reduce the financial risk associated to option trading, and for the very same reason the variance of our hmc scheme as compared to the previously discussed rnmc schemes. The explicit accounting of the hedging cost naturally converts the objective probability into the ‘risk-neutral’ one. This allows a consistent use of purely historical time series to price derivatives and obtain their residual risk. We believe that there are many extensions and applications of the scheme, for example to price interest rate derivatives with faithful historical models (such as the one proposed in ), and market hedging instruments. With some modifications and extra numerical cost, the method presented here could be used to deal with transaction costs, or with non quadratic risk measures (VaR hedging).

### Acknowledgments

We thank J.-P. Aguilar, A. Matacz and M. Bazant for interesting discussions.


text:
# LYMAN CONTINUUM EMISSION FROM GALAXIES AT 𝐳≃3.4 Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and NASA, and was made possible by the generous financial support of the W.M. Keck Foundation.

## 1 INTRODUCTION

The question of the origin of the UV radiation field that re-ionized the universe at $`z>6`$ and maintained the ionization of the intergalactic medium at lower redshifts has been addressed numerous times since it was first considered 20 years ago (Sargent et al. 1980). The relative contributions of QSOs and massive stars in galaxies to the ionizing photon budget at high redshifts has important implications for the spectral energy distribution of the ionizing background (Miralda-Escude & Ostriker 1990; Steidel & Sargent 1989) which in turn affects the equation of state of the diffuse intergalactic medium. Early consideration of the evolution of the diffuse intergalactic medium suggested that the observed QSOs would have an increasingly difficult time producing enough ionizing photons to explain the absence of an observed Gunn-Peterson (1965) effect at redshifts $`z\stackrel{>}{}3`$ (e.g., Donahue & Shull 1987, Shapiro & Giroux 1987), suggesting the need for additional sources of ionizing photons at high redshift.

The contribution of star-forming galaxies to the UV background was first explicitly considered by Bechtold et al. (1987), but at the time there was little relevant data. Progress in very deep galaxy surveys (e.g., Songaila, Cowie, & Lilly 1990, Steidel & Hamilton 1993) suggested that high redshift ($`z\stackrel{>}{}3`$) star-forming galaxies might well exist, but many uncertainties in their contribution to the ionizing background remained due to the lack of spectroscopic redshift identifications and the uncertain opacity of galaxies to their own Lyman continuum radiation.

The advent of the Keck telescopes in the mid-1990s made possible spectroscopy of large numbers of star forming galaxies at redshifts $`z3`$, where the strongest constraints on the properties of the intergalactic medium were available from QSO absorption line surveys, where the ability of known QSOs to maintain the ionization of the IGM became uncertain, and where the rest-frame Lyman limit of hydrogen could be observed from the ground. This population of UV-bright star forming galaxies (Steidel et al. 1996), generally referred to as Lyman break galaxies (LBGs), is selected photometrically by exploiting the continuum discontinuity near the rest frame Lyman limit due to the intrinsic spectral energy distribution of hot stars, Lyman continuum opacity from H I layers in the galaxy itself, and, most model-independently, the opacity of the intergalactic medium due to intervening H I (e.g., Madau 1995). While this large continuum break (only a portion of which is expected to be intrinsic to the galaxy) makes the discovery of LBGs extremely straightforward and statistically complete, it also makes the measurement of leaking Lyman continuum radiation exceedingly difficult, even with 10-m class telescopes.

The problem of intergalactic HI opacity can be minimized by observing the Lyman continuum in galaxies at smaller redshifts, but this requires observations in the vacuum UV from space. The most sensitive direct measurement to date, by Leitherer et al. (1995), used Hopkins Ultraviolet Telescope far-UV spectra of 4 nearby starburst galaxies to place constraints on the fraction of Lyman continuum photons escaping the galaxies, $`f_{\mathrm{esc}}<3`$% (but see Hurwitz, Jelinsky, & Dixon 1997). Similarly small upper limits for relatively nearby galaxies have been set by several groups using less direct methods. Nevertheless, it is not clear how relevant these limits are to the situation at high redshift, where galaxies which simultaneously possess very large UV luminosities and low UV extinction are much more common (cf. Adelberger & Steidel 2000).

This paper presents what we believe to be the first direct detection of Lyman continuum flux from galaxies, made possible by a high S/N composite spectrum formed from a relatively small subset of $`z3`$ LBGs that are at high enough redshift to bring the rest-frame Lyman limit into a region of good spectroscopic sensitivity for current instrumentation. This measurement, which we currently regard as plausible, but tentative, has significant implications for the nature of the UV radiation field at high redshifts. In the context of evaluating the significance, we also present some preliminary results on the contribution of faint QSOs to the $`z3`$ UV background.

## 2 DATA

The spectra selected for this study were chosen from among the 875 spectra of $`z3`$ LBGs we have obtained with the Keck telescopes and Low Resolution Imaging Spectrometer (Oke et al. 1995) between 1995 October and 1999 November. All of the spectra were obtained using the same instrumental configuration, a 300 lines/mm grating blazed at 5000 Å and a slit mask with 1$`^{\prime \prime }.`$4 slits, providing spectral resolution of $`12.5`$ Å in the observed frame. The typical total exposure time was 7200 seconds. In its current configuration, the total throughput of LRIS with this grating is $``$20% at 4000 Å and $`10`$% at 3800Å; in order to include only spectra with adequate sensitivity at the rest-frame Lyman limit, we have confined ourselves to the small subset of our spectra which have $`z3.300`$ and $`25.0`$ (81 objects). Some of these spectra did not reach the Lyman limit on the blue side, due to the geometry of the slit masks, and some were obtained under poor conditions, were contaminated by light from neighboring galaxies, or were affected by other obvious problems. Two dimensional sky-subtracted spectrograms of each object were inspected for such problems, and only the cleanest looking slits with no obvious sky subtraction problems were retained. The final sample consists of 29 spectra with $`z=3.40\pm 0.09`$; the average properties of the LBGs are summarized in Table 1. The important observed-frame spectral region for evaluating the Lyman continuum is 3880-4020 Å at the average redshift.

The composite spectrum of the LBGs was constructed by shifting each extracted, one-dimensional, flux calibrated spectrum into the rest frame and averaging after scaling to a common median in the rest wavelength range 1250–1500 Å. A simple rejection scheme was used in which equal numbers of positive and negative outliers were excluded from the average at each dispersion point (primarily in order to exclude sky subtraction residuals near bright night sky lines). The overall nature of the co-added spectrum was found to be insensitive to the number of points excluded from the average at each dispersion point. The final spectrum, shown in Figure 1, results from the rejection of 2 positive and 2 negative outliers from each dispersion point. The composite spectrum does not depend strongly on the adopted scheme for combining the individual spectra.

The composite spectrum clearly shows the net effects of the Lyman $`\alpha `$ forest shortward of the Lyman $`\alpha `$ emission line, reducing the flux by a factor of $`1.8`$ relative to the flux in the $`12501500`$ Å range. Also evident is the general diminution of the flux shortward of Lyman $`\beta `$ due to increasing numbers of intervening and intrinsic absorption features. Many stellar and interstellar absorption features can be recognized in the composite spectrum; the most prominent are labeled in Figure 1. Interestingly, averaging the Lyman $`\alpha `$ forest absorption over many sight-lines brings out stellar and interstellar lines which are normally masked by the forest in any single spectrum. Stellar C III $`\lambda 1176`$ and He II $`\lambda 1084`$ (cf. Walborn & Bohlin 1996; Fullerton et al. 2000), and interstellar Lyman $`\beta `$, $`\gamma `$, and $`\delta `$ are the clearest examples. The Lyman $`\alpha `$ equivalent width in the composite spectrum in Figure 1 is 19.7Å , placing it in the top 20-25% of LBGs in terms of Lyman $`\alpha `$ line strength (cf. Steidel et al. 2000). This large Lyman $`\alpha `$ equivalent width results from a selection bias in our Lyman-break sample towards bluer objects at higher redshifts; at $`z3.4`$ galaxies will satisfy our color-selection criteria only if their intrinsic broad-band colors lie in the bluest quartile of observed LBG UV colors, and a strong contribution from Lyman-$`\alpha `$ emission helps galaxies achieve the required blue broad-band color. See Steidel et al. (1999) for a more complete discussion. The impact of this selection effect on our present results is considered in §4 below.

The bottom panel of Figure 1 shows residual flux in the Lyman continuum region. We have chosen to evaluate the residual flux over the rest wavelength range 880-910 Å, for two main reasons: first, the incidence of intervening Lyman limit absorption systems (Sargent, Steidel, & Boksenberg 1989; Stengler-Larrea et al. 1995) at $`z3.4`$ means that broad-band evaluation of the decrement across the Lyman limit is not likely to provide a good estimate of the emergent spectrum from the galaxies—the mean free path of a Lyman continuum photon at $`z3`$ is only $`\mathrm{\Delta }z0.17`$ (Haardt & Madau 1996), corresponding to $`35`$ Å in the rest frame at $`z3.4`$. Second, the 880-910 Å range includes observed spectral regions that are still within the range for which LRIS provides adequate throughput for the entire range of redshifts included in the composite spectrum. While the significance per resolution element of the residual flux in the Lyman continuum is low, the net flux observed over the rest-frame 30 Å bandpass is formally significant at the $`4.8\sigma `$ level (see Table 2).

An obvious concern in gauging the reality of the residual flux, which corresponds to an average observed flux density of $`m_{\mathrm{AB}}=27.4`$ in the narrow wavelength range, is the extent to which systematic problems with sky subtraction might mimic residual Lyman continuum flux. We experimented extensively with the data, using a number of different tests. First, we confirmed that the signal is not contributed by a few outlier spectra by combining the 1-D spectra using different subsets with different outlier rejection algorithms. We checked for overall systematics in the 2-D sky subtraction by combining the 2-D background-subtracted slitlets after rebinning in the wavelength direction to a common rest wavelength interval, registering on the spatial position of the LBG, and then experimenting with various methods of residual background subtraction (i.e., second-pass subtraction of the background that might eliminate systematics that were too subtle to notice in any single slitlet). We then evaluated the residual flux in the appropriate wavelength range both on and off the spatial positions of the detected continuum longward of the Lyman limit; the results are consistent with the residual flux measured from the stacked 1-D spectra, although all of the additional image processing necessary to combine the spectra in 2-D could itself be a source of systematics errors. Probably the most direct test of sky subtraction systematics is to examine spectra, taken with the same instrumental configuration (all observed on slit masks as part of the general LBG survey), of similar objects that are known to have optically thick Lyman limits at the same observed wavelengths as the region of interest for the LBGs. In our spectroscopic LBG survey, we have observed 3 faint $`=22.0`$,22.2, 23.6) QSOs which have optically thick Lyman limit systems (confirmed subsequently with long-slit observations at higher dispersion) in the redshift range $`3.29z_{LLS}3.43`$. These spectra were shifted into the rest frame of the intervening Lyman limit system, scaled, and averaged to produce the composite spectrum shown in Figure 2. Figure 2 illustrates that optically thick Lyman limit systems result in a measured flux in the  Å interval consistent with zero, making less likely the possibility that scattered light in the instrument is responsible for the residual flux shortward of 912 Å observed in the composite LBG spectrum.

In summary, none of the tests that we have performed gives us reason to believe that the residual flux shortward of 912 Å in Figure 1 is an artifact of systematics. Nevertheless, as discussed below, it would be prudent to view the detection with some skepticism until confirmed by higher quality spectra of individual objects. This is very difficult with existing instrumentation, but we expect that LRIS-B, the blue channel for LRIS optimized for the 3100-5500 Å wavelength range and soon to be commissioned on the Keck I telescope, will be ideally suited to making improved measurements on the much larger numbers of bright LBGs in our sample at $`z3.0`$.

## 3 THE EMERGENT FAR-UV SPECTRUM OF LBGs

One of the complications to the interpretation of the flux measurement above is that the opacity in the galaxies’ rest-frame Lyman continuum (in fact, for all spectral regions shortward of Lyman $`\alpha `$ emission) has a significant contribution from intervening material not directly associated with the galaxies. This contribution must be estimated accurately in order to infer the emergent spectrum of the star-forming galaxies. A common way to accomplish this is with models of intergalactic opacity based on column density distributions of Lyman series lines and the incidence of higher column density systems that contribute significant Lyman continuum opacity (e.g., Lyman limit systems). The distribution functions, obtained from QSO absorption line surveys, are used to estimate the opacity of the IGM as a function of observed wavelength for objects at a given redshift (e.g., Madau 1995; Bershady, Charlton, & Geoffroy 1999). For the present purposes we have opted to use a more direct empirical method that bypasses the need to model the opacity via fits to observed distribution functions. We have formed a composite QSO spectrum by combining spectra obtained for surveys of intervening Lyman limit systems (see Sargent, Steidel, & Boksenberg 1989, Stengler-Larrea et al. 1995) of QSOs in the same range of redshift as the galaxies considered above. The spectra were chosen from among more than 100 QSOs in the redshift range $`2.75z4.10`$ all observed using the Palomar 200-inch telescope and Double Spectrograph; the spectral resolution in the $`31504700`$ Å range was 4Å. The composite spectrum, which is effectively an average over 15 lines of sight for QSOs with a mean redshift $`z_{\mathrm{QSO}}=3.47\pm 0.14`$, has been scaled to the composite LBG spectrum, rebinned to a similar resolution, and over-plotted on Figure 1. Note the similarity of the QSO and galaxy SED over the rest wavelength range 1050-1700 Å, and that the Lyman $`\alpha `$ forest decrements in the rest wavelength range 1050-1170 are equivalent for the QSO and galaxy spectra (the measured value of the $`D_A`$ parameter \[Oke & Korycansky 1982\] is $`0.46`$ for both).

The observed flux densities in the arbitrary units of Figure 1 are summarized in Table 2 for both the QSO and LBG spectra. To estimate the change in net opacity due to intergalactic material between 1100 Å (where it is obvious in both the galaxy and QSO spectra and is dominated by Lyman $`\alpha `$ line opacity from the forest) and rest-frame 900 Å, we adopt the intrinsic QSO spectral energy distribution shortward of 1050 Å suggested by Madau, Haardt, & Rees (1999), $`f_\nu \nu ^{1.8}`$. The observed change in the QSO continuum level between the  Å interval and the  Å interval is a factor of $`3.7\pm 0.2`$; the intrinsic component of this flux decrement is estimated from the power law assumption to be a factor of 1.5, so that the inferred flux decrement at 900 Å relative to that at 1100 Å due to intervening material is an additional factor of $`2.5`$. Using this factor to correct the observed ratio $`f[1100]/f[900]=9.4\pm 2.0`$ for the galaxy spectrum, and accounting for the Lyman $`\alpha `$ line blanketing in the  Å interval and a small difference in unabsorbed flux between 1500 and 1100 Å intrinsic to the galaxy SED (based on the mean color of the composite spectrum after correction of the observed $`G`$ band for Lyman $`\alpha `$ line blanketing), one obtains an intrinsic flux density ratio for the LBGs of $`f[1500]/f[900]=4.6\pm 1.0`$. It is important to note that the observed ratio is $`f[1500]/f[900]17.7`$, or a contrast of 3.1 magnitudes on the AB scale.

Because of the effects of intergalactic absorption on the observations, a positive measure of Lyman continuum flux becomes very rapidly more difficult with increasing redshift; the observation is probably impossible using present observational facilities by redshift $`z4.5`$, despite higher system throughputs for ground-based spectroscopy—the intervening opacity is approximately 3 times higher, the sky background is about 3 times brighter, and typical galaxies are two times fainter, than at $`z3`$.

## 4 IMPLICATIONS

The flux density ratio $`f[1500]/f[900]=4.6\pm 1.0`$ obtained above is unexpectedly large for star forming galaxies; it is similar to typical models of the UV SEDs of star forming galaxies (e.g., Bruzual & Charlot 1996) without any Lyman continuum self–absorption from H I in the galaxy. However, it must be realized that at present there are no direct measurements of the intrinsic amplitude of the stellar Lyman break and we have to rely exclusively on models. The model predictions depend sensitively on the initial mass function and stellar lifetimes at the very high mass end, and on the age of the star formation episode. The amplitude of the stellar Lyman limit for most spectral synthesis models (e.g., Bruzual & Charlot 1996; Leitherer et al. 1999) ranges from a factor of $`3`$ to $`5.5`$ over a plausible range of ages and initial mass functions. More recent models (e.g., Schaerer 1998), which include non-LTE effects, generally produce somewhat smaller break amplitudes than the LTE models that are used by the popular spectral synthesis programs, as small as a factor of $`2.5`$ for constant star formation over periods $`>10^7`$ years (Smith & Norris, private communication) . Adopting for the moment the assumption that the stellar energy distribution has an intrinsic Lyman discontinuity of a factor of 3, our result for $`f[1500]/f[900]`$ implies that the fraction of UV ionizing photons escaping the galaxy is $`f_{\mathrm{esc}}\stackrel{>}{}0.5`$.

Our definition of $`f_{\mathrm{esc}}`$ differs from a definition sometimes encountered elsewhere in the literature (e.g. Leitherer et al. 1995). By $`f_{\mathrm{esc}}`$ we mean the fraction of emitted 900Å photons that escapes the galaxy without being absorbed by interstellar material divided by the fraction of 1500Å photons that escapes. This is the quantity required to calculate the ionizing radiation emitted by Lyman-break galaxies from their well constrained 1500Åluminosity density. An alternate definition of $`f_{\mathrm{esc}}`$ is simply the fraction of emitted 900Å photons that escapes (e.g. Leitherer et al. 1995); in this definition the 1500Å normalization is omitted. Because perhaps only about 15–20% of 1500Å photons escape from typical LBGs without being absorbed by dust (e.g. Pettini et al. 1998, Adelberger & Steidel 2000), values of $`f_{\mathrm{esc}}`$ calculated under these two definitions will differ significantly.

The fraction of escaping Lyman continuum photons is likely to be highly variable from galaxy to galaxy, depending (at least) on the morphological details of the galaxy-scale gaseous outflows that appear to be a ubiquitous signature of both high redshift and low-redshift starburst galaxies (see Heckman 2000). As an example, consider the best far-UV spectrum of a galaxy (at any redshift) obtained to date, the gravitationally-lensed LBG at $`z=2.72`$, MS1512-cB58 (Pettini et al. 2000). In this case it would be surprising to find significant leakage of Lyman continuum photons, given the extremely large covering fraction of optically thick (in the Lyman-continuum) gas implied by the the observed damped interstellar Lyman $`\alpha `$ line and several lines of low-ionization metallic species that are black at line center. On the other hand, the interstellar lines in the spectrum of the composite $`z=3.4`$ LBG presented above are less than half the strength of those in cB58; it is possible that this may be due to a smaller outflow covering fraction along our line of sight (or a smaller average velocity spread—it is difficult to tell the difference in spectra of 12.5Å resolution), in which case a higher fraction of leaking Lyman continuum photons would be expected. As mentioned in §2, the galaxies chosen for the composite spectrum are drawn from the bluest quartile of the observed LBG population. While on average they have observed UV luminosities similar to the de-magnified luminosity of cB58, they are a priori more likely to be younger and/or less dusty. Thus, it may be that, because of the selection effect for the high redshift tail of the $`z3`$ LBG distribution to be drawn from the bluest objects, they may be those “most likely to succeed” in getting Lyman continuum photons out into the intergalactic medium.

Since there is now a published luminosity function for $`z3`$ LBGs to apparent magnitudes of $`R27`$, 2.5 magnitudes fainter than $`L^{}`$ (Steidel et al. 1999), it is straightforward to calculate the contribution of LBGs to the general metagalactic radiation field at $`z3`$ under the assumption that the composite spectrum is representative of the population as a whole. We use the measured ratio $`L(1500)/L(900)=4.6\pm 1.0`$ to convert the published LBG luminosity function (evaluated at a rest frame wavelength close to 1500 Å) into a distribution of Lyman continuum luminosities. This calculation is identical to that already carried out by Madau (2000), substituting a slightly different constant of proportionality. Integrating the Schechter (1976) function over all observed luminosities (down to 0.1$`L^{}`$) we obtain a co-moving emissivity at 1 Ryd of $`1.2\pm 0.3\times 10^{26}h`$ ergs s<sup>-1</sup> Hz<sup>-1</sup> Mpc<sup>-3</sup> (for an Einstein de-Sitter cosmology) at $`z3`$, exceeding (according to Madau, Haardt, & Rees 1999) the contribution from QSOs by a factor of about 5. Ignoring any contribution to the radiation field from “re-processed” recombination and 2-photon radiation (cf. Haardt & Madau 1996), and assuming an effective absorption distance of $`\mathrm{\Delta }z0.17`$ for a Lyman continuum photon at $`z3`$, the implied contribution of the LBGs to the metagalactic radiation field intensity is $`J_\nu (912)1.2\pm 0.3\times 10^{21}`$ ergs s<sup>-1</sup> cm<sup>-2</sup> Hz<sup>-1</sup> sr<sup>-1</sup>. This is very close to the value typically obtained from the QSO “proximity effect” (Bajtlik, Duncan, & Ostriker 1989, Scott et al. 2000, and references therein). Estimates of $`J_\nu (912)`$ from the proximity effect are subject to a variety of systematic uncertainties (see discussion in Rauch et al. 1997); independent upper limits have been placed on the value of $`J_\nu (912)`$ based on the lack of detection of fluorescent Lyman $`\alpha `$ emission from Lyman limit absorption clouds (Bunker et al. 1999), $`J_\nu (912)<2\times 10^{21}`$, but the validity of this limit is somewhat model-dependent. Given the inherent uncertainties in both our crude calculation and in the measurements of $`J_\nu (912)`$, we do not believe that there is currently reason to be concerned about over-producing the UV background.

There has been extensive speculation in the literature as to whether QSOs are sufficient to provide the implied value of $`J_\nu (912)`$ at high redshift (e.g., Donahue & Shull 1987, Shapiro & Giroux 1987, Haardt & Madau 1996, Scott et al. 2000, Meiksin & Madau 1993, Madau, Haardt, & Rees 1999). Most of the recent studies conclude that QSOs are marginally able to produce the proximity effect background at $`z3`$, but begin to have a serious problem by $`z4`$ due to the observed decreasing space density of bright QSOs (assumed to apply to faint QSOs as well) that has been observed by several groups (Kennefick et al. 1995, Schmidt et al. 1995, SDSS collaboration). At $`z3`$, we have perhaps the best data for testing whether the QSO luminosity function at apparent magnitudes fainter the $`R21`$ is consistent with the extrapolated luminosity functions used in computing the QSO contribution to the background <sup>1</sup><sup>1</sup>1Unlike the QSO contribution, the LBG contribution to the metagalactic UV radiation field is based on a luminosity function where the full range of luminosities has been observed.. For the LBG survey at $`z3`$ ($`0.3`$ square degrees) our observational selection of Lyman break objects should be at least as complete for QSOs as it is for galaxies, since no attempt has been made to remove stellar objects from the sample, the SEDs of the galaxies and QSOs are quite similar, and the spectroscopic identification of objects with broad emission lines is more straightforward than for typical LBG spectra. Adopting the QSO luminosity function of Pei (1995), we should have found a total of $`12`$ spectroscopically identified QSOs in the redshift range $`2.5\stackrel{<}{}z\stackrel{<}{}3.5`$ over the apparent magnitude range $`20.525.5`$, taking into account the spectroscopic completeness and effective survey volumes as a function of apparent magnitude and redshift as described by Steidel et al. (1999). Our sample of 875 spectroscopically identified objects contains a total of 8 QSOs with $`2.5z3.5`$, ranging in apparent magnitude from $`=20.624.8`$. While a more careful analysis of the observed QSO density is beyond the scope of this letter, the point is that the QSO contribution to the metagalactic background is certainly no higher than has been estimated in recent analyses (cf. Haardt & Madau 1996). Steeper faint end luminosity function slopes for QSOs at $`z3`$ are probably excluded by our survey data (cf. Haiman, Madau, & Loeb 1999).

## 5 SUMMARY AND DISCUSSION

We have presented observational evidence for significant Lyman continuum emission over the rest-frame wavelength range 880-910 Å in the composite spectrum of Lyman break galaxies at $`z3.4`$. After accounting for the effects of absorption due to intervening material in the spectral range of interest, we find an emergent flux density ratio $`f(1500)/f(900)=4.6\pm 1.0`$. This ratio has been used to convert the observed far-UV luminosity function of LBGs into the luminosity density of Lyman continuum photons at $`z3`$. Under the assumption that this spectral energy distribution is characteristic of the whole LBG population, then LBGs provide a factor of $``$ 5 times more hydrogen ionizing photons than the estimated contribution from QSOs at the same redshifts. We have emphasized that the LBGs comprising the composite spectrum are drawn from the bluest quartile of intrinsic far-UV colors, so that they may be more likely than typical LBGs to exhibit significant escaping Lyman continuum emission. However, adopting the steep $`\alpha =1.6`$ far-UV luminosity function observed for LBGs (Steidel et al. 1999), intrinsically faint objects contribute a substantial fraction of the UV luminosity density. It may be that these faint objects are bluer on average than the $`L^{}`$ objects which comprise most of the LBG spectroscopic samples (cf. Meurer, Heckman, & Calzetti 1999), in which case the composite spectrum may be a good representation of the objects which dominate the UV luminosity density. Even if the distribution of rest–frame colors is independent of UV luminosity (in which case the composite spectrum would be representative of only $`25`$% of the LBGs), the implication is that LBGs contribute at least as many ionizing photons as QSOs at $`z3`$. We have presented preliminary results on the QSO density at faint magnitudes which emphasize that, if anything, the QSO contribution has been slightly over-estimated.

A metagalactic radiation field with a contribution from young galaxies that is at least as large (at $`z3`$) as that from QSOs makes it much easier to understand how the intergalactic medium can remain highly ionized at redshifts $`z\stackrel{>}{}4`$, where the QSO contribution to the radiation field decreases significantly (cf. Madau, Haardt, & Rees 1999). While the UV luminosity function of LBGs is relatively poorly constrained at $`z4`$, and even more poorly constrained at $`z\stackrel{>}{}5`$, there are indications that the co-moving density of star forming galaxies changes much more slowly than that of the QSOs at these redshifts (cf. Steidel et al. 1999), making it quite plausible that the universe is re-ionized by stars rather than AGN, and that stars play at least as important a role as QSOs in maintaining the high ionization level at later epochs.

While there is much to discuss about the possible implications of a metagalactic ionizing radiation field dominated by massive stars at $`z3`$, we suggest that the results above be treated as preliminary until high quality observations of individual galaxies exhibiting clear evidence for Lyman continuum photon leakage become available. In particular, it will be very interesting to explore the variation of Lyman continuum leakage as a function of other galaxy properties (e.g, UV color, strength of stellar and interstellar absorption lines, UV or bolometric luminosity, etc.) as these are all factors of relevance for understanding the physics of rapidly star-forming galaxies at high redshifts, and their importance in affecting the hydrodynamics and chemical enrichment of the intergalactic medium. The limitation of the current data is largely instrumental: most spectrographs are blind in the near-UV. A new generation of near-UV optimized spectrographs on 8-10m telescopes should allow wholesale observations of the Lyman continuum region among the large existing samples of UV-bright galaxies at $`z3`$.

We would like to thank our collaborators in the LBG survey project, Mark Dickinson, Mauro Giavalisco, and Alice Shapley, who played important roles in making this work possible. Informative discussions with Raul Jimenez, Linda Smith, Danny Lennon, and Stephane Charlot are gratefully acknowledged. CCS and KLA have been supported by grant AST95-96229 from the U.S. National Science Foundation and by the David and Lucile Packard Foundation.