TRANSIT AND ECLIPSE ANALYSES OF THE EXOPLANET HD 149026b USING BLISS MAPPING

We observed secondary eclipses of HD 149026b at 3.6, 4.5, 5.8, and 8.0 μm with the Infrared Array Camera (IRAC; Fazio et al. 2004) and at 16 μm using the Infrared Spectrograph's (IRS; Houck et al. 2004) photometric blue peak-up array. The program also observed a primary transit at 8.0 μm. Including the four previously analyzed data sets labeled in Table 1, we present 14 observations spanning more than 3.5 years.

Our Photometry for Orbits, Eclipses, and Transits (POET) pipeline produces systematics-corrected light curves using Spitzer-supplied basic calibrated data, fits a multitude of models with a wide range of analytic forms for systematic effects, chooses the best-fit model, and assesses the uncertainty of each free parameter. Below, we describe each of these steps in detail.

We calculate the Julian date of each image at mid-exposure using the UTCS_OBS and FRAMTIME keywords in the Spitzer-supplied headers. Following Eastman et al. (2010), we convert dates to Barycentric Julian Dates in the Coordinated Universal Time standard (BJD_UTC) using the JPL Horizons system⁴ to interpolate Spitzer's position relative to our solar system's barycenter. Additionally, converting from UTC to the Barycentric Dynamical Time (TDB) standard addresses any discontinuities due to leap seconds. This paper reports both time standards to facilitate comparisons with previous work, which mostly does not apply the leap-second correction, and to ease the transition to the more-accurate standard.

The POET pipeline flags bad pixels (from energetic particle hits and other causes) by grouping sets of 64 frames and doing a two-iteration, 4σ rejection at each pixel location in the set. Stellar centers for photometry come from a Gaussian fit and 5× interpolated aperture photometry (H07 Supplementary Information, SI) produces the light curves. We test a broad range of aperture sizes in 0.25 pixel increments and omit frames with bad pixels in the photometry aperture. The background, subtracted before photometry, is an average of the good pixels within the specified annulus centered on the star in each frame.

Exoplanet characterization requires photometric stability well beyond Spitzer's design criteria (Fazio et al. 2004). Detector sensitivity models vary by channel and can have both temporal (detector ramp) and spatial (intrapixel variability) components. The main systematic effect at 3.6 and 4.5 μm is intrapixel sensitivity variations (Charbonneau et al. 2005), in which the photometry depends on the precise location of the stellar center within its pixel. We fit this systematic using the new BiLinearly-Interpolated Subpixel Sensitivity (BLISS) mapping technique described in Section 3. This technique maps the spatial sensitivity variations at high resolution.

The 5.8, 8.0, and 16 μm arrays primarily suffer from temporal variability, attributed to charge trapping (K09) in the 8.0 μm case. Weak spatial dependencies can also occur at these wavelengths (Stevenson et al. 2010), so we consider both systematics when determining our best-fit model. Typically, we omit the initial portion of each light curve from the model fit to avoid the worst of the ramp effect and to allow the telescope pointing and detector to stabilize. The clipping parameter, q, defines the number of unmodeled points from the start of a data set.

Pixelation is an infrequently discussed systematic error inherent to all array detectors. Sufficiently small stellar center motions between frames will not add or subtract pixels from an aperture, but these motions will cause the total flux within the aperture to vary. This means there is a (potentially large) range of stellar centers that utilizes the same set of aperture pixels, introducing a position dependence to the photometric sensitivity. We provide an illustrative example in Figure 1 and display the magnitude of this effect in Figure 2. A flux-conserving, subpixel image interpolation, combined with precise centering and applied before photometry, mitigates the pixelation effect by decreasing its range and amplitude. As demonstrated by our BLISS maps in Section 3, uninterpolated photometry exhibits strong sensitivity peaks at 1 pixel increments, while 2×- and 5×-interpolated pixels exhibit progressively weaker peaks at 0.5 and 0.2 pixel increments, respectively, in each spatial direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Pixelation is most apparent with small apertures placed on underresolved point-response functions (PRFs) such as Spitzer's. It may not be apparent in other situations, such as when the aperture contains almost all of the integrated PRF, when the centroid wander is small relative to the subpixel size, and when other noise sources dominate (e.g., systematics and variable PRFs). Increasing the aperture size lessens the pixelation effect by decreasing the fraction of uncaptured light outside the aperture, but doing so may decrease the signal-to-noise ratio (S/N) by increasing the amount of background noise included in the aperture. Thus, choosing the best aperture size may introduce this position-dependent systematic. We have determined the intrapixel variability at 5.8, 8.0, and 16 μm to be a pixelation effect, previously reported at 5.8 μm by Stevenson et al. (2010) and at 8.0 μm by Anderson et al. (2011).

There are several ways to correct pixelation. First, one could shift model PRFs to match each frame's precisely determined stellar center. Dividing the stellar flux by the PRF flux in the aperture should remove the effect, but requires a highly accurate model PRF. Second, using smaller subpixels could decrease the amplitude of the pixelation effect until it is insignificant relative to the noise, but there is a limit: interpolation can only approximate the information destroyed when photons fall into the detector's finite-sized pixels. Third, if the pointing is sufficiently consistent and compact, one could choose an interpolation factor that happens to place the flat portion of the pixelation response on the stellar centers (since the peaks in Figure 2 move with different interpolation factors). Fourth, with high-precision centering, one could use a series of images taken at slightly different positions to model the position sensitivity analytically or with pixel-mapping techniques such as BLISS, but one would first have to remove any time-dependent components from the light-curve model (see Section 2.5). The accuracy of these models depends on the centering and photometric precisions, the former being ∼0.01 pixels for 0.4 s IRAC subarray exposures of bright sources (see below and Stevenson et al. 2010). We test for pixelation in all data sets using BLISS mapping, which corrects the effect when it is significant.

The full light-curve model is

$$\begin{eqnarray} &&F(x, y, t) = F_{\rm s} E(t)R(t)M(x,y)V(\upsilon)P(p), \end{eqnarray} \tag{ 1 }$$

where F(x, y, t) is the measured flux centered at position (x, y) on the detector at time t, F_s is the (constant) system flux outside of secondary eclipse or primary transit, E(t) is the primary-transit or secondary-eclipse model component, R(t) is the time-dependent ramp model component, M(x, y) is the position-dependent intrapixel model component or sensitivity map, V(υ) is the visit sensitivity as a function of visit frame number υ, and P(p) is the flat-field correction at position p. Below, we discuss some of these components in more detail.

2.5.1. Eclipse and Transit Models

2.5.2. Ramp Models

We consider a multitude of ramp equations, R(t), all of which stem from three basic forms: exponential, logarithmic, and/or polynomial.

$$\begin{equation} R(t) = 1 {\pm } e^{-r_{0} t + r_{1} } \end{equation} \tag{ 2 }$$

$$\begin{equation} R(t) = 1 {\pm } e^{-r_{0} t + r_{1} } + r_{2} (t-0.5) \end{equation} \tag{ 3 }$$

$$\begin{equation} R(t) = 1 {\pm } e^{-r_{0} t + r_{1} } + r_{2} (t-0.5) + r_{3} (t-0.5)^{2} \end{equation} \tag{ 4 }$$

$$\begin{equation} R(t) = 1 {\pm } e^{-r_{0} t + r_{1} } {\pm } e^{-r_{4} t + r_{5} } \end{equation} \tag{ 5 }$$

$$\begin{equation} R(t) = 1 + r_{1} (t-0.5) + r_{6} \ln (t-t_{0}) \end{equation} \tag{ 6 }$$

$$\begin{equation} R(t) = 1 + r_{1} (t-0.5) + r_{2} (t-0.5)^{2} + r_{6} \ln (t-t_{0}) \end{equation} \tag{ 7 }$$

$$\begin{equation} R(t) = 1 + r_{6} \ln (t-t_{0}) + r_{7} [\ln (t-t_{0})]^{2} \end{equation} \tag{ 8 }$$

$$\begin{eqnarray} R(t) &=& 1 + r_{1} (t-0.5) + r_{2} (t-0.5)^{2}\nonumber\\ && + r_{6} \ln (t-t_{0})+ r_{7} [\ln (t-t_{0})]^{2} \end{eqnarray} \tag{ 9 }$$

$$\begin{equation} R(t) = 1 + r_{1} (t-0.5) \end{equation} \tag{ 10 }$$

$$\begin{equation} R(t) = 1 + r_{1} (t-0.5) + r_{2} (t-0.5)^{2}. \end{equation} \tag{ 11 }$$

The time, t, is in units of phase (days) for secondary-eclipse (primary-transit) events. We use "+" and "−" subscripts in Equations (2)–(5) to denote the corresponding functional form. For example, Equation (2₊) describes an exponentially decreasing, asymptotically constant ramp while Equation (2_–) describes an exponentially increasing, asymptotically constant ramp.

There is a physical interpretation applicable to the rising exponential ramps (Agol et al. 2010). Consider a population of charge traps due to an impurity in the detector's infrared material and a flux of photoelectrons through the material. The traps collect some fraction of electrons, releasing them randomly with some characteristic timescale that depends on the impurity. Bright sources saturate the traps, decreasing the fraction of captured electrons and raising the detected signal of a steady source according to the asymptotic rising exponential function in Equation (2_–). A double rising exponential (Equation (5_–)) approximates a rapidly saturating point spread function (PSF) core and slowly saturating wings (Knutson et al. 2008). It could also represent two impurities; the very short ramps of the HD149bs41 data set's visit sensitivity suggests the presence of an impurity that has a very short timescale and releases many of its electrons over the course of one cycle (∼30 minutes; see H07 Supplementary Figure 6). We tested for a common set of characteristic timescales in all of the rising exponential free parameters; however, even for the same planet in the same array, we did not achieve reasonable fits with all data sets. This may be due to inadvertent and inconsistent preflashing by the objects observed prior to our own observations. Despite its potential for providing a physical explanation for the ramps, the rising exponential does not always provide the best model according to the criteria in Section 2.6. This may be due to our fitting the final photometry and not the individual pixels' data and/or to the pointing instability producing unsteady illumination at the precision levels relevant here.

Zoom In Zoom Out Reset image size

Agol et al. (2010) advocate using a double rising exponential (Equation (5_–)) for 8.0 μm data due to its improved fit and smaller residuals, weaker dependence on aperture size, and less sensitivity in the clipping parameter. Similarly, Knutson et al. (2011) find that Equation (2_–) is sufficient for their 8.0 μm preflashed data sets, while Equation (5_–) is necessary for their non-preflashed data. We test all relevant ramp equations on data sets that exhibit time-dependent systematics and orthogonalize any correlated parameters that inhibit convergence (see Sections 2.5.3 and 2.7). Upon doing so, we find that we cannot corroborate the claims by Agol et al. (2010). Equation (5_–) should not be used for all 8.0 μm data sets because we typically find correlations between the eclipse depth and its ramp parameters that can double the latter's uncertainty relative to other models (see HD149bp42). Rather, we recommend a comprehensive examination of all relevant ramp equations before selecting the final model. See Sections 4 and 5 for discussion relevant to particular events.

2.5.3. Orthogonalization

The exponential model components have one difficulty: the Markov chain method used to assess the uncertainties does not converge to the posterior probability distribution, according to the criteria discussed in Section 2.7, even after tens of millions of iterations. The problem is a correlation between the exponential ramp parameters. Several recent analyses (e.g., K09; Stevenson et al. 2010; Campo et al. 2011) have "solved" the problem by fixing one of the exponential parameters. This is effectively profiling (slicing) rather than marginalizing (integrating) over the posterior parameter distribution, an approach that ignores correlated errors and may reduce the calculated error bar (see Press et al. 2007, Figure 15.6.3 and related text). In one case discussed below (HD149bp42), fixing parameters incorrectly decreased the calculated eclipse-depth uncertainty by 50%. There are two legitimate escapes from the problem. The first is to write a Monte Carlo sampler that explores the phase space more intelligently than the Metropolis random walk we and many others use. The second is to re-cast the equations in a form that eliminates the nonlinear correlation of the ramp parameters.

For this paper, we have expediently chosen the second approach. The version of Equation (2_–) presented here produces a more linear correlation between the r₀ and r₁ parameters than the original version in H07, whose exponent was −r(t − t₀). In some cases this converges quickly and the job is done. In others it still does not converge, so we run at least 10⁵ iterations to sample the posterior distribution and then rewrite the model with a change of variables that orthogonalizes the most correlated parameters. This method does not modify the number of free parameters, nor involve any interpolations or approximations. For the selected correlated parameters, the orthogonalization shifts the origin to the center of the distribution, divides each parameter by its standard deviation to give a uniform scale in all directions, and rotates the subspace to minimize correlations (see Figure 3). A routine that prepares for principal components analysis (PCA) finds the transformation matrix from the distribution sample (we do not actually perform the PCA). We rerun the Markov chain using the new equations until it converges according to the criteria given below. Then we transform the points back to the familiar equations to assess parameter uncertainties. A simple example of a two-parameter orthogonalization of Equation (2) is

$$\begin{equation} R(t) = 1 {\pm } e^{-c_{0} [t\cos (\theta)-\sin (\theta)]} e^{c_{1} [t\sin (\theta)+\cos (\theta)]} , \end{equation} \tag{ 12 }$$

where c₀ and c₁ are the ramp parameters in the rotated frame and θ is the rotation angle between coordinate systems. The arctangent of the slope of the best-fit line through the initial sample, projected into the r₀ − r₁ plane, determines θ. Such approximate orthogonalization of the posterior distribution has long been standard practice for improving Markov chain Monte Carlo (MCMC) performance (e.g., Hills & Smith 1992; Brooks 1998). Similar discussion can be found in Connolly et al. (1995), Pál (2008), and Cowan et al. (2009).

Zoom In Zoom Out Reset image size

In effect, this method is the same as the first, accomplished by rotating the data and using the original sampler rather than writing a new one. This method works best for linearly correlated parameters and achieves moderate success with more exotic correlations. Converting to curvilinear coordinates or implementing manifold learning algorithms from nonlinear dimensionality reduction may offer further improvements in convergence time for the extreme cases (Lee & Verleysen 2007).

2.5.4. Flat-field (Position) Sensitivity Models

2.5.5. Visit-sensitivity Model

With HD149bs41, Spitzer completed 12 cycles through the 9 nod positions mentioned above for a total of 108 visits. Each visit has a briefer and steeper ramp compared to the overall ramp, R(t). As discussed in their SI, H07 use a 12-knot spline to model the visit-sensitivity effect, V(υ). They fix the final three knots to unity and allow the remaining nine parameters to vary. We model the visit-sensitivity ramp using

$$\begin{equation} V(\upsilon) = \upsilon _{1} \cdot \ln (\upsilon -\upsilon _{0}) + 1, \end{equation} \tag{ 13 }$$

which is identical in form to the model component used by K09. The only difference is that K09 use time as the independent variable while we use visit frame number, υ. In either case, the independent variable resets to zero upon moving to a new nod position.

For each data set, we test different photometry aperture sizes, detector ramp model components, and intrapixel sensitivity model components looking for the best combination. One must be careful in assessing which model is the "best" because χ²-like comparisons must derive from the same data set and different photometric apertures produce different data sets for this purpose. We use the standard deviation of the normalized residuals (SDNR) when comparing models of the same analytic form to different data sets. Once we have identified the best aperture size, we use the Bayesian information criterion (BIC; Schwarz 1978; Liddle 2007; Stevenson et al. 2010; Campo et al. 2011) to compare models with different numbers of free parameters:

$$\begin{equation} {\rm BIC} = \sum \limits _{1<i<n} \frac{\epsilon _{i} ^{2} }{\sigma _{i} ^{2} } + k\ln n. \end{equation} \tag{ 14 }$$

Here, _i and σ_i are the residual and uncertainty of the ith data point, k is the total number of free parameters, and n is the number of data points used in the model fit. The Spitzer-supplied uncertainty frames are typically overestimated, so we scale σ_i such that the reduced χ², χ²_ν, equals unity for our best-fit model. Using the unscaled σ_i values in a least-squares minimizer improperly weights the data and results in a sub-optimal fit. When selecting between competing models for an observed data set, what matters is how the models compare in predicting the actually observed data, not how variable individual model fits may be among hypothetical realizations of the data. The BIC is defined so that the marginal likelihood ratio—the ratio of predictive probabilities for the observed data—is approximately e^0.5ΔBIC. Note that what matters is the difference in BIC values (and thus the difference of maximum likelihood or minimum χ² values), not the absolute values. Information criteria such as BIC may not apply when comparing fits with intrapixel maps to those with polynomial model components. See Appendix A for more discussion on the matter.

We explore phase space and estimate uncertainties using an MCMC routine following the Metropolis random walk algorithm. See Campo et al. (2011) for more discussion on MCMC. The POET pipeline can test any combination of systematic model components, computing the SDNR and BIC for each fit, and can model multiple events simultaneously while sharing parameters such as the eclipse midpoint and depth. Each MCMC run begins with a least-squares minimization, a rescaling of the Spitzer-supplied uncertainties, a second least-squares minimization using the new uncertainties and, finally, at least 10⁵ MCMC iterations to populate the posterior parameter distributions, from which we derive parameter uncertainties. We study parameter correlation plots and the posterior distribution to ensure that we have a reliable result, then publish them so that others may evaluate our work and compare to their own. We test for convergence every 10⁵ steps, terminating only when the Gelman & Rubin (1992) diagnostic for all free parameters has dropped to within 1% of unity using all four quarters of the chain. We also examine trace and autocorrelation plots of each parameter to confirm convergence visually. We estimate the effective sample size (ESS; Kass et al. 1998) and autocorrelation time, τ, for the ith free parameter as follows:

$$\begin{equation} {\rm ESS_{\it i} } = \frac{m}{\tau _{i} } = \frac{m}{1 + \sum \limits _{k=1} ^{k^{\prime } } \rho _{k} (\theta _{i})}, \end{equation} \tag{ 15 }$$

where m is the length of the MCMC chain, ρ_k(θ_i) is the autocorrelation of lag k for the free parameter θ_i, and k' is the cutoff point such that ρ_k < 0.01. We use the longest autocorrelation time from each event to determine the number of steps between effectively independent samples for thinning each MCMC chain. The distribution of thinned samples quantifies parameter uncertainties.

The change in pixel sensitivity with respect to stellar position on the detector is a well-known systematic with the Spitzer Space Telescope (Charbonneau et al. 2005; Knutson et al. 2008). This effect is particularly strong in IRAC's 3.6 and 4.5 μm channels but has also been seen at 5.8, 8.0, and 16 μm (Stevenson et al. 2010; Anderson et al. 2011). The position sensitivity in the latter channels is due to a pixelation effect (see Section 2 for a description) rather than an actual change in sensitivity over the pixel surface. In this section, we present a new technique for modeling these position-dependent systematics, called BLISS mapping.

The most common method for removing the intrapixel variability is to fit a polynomial in both spatial directions (Knutson et al. 2008). The polynomial order typically ranges from quadratic to sextic and may include cross terms. Other variations of this modeling method have applied multiple polynomials, one for each pixel quadrant, in order to find the best fit. Polynomial methods work reasonably well for data sets with small stellar position wander, resulting in a smoothly varying intrapixel sensitivity. An analysis becomes exceedingly complicated if the variation is not smooth or if strong correlations arise between parameters. These complications can increase uncertainty estimates in the best case scenario and, in the worst case, lead to incorrect results.

A new approach, pioneered by Ballard et al. (2010, hereafter B10), attempts to map the intrapixel variability on a subpixel-scale grid without assuming a specific functional form. For their particular light curve, they bin their flux and stellar positions into 20 s bins (∼145 points bin⁻¹) before computing a sensitivity correction for each binned point. Each correction considers a set of flux values that does not include in-transit frames or frames from the current binned position. This set of flux values is Gaussian-weighted in both spatial directions relative to the position of the binned point being corrected, summed to a single value, then normalized by dividing by the summed Gaussian weighting function. This method effectively models not only the large-scale features in their data, but also a smaller-scale "corrugation" effect that a low-order polynomial cannot remove.

Moving away from a polynomial model is an excellent concept; however, this particular implementation has some drawbacks that limit its scope. For instance, ignoring the points during secondary eclipse requires that the out-of-eclipse portion of the data set be significantly longer than the in-eclipse portion, which is atypical in primary-transit and secondary-eclipse observations. Also, the weighting function computes too slowly to be used in an MCMC routine and is even slow when using a minimizer. The calculation would be even slower if the data were not binned into relatively long, 20 s time intervals. Figure 4 shows 120 s of vertical pixel positions from HD149bs11, which has a 0.4 s exposure duration versus 0.1 s for B10. The stellar center can vary significantly over a 20 s interval, indicating that positional information is lost when binning over such timescales.

Zoom In Zoom Out Reset image size

The idea of using a spline to interpolate position-dependent systematics stems from observed fine-scale sensitivity variations in some of our data sets that cannot be modeled by a low-order polynomial (see the pixelation effect in Figures 2 and 5 for examples). We initially attempted to map the pixel surface using a bicubic spline because we wanted a smoothly varying model; however, this type of interpolation is prohibitively computationally expensive. A typical secondary-eclipse observation spans a 0.3 × 0.3 pixel region. Placing knots at 0.05 pixel intervals requires 49 free parameters. Polynomial models typically require less than 10 parameters. Adequately describing the fine-scale sensitivity variations requires a large number of knots, but varying all of these knot parameters at each step of an MCMC routine leads to extremely slow convergence. Our new BLISS mapping technique circumvents the problem of slow convergence by directly computing the knot parameters, rather then allowing them to vary freely. Thus, we can use >1000 knots to map the pixel surface at high resolution (see Figure 5).

Zoom In Zoom Out Reset image size

In seeking methods that compute faster than bicubic interpolation, we give up the constraint of differentiability. Nearest-neighbor interpolation (NNI) is simple, assigning each point the value of its nearest knot. Bilinear interpolation (BLI) is a straightforward calculation (see Equation (16)) that maintains continuity over the pixel surface, unlike NNI, and, given sufficiently precise centering, should be more accurate. Using BLI, the flux at a point (x, y) is

$$\begin{eqnarray} M(x,y) &=& F_{\rm IP} (x_{1} ,y_{1}) \frac{x_{2} -x}{x_{2} -x_{1} } \frac{y_{2} -y}{y_{2} -y_{1} } \nonumber\\ &+& F_{\rm IP} (x_{2} ,y_{1}) \frac{x-x_{1} }{x_{2} -x_{1} } \frac{y_{2} -y}{y_{2} -y_{1} } \nonumber\\ &+& F_{\rm IP} (x_{1} ,y_{2}) \frac{x_{2} -x}{x_{2} -x_{1} } \frac{y-y_{1} }{y_{2} -y_{1} } \nonumber \\ &+& F_{\rm IP} (x_{2} ,y_{2}) \frac{x-x_{1} }{x_{2} -x_{1} } \frac{y-y_{1} }{y_{2} -y_{1} }. \end{eqnarray} \tag{ 16 }$$

This is a distance-weighted average of the flux of the four nearest knots, F_IP(x_i, y_j), where i and j are horizontal and vertical indices for a rectangular grid of knots. This method computes faster than bicubic interpolation and may achieve comparable smoothness within the errors with less computing time simply by increasing the number of knots (see Figure 2).

We create a rectangular grid of knots that spans the range of centers in x and y. Each point in the data set associates with its nearest knot. For BLI, we compute the distances from each point to its four nearest knots, for Equation (16). If one or more knots in Equation (16) does not have any assigned points, we use NNI there instead, or the calculation would fail. This usually only occurs near the boundary of the grid of knots. We precompute the knot associations and distances prior to initiating the MCMC as they remain constant from iteration to iteration.

We do not treat the knots as MCMC jump parameters. Rather, we step all other free parameters from Equation (1), generate a new model using these new jump parameters, then divide the observed flux by the new model (F_IP(x, y) = F_obs/F_sE(t)R(t)V(ν)P(p)). Hypothetically, the residuals of F_IP(x, y) contain only position-dependent flux variations. The flux value of a particular knot is the mean of F_IP(x, y) for the points associated with that knot. We also tried median and weighted average knot values but the results did not improve and these calculations are much slower. Next, we generate the sensitivity map (M[x, y], Figure 6) by interpolating the flux from the knots to all of the observed points using BLI and/or NNI. We tested various weighted smoothing functions when generating the sensitivity maps but, again, there was no improvement in the results. Finally, M(x, y) enters Equation (1) for comparison with the observed flux to obtain an estimate of the goodness of fit and determine the MCMC acceptance probability. This process repeats for each step of the MCMC routine or minimizer.

The accuracy of the BLISS mapping technique depends critically on the bin size, or resolution in position space; however, there is a tradeoff between bin size and speed. Decreasing the bin size requires more knots and runs slower but may be necessary to adequately resolve sensitivity changes on the pixel surface. There is, however, a practical limit to how small the bins can be. A bin for every measurement will always produce a perfect fit, resulting in a negative number of degrees of freedom and leaving the eclipse parameters unconstrained. Bin sizes must be small enough to resolve real, small-scale variations on the pixel surface but large enough to mix in- and out-of-eclipse points. This mixture helps to minimize correlations between the eclipse parameters and the knots in the sensitivity map.

To establish the optimal bin size, we consider a range of bin sizes for both BLI and NNI using the 3.6 μm data set (HD149bs11), which has the strongest position-dependent systematic. We draw several conclusions from the results in Table 2. First, the SDNR (our measure of goodness of fit) decreases with decreasing bin size, indicating a better fit. Unfortunately, the SDNR decreases indefinitely, so it cannot constrain the minimum bin size. Second, BLI fits the data better than NNI for bin sizes greater than 0.015 pixels. The opposite is true for smaller bin sizes, which is counterintuitive because BLI should always outperform NNI, assuming no uncertainty in the position. Thus, the bin size at which NNI outperforms BLI is indicative of the centering precision for a particular data set. We estimate the precision for our analysis of HD149bs11 to be 0.009 pixels in x and 0.007 pixels in y by calculating the root mean squared (rms) frame-to-frame position difference. This agrees well with Figure 4 and is consistent with the crossover bin size of 0.015 pixels, where NNI outperforms BLI. Last, we place an upper limit on the bin size by noting that the eclipse depths become inconsistent with each other for pixel sizes ⩾0.050 using BLI and ⩾0.020 using NNI.

We conclude that, whenever possible, BLI should be used with a bin size that is independent of the eclipse depth and has a lower SDNR than NNI. We have found cases where BLI is never better than NNI. In those instances, the position dependence is so weak that the intrapixel model component is unnecessary. A better fit with BLI, compared to NNI, is thus a good indicator that a position-dependence systematic is present in a given data set. For the test cases shown in Table 2, using BLI with a bin size between 0.015 and 0.020 pixels is recommended based on our criteria.

Precise centering is important for this method because imprecision limits the smallest meaningful bin size. Our preliminary work (see the SI of Stevenson et al. 2010) indicates that the Gaussian and least-asymmetry centering methods are better than the center of light; additional work is in preparation.

To compare the BLISS mapping technique with other intrapixel methods, we fit six different intrapixel models to the HD149bs11 data set. These models are quadratic, cubic, and sextic polynomials (including lower-order cross terms), B10's new weighted sensitivity function (fixing σ_x to 0.021 and σ_y to 0.0079), BLI, and NNI. The eclipse depth in B10's model is slightly shallower than the other models, which are all well within 1σ of one another, but the uncertainties are essentially identical (see Table 3). The BLISS models show significant improvement in SDNR compared to the others. We cannot use the BIC to compare the BLISS model components with the polynomial model components for the reasons discussed in Appendix A.

BLISS mapping represents a substantial improvements over polynomial model components because BLI and NNI can model real structure (such as pixelation) that cannot be modeled with low-order polynomials, they encounter fewer correlations between free parameters, and require fewer iterations to assess parameter uncertainties.

At each 0.25 pixel increment in photometry aperture size, we model the time-dependent systematic with the ramps listed in Equations (2)–(11). An aperture size of 3.5 pixels produces the lowest SDNR values for all ramps and for all transit events. We estimate the background flux using an annulus from 7 to 15 pixels centered on the star. We follow the method described in Section 3.3 when determining the optimal bin sizes of the BLISS maps.

4.1.1. HD149bp41

4.1.2. HD149bp42

4.1.3. HD149bp43

We perform two joint-model fits, each requiring less than 2 × 10⁶ iterations to estimate uncertainties. The first considers only the three transits analyzed here while the second also considers the more precise NICMOS data from Carter et al. (2009) by placing priors on i and a/R_⋆. Both fits in Table 5 are consistent with previous results from Knutson et al. (2010; R_p/R_⋆ = 0.0522 ± 0.0008) and Carter et al. (2009; R_p/R_⋆ = 0.0519 ± 0.0008) and have improved estimates of the radius ratio. The uncertainties in the duration and ingress/egress times for the independent fit are significantly larger than those from the fit with Carter et al. (2009) priors.

For this data set, we find that BLI outperforms NNI down to a bin size of 0.015 pixels when we exclude bins with less than four measurements. Bins with fewer data points are insufficiently sampled to compute a reliable mean flux for the knot value. The linear ramp Equation (10) fits best. The posterior parameter distributions (histograms) and parameter correlation plots (including knot values in the BLISS map) are in Figures 11–13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (13 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (13 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Using the strategy described in Section 3.3, BLI achieves a better fit than NNI with bin sizes of 0.012 pixels along x and 0.006 pixels along y. Additionally, we fit the intrapixel sensitivity with three different polynomial models ranging between second and sixth order. After clipping the first 5000 points (q = 5000), the eclipse depths using various ramp and intrapixel model components are consistent (see Table 6); however, the SDNR clearly favors BLI and the BIC favors Equation (2₊) to model the systematics. To minimize the convergence time in our MCMC chains, we orthogonalize the eclipse depth, system flux, and both ramp parameters (r₀ and r₁). All of the model fits exhibit, to various degrees, bimodal distributions in the eclipse-midpoint histograms. The lesser peak occurs at a phase of 0.497 and is likely a result of the model trying to fit the eclipse egress to the points from phases of 0.514 to 0.520, which are consistently above the secondary eclipse by 1σ–2σ (see Figure 10, center panel).

Using the processes described below, we choose the best-fit model for each event, then perform a joint fit with a single eclipse-depth parameter.

5.3.1. HD149bs31

5.3.2. HD149bs32

5.3.3. HD149bs33

5.3.4. 5.8 μm Joint Fit

To improve S/N on the eclipse depth, we share this parameter in a joint fit of all three data sets. We retain individual eclipse-midpoint times for the subsequent orbital analysis. The MCMC chain converged after 6 × 10⁵ iterations. The combined light curve in Figure 14 illustrates the improvement when compared to the three 5.8 μm light curves in Figure 10. The best simultaneous fit favors an eclipse depth of 0.044% ± 0.010%.

Zoom In Zoom Out Reset image size

Similarly to the fits at 5.8 μm, we fit models to each data set individually then use the best models in a 1.1 × 10⁶ iteration joint fit that shares a common eclipse depth.

5.4.1. HD149bs41

The significant improvements in our pipeline since the original analysis by H07 warrant a new analysis of HD149bs41. We follow all of our current techniques described in Section 2 and test all of the listed ramp model components. As with H07 and K09, we find that Equation (2_–) best describes the overall ramp. The smaller ramps associated with each telescope movement are best described by Equation (13), according to BIC; however, we also present H07's 12-point spline for comparison (see Table 7). Each model employs a constant-flux offset at each of the nine nod positions, of which eight are free parameters as described in Section 2.

Table 7. HD149bs41—Comparing Model Fits

R(t)	M(x, y)	Visit Sensitivity	SDNR	ΔBIC	Ecl. Depth
					(%)
2–	...	Equation (13)	0.0084575	0	0.049
2–	...	12-pt Spline	0.0084559	59	0.049
2–	9 ×Linear	Equation (13)	0.0084510	126	0.050
2–	9 ×Linear	12-pt Spline	0.0084494	184	0.050

Note. Bold indicates the model selected for the joint fit.

Download table as:  ASCIITypeset image

Due to the nodding motion with this particular data set, BLI and NNI are inappropriate models to use. We can see in Figure 15, which illustrates 1 of the 9 nod positions, that the pixel position is slightly different for each of the 12 visits to this position. This behavior introduces a strong time dependence in the position sensitivity correction that cannot be disentangled. Our best attempt to correct for the position sensitivity uses a linear correction in two dimensions for each of the nine nod positions (9 ×Linear). Table 7 compares the four best model combinations. Compared to K09, our flux offsets are multiplicative rather than additive (see Equation (1)) and our final model does not fix either of the ramp parameters. As a result, our eclipse-depth uncertainty is larger (0.049% ± 0.016%).

Zoom In Zoom Out Reset image size

5.4.2. HD149bs42

5.4.3. HD149bs43

5.4.4. HD149bs44

Each event consists of 1050 exposures, divided serially into three 350-image sequences at non-overlapping stellar centers on the detector. The second sequence (p = 1) is completely within the secondary eclipse and, because of the free position offset parameter, has no impact on the eclipse depth (but still helps to model the systematics). Unfortunately, this means that the eclipse depth is completely determined by the small number of points after ingress in the first sequence and before egress in the last sequence.

To avoid residual effects from potentially bright previous targets, neither of the two 16 μm observations was positioned at the center of the array. As a result, the outer radius of the sky annulus is limited to 11 and 12 pixels for HD149bs51 and HD149bs52, respectively, to avoid the flux falloff at the edges of the blue peak-up array. We apply both aperture and optimal photometry (Horne 1986; Deming et al. 2005); the former produces cleaner results for both data sets.

Figure 16 displays the best-fit eclipse depths and corresponding SDNR values versus aperture size for four competing models in each data set. The first applies no intrapixel correction, the second uses three linear components (one at each position), the third uses BLI, and the final model uses NNI. The first two models also apply a position offset to account for the differences in pixel sensitivity. This offset is redundant for BLI and NNI.

Zoom In Zoom Out Reset image size

Although these models find similar (<1σ) eclipse depths at the best aperture size of 2.0 pixels, this is not the case for other aperture sizes. Models using BLI or NNI produce the most consistent eclipse depths and result in lower SDNR values (see Tables 10 and 11). We test a range of bin sizes for BLISS mapping but find that NNI consistently outperforms BLI. We take this as evidence that there are insufficient data points in most bins to model the weak position dependence accurately. As such, we consider only the last two models and select no intrapixel model for HD149bs51 and the 3  × linear model for HD149bs52 to perform the joint fit. For improved convergence in our MCMC chains, we orthogonalize the system fluxes and both ramp terms in Equations (11) and (2_–) for HD149bs51 and HD149bs52, respectively. Due to the weak eclipse signal, the midpoint is fixed to a phase of 0.5 for the individual fits; the joint fit achieves a sufficient S/N to share the eclipse midpoint as a free parameter. The final model fit uses 1.6 × 10⁶ iterations and is in Appendix B.

We find that the eclipse depths for both data sets vary with the choice of eclipse midpoint. Using the midpoint from our joint model (0.5015, see Appendix B), we find that the best-fit eclipse depth differs by up to 1σ from the values listed in Tables 10 and 11, where the midpoint is fixed to 0.5. We also test joint fits without the data at p = 1 and find a comparable joint-fit eclipse depth of 0.099% ± 0.035%. The small change is likely the result of weak correlations with the ramp parameters. We include all positions in the final fit as this results in a lower SDNR and a small improvement in the eclipse-depth uncertainty.

In Section 3, we show that BLISS mapping improves the SDNR of models that include an intrapixel sensitivity term. Elsewhere in this paper and in our earlier work (Campo et al. 2011) we use the BIC to compare models of other systematic error components (e.g., rival parametric ramp models), but we did not use the BIC to compare intrapixel sensitivity models. In this Appendix, we discuss the approximations underlying the BIC to elucidate when it is useful, and in particular, to explain why we do not use it (or similar statistical criteria) for comparison of intrapixel sensitivity models.

The BIC provides an asymptotic approximation to quantities that may be used for Bayesian quantification of model uncertainty. The most simple use of the BIC is to approximate Bayesian explanatory model selection. A basic distinction among model selection criteria is between explanatory criteria (that seek the model that best describes the processes that produced the available data) and predictive criteria (that seek the model that will make the most accurate predictions of future data based on the available data). For example, the well-known Akaike information criterion (AIC = χ² + 2k) is a predictive criterion. Although different from the BIC in rationale, its derivation invokes similar asymptotic approximations to those we describe here for the BIC, and we consider it to be similarly hampered for comparing large models. To illuminate the nature of the approximations underlying the BIC, we sketch its derivation here.

Consider a set of models {M_i} for observed data D. Let θ_i denote the parameters of model i, with k_i dimensions. The (Bayesian) posterior probability for model i is proportional to the product of a prior probability for the model, and the model's marginal likelihood,

$$\begin{equation} \mathcal {M}_i = \int d\theta _i\; \pi _i(\theta _i) \mathcal {L}_i(\theta _i), \end{equation} \tag{ A1 }$$

where π_i(θ_i) ≡ p(θ_i|M_i) is the prior probability density function (PDF) for the model's parameters and $\mathcal {L}_i(\theta _i) \equiv p(D|\theta _i, M_i)$ is the likelihood function, the probability for the data presuming the model holds with parameters θ_i (the likelihood function is proportional to exp [ − χ²(θ_i)/2] for our models). As its name suggests, the role of the marginal likelihood in quantifying model uncertainty is completely analogous to the role of the more familiar likelihood function in quantifying parameter uncertainty (within a particular model). Note that $\mathcal {M}_i$ is an average of the likelihood function for model i, not a maximum: in Bayesian inference, the weight of the evidence for a model is given by the typical value of the likelihood function for its parameters, not the optimum (largest) value.

For nonlinear models with more than a few dimensions, calculation of the integral in Equation (A1) is not feasible using standard quadrature methods. The development of algorithms for accurate calculation of marginal likelihoods is an active research area, and existing algorithms are typically problem-specific and computationally expensive (see Clyde et al. 2007 for a recent review targeting astronomers).

The BIC approximates $-2\ln \mathcal {M}$ for straightforward models in the limit of voluminous data, i.e., asymptotically (we drop the model index, i, when referring to a generic model, to simplify notation). In this limit, the likelihood function $\mathcal {L}(\theta)$ will be strongly peaked at the maximum likelihood estimate, $\hat{\theta }$, and it will be much more strongly concentrated than the prior; Figure 19 depicts the situation. As a first step in approximating the integral in Equation (A1), we evaluate the prior at $\hat{\theta }$ and pull it out of the integral, so

$$\begin{equation} \mathcal {M}\approx \pi (\hat{\theta }) \int d\theta \; \mathcal {L}(\theta). \end{equation} \tag{ A2 }$$

The prior PDF has dimensions of [1/θ], and we can express $\pi (\hat{\theta })$ as the inverse of a local prior scale, Δθ (for a normalized flat prior spanning a range Δθ, we have $\pi (\hat{\theta }) = 1/\Delta \theta$ exactly). We next approximate the integral of the likelihood function in Equation (A2) as the product of its height (the maximum likelihood value, $\mathcal {L}(\hat{\theta })$) and a characteristic width, δθ (describing the posterior uncertainty in θ). This gives

$$\begin{equation} \mathcal {M}\approx \mathcal {L}(\hat{\theta }) \frac{\delta \theta }{\Delta \theta }. \end{equation} \tag{ A3 }$$

With suitable definitions of the prior and posterior scales, we can make this equation exact (e.g., for a 1D model with a Gaussian likelihood function of standard deviation σ, and a flat prior with range Δθ ≫ σ, as long as $\hat{\theta }$ is not near a boundary of the prior range, then setting $\delta \theta = \sigma \sqrt{2\pi }$ ≈1.06 times the full width at half-maximum makes the equation accurate). The factor multiplying the maximum likelihood is sometimes called the Ockham factor and will be ⩽1; it quantifies how much of the parameter space of the model is wasted, in the sense of including parameter values that are ruled out by the data. Note that for dimension k > 1, these scales are volumes, i.e., products of scales in each dimension.

Zoom In Zoom Out Reset image size

Asymptotically, for a simple model of fixed dimension, we expect the uncertainty for each parameter to scale like $1/\sqrt{N}$ for sample size N. So for a model of dimension k, we expect δθ to eventually decrease proportional to N^−k/2. Let N_a be the sample size where the asymptotic behavior kicks in, and δθ_a be the typical scale of the uncertainties at that sample size. Then we expect

$$\begin{equation} \mathcal {M}\approx \mathcal {L}(\hat{\theta }) \frac{\delta \theta _a}{\Delta \theta } \left(\frac{N_a}{N}\right)^{k/2}. \end{equation} \tag{ A4 }$$

In the simple case of estimating linear parameters for data with additive Gaussian noise of fixed noise variance, we expect asymptotic behavior right away, so N_a = 1. Now take the logarithm and group terms according to their dependence on N:

$$\begin{equation} \ln \mathcal {M}\approx \ln \mathcal {L}(\hat{\theta }) - \frac{k}{2}\ln N + \ln \left(\frac{N_a^{k/2} \delta \theta _a}{\Delta \theta }\right). \end{equation} \tag{ A5 }$$

The first term may have a nontrivial N dependence; the last term is constant with respect to N. Schwarz (1978) derived a more precise expression like this, explicitly calculating the first term in the case of linear models with sampling distributions in the exponential family (which includes, e.g., normal, Poisson, and multinomial distributions). In that case the maximum log-likelihood term is expected to grow increasingly negative, roughly proportionally to N.

The BIC keeps the N-dependent terms in Equation (A5) and multiplies by −2;

$$\begin{equation} \mbox{BIC} = -2\ln \mathcal {L}(\hat{\theta }) + k\ln N = \chi ^2 + k\ln N. \end{equation} \tag{ A6 }$$

It is attributed to Schwarz (and sometimes called the Schwarz criterion), but notably, Schwarz did not drop the N-independent term in his approximation to $\ln \mathcal {M}$, although he termed it a "residual" with respect to the N-dependent terms.

If the models under consideration are considered equally probable a priori, the most probable model is the one with the largest marginal likelihood. BIC-based model selection uses the BIC to approximate the logarithm of the marginal likelihood, choosing the model with the smallest value of the BIC. The derivation sketched above provides some insight regarding when we might expect this procedure to identify the highest likelihood model. There are two main considerations.

First, the BIC is an asymptotic criterion. Its accuracy requires sample sizes large enough so that the parameter uncertainties are decreasing at the O(N^−k/2) rate. For complex models with many parameters, this is not simply a matter of sample sizes being "large enough." For some models—e.g., nonparametric models—the number of parameters may grow with sample size (explicitly or effectively); for others, some parameters may be sensitive to only a subset of the data. For example, the BLISS model uses a piecewise linear intrapixel sensitivity map, so a particular coefficient is determined by only a subset of the image pixels. In these cases, a BIC-like criterion may be valid, but with the kln N term replaced with a term that more accurately describes the asymptotic behavior of parameter uncertainties. Determining the form of such a term can be subtle (see Kass & Raftery 1995, Section 4.2). In these settings, it may take very large sample sizes to reach asymptotic behavior.

Second, the BIC drops a constant (in N) term from the logarithm of the marginal likelihood—Schwarz's "residual." That is, it drops a multiplicative factor from the estimated marginal likelihood. This factor depends on the prior volume, Δθ. For models with many similar degrees of freedom, like the various intrapixel sensitivity models, the prior volume is the product of the ranges of many variables. It can vary sensitively with the choice of a priori scale per parameter, and if the competing models have different types of parameters, the omitted residual terms may be very different from one model to another. The difference between the residual terms can be large when the models are large. As a result, the change in the BIC between two large models cannot be relied upon for identifying the model with the larger marginal likelihood.

Kass & Wasserman (1995, hereafter KW95) have examined the role of the residual term in the asymptotic approximation of the log marginal likelihood, arguing that for some problems there may be a reasonable argument for it to be negligible. Note that the last term in Equation (A5) will vanish if

$$\begin{equation} \Delta \theta = N_a^{k/2} \delta \theta _a. \end{equation} \tag{ A7 }$$

When the asymptotic sample size N_a = 1 (e.g., for linear models, like estimating the mean of a normal distribution with known variance), this requirement corresponds to having the prior range equal to the width of the likelihood for a single-sample data set. For N_a > 1, the N^k/2_a factor scales the δθ_a likelihood volume to what the single-sample volume would be if the model were asymptotic starting with N = 1. Thus KW95 dubbed a prior satisfying Equation (A7) a unit information prior, i.e., a prior that is as informative as a single datum.⁵ To the extent that one could consider the uncertainty scale associated with a single measurement to reflect prior uncertainty, the BIC may be an accurate approximation to $\ln \mathcal {M}$.

However, even when a unit information prior scale appears reasonable, for large models, even a small variation from this scale for the prior range of each parameter could produce a large net residual term. We thus do not consider the unit information prior argument to provide a sound justification for using the BIC to compare large models.

For these reasons, in our work we limit use of the BIC to comparing small models, or large models that are nested, so that rival models share the vast majority of parameters. For example, we rely on the BIC to compare different ramp models that share a common BLISS map model, but we do not consider the BIC to be valid for comparing, say, a model using a BLISS map to a model using B10's intrapixel variability correction, or to polynomial intrapixel models.

The residual issue, in part, reflects an inherent weakness of marginal-likelihood-based model comparisons with large models. Small changes in the per-parameter prior scale for such models can lead to large changes in marginal likelihoods, even with an accurate numerical calculation of the marginal likelihoods. This motivates developing a way to allow the scale to adapt to the data. In a Bayesian framework, this can be accomplished using a hierarchical model to implement regularization. One considers the prior range to be an adjustable parameter itself that is learned from the data. This can be done in a manner that essentially lets the data determine the effective number of free parameters (the total number of knots) required for the intrapixel map. We would then compare our best fit using BLISS mapping to those obtained using other intrapixel models and select the appropriate model component. Such calculations are beyond the scope of this paper, but would be a productive avenue for future research.

Table 14 displays the best-fit parameter values and uncertainties from our joint transit light-curve fit. Tables 15 and 16 display the best-fit parameter values and uncertainties from our joint eclipse light-curve fits. The tables also contain information to help evaluate the goodness of fit.

Table 14. Best-fit Joint Transit Light-curve Parameters

Parameter	HD149bp41	HD149bp42	HD149bp43
Wavelength (μm)	8.0	8.0	8.0
Array position ($\bar{x}$, pixel)	14.93	14.96	15.14
Array position ($\bar{y}$, pixel)	15.14	14.56	14.47
Position consistencya (δx, pixel)	0.013	0.011	0.011
Position consistencya (δy, pixel)	0.012	0.013	0.013
Aperture size (pixel)	3.5	3.5	3.5
Sky annulus inner radius (pixel)	7.0	7.0	7.0
Sky annulus outer radius (pixel)	15.0	15.0	15.0
System flux Fs (μJy)	117229 ± 14	117460 ± 60	117356 ± 13
Transit midpointb (MJDUTC)	4327.3720 ± 0.0005	4356.1315 ± 0.0005	4597.7067 ± 0.0005
Transit midpointb (MJDTDB)	4327.3727 ± 0.0005	4356.1323 ± 0.0005	4597.7075 ± 0.0005
Rp/R⋆	0.0518 ± 0.0006	0.0518 ± 0.0006	0.0518 ± 0.0006
cos i	0.095 ± 0.009	0.095 ± 0.009	0.095 ± 0.009
a/R⋆	5.98 ± 0.17	5.98 ± 0.17	5.98 ± 0.17
Ramp equation (R(t))	3–	8	2–
Ramp, r0	24 ± 4	0	18 ± 2
Ramp, r1	−0.6 ± 0.7	0	4.1 ± 1.4
Ramp, r2	0.0110 ± 0.0015	0	0
Ramp, r6	0	0.0009 ± 0.0005	0
Ramp, r7	0	−0.00048 ± 0.00011	0
Ramp, t0	0	−0.0248 ± 0.0011	0
BLISS map (M(x, y))	Yes	Yes	Yes
Min. number of points per bin	4	4	4
Total frames	67008	54080	70000
Rejected frames (%)	0.714840	0.377219	0.504286
Frames usedc	61520	53865	68646
Free parameters	8	5	4
AIC value	184044	184044	184044
BIC value	184216	184216	184216
SDNR	0.0083440	0.0083556	0.0083691
Uncertainty scaling factor	0.818677	0.818525	0.821463
Photon-limited S/N (%)	83.1	83.1	82.6

Notes. ^arms frame-to-frame position difference. ^bMJD = BJD−2450,000. ^cWe exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for bad pixels in the photometry aperture.

Download table as:  ASCIITypeset image

Table 15. Best-fit Joint Eclipse Light-curve Parameters

Parameter	HD149bs11	HD149bs21	HD149bs31	HD149bs32	HD149bs33
Wavelength (μm)	3.6	4.5	5.8	5.8	5.8
Array position ($\bar{x}$, pixel)	14.58	14.57	14.50	14.42	14.78
Array position ($\bar{y}$, pixel)	15.66	14.99	14.37	14.17	14.67
Position consistencya (δx, pixel)	0.009	0.009	0.020	0.020	0.015
Position consistencya (δy, pixel)	0.007	0.004	0.018	0.013	0.017
Aperture size (pixel)	3.75	2.75	2.75	2.75	2.75
Sky annulus inner radius (pixel)	7.0	7.0	7.0	7.0	7.0
Sky annulus outer radius (pixel)	15.0	15.0	15.0	15.0	15.0
System fluxb Fs (μJy)	528456 ± 10	334037 ± 40	216010 ± 70	214102 ± 40	212621 ± 80
Eclipse depth (%)	0.040 ± 0.003	0.034 ± 0.006	0.044 ± 0.010	0.044 ± 0.010	0.044 ± 0.010
Brightness temperature (K)	2000 ± 60	1650 ± 120	1600 ± 200	1600 ± 200	1600 ± 200
Eclipse midpointc (orbits)	0.5003 ± 0.0004	0.4994 ± 0.0007	0.502 ± 0.004	0.501+0.002− 0.005	0.500+0.003− 0.005
Eclipse midpointd (MJDUTC)	4535.8756 ± 0.0010	4596.2669 ± 0.0019	4325.941 ± 0.013	4633.658 ± 0.010	4903.989 ± 0.012
Eclipse midpointd (MJDTDB)	4535.8764 ± 0.0010	4596.2676 ± 0.0019	4325.942 ± 0.013	4633.659 ± 0.010	4903.990 ± 0.012
Eclipse duration (t4–1, hr)	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02
Ingress/egress time (t2–1, hr)	0.230 ± 0.011	0.234 ± 0.012	0.234 ± 0.012	0.234 ± 0.012	0.234 ± 0.012
Ramp equation (R(t))	10	2+	2+	2+	2+
Ramp, r0	0	29 ± 9	30 ± 6	42 ± 5	26 ± 5
Ramp, r1	0	7 ± 4	8 ± 3	14 ± 2	6.4 ± 2.0
Ramp, r2	−0.0038 ± 0.0008	0	0	0	0
BLISS map (M(x, y))	Yes	Yes	Yes	No	No
Min. number of points per bin	4	5	5	...	...
Total frames	54080	54080	54080	54080	54080
Rejected frames (%)	0.308802	0.160873	0.268121	0.277367	0.312500
Frames usede	50769	48769	50919	50930	53911
Free parameters	6	7	7	4	4
AIC value	50775	48776	155775	155775	155775
BIC value	50828	48837	155924	155924	155924
SDNR	0.00280365	0.00388070	0.0120115	0.0119551	0.0119188
Uncertainty scaling factor	0.401020	0.489317	0.914694	1.02508	1.02476
Photon-limited S/N (%)	93.9	89.9	73.6	73.7	74.3

Notes. ^arms frame-to-frame position difference. ^bWe multiply the HD149bs32 and HD149bs33 measured system fluxes by 0.968 to correct for an IRAC flux conversion issue in the S18.18 pipeline. ^cBased on the period and ephemeris time given by K09. ^dMJD = BJD−2450,000. ^eWe exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for bad pixels in the photometry aperture.

Download table as:  ASCIITypeset image

Table 16. Best-fit Joint Eclipse Light-curve Parameters

Parameter	HD149bs41	HD149bs42	HD149bs43	HD149bs44	HD149bs51	HD149bs52
Wavelength (μm)	8.0	8.0	8.0	8.0	16	16
Array position ($\bar{x}$, pixel)	15.06	14.40	15.04	14.64	26.07	23.95
Array position ($\bar{y}$, pixel)	14.45	15.09	14.13	15.08	22.16	20.65
Position consistencya (δx, pixel)	0.091	0.013	0.011	0.012	0.016	0.018
Position consistencya (δy, pixel)	0.077	0.011	0.012	0.011	0.021	0.018
Aperture size (pixel)	4.0	3.5	3.75	3.5	2.0	2.0
Sky annulus inner radius (pixel)	7.0	7.0	7.0	7.0	6.0	6.0
Sky annulus outer radius (pixel)	15.0	15.0	15.0	15.0	11.0	12.0
System fluxb Fs (μJy)	117190 ± 100	116960 ± 10	117818 ± 6	118365 ± 70	18618 ± 6	18805 ± 11
Eclipse depth (%)	0.052 ± 0.006	0.052 ± 0.006	0.052 ± 0.006	0.052 ± 0.006	0.085 ± 0.032	0.085 ± 0.032
Brightness temperature (K)	1650 ± 110	1650 ± 110	1650 ± 110	1650 ± 110	1800 ± 600	1800 ± 600
Eclipse midpointc (phase)	0.5002 ± 0.0007	0.5010 ± 0.0009	0.4961 ± 0.0012	0.4988 ± 0.0006	0.5013 ± 0.0014	0.5013 ± 0.0014
Eclipse midpointd (MJDUTC)	3606.962 ± 0.002	4567.513 ± 0.003	4599.133 ± 0.003	4912.613 ± 0.002	4317.311 ± 0.004	4343.194 ± 0.004
Eclipse midpointd (MJDTDB)	3606.963 ± 0.002	4567.513 ± 0.003	4599.134 ± 0.003	4912.614 ± 0.002	4317.311 ± 0.004	4343.194 ± 0.004
Eclipse duration (t4–1, hr)	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02	3.29 ± 0.02
Ingress/egress time (t2–1, hr)	0.234 ± 0.012	0.234 ± 0.012	0.234 ± 0.012	0.234 ± 0.012	0.234 ± 0.012	0.234 ± 0.012
Ramp equation (R(t))	2–	11	10	2–	11	2–
Ramp, r0	15.4 ± 1.2	0	0	14 ± 8	0	62 ± 10
Ramp, r1	3.1 ± 0.5	0	0	0 ± 4	0	24 ± 4
Ramp, r2	0	0.083 ± 0.002	−0.0014 ± 0.0016	0	0.413 ± 0.012	0
Ramp, r3	0	−0.67 ± 0.11	0	0	−5.1 ± 0.2	0
BLISS map (M(x, y))	No	Yes	Yes	Yes	No	No
Min. number of points per bin	⋅⋅⋅	8	4	10	⋅⋅⋅	⋅⋅⋅
Total frames	44352	54080	60500	54080	1050	1050
Rejected frames (%)	0.47574	0.432692	0.634711	0.488166	0.285714	0.571429
Frames usede	44041	48801	60100	46274	1047	1044
Free parameters	15	4	3	4	10	12
AIC value	194243	194243	194243	194245	2113	2113
BIC value	194518	194518	194518	194540	2237	2237
SDNR	0.00845740	0.00833381	0.00837245	0.00847674	0.00311603	0.00337752
Uncertainty scaling factor	0.816657	0.813220	0.819215	0.981709	0.374578	0.402919
Photon-limited S/N (%)	82.3	83.1	82.4	81.6	27.2	24.9

Notes. ^arms frame-to-frame position difference. ^bWe multiply the HD149bs44 measured system flux by 0.973 to correct for an IRAC flux conversion issue in the S18.18 pipeline. ^cBased on the period and ephemeris time given by K09. ^dMJD = BJD−2450,000. ^eWe exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for bad pixels in the photometry aperture.

Download table as:  ASCIITypeset images: 1 2