Abstract
We present a new Eclipse Mapping method for reconstructing the properties of accretion discs in eclipsing catatclysmic variables. Physical Parameter Eclipse Mapping reconstructs spatially resolved maps of physical parameters such as temperatures and surface densities by fitting observed eclipse light curves at several wavelengths with the aid of a model describing the dependence of the disc spectrum on the physical parameters. Physical Parameter Eclipse Mapping has clear advantages over a classical approach for many applications. Most importantly, all of the available spectrophotometric data can be fitted simultaneously. The only drawback is the assumption of a physical model prior to the reconstructions. However, the method is very flexible in the choice of the spectral model, since the method itself does not depend on the choice. Furthermore, physical considerations, observations and experience will help determine the model to be used.
The analysis with synthetic data shows that reliable parameter distributions can be obtained if the reconstruction behaviour of the model used is well understood. This means that the reconstructions are only as good as the model allows. In ambiguous parameter regions the method can (by definition) not give us the correct answer. However, in ambiguous-free zones, we show that the method indeed gives us the correct answer.
1 Introduction
Accretion discs are common objects in the Universe. They appear in systems as different as pre-main-sequence stars and active galactic nuclei, indicating a variation in the magnitude of the relevant properties by many orders of magnitude, and a very general mechanism of transforming energy into radiation. Unfortunately, the underlying viscous mechanisms which power accretion discs are still not fully understood.
The ideal laboratories for investigating the physical processes in such discs are eclipsing cataclysmic variables, close interacting binaries in which a main-sequence star (hereafter the ‘secondary’) loses mass to a white dwarf (the ‘primary’) through an accretion disc. At certain phases, the secondary partly occults the disc, creating a characteristic light curve containing information about the spatial structure of the disc. The classical Eclipse Mapping method invented by Horne (1985) takes advantage of this fact, and attempts to reconstruct the intensities in the disc at a particular wavelength. The ambiguities inherent in transforming one-dimensional light curve data into a two-dimensional intensity map can be avoided by using, e.g., a maximum entropy algorithm (hereafter ‘MEM’) such as that coded by Skilling & Bryan (1984).
Various authors have produced intensity maps from multicolour high-speed eclipse photometry, using either a version of Horne's original algorithm (e.g. Horne & Cook 1985; Rutten, van Paradijs & Tinbergen 1992a; Wood, Horne & Vennes 1992; Bobinger et al. 1997) or somewhat different algorithms (e.g. Baptista & Steiner 1993; Spruit 1994). However, the principles in all these analyses are generally the same: the intensity distribution is reconstructed with only three assumptions: (a) the geometry of the secondary is known; (b) the geometry of the accretion disc is known; and (c) the emission of the disc is independent of the viewing angle and time. The first assumption is wellfulfilled in semicontact binaries like CVs. The disc is either assumed to be restricted to the orbital plane ( flat disc; Horne & Cook 1985; Baptista & Steiner 1991) or to have a finite opening angle (Robinson et al. 1995; Bobinger et al. 1997). The latter is physically motivated, but does not necessarily improve the resulting distributions dramatically (Rutten 1998). The independence of the disc emission on time is critical, because photometric ‘flickering’ in the disc occurs on shorter time-scales than the duration of the eclipse, making it necessary to average a number of light curves. However, the emission can usually be considered independent of the viewing angle. Only if we use light curves in narrow passbands in the emission lines, shearing effects cannot be neglected (Horne 1995) and the emission from the disc depends on the position of the observer (i.e., for a given system on the phase of the orbit).
Both opaque and transparent discs are predicted for cataclysmic variables on the basis of disc instability theories, which account for the outbursts of dwarf novae by using the ionization of hydrogen in the disc to modulate the viscosity and hence the mass transfer rate M through the disc. In high-M discs, in nova-like variables and in dwarf novae during outbursts, we should find a hot (T>;104 K), optically thick, high-viscosity disc. In contrast, the low-viscosity disc in dwarf novae during quiescence should be cool (Teff∼3000 K) or optically thin. The transition downward to the cool state occurs at a critical surface density Σ that increases linearly with radius R. This prediction may—in principle—be tested if we can determine the radial profiles of T and Σ for real discs by using eclipse observations.
The classical Eclipse Mapping method (e.g. Rutten et al. 1992b; Wood et al. 1992; Baptista, Steiner & Horne 1996) showed that the disc are largly circular (however, minor deviations would not be reconstructed due to the assumption of partial azimuthal symmetry), and disc radii can be derived (e.g., of decline from outburst; Bobinger et al. 1997), especially if one applies the method to infrared (IR) data (the outer cooler parts of the disc are more sensitive for the IR wavelengths (e.g. Froning et al. 1999). Usually, authors also derive brightness temperature distributions and compare them to theoretical effective temperature distributions. To allow for optically thin emission, Wood et al. (1992) applied an optically thin hydrogen model model to the derived U,B,V and R maps of the dwarf nova HT Cas.
However, the major weaknesses of classical Eclipse Mapping stem from the fact that individual, independent intensity maps for each light curve are produced. Every single reconstruction is individually smeared, leading to a deterioration of any derived spatially resolved physical parameter. The reason for performing the reconstructions is to extract physical information about the spatial and spectral properties of the disc. The goal of Physical Parameter Eclipse Mapping (hereafter ‘PPEM’) is to infer the physical conditions present in real accretion discs by using all of the available observations simultaneously in the fitting procedure.
In PPEM, the physical model has to be chosen before a fit can be made. In contrast, classical Eclipse Mapping does not require any such assumption. However, in order to achieve the same results, a physical model often has to be compared with the classically reconstructed intensity distributions. Thus, when making physical comparisons between observations and models, there are no advantages to using the classical method, especially since the method does not depend on the physical model used and one can therefore be very flexible in the choice of the model. A clear exception to this rule applies to situations where the use of a theoretical model is undesirable—e.g., when searching for previously unknown or unexpected features in classical maps with high spectral resolution (e.g. Baptista et al. 1998).
2 The Method
Our approach is very similar to the classical Eclipse Mapping algorithm described in Horne (1985) or Baptista & Steiner (1993). The main formal difference is the replacement of fitted two-dimensional intensity distributions by distributions in physical parameters. This demands the explicit choice of a model describing the intensity from the disc as a function of the physical parameters to be mapped. The model can be as simple as a blackbody spectrum, or, in principle, as complicated as theoretical non-LTE stellar atmosphere models including UV illumination from the central disc.
The choice of the model is necessarily based on prior knowledge of the physical properties of the disc (e.g., whether it is optically thick or thin), the quality and wavelength coverage of the data, and our often poor understanding of the physical processes in accretion discs (e.g., the use of a viscosity prescription in atmosphere models). Since the use of a model to fit observational data is often dangerous—the agreement between model predictions and data is no garantee that the model is correct—the new method requires one not only to start with the simplest conceivable model and smallest number of fitted parameters, but to understand the systematic effects introduced by the use of a certain model in order to estimate the reliability of the reconstructed parameter values. Thus we have purposefully chosen the term ‘Physical Parameter Mapping’ rather than ‘Physical Quantity Mapping’ to express the inherent danger of parameter fitting.
In χ2 fitting of parametrized models, χ2 is minimized by adjusting the model parameters, and the number of parameters determines the complexity of the model. Here we are fitting a whole map of parameters. Every pixel is potentially a new parameter, but they are not all independent. In MEM fitting, we minimize χ2−αS, where S is the entropy of the map, and α is a mixing parameter which controls how many independent parameters are being used in the fit. As we decrease α, the fit improves but the effective number of parameters rises. We set α so that the fit is ‘good enough’, i.e., χ2/N∼1, where N is the number of degrees of freedom.
The changes to the classical method are described in the following sections.
2.1 The entropy
Our definition of the entropy is essentially the same as in Horne (1985), and we also use Skilling & Bryan's (1984) MEMSYS algorithm. The entropy of a map is defined relative to a map of default values, the so-called default map. We have the freedom to select the default map, and thereby to define what we regard as a ‘simple’ map. Our usual procedure is to let default values be averages over nearby pixels, and the entropy thereby ‘steers’ the solution toward a smooth map, each pixel value resembling its neighbours. The steering is determined by the nature and magnitude of symmetry introduced to the reconstruction algorithm. The (Keplerian) orbiting material suggests that some amount of azimuthal symmetry is sensible. This default map is therefore used to steer the map into a certain symmetric direction, as long as it is compatible with the observed data. Thus the entropy is measured relatively rather than absolutely (e.g., in comparison to a ‘flat’ set of parameters). As in almost all astrophysical tomographic applications, the user thus has considerable freedom in choosing which default map—and hence final result—is appropriate.
For our purpose in the definition of the entropy, the intensities (as in the classical Eclipse Mapping) are replaced by the parameter values in both the map Pi,j and the default Di,j map (where i denotes the type of parameter, and j the number of pixel) to calculate the entropy S:
This way, the entropy still expresses the difference of the map compared with the default map. For the symmetry, as introduced through the default map, we chose small amounts of radial and azimuthal smearing.
The MEMSYS algorithm is tailored for linear models in which the predicted data values are linear functions of values in the map. PPEM model spectra will almost always be very non-linear functions of the fitted parameters. Our method works due to the fact that only short steps using the local response matrix connecting the observed data with the parameters of the model are taken. Not unexpectedly, PPEM solutions converge much more slowly than those obtained in classical Eclipse Mapping, and considerable care must be made to ensure that a globally optimal solution is obtained.
2.2 The central star
In quiescent dwarf novae, the eclipse light curves clearly show the eclipse of the white dwarf by the late-type secondary dwarf star. In classical Eclipse Mapping this eclipse is removed before an intensity map is reconstructed. In our new algorithm, we include the white dwarf in the tomogaphic analysis as a spherical object at the centre of the disc. Hereby, we assume a radius of the white dwarf with no limb darkening (see below) as derived from the white dwarf eclipse width (e.g. Wood et al. 1989, 1992). ‘Spherical’ is here meant in the sense that the region hidden from sight by the white dwarf is calculated as the shadow of a spherical object, and the visibility of the white dwarf due to the eclipse of the secondary is calculated using a spherical surface.
As was already suggested by Wlodarczyk (1986), the eclipse of the inner back part of the disc by the central star can be an important effect, depending on the size of the white dwarf, the brightness of the inner disc as well as the inclination angle of the system. For example, the system HT Cas has a large and bright white dwarf and a high inclination (81°) (Horne, Wood & Stiening 1991), so the effect is strong and affects a relatively large area behind the white dwarf (the length of the white dwarf shadow is RWtani=0.164RL1). In contrast, the length of the shadow of IP Peg's white dwarf is short (RWtani=0.053RL1) (Marsh 1988), and the relative contribution of the white dwarf to the total light much less.
Our program includes the occultation of the pixels behind the star and the shadowing by the secondary star of pixels on the white dwarf's spherical surface. Currently, any geometrical or mutual irradiation effects of, e.g., the boundary layer are ignored, and the white dwarf radius is not adjusted in the course of the program. At the moment, the white dwarf is effectively treated as a single pixel, so that intensity variations on the surface of the white dwarf, like limb darkening, are not included. A spatially resolved white dwarf surface could be used if the phase resolution of the observations is large enough that the shadow on the secondary moves across the surface of the white dwarf over the course of several data points. This number is dependent on the inclination angle i and the relative radius of the white dwarf RW/a, where a is the binary separation. The phase width Δφ corresponding to NW pixels on the visible side of the white dwarf is:
During the reconstruction, the entropy is calculated as a sum of contributions from all pixels. Here, we exclude the ‘pixel’ of the white dwarf, because the entropy constraint steers the reconstruction towards a smooth map, where neighbouring pixel assume similar values. This way, we allow an optically thick white dwarf adjacent to an optically thin inner accretion disc.
In order to treat the white dwarf correctly, spectra for the white dwarf have been included for the temperature range of 10 000 to 30 000 K. Outside this range, blackbody spectra are used.
2.3 The disc geometry
Following the dominant symmetry in the disc, we use a polar rather than the usual Cartesian grid. One can still allow for equal pixel areas by choosing the number of pixels in the following way: if the number of pixels in the central circle is Na, the number of pixels in the following equidistant rings are 3Na,5Na,7Na, etc. The number of pixels in the map N is then determined by N=Nr2Na, where Nr is the number of rings. The area of each pixel is (fRL1)2π/N, with fRL1 the radius of the circular region in which the parameters are reconstructed.
If the central star is included, the number of pixel in the map increases accordingly by the number of pixels NW used for it to N+NW.
2.4 The uneclipsed component
Following Baptista & Steiner (1993), we allow for an unspecified uneclipsed component. This may include contributions from the secondary, disc regions which are never eclipsed, winds, etc. No assumption about the nature of this component is made. Unlike in Baptista & Steiner's classical Eclipse Mapping solution, the uneclipsed component is treated as a phase-constant flux contribution in the light curve. No pre-determined correlation between contributions in different passbands is assumed. Instead, after the reconstruction has finished, these reconstructed fluxes can be analysed: e.g., a strong IR flux originates very likely on (the back side of) the secondary. This procedure avoids systematic errors, because a real contribution from, e.g., the back side of the secondary would otherwise have to be reconstructed in a never-eclipsed part of the disc region, or, more likely, no good fit to the data can be achieved.
In classical Eclipse Mapping this uneclipsed component is determined by the entropy, since if it exists but is not considered in the mapping procedure, it will produce additional intensity in the back part of the disc which are least eclipsed. Therefore the constraint on the uneclipsed component is not very strong. In PPEM we force the disc to emit in the fashion of a given spectral model. Therefore any additional component, which is not fitting to the spectral emission and at the same time independent of phase, will be reconstructed in the uneclipsed component. In case we choose a wrong model, the uneclipsed component may be over-estimated, but it cannot cover all emission not fitting to the observations, since it does not change the eclipse depth. For example, a real disc emits a very hot, eclipsed component, and we choose a model which forces the emission in U to be always similar to the emission in B (optically thin versus optically thick). Then the eclipse depth in the U light curve will not be well reproduced, no matter how much flux will be reconstructed in the uneclipsed component.
The flux at mid-eclipse is the sum of the uneclipsed component and the contribution to the specific passband from uneclipsed parts of the disc, i.e., at the edges. Usually, the infrared eclipse light curves determine the size of the disc and the temperature at the edge. In contrast, the UV light curves tend to be more sensitive to the inner regions of the disc.
2.5 Model spectra
PPEM can handle any physical model which provides a relationship between given physical parameters and the spectral intensity distribution, which we average over a bandpass and multiply by a solid angle to obtain the flux of each disc surface element. If a spectral model much more complex than a blackbody is used, it becomes impractical to compute the spectra in real time. Thus the program provides the option to use a precalculated grid of model spectra. The disadvantage is that grid-interpolations may cause problems in regimes of steep gradients (e.g., the Balmer jump). This can be minimized by using a flexible grid with higher resolutions where it is needed.
2.6 Reliability
In order to judge how powerful a given model really is, it must be tested on artificial data. We therefore create a realistic disc map for each parameter, calculate the observed or here rather original light curve, and apply the PPEM code in an attempt to reproduce the parameter distribution(s) by fitting the original light curves. Comparisons with the original parameter distributions together with a study of the model (particularly important for multiple-parameter versions) lets us assess the reliability of the method, its power and—most importantly—its limits.
Because of forthcoming analyses of real data, we adopt the geometry of real systems, i.e., that of IP Peg (inclination angle i=79°.3, mass ratio q=M2/M1=0.587; Marsh 1988) for the one-parameter, optically thick case (Section 3), and those of HT Cas (i=81, q=0.15; Horne et al. 1991) for the two-parameter, optically thin case (Section 4).
3 One-Parameter Mapping: Temperature
The simplest model spectrum one can think of is that of a blackbody, Bν(T), leading to a single-parameter version of the new method: Temperature Mapping (Vrielmann et al. 1996).
3.1 Application to synthetic data
The original and reconstructed light curves are displayed in Fig. 1. The fit to the data is very good, with a χ2 of 1.0. The temperature distribution is displayed as a radial profile in Fig. 2, and as grey-scale plots in Fig. 3. The input temperature distribution is very well reproduced: the radial temperature profile follows the original and shows almost no azimuthal variation. Only minor deviations can be found; these can be attributed to not unexpected systematic effects, whose origin we will now discuss in some detail.
The test eclipse light curves in U, B, V, R, J, H and K (dots with error bars, S/N = 100) and the fitted eclipse light curve of the reconstructed blackbody disc map (dotted-line light curves are offset by 0.4 mJy each). The residuals between the input and the model light curve are displayed at the bottom (χ2=1.0).
The test eclipse light curves in U, B, V, R, J, H and K (dots with error bars, S/N = 100) and the fitted eclipse light curve of the reconstructed blackbody disc map (dotted-line light curves are offset by 0.4 mJy each). The residuals between the input and the model light curve are displayed at the bottom (χ2=1.0).
The temperature distribution of the blackbody test map (solid line) and the reconstructed one (circles) by fitting all seven light curves simultaneously.
The temperature distribution of the blackbody test map (solid line) and the reconstructed one (circles) by fitting all seven light curves simultaneously.
Grey-scale maps of the original (top) and reconstructed (bottom) temperature distributions. The line surrounding the disc is the Roche lobe of the primary; the system is rotating counter clockwise. Theoretical accretion stream paths for mass ratios 0.587 ± 50 per cent are plotted as well.
Grey-scale maps of the original (top) and reconstructed (bottom) temperature distributions. The line surrounding the disc is the Roche lobe of the primary; the system is rotating counter clockwise. Theoretical accretion stream paths for mass ratios 0.587 ± 50 per cent are plotted as well.
Beyond the outer rim of the input disc, the reconstructed temperature distribution displays a residual temperature of about 800 K. This temperature is undetectable at the wavelengths used: the derived parameter distributions can be trusted only in regions where the intensity calculated from the current parameter set is significant in at least one of the wavelength bands. Otherwise, the maximum entropy criterion forces a value from the default map. Thus the distribution is as smooth as long as it is compatible with the observed data.
The temperatures of the white dwarf and the central ring fail to reach exactly the original values by about 5 per cent, while the temperature in the second ring is slightly increased by a similar amount, i.e., the temperature gradient is slightly less steep in the reconstruction compared to the original one. This is also caused by the MEM, which prefers smoother gradients as long as they are compatible with the data.
In real data with a lower S/N than we used in the simultations and additional ‘noise’ due to flickering this effect may even be stronger. Small-scale structures in the temperature distribution may be smeared out slightly, but not as strongly as in classical Eclipse Mapping (see Section 3.4), since the colour information provided by the use of light curves at different wavelengths restricts the temperatures to a certain range. The eclipse shape and, in particular, the eclipse width hereby force the smoothing of the reconstruction to stay within a limited spatial region. However, while analysing reconstructed maps of real data (Vrielmann 1997; Vrielmann, Hessman & Horne in preparation), we have to keep this effect in mind.
3.2 Distance effects
The previous fit was made assuming that the distance to the system is known beforehand. If this is not the case, then the distance must be fitted as another parameter (though, of course, a parameter with very different properties from those contained in the maps, justifying the continued use of the term ‘one-parameter map’). If the wrong distance is assumed, the solid angles of the disc pixels change. Since both the overall flux levels and the colours of the disc spectrum have to be fitted, the choice of an incorrect distance forces the fit to modify artificially both the effective size of the disc and the temperature distribution to stay within the required constraints. Fig. 4 shows the result of such an incorrect choice, where the distance was assumed to be too large by 20 per cent. Although we aimed only at a χ2 of 4, the following effects are already obvious: (a) the size of the reconstructed disc is significantly larger in order to reproduce the observed flux; (b) the average temperature is larger than the original one; and (c) the noise in the temperature distribution is much larger than for the reconstruction with the correct distance.1 (Aiming at a lower χ2 would have lowered the temperature of the white dwarf, but even increased the above mentioned effects. Since the fitting is very slow, we stopped at this stage.) While results (a) and (b) give no criteria for the correct distance (in real systems we do not know the correct disc size and the actual temperatures), the last one (c) can easily serve as a distance indicator by choosing the distance according to the smallest noise level in the reconstruction, expressed by the maximal entropy. Hereby, real structures like the bright spot should be distinguishable from the pure scatter.
3.3 Model effects
To test our optically thick model against choosing the wrong model, we created an artificial accretion disc which emits optically thin accoFor this disc we calcrding to the model given in Section 6. While the temperature decreased from 24 000 to 3000 K, the surface density increased from 0.01 to 0.06 g cm−2. ulated the light curve in the mean wavelength of the UBVRJHK passbands, and analysed this set of light curves with the Temperature Mapping algorithm.
As expected, the algorithm has most difficulties reproducing the light curve in the ultraviolet wavelength—in this case, the only passband beyond the Balmer jump. It tries to reproduce the deep eclipse in U by increasing the white dwarf temperature to unreasonably high values, but fails. The reduced χ2 of the fit in this case cannot get lower than about 600. Furthermore, the reconstructed discs show strong non-axisymmetric structure, similar, though not the same as, ingress and egress arcs. Additionally, the reconstructed light curve fails to reproduce the smoother ingress and egress.
The tests also showed that the closer we get to an optically thick model, the better the fit and the better the reconstruction of the temperature. As shown in Section 4.1, the location of the transition from an optically thin to thick region depends on the temperature and the radius in the disc. If the disc is partially optically thick, these regions can be reproduced with the Temperature Mapping analysis. However, the presence of truly optically thin regions leads to (a) too much structure (pixel-to-pixel variations and ingress and egress arcs in the distribution), (b) an unreasonably hot white dwarf, and (c) poor fits to the observed data. These three criteria can be used to determine if the chosen model is adequate.
3.4 Comparison with the classical approach
The classical Eclipse Mapping method yields intensity distribution reconstructions of the accretion disc in a given passband. Usually, these intensity maps are derived independently from the eclipse data at each wavelength, and then presented as a set of brightness temperature maps, where the brightness temperature is the temperature of a blackbody that has the same intensity in this passband (e.g. Rutten et al. 1992a).
For the purpose of comparing the results obtained by PPEM and classical Eclipse Mapping, a single, unique temperature distribution can be obtained from the latter in the following manner: after the intensity distributions are calculated, a Planck spectrum can easily be fitted to the flux data for each pixel, yielding a fitted temperature map (e.g. Rutten et al. 1993; compare Horne, Wood & Stiening's 1991 TΣ-model used to fit the intensity maps for HT Cas).
This way, we constructed a single fitted temperature map using the classical Eclipse Mapping method on similar test data, and compared the results with the input temperature distribution (Fig. 5). Each intensity map was weighted equally. The resulting temperature distribution qualitatively resembles the original distribution, but the reconstruction is strongly smeared out, both azimuthally and radially, and does not reproduce, e.g., the steep gradient at the edge or the centre of the disc. The values of the temperature in the inner disc assume much lower values of around 15 000 K, where the original reaches a maximum of twice this value.
The temperature distribution of the blackbody test map (solid line) and that obtained by fitting all seven light curves independently, in a classical Eclipse Mapping approach (circles). The fits to each of the seven light curves (cf. Fig. 1) reached a χ2 of 1.0. (The temperature of the white dwarf was reconstructed to a temperature of 4025 K.)
The temperature distribution of the blackbody test map (solid line) and that obtained by fitting all seven light curves independently, in a classical Eclipse Mapping approach (circles). The fits to each of the seven light curves (cf. Fig. 1) reached a χ2 of 1.0. (The temperature of the white dwarf was reconstructed to a temperature of 4025 K.)
This result is easily explained: a single light curve does not force any single intensity distribution to have steep gradients. As long as in certain regions the intensity is low enough to have no significant contribution to the light curve the distribution is determined mainly by the entropy constraint. Combining the maps obtained independently from the lightcurves at several wavelengths does not yield a steeper gradients (e.g., in the inner disc or the edge). Especially in the outer regions, where the entropy constraint is most important, the fits to the pixel spectra are very poor. To improve the fits, the intensity distributions with larger contributions could be given larger weights in the fit. However, this method would still not have the power of our direct approach, since the information can only be recovered if a single map is fitted simultaneously to the light curves at all the wavelengths.
The comparison of this reconstruction (Fig. 5) with that of the PPEM method (Fig. 2) shows a clear difference: while reaching the same reduced χ2, the PPEM reconstruction follows the original distribution in most parts of the disc more closely. The general shape, the steep gradient in the inner disc, and the edge are reproduced, whereas the reconstuction with the classical method only qualitatively follows the original distribution.
The classical method's advantage is that it makes no specific assumption about the local physics, allowing the spectral intensities to deviate considerably from a blackbody, i.e., allowing the brightness temperatures to be different at each wavelength. However, the disadvantage is that the maps at each wavelength are differently affected by the smoothing effect of the entropy. The smoothing is also more pronounced at wavelengths with low S/N eclipse data. This wavelength-dependent redistribution of intensity on the disc surface distorts and complicates the subsequent interpretation of the maps at different wavelengths.
The PPEM approach uses a specific model to fit all wavelengths simultaneously. This reduces the smoothing effects because all of the data are used, so features detected in the light curve at even one wavelength are sufficent to override the entropy's efforts to smooth the map.
4 Two-Parameter Mapping: A Uniform Slab
To illustrate a two-parameter model, we will use the spectrum of an LTE, pure hydrogen uniform slab with no line emission (Vrielmann 1997; Vrielmann et al. 1997). Here, in addition to the pure temperature dependence of the optically thick model, the optical depth of the material is important:
For τ/cosi≫1, equation (4) transforms into the optically thick case, Iν≈Bν(T). For the optically thin case τ/cosi≪1, equation (4) reduces to Iν(T)≈τ/cosi×Bν(T).
In this version of PPEM, we adopt T and Σ as the two physical parameters for which we seek reconstructed maps (‘T−Σ Mapping’). Other choices are possible: for example, we could have chosen Teff and ρ, and even considered the Gaussian density distribution of an isothermal disc. In any case, MEM forces the maps to take on positive values only, thereby automatically eliminating non-physical solutions.
4.1 A study of the model
Wood et al. (1992) have suggested that a two-parameter model describing the radiation in terms of temperature and (surface) density might contain ambiguities. Thus it is very important to understand the reconstruction behaviour of our chosen spectral model.
To evaluate the ability of U,B,V and R data to constrain the parameters T and Σ (for a forthcoming application to real data), we calculated the spectrum Iλ(Ti,Σi) for a grid of temperatures Ti and surface densities Σi. To stay close to the real data, we used the parameters of HTCas, a mean disc radius of R=0.32RL1=1×1010 cm, and the mean wavelengths of the U,B,V and R passbands (3550, 4370, 5340 and 6520 Å). Fig. 6 shows a few of these spectra for a temperature range of 3000 to 20 000K and a surface density range of 10−3 to 10−1 g cm−2. The given wavelengths cover one point in the Balmer continuum, and three in the Paschen continuum. As one progresses from lower to higher surface densities, the most prominent change is the decrease of the Balmer jump: it vanishes completely as the LTE slab becomes opaque. The rapid transition to the optically thick state is due to the ionization of hydrogen, which drives the opacity up. At the same time, the magnitude of the emitted spectrum decreases with decreasing surface density, since the material becomes more and more transparent. The disc scaleheight H increases as T1/2. Since ρ=Σ/2H, the optically thick limit moves slightly towards higher values of the surface density as one increases the temperature. For low temperatures, the temperature dependence itself appears as a changing shape of the continua (resembling the blackbody distribution), while for higher temperatures and low surface densities, only the magnitude of the emitted intensity decreases with increasing temperature (the maximum of the emission moves towards the blue spectral range, reflecting the Wien limit for optically thick emission).
Uniform-slab spectra for different combinations of temperature and surface density. The temperature increases from bottom to top (T= 3000, 6000, 12000 and 20 000 K) and the surface density from left to right (Σ=1×10−3, 5×10−3, 2×10−2 and 1×10−1 g cm−2). The points give spectra for 1 per cent variation in the parameters values to illustrate the sensitiveness of the model (different symbols only for different temperatures).
Uniform-slab spectra for different combinations of temperature and surface density. The temperature increases from bottom to top (T= 3000, 6000, 12000 and 20 000 K) and the surface density from left to right (Σ=1×10−3, 5×10−3, 2×10−2 and 1×10−1 g cm−2). The points give spectra for 1 per cent variation in the parameters values to illustrate the sensitiveness of the model (different symbols only for different temperatures).
For each spectrum in the grid we choose a region ΔTi,ΔΣi around the given pair (Ti,Σi) in which we set up a subgrid (Ti,j,Σi,j) to calculate spectra and determined the difference between the test spectra Iλ(Ti,j,Σi,j) and the original Iλ(Ti,Σi). The difference is expressed as the χ2 between the spectrum in the subgrid and the original spectrum, where we assume an uncertainty of 1 per cent of the maximum value in each spectrum. Fig. 7 shows contour lines for the χ2-levels, χ2=1,2 and 3 for each original spectrum. In the region of high surface densities (above ∼5×10−2 g cm−2), the spectra proceed into the optically thick limit and become independent of the surface density—expressed by the long χ2 contours. The elongated, slightly inclined χ2-levels in the region of high temperatures and low densities reflect the similarity in the decrease in the intensity for increasing temperatures and decreasing surface densities without changing the shape of the spectrum. A similar situation occurs in the low-temperature region, especially at intermediate surface densities. Here, the intensity increases with increasing temperature, leading to banana-shaped contours in an almost perpendicular direction. In this region, the emission in the short wavelength range is so low that the Balmer jump disappears; however, the intensity is still dependent on the surface density. The inclined shapes of the χ2-levels means that we will find a well-defined parameter (e.g., T) for a given value of the other parameter (e.g., Σ). However, since usually both parameters are to be fitted, ambiguities will occur in the reconstructions in such disc regions, and especially if the quality of the data is low. Therefore we expect that the resulting parameter distributions in regions where the parameter pair assumes such ambiguous combination of values will depend on the steering of the MEM and at least slightly on the initial distribution. This makes the choice of a reasonable initial map very important.
χ2 contours for different parameter combinations of temperature and surface density for a disc radius R=0.32RL1=1×1010 cm. See text for details.
χ2 contours for different parameter combinations of temperature and surface density for a disc radius R=0.32RL1=1×1010 cm. See text for details.
Between 5000 and 6000K for Σ=10−3 g cm2, and at ∼ 7000K for Σ=10−1 g cm2, the Balmer jump appears in the spectrum, giving good sensitivity to both the surface density and the temperature. This leads to very well-defined parameter pairs (T,Σ) reflected by the narrow χ2-levels. Here we expect PPEM to find a unique solution. A model including emission lines should lead to a wider ambiguity-free zone, because the lines have a different dependence on the temperature. However, their opacity is larger than in the continuum, which means that they become optically thick at lower surface densities. Still, Horne & Hubeny (1998) therefore find a very similar distribution for the χ2-levels using Hubeny's disc model atomspheres (Hubeny 1991).
This study shows the importance of understanding the reconstruction behaviour of the model with which the data are to be fitted. Only when we know the behaviour of the parameters can the resulting maps be interpreted with confidence. We have to keep in mind that the PPEM method cannot give us more information about the accretion disc than the used model allows. Nevertheless, as we will show in the next section, in the regions where at least one parameter is welldefined, the method still yields very good partial-reconstructions.
4.2 Application to synthetic data
Let us consider an artificial accretion disc with a temperature and surface density distribution similar to that expected for quiescent dwarf novae and covering a range of probable parameters. For this test case, the white dwarf is assumed to emit blackbody radiation, i.e., its effective surface density was set to Σ=1 g cm−2. The white dwarf temperature was set to a relative low value of 6000 K, in order to test if the method finds such a low white dwarf temperature if the surrounding disc is much hotter. The artificial disc gas temperature is assumed to follow a theoretical steady state effective temperature distibution (equation 3). This kind of distribution would only be expected for the effective temperature in a steady state disc, which might not reflect the distribution of the gas temperature in the optically thin case. However, this is only a test of the method, and it will show how far the method is capable of reproducing power-law distributions covering temperatures from 5000 K to almost 19 000 K. The surface density is chosen to cover the range from optically thin to thick emission, crossing the unambiguous T−Σ region at a radius of ∼0.2RL1. A linear increase in surface density with radius was chosen according to theoretically expected distributions (e.g. Ludwig, Meyer-Hofmeister & Ritter 1994).
From these parameter distributions, we calculated the corresponding light curves for the average wavelengths of the U,B,V and R passbands (3550, 4370, 5340 and 6520 Å). These light curves were used to derive reconstructions of the parameter distributions. For the reconstruction, the distance was assumed to be known.
In Fig. 8 the light curves are displayed with the resulting fits. Fig. 9 shows the reconstructed parameter distributions as a comparision with the original distribution. While the temperature distribution follows the original quite well in the intermediate disc regions, the surface density distribution is best reconstructed in the inner part of the disc.
The test eclipse light curves in U,B,V and R (dots with error bars, S/N = 100) and the fitted eclipse light curve of the reconstructed blackbody disc map (dotted line, light curves are offset by 0.1 mJy each). In the bottom part the residuals between the input and the model light curve are displayed. The fit is excellent (χ2=1.0).
The test eclipse light curves in U,B,V and R (dots with error bars, S/N = 100) and the fitted eclipse light curve of the reconstructed blackbody disc map (dotted line, light curves are offset by 0.1 mJy each). In the bottom part the residuals between the input and the model light curve are displayed. The fit is excellent (χ2=1.0).
The temperature and surface density distributions of the test map (solid lines) and the reconstructed ones (circles) by fitting all four light curves simultaneously. The value of white dwarf's surface density is not displayed to allow the examination of the behaviour in the disc. It was reconstucted to an optically thick value indistinguishable from the original.
The temperature and surface density distributions of the test map (solid lines) and the reconstructed ones (circles) by fitting all four light curves simultaneously. The value of white dwarf's surface density is not displayed to allow the examination of the behaviour in the disc. It was reconstucted to an optically thick value indistinguishable from the original.
Although the surface density of the white dwarf could have been held fixed for physical reasons, it was reconstruced in the same way as the parameter values in other pixel for a consistency test of the method.
4.3 Discussion
Considering the results of the reconstruction behaviour in Section 4.1, we can now interpret the reconstructed parameter distributions. Plots of spectra from four different regions of the reconstructed disc as displayed in Fig. 10 and the associated intensity distributions are shown in Fig. 11. The following features are apparent.
Spectra for four different regions in the reconstructed temperature and surface density distributions corresponding to the radii 0.1RL1, 0.25RL1, 0.4RL1 and 0.6RL1 as indicated. The points give spectra for 1 per cent variation in the parameters values to illustrate the sensitiveness of the model (different symbols only for different temperatures).
Spectra for four different regions in the reconstructed temperature and surface density distributions corresponding to the radii 0.1RL1, 0.25RL1, 0.4RL1 and 0.6RL1 as indicated. The points give spectra for 1 per cent variation in the parameters values to illustrate the sensitiveness of the model (different symbols only for different temperatures).
Intensity distributions for the four wavelengths corresponding to U,B,V and R as indicated. In each case the solid line represents the original distribution, while the circles are calculated from the reconstructed maps. The scatter in the intensities shows how sensitive it is for slightly deviating parameter combinations.
Intensity distributions for the four wavelengths corresponding to U,B,V and R as indicated. In each case the solid line represents the original distribution, while the circles are calculated from the reconstructed maps. The scatter in the intensities shows how sensitive it is for slightly deviating parameter combinations.
As expected, the temperature and surface density of the white dwarf are well constrained and fairly well reproduced (reconstructed temperature Trec=4831 K versus original Torg = 6000 K, and Σrec=Σorg=1.00).
At intermediate radii, the temperature distribution of the disc is very well reproduced. The low temperatures of about 3000 K in the outer regions (>;0.5RL1) contribute only very little to the total light (maximal in R with about 3 per cent of the intensity in the inner regions with a temperature of 12 000K; see Figs 10 and 11), which means that the entropy constraint determines the profile. In the inner disc regions, the temperatures are much lower than expected, but—with values of about 13 000 K—they lie in the region of the T−Σχ2-plane where the temperature is not very well constrained (see Fig. 7).
The reconstructed surface density distribution follows the original one only in the inner disc up to about 0.25RL1. Further out, it is significantly smoother and assumes lower values than the original. As seen in Fig. 10, the spectrum in this parameter region (5000 K, 0.028 g cm−2) is mainly determined by the flux in R. The intensity distribution (Fig. 11) shows that the contribution in R to the intensity in this disc region is only about 3 per cent. In this case, this ambiguity could be resolved by including further light curves in the infrared.
As mentioned above, the distance was assumed to be known. In practice, the distance can be pinned down by the need for the physical model to fit the spectrum over a large range of physical conditions simultaneously. Furthermore, IR observations of the secondary and/or UV observations of the white dwarf also constrain the distance, and values derived by those methods can serve as trial distances. Just as in Temperature Mapping, the assumption of a wrong value for the distance will degrade the quality of the fits to the observed data by introducing noise into the parameter distributions.
This analysis makes clear that—in principle—the temperature and surface density of optically thin accretion discs can be adequately reconstructed using two-parameter PPEM over a scientifically interesting range of radii. As the application to artificial data show, it is essential to study the used model in terms of ambiguities of the dependent physical parameters. However, the Physical Parameter Eclipse Mapping method cannot resolve ambiguities present in the model used. The reconstruction behaviour of a given model must be studied over a wide range of parameter values in order to determine where the parameters can be trusted (see Section 4.1). It is also advisable to use as wide a spectral range as possible for this kind of analysis, since each wavelength adds a different constraint on the physical parameters.
4.4 Model effects
Also in this case we performed a cross-examination according to the model used. In the same way as described in Section 3.3, we created light curves from an artificial optically thick, steady state disc (for a mass accretion rate of M=1017 g s−1), and analysed them with the optically thin version of the PPEM.
Since the optically thin model transforms to a optically thick model at a certain surface density (depending on temperature and radius), we should always be able to reconstruct the accretion disc. Even if we start with a completely and extremely optically thin emitting disc, we reconstructed a optically thick emitting disc with the surface density just beyond the transition zone from optically thin to thick. The effective temperature distribution follows mainly that of a steady state disc of −3.5m=1017g/s already for a rather poor fit to the data with χ2 of about 20. The fitting, however, is very slow, using several thousands of iterations compared to several hundreds if we use the right model. Thus the PPEM proves to be a reliable and stable method.
5 Conclusions
We have presented the method of Physical Parameter Eclipse Mapping and have described one-parameter (‘Temperature Mapping’) and two-parameter (‘T−Σ Mapping’) versions capable of being applied to real spectrophotometric data of cataclysmic variable systems. The physical models we have used so far are relatively simple, but they illustrate both the power and the limits of the new method.
Detailed studies of the multiple-parameter model spectra are necessary before results can be interpreted. Depending upon the model, there will be regions of parameter space—and hence spatial regions within the disc—for which one or more parameters may not be well constrained. Sometimes, but not always, the addition of more data at different wavelengths can resolve the inherent ambiguities.
Physical Parameter Eclipse Mapping is no panacea in our search for more detailed information on the physical properties of and processes occurring within accretion discs. However, by simultaneously using all of the available spectrophotometric information, it does hold the promise of providing more detailed and physically interpretable insights into these very intriguing astrophysical systems.
Acknowledgments
We thank R. Rutten, H. Spruit, and especially R. Baptista for many insightful discussions and copies of their classical Eclipse Mapping codes. SV was supported by the European Commission through a TMR fellowship. Furthermore, we thank the referee for helpful comments.
