FUSE EM Cyg

Data Reduction:  Time-averaged Spectrum

EM Cyg was observed with FUSE 2002 September 5.  The time-tag spectra were written as four exposures with a total exposure time of 8911 sec (7570 sec at night):

    Exposure                                           UT start      Texp            Phi
c01001010011alif4ttagfcal.fit[0]    11:34:11    1479 sec    0.415 - 0.474
c01001010021alif4ttagfcal.fit[0]    12:54:33    2651 sec    0.607- 0.712
c01001010031alif4ttagfcal.fit[0]    14:34:31    2647 sec    0.846 - 0.951
c01001010041alif4ttagfcal.fit[0]    16:22:59    2134 sec    0.104 - 0.189

The archive spectra were reduced using CALFUSE V.2.2.1.  I examined the spectra in fuse_register and there were no obvious problems during data acquisition. There is an error in this version of CALFUSE that affects the flux calibration of the spectra, so we need to re-reduce the spectra with CALFUSE V.2.3.  I created four concatenated files using ttag_combine.  Their names are c010010001attagfraw.fit, etc.  I examined these images using fuse_scan and didn't see any anomalies. CALFUSE does not appear to run properly on ttag_combined exposures when jitter files are present:  the jitter files cannot be combined with ttag_combine, and when CALFUSE only finds the first jitter file, it rejects the rest of the time windows in the spectrum.  [See email from Van Dixon for work-around.]  I reduced the 1a exposures individually to check whether any jitter correction is performed:  it isn't.  As a result, I re-ran CALFUSE on the combined spectrum with jitter correction turned off.

1A:  Ran fine, no rejection for bursts.  44 sec were removed from all 4 channels for a SAA incursion.
1B:  The centroid of the SiC1B spectrum was outside the allowed shift range, so Calfuse moved the extraction window back to its nominal poisition and consequently missed the spectrum. I upped the EMAX_SIC parameter in parm1b009.fit to 35 pixels and re-ran the pipeline.
2A:  Looks good.
2B:  Ditto.

Combined segments in fuse_register.  Masked out LiF1B longward of 1135 A and LiF2B shortward of 995.  There was no need to rescale any of the segments to correct for drift.  I cross-correlated the segments to line them up in overlapping wavelength regions.  The time-averaged spectrum, weighted using the errors, is  emcyg.fits.  I then used IRAF trebin to bin the table to 0.1 A dispersion and wrote the result to emcyg01.tab.  The ascii spectrum of that is in emcyg01.dat.  The S/N in the 0.1 A spectrum is not so hot, so I also created spectra binned to 0.5 A, emcyg05.tab and emcyg05.dat.

Properties of the time-averaged spectrum

The IUE quiescent spectrum (swp08088) shows a continuum of 2 to 3e-14.  CIV and SiIV are in emission;  NV and SiIII 1300 are in absorption.  In outburst, the continuum is 5 to 6e-13 at 1200 A.  All lines are in absorption except CIV and HeII 1640.

The FUSE spectrum has a continuum level around 2.5e-14 and is richly populated with spectral features.  There are absorption features with FWHM around 500 km/s.  CIII 977 and, I think, 1176 appear to be in emission (the latter with an absorption core) but the other lines appear to be in absorption.  It is difficult to be definitive on this point because there are so many lines that continuum placement is uncertain.

Lines in the FUV spectrum:

SVI 933 and 944
HI 972 ?
CIII 977 and 1175
OVI 1032 and 1038
SIV 1062 and (less clearly) 1072
He II 1085 (and/or NII 1084)
Si III 1108,1110,1113
PV 1117/1128 and SiIV 1122/1128 are tentative


Data Reduction:  Time-Sampled Spectra

The FUSE proposal states that a sufficient S/N (14-22) can be obtained in combined spectra of 60 sec duration and 5 A binning.  That's 148 spectra.  I don't think we are going to get the S/N in the lines than anticipated since they aren't so strong in emission, but I'll go ahead and extract spectra in 60 sec increments and bin up if necessary.  I wrote an IRAF script, docal.cl that hedits the scrn1a013.fit file for successive time increments, calls cal_fuse_exp_ttag.csh externally, deletes unwanted aperture spectra, and renames the output spectra and trailer files.  I turned off burst and jitter correction in CALFUSE, as that is often squirrelly in short time integrations.  I will also keep an eye on the deadtime correction, which was an issue in the UX UMa 60 spectra.

The output spectra were written to ascii files as lif1a_001.dat, etc. These can then be binned up in time or orbital phase and used for lightcurve construction.  I also created a file, all, that gives the exposure names, times relative to the start of the observation, and the mid orbital phase of each 60 sec exposure.  The deadtime correction was 1.003 for all of the exposures-- a-ok.

Checking the reliability of the 60 sec reduction:  I compared a time-averaged spectrum of the 60 sec data with a time average of the 8 segments from the ttag_combined spectrum.  The two are identical at longer wavelengths, but the 60 sec combined spectrum is fainter <980 A.   I tried removing first SiC1B, then SiC2A.  The spectrum sans SiC1B is good agreement with the overall spectrum, while the spectrum without SiC2A is too faint.  It looks like the response of SiC1B drops off rapidly below 980 A. Oops:  the problem centroiding SiC1B that we had above was repeated here.  I had changed the EMAX_SIC back to 30 pixels;   re-ran.

FYI, (after fixing a masking error in the fuse_register combined spectrum) the time-averaged spectrum created by fuse_register is virtually identical to one I made using my usual perl scripts (although I didn't weight by the errors in the latter) .

I used the 60 sec spectra and script shift.pl to create a time-averaged spectrum with the orbital motion of the primary removed.  The spectrum is emcyg_K1.dat/.fits and I put it in Fcal.  It's binned to 0.5 A with no weighting and the standard wavelength V.14 masking.  There is no improvement in the line profile sharpness that points to an obvious identification of the lines with the WD K1.

Time-Series Analysis

We're interested in variability in the lines with time, but it's a good idea to take a look at orbital variability too, to make sure we aren't confusing the two.  Actually, this won't  work because we only cover each orbital phase once.  So, time-series variability it is, although I don't know how we can address questions concerning orbital variability too.  Since each original expousure covered a reasonably narrow range of orbital phase, we can compare each of them first and then look within each exposure for short-term variability.

I made up spectra in different combinations of time and spectral resolution:

timebin2700/      Original 4 exposures
timebin300/        300 sec and 2.5 A resn
timebin60/           60 sec and 5 A resn

Of the original 4 exposures:  the 1st (0.41 - 0.47) is about 2x brighter than the 2nd (0.61 - 0.71).  The 3rd (0.85 - 0.95) is very similar to the second and at the same flux levels and the 4th (0.10 - 0.19) is similar to the first, though slightly (20% or so) fainter. There are no dramatic changes in the shapes of the spectrum between exposures.  The absorption line EWs appear to increase with the continuum flux (definitely in CIII 1175, less obviously in the other lines) but the emission line strengths do not go up when the continuum rises.  Indeed, ratio plots show CIII 977 lower than the surrounding regions, suggesting that the emission line gets weaker as the continuum goes  up.  With the data we have, there isn't enough information to determine if these represent orbital changes or secular, so we will have to be cautious in interpreting variability between exposures.

Interpretation is going to be problematic because it isn't clear where the lines are and whether we have lines in emission or absorption predominantly.  Ugh.  Anyway, I am going to make light curves incorporating CIII 977 and 1175 , OVI 1032 and 1038 and pick nearby regions for comparison.

CIII 977: 975 - 980 (two bins)
nearby:  995 - 1000

OVI:  1030 - 1040
nearby:  1050 - 1060

CIII 1175:  1170 - 1180
nearby:  1155 - 1165

OVI and CIII 1175 are wider b/c they are blends (and to avoid the 972 airglow near CIII 977).  Script timecurve.pl makes average flux light curves from the 60 sec spectra in timebin60.  The light curves are in lcurves.
 

 Here is a summary with figures of the initial analysis of the spectra (08/21/03).

May 10, 2004

After talking with Bill Welsh last week about the FUSE data on EM Cyg, we've decided to publish a short data paper. For the light curves, we decided to re-do the time variability analysis, recentering the pseudo continuum in low areas a bit further from the lines of interest and to plot continuum-subtracted light curves for ease of interpretation.  I will also send the new light curves to Bill so he can search for correlated, lagged variability. 

Light curves:  We'll keep the line bins as given above.  The pseudo-continua will move to the following:

CIII 977:  965 and 1010-1020
OVI:   1010-1020 and 1065
CIII 1175:  1110 and 1130

I created a new directory, newlcurves, in which I put the old line light curves and the new continuum light curves.  The sm.plot scripts average the continuum regions, subtract them and plot the line and continuum flux light curves.  The sm.cmp scripts plot line flux vs. continuum flux and print out the light curves.

Disk and white dwarf model fits

System parameters

North et al. 2000, MN, 313, 383:  discovered that a 2nd K star contributes to the light in EM Cyg (in addition to the donor star).  They recalculate the system parameters and get q=0.88, i=67, M1=1.12 Msolar
Winter and Sion 2001, ApJ, 582, 352:  they model the IUE spectrum of EM Cyg in quiescence (swp08088) and find that the disk is the principal emitter based on combined disk-wd fits.  They give a range of distances for EM Cyg of 59 - 411 pc with 350 pc being the best current estimate.

Other data sets

STIS data in outburst.  HUT data in outburst.  The last spectrum is on decline but the flux is still 1e-13 in the FUV.  Most IUE spectra are in outburst but two SWP spectra -- 17591 and 08088 -- are not.  The latter is fainter (3e-14 vs. 7e-14 at 1350 A) but both are consistent with FUSE spectrum in the narrow region of overlap. 

I printed out the fluxes and errors (the latter requires that tprint be used on the *mxlo.fits files) of the two quiescent IUE spectra and averaged them with error propagation using a perl script.  The mean spectrum is emcyg_iue.dat, which I concatenated with the FUSE emcyg05.dat spectrum (using the FUSE data i the overlap region) to make /ghola/data/froning/Progs/kslfit/emcyg/emcyg.dat for the model fits.  I also made emcyg2.dat, in which the FUSE spectra are also at 1.7 A so that there is equal weighting between the two spectral regions in the fit.

White dwarf model fits

I first took a look at pure WD model fits to the spectra, assuming log g = 8 (about right) to start.  I had to modify NWAVES in kslfit upward again so it would read in model spectra covering 60,000 grid points between 900 and 4000 A.  The stellar spectra for WDs and (later) disk models are in /ghola/data/froning/stargrid/uvdisk and are smoothed to 0.5A.  Models smoothed to 1.7 A are in uvdisk2. 

For the models, Rwd = D sqrt(N/4*pi), where N is the model normalization to the observed spectrum.  A 1.12 Msolar WD should have a radius ~ [1 - Mwd / Mch]**3/5 where Mch = 1.44 Msolar (Nauenberg 1972).  This gives a WD radius in EM Cyg of about 4e8 cm.

The best fit for the FUSE range alone gives:  N=2.141e-24, T=46,380 K.  Assuming D=350 pc, we get Rwd=4.33e8 cm.  Assuming Rwd=4e8 cm, we get D=323 pc.  Pretty consistent with previous assumptions above, but the model doesn't fit the shape of the FUV continuum very well and is too steep to match the observed UV flux.

The best fit for the FUSE plus IUE range (both binned to 1.7 A) gives:  N=1.610e-22, T=19,620 K.  Assuming D, we get Rwd=3.76e9 cm.  Assuming Rwd, we get D=37 pc.  This WD is either oversized for its mass or is located much closer to us than we think.  Not likely.  The fit looks slightly better for this case, but still isn't very good.

I also tried varying T and log g:

The best fit for FUSE only gives:  N=5.583e-23, T=29140, log g=2.25 (?).  Rwd= 2.2e9 cm or D=63 pc. Not a great fit.

The best fit to FUSE+IUE gives:  N=1.067e-21, T=12290, log g=3.92.  Not a bad looking fit, except that it gives Rwd=9.68e9 cm or D=14 pc.  The log g implies one odd WD, too, given its mass.

Disk models

I have to up NWAVES in cv_disk.c and 3 or 4 other programs in that directory for program cv to run on the larger model stellar atmospheres I'm using.  The disk models are in /ghola/data/froning/Progs/cvmodels/emcyg.  The "disk only" models actually have a 10,000 K WD in them as well as I can't get Knox's program to run without them properly.  Anyway, that WD contributes at most 6% of the flux in the FUV range.

Disk-only, full range, gives:  norm=13.4, log mdot = 15.0 g/s (the min), chi-sq 5.577.  The fit is similar to that of the funky WD with low T, log g above.   The normalization is an order of magnitude off for the 350 pc distance.  The normalization is around 1 for log mdot = 16.2 but the slope on a disk of that accretion rate is too steep to fit the full UV+FUV range.  A fit to the FUSE region alone gives N=0.13 and log mdot=16.83 but is of course too steep to account for the UV continuum as well as the FUV (the disk model provides about 1/3 of the observed flux at 2000 A).

Disk plus WD models

Finally, I combined the disk and the WD.  Program cv can do this but my large stellar models are making that program rather slow.  I'm trying to run disk plus WD in kslfit which is also slow but not as slow.

Disk 1:  disk plus WD, 920-1975 A

WD: 2.01e-22, 17,790 K
disk:  0.008, log mdot=18.48
E(B-V)=0.025

The parameters are not terribly realisitic:  the WD is about an order of magnitude larger in radius than it should be and the disk is a much too high a rate for quiescence and has to be moved further away to match the observed flux.  The shape of the combined spectrum is roughly correct, though a bit too bright at 1100 A.

Disk 2:  Restrict mdot to 16.2 in the log or less

WD:  1.495e-22, 17,080 K
disk:  4.39, 15.93
E(B-V) = 0.002
chi-sq=5.3

I didn't record the chi-sq above, but the main point here is that a similar fit can be obtained with more realisitic parameters.  The disk mdot is very similar to that found by Winter & Sion, who obtained 5e-11 Msolar/yr, or 15.5 g/s in the log.  Our fit isn't really like theirs, though, as the WD here contributes as much flux in the IUE range as the disk, while it  contributes <8% in their models.  The WD is also too big here again.

The last thing to try is to synthesize the WD and disk spectra together so that they are at the same distance.  I've re-run cv with the disk at 20,000K.  The original, "cold" WD models with T=10,000K are now in subdirectory wd10k.

Disk 3:  WD temperature set to 20,000 K.  Model:  N=12.65, log mdot=15.00 (the minimum), chi-sq=5.641.  This model looks like the disk-only fit above.  When the WD and disk are treated self-consistently, a cool WD is not a significant contributer to the UV flux.  This fit is acceptable.  It's a bit too bright around 1100 A, but otherwise is a good rough approximation of the shape of the UV spectrum.  In fact, just multiplying the normalization by about 2/3 may bring the model in better line with the spectral flux.

Disk 4:  WD temp set to 30,000 K.  Model:  N=4.619, mdot = 15.4, chi-sq=7.38.  Too blue with a hot WD.

Basically, the disk-only fit is the best;  adding a WD just creates excess blue flux above that observed. I've added some lower m-dot models to the spectra with the 10,000 K WD to see what the best fit is when mdot isn't bounded by the limits of the models.

Disk 5:  WD 10,000 K (basically off).  Model:  N=46.6, log mdot=14.2, chisq = 5.406.  This is basically the best fit--a low accretion rate disk and a cool WD.  Screening out the strongest emission lines has a minor effect on the fit (slightly lower mdot, slightly higher N, basically the same spectrum).  I also checked on using a 20,000 K WD but the fit is still worse in that case.

This fit still has a large normalization, indicating that a disk with that mdot doesn't provide enough flux at 350 pc to match the observed flux.

Disk 6:  WD 10,000 K plus disk fit only to the FUSE range (HI 972 airglow, CIII 977 and Ly beta airglow screened out):  Model:  N=0.27, mdot=16.45, chisq=0.55 (479 pts).  This fit is a reasonable fit to the shape of the FUV continuum, though it looks too high for <1000 A, as the peaks in those areas are probably lines.  The normalization isn't bad, either, though the disk is too close and the accretion rate a bit high for quiescence.  The IUE data would have to be fainter by a factor of 2 - 3 to match the model.

Disk 7:  Same as above, but I restricted possible disk models to 15.2 or less.  Model:  N=5.84, log mdot=15.2 (the max), chisq=0.75.  The normalization is about right.  The lower accretion rate models drop off below 1000 A faster than the observed spectrum, though.  The IUE spectrum is only 1-2x brighter than the model in this case.

06/28/04

High time resolution light curves of EM Cyg. 

I used tprint to get the time-tag photon list of each channel:  these files are called c01001010001attagfraw.dat, etc.  They have time (from the start of the first exposure), x, y (positions of the events), and pulse heights for each photon event.  I edited out by hand all events at the start of each exposure, because they often cover more or less than a second of time and we found in the UX UMa data that they had biased count rates as a result. 

To decide on which photon events to use, I looked at the raw images in fuse_scan:

For LiF1A, approximately 990 -- 1080 A, excluding Ly beta airglow at 1025:  1460 <= x <= 6640 and 6800 <= x <= 15200; 520 <= y <= 620 for the spectrum and 400 <= y <= 500 for the background sample (away from all apertures)

For SiC1A, approximately 1000 -- 1090 A, excluding Ly beta airglow:   1000 <= x <= 11500 and 11650 <= x <= 15200; 60 <= y <= 160 for the spectrum and 400 <= y <= 500 for the b.g.

For LiF1B, approximately 1090 -- 1190 A but worm renders all data >1135 A unreliable:  1065 <= x <= 7200;  500 <= y <= 600 object, 400 <= y <= 500 b.g.

For SiC1B, 900 - 990 A:  1100 <= x <= 15175; 50 <= y <= 150 object, 400 <= y <= 500 b.g.

For LiF2A, 1090 - 1180:  1000 <= x <= 15300;  650 <= y <= 750 object, 100 <= y <= 200 b.g.

For SiC2A, 920 - 1010:  1500 <= x <= 16000; 350 <= y <= 450 object, 100 <= y <= 200 b.g.

For LiF2B, 980 - 1080 and excluding Ly beta airglow:  1150 <= x <= 8350 and 8500 <= x <= 14500; 650 <= y <= 750 object, 800 <= y <= 900 b.g.

For SiC2B, 1020 - 1100 and excluding Ly beta:  1200 <= x<= 2550 and 2650 <= x <= 15400; 450 <= y <= 550 object, 800 <= y <= 900 b.g.

Let's start with LiF1A, as this is the guide channel for FUSE and is the most reliable in relative count rates as a result.  I wrote a quick script, mklc.pl, to add up all events in the spectrum window in each 1 sec time bin and subtract all events in the same time bin from the background window.  I also screened for Ly beta airglow and for pulse heights between 2 and 31 only.  The output file is lif1a_lc.dat. 

I also made up LiF1B and sent both off to Bill for power spectra analysis.