The provenance of Late Ediacaran and Early Ordovician siliciclastic rocks in the Southwest Central Iberian Zone: Constraints from detrital zircon data on northern Gondwana margin evolution during the late Neoproterozoic
Pereira M.F.,University of Evora |
Chichorro M.,New University of Lisbon |
Sola A.R.,National Laboratory of Energy and Geology |
Medina J.,University of Aveiro |
Precambrian Research | Year: 2012
U-Pb geochronology of detrital zircon from Late Ediacaran (Beiras Group greywackes) and Early Ordovician (Sarnelhas arkosic quartzites and Armorican quartzites of Penacova) sedimentary rocks of the southwest Central Iberian Zone (SW CIZ) constrain the evolution of northern Gondwana active-passive margin transition. The LA-ICP-MS U-Pb data set (375 detrital zircons with 90-110% concordant ages) is dominated by Neoproterozoic ages (75% for the greywakes and 60% for the quartzites), among which the main age cluster (more significant for Beiras Group greywackes) is Cryogenian (c.840-750. Ma), while a few Mesoproterozoic and Tonian ages are also present (percentages <8%). These two features, and the predominance of Cryogenian ages over Ediacaran ages, distinguish the Beiras Group greywackes (SW CIZ) from the time-equivalent Serie Negra (Ossa-Morena Zone - OMZ), with which they are in inferred contact. The age spectra of the Beiras Group greywackes also reveal three major episodes of zircon crystallisation in the source area during the Neoproterozoic that are probably associated with a long-lived system of magmatism that developed either along or in the vicinity of the northern Gondwana margin at: (1) c. 850-700. Ma - Pan-African suture (not well represented in OMZ); (2) c. 700-635. Ma - early Cadomian arc; and (3) c. 635-545. Ma - late Cadomian arc. Comparison of Neoproterozoic ages and those of the Paleoproterozoic (c. 2-1.8. Ga) and Archean (mainly Neoarchean - 2.8-2.6. Ga, but also older) in the Beiras Group greywackes with U-Pb ages of Cadomian correlatives shows that: (1) SW CIZ, OMZ, Saxo-Thuringian Zone, North Armorican Cadomian Belt and Anti-Atlas) evolved together during the formation of back-arc basins on the northern Gondwana active margin and (2) all recorded synorogenic basins that were filled during the Ediacaran by detritus resulting from erosion of the West African craton, the Pan-African suture and a long-lived Cadomian magmatic arc. Differences in detrital zircon age populations in the greywackes of the Beiras Group (SW CIZ Cadomian basement) and the Serie Negra (OMZ Cadomian basement) are also observed in their respective overlying Early Ordovician quartzites. Since both these SW Iberia Cadomian basements evolved together along the active margin of Gondwana (but sufficiently separated to account for the differences in their detrital zircon content), this continuation of differing zircon populations into the Early Ordovician suggests that the inferred contact presently juxtaposing the Beiras Group and the Serie Negra is not pre-Early Ordovician and so is unlikely to demonstrate a Cadomian suture. © 2011 Elsevier B.V. Source
We observed the transit of TRAPPIST-1c followed 12 minutes later by the transit of TRAPPIST-1b on 4 May 2016. Observations were conducted using the HST/WFC3 infrared G141 grism (1.1–1.7 μm) in round-trip scanning mode10. Using the round-trip scanning mode involves exposing the telescope during an initial forward slew in the cross-dispersion direction, and exposing during an equivalent slew in the reverse direction (details on the trade-offs behind round-trip scanning are below). Scans were conducted at a rate of ∼0.236 pixels per second, with a final spatial scan covering ∼26.4 pixels in the cross-dispersion direction on the detector. We use the IMA output files from the CalWF3 pipeline, which have been calibrated using flat fields and bias subtraction. We applied two different extraction techniques which lead to the same conclusions. The first technique extracts the flux for TRAPPIST-1 from each exposure by taking the difference between successive non-destructive reads. A top-hat filter27 is then applied around the target spectrum, measured ±18 pixels from the centre of the TRAPPIST-1 scan, and sets all external pixels to zero. Next, the images are reconstructed by adding the individual reads per exposure back together. Using the reconstructed images, we extracted the spectra with an aperture of 31 pixels around the computed centring profile for both forward and reverse scan observations. The centring profile is calculated on the basis of the pixel flux boundaries of each exposure, which was found to be fully consistent across the spectrum for both scan directions. The second technique uses the final science image for each exposure and determines for each frame the centroid of the spectrum in a box 28 pixels by 136 pixels, which corresponds to the dimensions of the irradiated region of WFC3’s detector for our present observations. It then extracts the flux for 120 apertures of sizes ranging along the dispersion direction from 24 pixels to 38 pixels (with 1-pixel increments), and along the cross-dispersion direction from 120 pixels to 176 pixels (with 8-pixel increments)—we found the SDNR to be mostly insensitive to the aperture size along the dispersion direction. The best aperture was selected via minimization of the SDNR of the white-light-curve best fit, which is minimum for an aperture of 32 pixels by 157 pixels. Both techniques subtract the background for each frame by selecting a region well away from the target spectrum, calculating the median flux, and cleaning cosmic-ray detections with a customized procedure28. Our observations present three cosmic-ray detections that were not flagged by the CalWF3 pipeline. The exposure times were converted from Julian date in universal time (jd ) to the barycentric Julian date in the barycentric dynamical time (bjd ) system29. Both extraction methods result in the same relative flux measurements from the star and SDNR (~240 p.p.m. in the white-light curve), as the build-up of flux over successive reads is stable. We elected to obtain our observations using the round-trip scan mode in order to increase the integration efficiency compared with the standard forward scan mode. We note that, owing to slight differences in scan length/position and to the way in which the detector is read out (that is, if the direction of the scan is in the same direction as the column readout, then the integration time will be marginally longer than if the reverse were true10), round-trip scan mode results in measurable differences in the total flux of the forward scan exposures compared with the reverse scan exposures. This effect has been seen for previous WFC3 observations14, 30 in round-trip mode, and has been corrected for in two main ways. The first method involves splitting the data into two sets, one for forward scan exposures and one for reverse scan exposures, effectively halving the number of exposures per light curve, but doubling the number of light curves obtained. Each of these data sets is then analysed separately and the results combined at the end14. The second method uses the median of each scan direction to normalize the two light curves, which are then recombined and normalized before the light-curve analysis to obtain the transit parameters30. In the TRAPPIST-1 data, we measure a ~0.1% difference in flux level between the two scans. Because of the limited phase coverage of the combined transits, to retain the most information about the combined and separate effects of each planet (the transit of TRAPPIST-1c followed by that of TRAPPIST-1b), we cannot apply the first method. However, by applying the second method we found significant remaining structure in the residuals, suggesting that the correction is only partial. Previous observations using the round-trip scan30 show that the offset between the light curves obtained with each scan varies significantly from orbit to orbit, suggesting that correcting via a median combine across visits is not optimal. In addition, the total flux is affected asymmetrically by other instrumental systematics—for example, the detector ramp consistently yields a first measurement in the forward direction that is significantly lower than average—thus biasing the median combine. Therefore, we corrected for the flux offset induced by the round-trip scan mode on the basis of the offset in the residuals for each HST orbit individually. To do so, we estimate in our forward model the ‘intermediate residuals’, based on the data corrected for the transit model and the instrumental systematics. For each orbit, we estimate the mean of these residuals for each scan direction (m and m , for the mean of the residuals of the forward-scan exposures and the reverse-scan exposures, respectively). The ratio of the fluxes measured in reverse-scan exposures to the shared baseline level is 1 + m ; the ratio is 1 + m for forward-scan exposures. We therefore correct for their offsets by dividing each set of exposures by their respective ratio. We first analysed the white-light curve by summing the flux across all wavelengths. We fitted the transits of TRAPPIST-1b and TRAPPIST-1c by using the transit model of ref. 17, while correcting for instrumental systematics. We followed the standard procedure for analysing HST/WFC3 data by fixing the planets’ orbital configurations—all but the orbital inclinations, which are currently poorly constrained for TRAPPIST-1’s planets—to the ones reported in the discovery report3, while determining the transit times and depths. We used priors on the band-integrated limb-darkening coefficients (LDCs) derived from the PHOENIX model intensity spectra15, and on the planets’ orbital inclinations—these parameters being potentially correlated with the transit depth estimates—to adequately account for our present state of knowledge on TRAPPIST-1. We used different analysis methods to confirm the robustness of our conclusions. The first method uses a least-squares minimization fitting (L–M) implementation12 to investigate a large sample of systematic models—which include corrections in time, HST orbital phase, and positional shifts in wavelength on the detector—and marginalize over all possible combinations to obtain the transit parameters. The L–M implementation fits the light curves for each systematic model and approximates the evidence-based weight of each systematic model using the Akaike information criterion31. It does so while keeping the LDCs fixed to the best estimates presented below, and the orbital inclinations fixed to the estimates from ref. 3. The highest weighted systematic models include linear corrections in time, as well as linear corrections in HST orbital phase or in the shift in wavelength position over the course of the visit. Therefore, using marginalization across a grid of stochastic models allows us to account for all tested combinations of systematics and to obtain robust transit depths for both planets, separately and in combination. For this data set, the evidence-based weight approximated for each of the systematic models applied to the data indicates that all of the systematic models fit equally well to the data, and that no one systematic model contributes to the majority of the corrections required to obtain the precision presented (Extended Data Fig. 1). In other words, instrumental systematics affect our observations only marginally. We carried out independent analyses of the data by using adaptive Markov chain Monte Carlo (MCMC) implementations32, 33. For each HST light curve, the transit models17 of TRAPPIST-1b and TRAPPIST-1c are multiplied by baseline models that account for the visit-long trend observed in WFC3 light curves, WFC3’s ramp, and the ‘HST breathing’ effect12. For these analyses, priors are used for the LDCs and the orbital inclinations. We find that the visit-long trend is adequately accounted for with a linear function of time, the ramp with a single exponential in time, and the breathing with a second-order polynomial in HST’s orbital phase. More-complex baseline models were tested and gave consistent results, as revealed by the marginalization study. We calculated the transmission spectrum by fitting the transit depth of TRAPPIST-1b and TRAPPIST-1c simultaneously in each spectroscopic light curve. We divided the spectral range between 1.15 μm and 1.7 μm into 11 equal bins of Δλ = 0.05 μm. We applied again the two techniques described above to analyse each spectroscopic light curve, resulting in the combined and independent transmission spectra of TRAPPIST-1b and TRAPPIST-1c. An L–M implementation12 and the adaptive MCMC implementations produced consistent results for each stage of the analysis. We determined limb-darkening coefficients by fitting theoretical specific intensity spectra (I) downloaded from the Göttingen spectral library (http://phoenix.astro.physik.uni-goettingen.de/?page_id=73), which is described in ref. 15. The intensity spectra are provided on a wavelength grid with 1-Å cadence for 78 μ values, where μ is the cosine of the angle between an outward radial vector and the direction towards the observer at a point on the stellar surface. We integrated I over one broad and 11 narrow wavelength intervals, used in our analysis of the transit light curve. We divided I for each wavelength interval by I , the value of I at the centre of the stellar disc (where μ = 1). Because the PHOENIX code calculates specific intensity spectra in spherical geometry, the PHOENIX μ grid extends above the stellar limb relevant to exoplanet transit calculations. When fitting limb-darkening functions, PHOENIX μ values should be scaled to yield μ′ = 0 at the stellar radius34. We define μ′ = (μ − μ )/(1 − μ ), where I/I = 0.01 at μ = μ . The value of μ is a function of wavelength. We then fitted two commonly used functional forms for limb darkening18: When fitting, we ignored points with μ′ < 0.05. Extended Data Fig. 2 shows the limb-darkening fits for the 12 wavelength intervals in our transit light curve analysis. We calculated fits for four stellar models with effective temperatures of 2,500 K and 2,600 K and logarithmic surface gravities of 5.0 and 5.5. We then linearly interpolated the limb-darkening coefficients to an effective temperature of 2,550 K and gravity 5.22, appropriate for TRAPPIST-1 (ref. 3). We simulated the theoretical spectra for TRAPPIST-1b and TRAPPIST-1c using the model introduced in ref. 19. We used atmospheric temperatures equal to the planets’ equilibrium temperatures assuming a Bond albedo of 0.3 (these temperatures being 366 K for TRAPPIST-1b and 315 K for TRAPPIST-1c). The use of isothermal temperature profiles set at the equilibrium temperatures is conservative, as it does not account for possible additional heat sources or temperature inversion and results in a possible underevaluation of the atmospheric scale height. Our assumption regarding the temperature profiles does not affect our conclusion; variations of 50 K (that is, ∼15%) in the atmospheric temperature modify the amplitude of the transmission spectra by up to ∼15%, because at first order their amplitudes scale with the temperature. The planetary masses being unconstrained, we conservatively use a mass of 0.95M and 0.85M for TRAPPIST-1b and TRAPPIST-1c respectively—the maximum masses that would allow them to possess hydrogen/helium envelopes greater than 0.1% of their total masses given their radii20. We use the atmospheric compositions of the ‘mini-Neptune’ and ‘Halley world’ models introduced in ref. 35 to simulate the hydrogen-dominated and water-dominated atmospheres, respectively. We simulated the effect of optically thick cloud or haze at a given pressure level by setting to zero the transmittance of atmospheric layers with a higher pressure. The feature at 1.4 μm arises from water absorption; the feature at 1.15 μm for the water-dominated atmosphere arises from methane absorption. We compared the transmission spectra, allowing for a vertical offset to account for our a priori ignorance of the optically thick radius by setting the mean of each spectrum to zero. The significance of the deviation of each transmission spectrum from the WFC3 measurements is shown in Fig. 3. Significance levels less than 3σ mean that the data are consistent with that model within the reported errors. We rule out the presence of a cloud-free hydrogen-dominated atmosphere for either planet at the 10σ level through the combined transmission spectrum (and at a lesser 7σ level through their individual spectra). The measurements are consistent with volatile (for example, water)-rich atmospheres or hydrogen-dominated atmospheres with optically thick clouds or hazes located at larger pressures than 10 mbar. Conversion of the ut times for the photometric measurements to the bjd system was carried out using the online program created by J. Eastman and distributed at http://astroutils.astronomy.ohio-state.edu/time/utc2bjd.html. We have opted not to make available the codes used for data extraction as they are an important part of the researchers’ toolkits. For the same reason, we have opted not to make available all but one of the codes used for data analysis. The MCMC software used by M.G. to analyse independently the photometric data is a custom Fortran 90 code that can be obtained upon request. The custom IDL code used to determine limb-darkening coefficients can be obtained upon request.
Initial optical spectroscopic observations of M81 ULS-1 were carried out with Keck/LRIS on 13 April and 17 April 2010, revealing broad Balmer emission lines as if from an accretion disk21. The blueshifted Hα emission line (Hα−) is shown at 5,530 Å in the spectra, with a shift in the line centre of 10 Å ± 2 Å between those two observations (that is, two epochs separated by four nights). M81 ULS-1 was later observed using GTC/OSIRIS on 8 April, 7 May and 8 May 2015, masked with the 0.6″ slit followed by the R1000R grating, which yields a resolution of about 7 Å. The spectra were reduced in a standard way with IRAF (Image Reduction and Analysis Facility) software (http://iraf.noao.edu). After bias subtraction and flat correction, dispersion correction was carried out on the basis of the line lists given in the OSIRIS manual (http://www.gtc.iac.es/instruments/osiris/). Raw spectra were then extracted with an aperture size of 1″, and a standard star taken at each night was used to make the flux calibration. On 22 April 2015, another observation of M81 ULS-1 was carried out using Keck/LRIS with the 1.0″ slit. The light was split with a beam dichroic of 6,800 Å to the blue and red sides, followed by using the 300/5,000 and 400/8,500 gratings, which yields a resolution of ~8 Å. The spectrum was reduced with the IDL (Interactive Data Language) pipeline designed for the W. M. Keck Observatory. Extended Data Table 1 lists the basic information obtained from the 2010 and 2015 observations. Both Hα and Hα− emission lines are detected in all the spectra, and their line properties are calculated from Gaussian line profile fitting. Extended Data Table 2 lists the central wavelength, the full width at half-maximum (FWHM) and the equivalent width for each fitted emission line. The observed Hα− central wavelengths, λ , correspond to projected velocities, v , from −0.17c to −0.14c in these observations, given by . In the case of SS 433, the equivalent widths of Hα emission lines are tightly correlated with the phases of the precession, and those of the Hα− emission follow a similar trend but with a phase delay31. We use the power of the emission lines, calculated from the area of the Gaussian fitting, as representative of the emission intensity, because the observed continuum from M81 ULS-1 varies markedly between observations. Extended Data Fig. 1 shows that the power of the Hα emission lines from the accretion disk is positively correlated with that of the Hα− emission lines, suggesting a link between the accretion and the jet. The variations in the power of the emission lines are asymmetrical, with smooth rises and steeper declines around 7 May 2015 (Extended Data Fig. 2)—similar to the variations seen in the SS 433 emission lines31. The rate at which the projected Hα− velocity changes seems to be slower during 2015 than it was during 2010. The rate of change in 2015 was roughly 0.8 Å per day, whereas it was 2.6 Å per day in 2010; if the velocity shift is due to precession, this difference may be explained naturally, because the 2015 observations are sampling a different part of the precession cycle. We can estimate a minimum likely precession period by assuming that the turning point of the precession cycle occurred at around the time of the observations of 7 and 8 May 2015 (see Extended Data Figs 2, 3). If so, then, after the wavelength of the emission lines reached the maximum on 8 May (Extended Data Fig. 2), the Hα− emission probably turned back to the short wavelength with the rate of roughly 0.8 Å per day, indicating that the half-precession period must be longer than 30 days. There is a 115-Å gap between the observations of 13 April 2010 and 8 April 2015 (Extended Data Fig. 3), and if we assume that the maximum rate at which the wavelength decreases is 2.6 Å per day, then the Hα− line needs 44 days to move by the required amount. Therefore, the half-precession period is probably longer than 74 days, and lower limit of the precession period is about 148 days. More time-resolved spectra are needed in order to derive an accurate period and to characterize further the apparent precession of the jets. Given the existence of blueshifted Hα− emission lines from the approaching jets, redshifted Hα emission lines (Hα+) would be expected from receding jets, albeit with much lower intensities (because of Doppler boosting effects). Assuming symmetrical and steady jets, the boosting factors (D) for the lines emitted from the approaching and the receding jets are given by and , respectively. The total flux of a blueshifted or redshifted line in the observer frame is boosted by a factor of D3. The expected central wavelengths of the two lines are given by , and , and the corresponding redshifts are and . D and z values for Hα− in all observations are listed in Extended Data Table 3. We have searched for the redshifted Hα− emission lines in all observations. A weak emission line feature was detected at ~3σ at around 7,524 Å (Fig. 1), roughly symmetrical to Hα− at 5,648 Å, in one of the GTC exposures during the night of 8 April 2015. If this marginal detection were the redshifted Hα− line, its boosting factor would be , and the ratio of the received total flux of the blueshifted line to that of the redshifted line should be , which is roughly consistent with the observed flux ( ). However, the observed wavelength is not consistent with the expected Hα− wavelength given the blueshifted Hα− at 5,648 Å, that is, . Assuming the extreme case, θ = 0º, then we have β = 0.1491. If the receding jet has the same velocity, then the expected central wavelength of Hα+ should be 7,626.4 Å, which is about 104 Å larger than the detected line. If we assume that θ = 10°/20°/30°, then β = 0.152/0.160/0.177, and the expected central wavelength of Hα+ is 7,632 Å/7,650 Å/7,688 Å, which is about 108 Å/126 Å/164 Å larger than the detected line. The discrepancy becomes larger for larger inclination angles. This casts doubt on the identification of the 7,524-Å line feature as Hα+, unless the jets are asymmetrical or fast-changing. We may have not detected the redshifted Hα+, but the non-detection is not surprising given the Doppler boosting effects, and other realistic explanations. For example, the receding jets may be blocked by the optically thick outflows if this system is a supercritically accreting black-hole system, as described in the text. No candidate Hα+ emission lines were detected at all in the 7 and 8 May GTC observations, or in the Keck spectrum. Even if the 8 April line were a true Hα+ emission line, this non-detection would not be surprising, given the lower equivalent widths of Hα− on 7 and 8 May, and the relatively lower sensitivity in the red channel of LRIS. There have been 19 Chandra/ACIS observations of the nuclear region of M81, where ULS-1 resides. All of these observations were derived from the Chandra archive and analysed uniformly with CIAO 4.7 software tools (http://cxc.harvard.edu/ciao/). Point sources were detected with WAVDETECT on the individual Chandra images. As listed in Extended Data Table 4, the photon counts were extracted from the source ellipses enclosing 95% of the total photons as reported by WAVDETECT, which was run with scales of 1″, 2″, 4″ and 8″ in the 0.3 to 8.0 keV band. The spectra in the high-flux states (>10 counts per kilosecond) were fitted by absorbed blackbody models, with the spectral parameters presented in Extended Data Table 4, all of which show that M81 ULS-1 has been persistently supersoft in these observations. In addition, the spectra in the high-flux and low-flux states were added together into combined high- and low-state spectra, and were also fitted in the band 0.3–8.0 keV. Using the fitted absorbed blackbody model, we calculated the 0.3–8.0 keV flux, the 0.3–8.0 keV luminosity and the bolometric luminosity with the distance of 3.63 megaparsecs for M81 (ref. 32). As plotted in Extended Data Fig. 4, M81 ULS-1 displays a soft excess below 0.3 keV as compared to the best-fit model for 0.3–8.0 keV. However, considering that the response matrix for Chandra observations is not well calibrated below 0.3 keV, we refrain from interpreting this soft excess. Nonetheless, it is clear that M81 ULS-1 has very different spectral properties from the other known microquasars. Moreover, these uncertainties in calibration below 0.3 keV might merely make the intrinsic spectral differences between M81 ULS-1 and the other known microquasars even larger (that is, the energy distribution from M81 ULS-1 might be even softer than observed). The optical spectra were reduced with IRAF, available at http://iraf.noao.edu/. All of the emission lines in Extended Data Table 2 were fitted with the curve-fitting toolbox based on Matlab (http://www.mathworks.com/help/curvefit/index.html). The Chandra archive data were analysed with CIAO 4.7, which can be downloaded from http://cxc.harvard.edu/ciao/download/.
We observed LkCa 15 using non-redundant masking (NRM)11 with LBTI/LMIRCam30, 31 in December 2014 and February 2015. NRM transforms a conventional telescope into an interferometric array through the use of a pupil-plane mask, providing better point-spread function (PSF) characterization and resolving angular scales even within λ/D. We used LMIRCam’s 12-hole mask in single-aperture mode, yielding 1.4–7.0 m baselines. We broke up the observations into “visits,” each consisting of an identical set of integrations on LkCa 15 and an unresolved calibrator star (see Extended Data Table 1). We used three calibrators to lessen the possibility of contamination by a binary calibrator, and included one calibrator, GM Aur, from those observed previously at Keck8. We let the sky rotate throughout the observations, facilitating calibration of quasi-static speckles. At Ks and L′ we observed LkCa 15 at parallactic angles between −37° and 65°, and −65° and 65°, respectively. The NRM images show the interference fringes formed by the mask, the Fourier transform of which yields complex visibilities. Sampling the complex visibilities, we calculated squared visibilities (power versus baseline) and closure phases (sums of phases around three baselines forming a triangle). Closure phases eliminate atmospheric phase errors, leaving behind the sum of the intrinsic source phases. The LBT raw closure phase scatter was ~3°, significantly lower than the uncertainties in previous NRM observations8 (~4°). For each closing triangle, we fitted a polynomial to all calibrator closure phases as a function of time. We sampled the polynomial at the time of each target observation and subtracted it from each target closure phase. We calibrated using a variety of functions; of these, polynomials up to 2nd order in time provided the lowest-scatter calibrated closure phases, with standard deviations of 1.7° at Ks and 1.9° at L′. We calibrated the squared visibilities similarly, dividing by the calibrator rather than subtracting. We calibrated the mask baselines using the observed power spectra and knowledge of the filter bandpass and plate scale32. We fitted models directly to kernel phases33, 34, linearly independent combinations of closure phases, to search for companions. We modelled the central star as a delta function and each companion as another delta function located a distance s from the star, at position angle PA, and with contrast Δ. We sampled the synthetic complex visibilities at the same baselines and parallactic angles as the data, and performed a grid fit, using a Δχ2 to determine our parameter confidence intervals. Due to a known degeneracy between companion separation and contrast35, brighter companions at smaller separations provide equally good fits as those fainter and farther out. We thus performed fits to individual wavelengths to verify that the positions of b and c agreed across wavelength, then calculated a best fit where the companions coincided at Ks and L′ (see Table 1). The model grids in this study required ~50,000 CPU hours to generate, but were computed in a reasonable amount of time using the University of Arizona’s El Gato supercomputer. We also reconstructed images from the closure phases. To produce cleaner images, we re-calibrated the closure phases towards the best-fit Ks + L′ model using an optimized calibrator weighting technique applied in previous NRM studies8. This calibration is similar to the locally optimized calibration of images (LOCI)36 technique applied in direct imaging. Since this scheme can remove signal and underestimate errors, we applied it only to produce images (see Extended Data Fig. 1), using the polynomial calibration to estimate companion parameters. As a consistency check, we reconstructed images using both the BiSpectrum Maximum Entropy Method (BSMEM)37 and the Monte-Carlo MArkov Chain IMager algorithm (MACIM)38. The companion positions agree between the two algorithms, although BSMEM produces more compact emission towards each component. BSMEM has been shown to provide the most faithful image reconstruction of any available algorithms in blind tests39. Orbital fitting required parameter errors for the previously published8 Keck observations and the LBT observations to be consistently estimated. The published errors for the 2009–2010 companion parameters were generated using the nonlinear algorithm MPFIT40. While nonlinear fitters are often employed for computational expediency, the Levenberg–Marquardt algorithm can easily fall into a local minimum and underestimate parameter errors. The LBT grid χ2 surfaces show local minima for both two- and three-companion fits, rendering MPFIT unreliable unless the starting parameters were very close to the global minimum. We compared MPFIT and grid-based parameter errors for the LBT data, and found that MPFIT significantly underestimated the errors (Fig. 2). To create a “typical” error bar for each Keck companion, we estimated the error bar dependence on contrast using the LBT fits. Errors increased with decreasing companion flux, which we parameterized as a square root dependence. For a given Keck companion we thus scaled our LBT errors by the square root of the LBT-to-Keck flux ratio. We inflated the Keck error bars by a factor of 1.3, the ratio of the uncalibrated closure phase scatter in the Keck data (~4°) to that for the LBT data (~3°). We capped the separation upper limits at 3λ/D, where D is Keck’s telescope diameter, 10 m, since the largest LBT upper limit was at nearly 3λ/D, and companions at those distances are no longer subject to the separation-contrast degeneracy. We observed LkCa 15 on November 15 and 22, 2014, as part of the Giant Accreting Protoplanet Survey (GAPplanetS), a visible-wavelength survey of bright transition disks. GAPplanetS stars are imaged simultaneously in Hα (0.656 μm, Δλ = 6 nm) and the nearby stellar continuum (0.642 μm, Δλ = 6 nm) with the 585-actuator Magellan Adaptive Optics systems SDI camera12, 41, 42. The continuum channel provides a sensitive, simultaneous probe of the stellar PSF, allowing for effective removal of residual starlight and isolation of Hα emitting sources12, 43. The observations used new single-substrate narrowband Hα and continuum filters, a significant improvement over the previous VisAO SDI filters, which suffered from ghost images12. Seeing during the November 15 observations was better than the site median (0.56 ± 0.06″), winds were low (3.6 ± 0.9 mph), and humidity was typical of the season (37.0 ± 2.8%). Strehl ratio was low (<10%), and difficult to measure meaningfully. We characterized image quality using the stellar full-width at half-maximum (FWHM), 0.07″ (at 0.65 μm over 30 s integrations), a significant improvement over the seeing. We collected 316 30-s closed-AO-loop images, with a total of 156 min of integration time and 48.6° of sky rotation. We selected the 149 LkCa 15 images with the lowest residual wavefront error (~50%), equivalent to 74.5 min of exposure time. This image subset had 47.6° of sky rotation, with the rotational space well sampled. The November 22 data were not of sufficient quality to recover LkCa 15 b, due to lower sky rotation (27.0°), shorter total integration (91 min), and shallower individual exposures (15 s). Injected positive planets with the same separation as LkCa 15 b were only recoverable with SNR > 3 at contrasts > 5 × 10−2 (nearly an order of magnitude brighter than the measured November 15 LkCa 15 b contrast). For this reason, we discuss only the November 15 data set in the rest of the paper. Images were first bias-subtracted, registered, and aligned via cross-correlation. The Hα flat field image showed very little non-uniformity across the field (<1%), so a flat field was not applied. We masked CCD dust spots visible in the flat field wherever they affected the image throughput by more than 2%. We processed the aligned data using angular differential imaging (ADI44), comparing the “classical” method of using a single median PSF for all images (cADI45) to the Karhunen–Loeve image processing (KLIP46) algorithm, which calculates a least-squares optimum PSF for each image. LkCa 15 b was detected in the Hα channel via both methods, as shown in Extended Data Fig. 2. The planet was not detected in continuum with either method, so continuum images were used as a probe of PSF residuals and scattered light emission from the inner disk. Subtraction of the processed continuum images from the Hα images (“ASDI”) left behind only true Hα emission12. We constructed the stellar PSF by median combining images in 0.5° rotational bins and then median combining again to produce a PSF evenly sampled in rotational space. We subtracted the stellar PSF from the individual images, rotated them to a common on-sky orientation and combined them. Given the small separation between LkCa 15 A and b, the planet moved by only 1.5 FWHM over the course of the observations, resulting in self-subtraction and decreasing the FWHM of the processed planet PSF to 4–5 pixels in azimuth. KLIP reductions were carried out using a well-tested custom IDL code47. To optimize reduction parameters, we maximized the signal to noise of injected planets (with the same separation and contrast as LkCa 15 b) inserted after using a negative planet to erase the LkCa 15 b signal. Planets were placed at position angles distant from known artefacts, and east or west of the star to avoid the noisier north/south region of the PSF, corresponding to the wind direction during the observations. To limit self-subtraction, the library from which KLIP builds the stellar PSF is limited to images where a planet would have rotated away from its original position. We explored the size of this exclusion region (“rotational mask”) systematically through fake planet injection, and found that a 5° mask (~1 pixel at r = 11 pixels) produced the highest signal-to-noise recoveries of injected planets. Given the stellar FWHM of 0.07, this resulted in azimuthal self-subtraction, with a processed planetary PSF of 2 pixels in azimuth. Noise in the KLIP processed images was mostly Gaussian when images were divided into several independently-optimized radial zones, indicating efficient removal of speckles. Dividing these zones azimuthally provided no additional advantage, and the final KLIP reductions shown in Extended Data Fig. 2 reflect a PSF divided into 50-pixel (0.4″) annuli. Removal of the median PSF radial profile for the entire image set aided significantly in attenuating the stellar halo, improving the ability of the KLIP algorithm to match residual speckles and enhancing contrast close to the star. We estimated photometry and astrometry by minimizing residuals after injecting a negative planet at the location of LkCa 15 b. The cube of registered and aligned Hα channel images was scaled by the chosen contrast value, multiplied by −1, and injected into the raw images before KLIP processing. Using the full Hα image cube rather than its median combination simulated variability of the PSF between images. We generated error bars by injecting false positive planets with similar separations and contrasts to LkCa 15 b after using a negative planet to eliminate the true signal. Planets were placed at position angles away from the wind direction, and spaced by at least 75% of the measured stellar FWHM. We computed the centroid and peak pixel using a 5-pixel aperture around each planet, and assigned the standard deviations in recovered flux and position as our 1σ photometric and astrometric uncertainties, respectively (see Extended Data Table 2 and Extended Data Fig. 3). To create signal-to-noise ratio (SNR) maps, we calculated a radial noise profile using the standard deviation of 1-pixel-wide annuli and divided it into the raw images. In the raw maps, LkCa 15 b has SNR ≈ 3–4. Smoothing by a Gaussian with a 2-pixel FWHM maximized the SNR of injected fake planets, so we applied this smoothing to the final science images, resulting in peak SNRs of 4.4 and 6.8 in the KLIP Hα and ASDI images, respectively. However, directly-imaged exoplanets at small separations suffer from small number statistical effects48. The underlying speckle distribution is difficult to probe given the small number of independently sampled noise regions. In an annulus at the distance of LkCa 15 b (1.3 FWHM), seven noise regions exist, leading to corrected48 SNRs of 4.1 and 6.4 for the Hα and ASDI images, respectively. The ASDI detection corresponds to a false positive probability of 3 × 10−4 using the Student’s t-distribution with 6 degrees of freedom. Comparing the LkCa 15 b SNR to the distribution of values in the ASDI SNR map (Extended Data Fig. 4), shows that it is a clear outlier. Comparison of the peak pixel in an aperture centred on b to those in the surrounding noise apertures (Extended Data Figs 4, 5) further demonstrates b’s statistical significance. In addition to the high SNR, low false positive fraction, and the statistics presented in Extended Data Fig. 4, the Hα detection is significant because it occurs at the same location as the independent LBT detection. This further reduces the probability of a false positive detection in the MagAO data, since speckles have no preferred location. Neither the existence of LkCa 15 b nor its derived parameters are dependent on our choice to include only the top 50% of raw images. The planet appears at the same location and with the same approximate brightness when processing all 316 images, as well as only the top 25% of images. An excess with SNR >3 appears at LkCa 15 b’s location with a wide range of KLIP zone geometries and rotational masks, when any number of KL modes from 2 to 100+ are removed, and whether or not the median radial profile of the PSF is subtracted before processing. We used simulated planet detections to place an upper limit on LkCa 15 b’s continuum flux. We injected planets into the raw continuum channel images with a range of contrasts and at positions near LkCa 15 b. We then measured the SNR of each simulated detection to determine the confidence at which we could detect a given contrast. As above, we apply a small number statistical correction48 to the SNR of each recovered planet. The simulations suggest that we would have detected an excess with a corrected SNR of 3 (false positive fraction of 10−2) for a continuum source with contrast greater than 5 × 10−3. Since LkCa 15 A is 1.8 times brighter at Hα than continuum, this corresponds to an Hα-to-continuum-flux ratio lower limit of 2.7. We established limits on the LkCa 15 c Hα contrast using false planet injections, first using a negative planet to eliminate the LkCa 15 b signal. We injected planets with a range of contrasts at positions sampling the LBT error ellipse for LkCa 15 c. At position angles between −40 and −52°, several 2–2.5σ peaks near c’s location boost the SNRs for recovered planets. Here, we can detect contrasts down to 2 × 10−3 with corrected48 SNRs of 3 (false positive fraction of 2 × 10−2). Position angles greater than −40° approach a noisier region of the PSF, leading to decreased sensitivity; here contrasts of 6 × 10−3 are required for adjusted signal-to-noise ratios of 3. We cannot reject or confirm accretion onto LkCa 15 c below 6 × 10−3 contrast (ΔHα = 5.6 mag) with the current data set. This improves upon previous spectro-astrometric observations, which yielded a contrast limit of ΔHα = 3.4 mag for a position angle near LkCa 15 c10. We ran a series of orbit integrations to demonstrate that stable solutions exist for b, c, and d at separations within the 1σ semimajor axis error bars (see Extended Data Figs 6, 7). We used the publicly available Swifter package49–specifically, the symplectic integrator, SyMBA50, which switches to a Burlisch-Stoer algorithm for planetary close approaches. We also ran comparison integrations with the Gauss Radau 15th order integrator and found comparable results, with minimum energy conservation of 1 part in 107 over a 10 Myr integration. We required all orbits to be nearly co-planar, with a random scatter <1°, and assigned each planet a random eccentricity below 10−4. To assess stability we integrated three different random phase combinations for 10 Myr. We found stable three body solutions out to 1–2 Myr with semi-major axes of a = 12.7 au, a = 18.6 au, a = 24.7 au. To ensure stability out to 10 Myr with orbits in the 1σ errors requires that all planets be ≤0.5 M . A wider range of orbits is allowed if d’s mass is decreased further. These constraints are in line with previous large numerical studies of equally spaced (in R ), equal-mass planets51. Note for planets b and c, there are possible resonant configurations within the predicted period ranges, which would admit somewhat higher masses.
We observed four phase curves of 55 Cancri e with the Spitzer Space Telescope in the IRAC/4.5-μm channel as part of our program ID 90208. Because of downlink constraints, these four phase-curve observations were split into eight separate observations (or Astronomical Observation Requests, AORs) each lasting half an orbit of 55 Cancri e. Details of each AOR are provided in Extended Data Table 1. The corresponding data can be accessed from the Spitzer Heritage Archive (http://sha.ipac.caltech.edu). All AORs were acquired in stare mode using a constant exposure time (0.02 s). All our data were obtained using the Pointing Calibration and Reference Sensor (PCRS) peak-up mode, which allows the observer to place the target on a precise location on the detector to mitigate the intra-pixel sensitivity variations. This observing mode increases the pointing stability and reduces the level of correlated noise in the data by a factor of 2–3 (ref. 31). AOR 48072960 experienced a 30-min interruption during data acquisition, which forces us to treat both parts of that AOR separately in the rest of this section. We do not retain the 30-min-long PCRS sequences in our analysis because the motion of the star on the detector yields large correlated noise in these data sets. Our reduction uses the basic calibrated data (BCD) that are downloaded from the Spitzer Archive. The BCD are Flexible Image Transfer System (FITS) data cubes consisting of 64 frames of 32 × 32 pixels each. Our data reduction code reads each frame, converts fluxes from the Spitzer units of specific intensity (MJy sr−1) to photon counts, and transforms the data time-stamps from BJD to BJD using existing procedures32. We did not deem it necessary to discard specific subarray frames. During the reduction process, we compute the flux, position and FWHM in each of the 64 frames of each data cube; the frames for which any of these parameters differ from the median by more than 5σ are discarded. The centroid position on the detector is determined by fitting a Gaussian to the marginal x, y distributions using the GCNTRD procedure of the IDL Astronomy User’s Library33. We also fit a two-dimensional Gaussian to the stellar PRF following previous studies34. We find that determining the centroid position using GCNTRD results in a smaller dispersion of the fitted residuals by 10% to 15% across our data set, in agreement with other Spitzer analyses35. We then perform aperture photometry for each data set using a modified version of the APER procedure using aperture sizes of ranging from 2.2 to 4.4 pixels in 0.2-pixel intervals. We choose the optimal aperture size on the basis of minimizing for each AOR, where r.m.s. is the root-mean-square of the photometric time series and β is the red-noise contribution8. The red noise is assessed over 60-min timescales because shorter timescales are irrelevant for the phase-curve signal whose periodicity is 18 h. We measure the background contribution on each frame using an annulus located 10 to 14 pixels from the centroid position. Our code also determines the FWHM of the PRF along the x and y axes. We use a moving average based on forty consecutive frames to discard data points that differ from the median by more than 5σ in background, (x, y) position or FWHM. We find that, on average, 0.06% of the data points are discarded. The resulting time series are combined into 30-s bins to speed up the analysis; this binning has been shown to have no influence on the values or uncertainties of the system parameters7. We show the optimal aperture size, corresponding r.m.s. and β for each data set in Extended Data Table 1. Intra-pixel sensitivity correction. We use an implementation of the BLISS (BiLinearly-Interpolated Sub-pixel Sensitivity)9 to account for the intra-pixel sensitivity variations, as was similarly used in a previous study using the same data set7. The BLISS algorithm uses a bilinear interpolation of the measured fluxes to build a pixel-sensitivity map. The data are thus self-calibrated. Our implementation of this algorithm is included in the Markov chain Monte Carlo (MCMC) framework presented in ref. 8. The improvement introduced by a pixel-mapping technique such as BLISS requires that the stellar centroid remains in a relatively confined area on the detector, which warrants an efficient sampling of the x–y region and, hence, an accurate pixel map. In our implementation of the method, we build a sub-pixel mesh of n2 grid points, evenly distributed along the x and y axes. The BLISS algorithm is applied at each step of the MCMC fit. The number of grid points is determined at the beginning of the MCMC by ensuring that at least five valid photometric measurements are located in each mesh box. Similar to two recent studies7, 10, we find that a further reduction of the level of correlated noise in the photometry is achieved by the inclusion of the FWHM of the PRF along the x and y axes as extra parameters in the baseline model. The PRF evolves with time and its properties are not accounted for by the BLISS algorithm. We thus combine the BLISS algorithm with a linear function of the FWHM of the PRF along the x and y axes. In addition, the baseline model for each AOR includes a flux constant. We find that including a model of the FWHM of the PRF decreases the Bayesian Information Criterion (BIC)36 by ΔBIC = 591. We show the raw data sets with the best-fit instrumental + astrophysical model superimposed in red in Extended Data Figs 1, 2, 3. The corrected photometry is shown in Extended Data Figs 4, 5, 6. The phase-curve modulation is clearly noticeable in each AOR. The behaviour of the photometric r.m.s. as a function of binning is shown for each data set in Extended Data Fig. 7. Model comparison. In our first MCMC analysis, to model the variation in the infrared emission of the planet we use F = F + Tr + Oc (in which F is the observed flux, F is the phase modulation driven by the planet, Tr is the transit model and Oc is the occulation model), and a Lambertian37 functional form for F in which A is the phase amplitude and where θ is the phase-curve offset, ϕ is the phase angle, i is the orbital inclination of the planet, P the orbital period, T the transit centre and t is time. The transit- (Tr) and occultation- (Oc) light-curve model MA (ref. 38) are summarized as in which and are the transit and occultation depths, respectively, b is the impact parameter, M is the stellar mass, and c = 2u + u and c = u − 2u are the limb-darkening linear combinations, with u and u the quadratic coefficients obtained from theoretical tables39 using published stellar parameters40. We also experimented using a sinusoid for the phase variation: F = A cos(ϕ + θ ). The fit using a sinusoid results in an amplitude A = 218 ± 50 p.p.m. and an offset value θ = 68° ± 24° east of the substellar point, in agreement with our results using a Lambertian functional form (A = 197 ± 34 p.p.m. and θ = 41° ± 12°). The Lambertian sphere model provides a better fit to the data than does the sinusoid model, with ΔBIC = 11. We also perform another MCMC analysis with no phase-curve model, hence removing two degrees of freedom (phase amplitude and phase offset). We find ΔBIC = 21 in favour of the model including the phase-curve model. We also run an MCMC fit that includes the phase amplitude, but not the phase offset. We find that this fit produces only a marginal χ2 improvement over the MCMC fit with the no-phase-curve model, but this improvement is penalized by the extra degree of freedom according to the BIC. We indeed obtain a ΔBIC = 25 in favour of the model including the phase-curve offset. Altogether, this model comparison confirms that a phase-curve model that includes a phase offset is the favoured functional form according to the BIC. We conduct two additional analyses of our entire data set to assess the robustness of our initial detection that used the BLISS mapping technique. In these two analyses, we use different approaches to (1) model the intra-pixel sensitivity of the detector and (2) change the input data format. In the first analysis, we use a simple polynomial detrending approach with a functional form that includes only the centroid position (fourth-order) and FWHM (first-order). We experimented with different polynomial orders (from one to four) for these two parameters and found that this combination globally minimizes the BIC. Each AOR has its own set of baseline coefficients. As for the BLISS mapping, the polynomial detrending is included in the MCMC fit so the baseline model and the system parameters are adjusted simultaneously to efficiently propagate the uncertainties to the final parameters. We find a level of correlated noise in the data that is only slightly larger (about 10%) than that obtained with the BLISS mapping technique. Using this method we find a phase-curve minimum of 36 ± 41 p.p.m., a maximum of 187 ± 41 p.p.m. and an offset of 50° ± 13° east of the substellar point; using the BLISS mapping, the corresponding values are 47 ± 34 p.p.m., 197 ± 34 p.p.m. and 41° ± 12° east. As previously shown10, the addition of the FWHM of the PRF in the baseline model substantially improves a fit based on only a centroid position, and, most importantly, it enables an acceptable fit to 8-h time series. In the second analysis, we aim to assess whether the phase-curve signal persists when we split our input data. All our AORs have durations of nearly 9 h, and we elect to split each of them in two to reduce the duration of each individual data set to 4.5 h. The functional forms of the baseline models are the same as for the analysis using the unsegmented input data, described above. In this additional test, we find a phase-curve minimum of 51 ± 51 p.p.m., a maximum of 216 ± 51 p.p.m. and an offset of 54° ± 16° east of the substellar point. These results are in good agreement with our main analysis. The uncertainties in the phase-curve parameters are larger in this case because of the time-series segmentation, which does not constrain the baseline coefficients as effectively as for longer data sets. The phase curves obtained from these additional analyses are shown in Extended Data Fig. 8. We finally note that the phase-curve peak is located close to the start of half of our observations and towards the end of the other half data sets, which was necessary owing to Spitzer downlink limitations. We deem this pattern purely coincidental for two reasons. First, if our reported phase curve was due to uncorrected systematics, then it would be unlikely that the systematics would produce an upward trend in half of the data and a downward trend in the other half. These data sets are independent and there is no relationship between those obtained from transit to occultation and those obtained from occultation to transit. There is also no correlation with the centroid position on the detector. Second, if the phase-peak offset occurred after or before this discontinuity, then it would have been clearly detected in the continuous parts of our data set; however, only gradual slopes are seen in both data sets. A comparison with data obtained in the same year with the Microvariability and Oscillations of STars (MOST) satellite (D. Dragomir, personal communication) shows an agreement in the phase-curve amplitude and offset values derived from both facilities. The key features of the phase curve of 55 Cancri e translate directly into constraints on maps41, 42 assuming a tidally locked planet on a circular orbit. A planetary phase curve F /F (F is the flux from the star) measures the planetary hemisphere-averaged relative brightness 〈I 〉/〈I 〉 as in which α is the orbital phase, R is the planetary radius and R is the stellar radius. The longitudinal mapping technique used here19 aims to mitigate the degeneracy between the distribution of the planetary thermal brightness and the system parameters. This part of the analysis is independent of the light-curve analysis presented above. Therefore, here we fix the system parameters to those derived from a previous study7, which is based on the entire 55 Cancri e Spitzer data set. Using this prior information for the purpose of longitudinal mapping is adequate because the degeneracy between the planetary brightness distribution and the system parameters is only relevant in the context of eclipse mapping19. We follow the same approach as for Kepler-7b (ref. 20) and use two families of models, similar to the ‘beach-ball models’ introduced in ref. 21: one using n longitudinal bands with fixed positions on the dayside, and another using longitudinal bands whose positions and widths are jump parameters in the MCMC fit. We choose a three-fixed-band model and one-free-band model to extract both the longitudinal dependence of the dayside brightness of 55 Cancri e and the extent of its ‘bright’ area. Increasing n to five yields a larger BIC than for n = 3. For both models, we compute the amplitude of each band from their simulated light curve using a perturbed singular-value decomposition method. The one-free-band model (Fig. 2, left) yields a uniformly bright longitudinal area extending from 5° ± 18° west to 85° ± 18° east with a relative brightness of 0.72 ± 0.18, compared to a brightness of 0.15 ± 0.05 for the rest of the planet. The three-fixed-band model yields bands of relative brightness decreasing from the west to the east: <0.21 (3σ upper limit), 0.58 ± 0.15 and 0.74 ± 0.15, compared to the nightside contribution of 0.17 ± 0.06. Variability in the thermal emission of 55 Cancri e between 2012 and 2013, has previously been determined from occultation measurements7. Several tests regarding the robustness of the variability pattern were conducted, including three different analyses that used BLISS mapping, polynomial detrending and a pixel-level decorrelation method11. These three approaches confirmed the variability of the thermal emission of the planet between 2012 and 2013 with similar uncertainties. Therefore, we consider it very likely that the emission of the planet is varying, but on timescales that are substantially longer than the timescale of the 2013 observations alone (a month) used here. No variability is reported in the 2013 data alone7. These factors justify our combining of the 2013 observations and our use of a single phase-curve model. Furthermore, we detect the phase-curve shape in all individual data sets in addition to the combined phase-folded time series. This strengthens our conclusion that it is unlikely that stellar variability would cause the combined phase-curve shape from individual stellar events taken at different times over the month of observations. We use an observed infrared spectrum of 55 Cancri e (ref. 43) to compute the brightness temperatures in the IRAC 4.5-μm bandpass from the F /F values derived from the MCMC fits. If an atmosphere was present, then the large temperature contrast between the dayside and nightside hemispheres suggests that the radiative cooling time (t ) is less than the dynamical time scale (t ), resulting in a poor redistribution of heat from the dayside to the nightside. This sets a constraint on the mean molecular weight, which we may estimate. The zonal velocity is , in which is the specific gas constant, ΔT = 1,460 K is the temperature difference between the hemispheres and we have ignored an order-unity correction factor associated with the pressure difference between the hemispheres44. If we enforce t < t , then for the mean molecular weight μ we obtain in which is the universal gas constant, T = 2,700 K is the dayside temperature, σ is the Stefan–Boltzmann constant, R = 1.91R is the planetary radius (R is the Earth radius), κ = 2/7 is the adiabatic coefficient, g = 103.33 cm s−2 is the surface gravity and P is the dayside pressure—the only unknown parameter in this expression. If we set P = 1 bar, then μ > 9. This estimate further suggests that a hydrogen-dominated atmosphere is unlikely, and sets a lower limit on the mean molecular weight. It is unlikely that 55 Cancri e is harbouring a thick atmosphere, owing to its proximity to its star. If we assume energy-limited escape45, then the atmosphere needs to have sufficient mass to survive for the stellar age, which translates into a lower limit on the required surface pressure in which is the X-ray luminosity of the star, t is the stellar age, G is Newton’s gravitational constant, M = 8.08M is the planetary mass (M is the Earth mass) and a = 0.01544 au is the orbital semi-major axis. If (ref. 24) and t = 8 Gyr, then P > 31 kbar; in other words, the surface pressure of 55 Cancri e needs to be larger than 31 kbar to survive atmospheric escape over the stellar lifetime. Despite the uncertainties associated with estimating the mass loss due to atmospheric escape, this estimate is conservative because the star probably emitted higher X-ray luminosities in the past. Our suggestion of an atmosphereless 55 Cancri e is consistent with the trends predicted for super-Earths45. Hence, it is unlikely that the large infrared peak offset is due to an atmosphere rich with volatiles. It is more likely that the infrared phase curve of 55 Cancri e is probing non-uniformities associated with its molten rocky surface. Because 55 Cancri is a multi-planet system, the eccentricity and obliquity of 55 Cancri e are excited, owing to the presence of the outer planets. This creates a tidal heat flux that is responsible, in part, for the thermal emission of the planet. To evaluate the contribution of the tidal heat flux to the measured thermal emission, we investigate the possible values of the eccentricity and obliquity of 55 Cancri e for different tidal dissipation using N-body simulations (using Mercury-T; ref. 46). We use the orbital elements and masses for the four outer planets30 and the most recent values7 for the mass, radius and orbital semi-major axis of 55 Cancri e. We find that the obliquity of 55 Cancri e is very low (<1°) and that the eccentricity is about 10−3 for the eight orders of magnitude (10−5–10 times the dissipation of Earth σ ) we consider for the tidal dissipation of 55 Cancri e. The corresponding tidal heat flux ϕ or tidal temperature (ϕ /σ)1/4 increase with the dissipation in 55 Cancri e, from 10−3 W m−2 (a few kelvin) to 106 W m−2 (about 2,000 K). We calculate the occultation depth at 4.5-μm for a range of eccentricities and albedos (0.0–1.0) to enable a comparison with the output of the dynamical simulations (Extended Data Fig. 9). We find that a combination of large dissipation (10σ ), eccentricity and obliquity can explain the level of thermal emission observed in 2013; however, these solutions do not allow us to reproduce the nightside temperature. In our configuration (no heat re-distribution and assuming an isotropic tidal heat flux), tides do not match our measurements, so an additional heat source is probably responsible for at least part of the large planetary thermal emission observed in 2013. The code used to perform the aperture photometry on the Spitzer data sets presented here is publicly available from the IDL Astronomy User’s Library at http://idlastro.gsfc.nasa.gov. We have opted not to make the MCMC code available, but the corrected photometry for each data set is available online as Source Data.