A search for molecules in the atmosphere of HD 189733b (original) (raw)

Journal Article

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

Search for other works by this author on:

Accepted:

02 September 2009

Published:

17 December 2009

Cite

J. R. Barnes, Travis S. Barman, H. R. A. Jones, R. J. Barber, Brad M. S. Hansen, L. Prato, E. L. Rice, C. J. Leigh, A. Collier Cameron, D. J. Pinfield, A search for molecules in the atmosphere of HD 189733b, Monthly Notices of the Royal Astronomical Society, Volume 401, Issue 1, January 2010, Pages 445–454, https://doi.org/10.1111/j.1365-2966.2009.15654.x
Close

Navbar Search Filter Mobile Enter search term Search

Abstract

We use signal enhancement techniques and a matched filter analysis to search for the _K_-band spectroscopic absorption signature of the close orbiting extrasolar giant planet, HD 189733b. With time-series observations taken with the Near Infrared Spectrometer (NIRSPEC) at Keck II, we investigate the relative abundances of H2O and carbon bearing molecules, which have now been identified in the dayside spectrum of HD 189733b. We detect a candidate planet signature with a low level of significance, close to the ∼153 km s−1 velocity amplitude of HD 189733b. However, some systematic variations, mainly due to imperfect telluric line removal, remain in the residual spectral time series in which we search for the planetary signal. Using principal components analysis, the effects of this pattern noise may be reduced. Since a balance between the optimum systematic noise removal and minimum planetary signal attenuation must be struck, we find that residuals, which are able to give rise to candidate planet signatures, remain. The robustness of our candidate signature is therefore assessed, enabling us to conclude that it is not possible to confirm the presence of any planetary signal which appears at _F_p/F* contrasts deeper than the 95.4 per cent confidence level. Our search does not enable us to detect the planet at a contrast ratio of _F_p/F*= 1/1920 with 99.9 per cent confidence.

Finally, we investigate the effect of model uncertainties on our ability to reliably recover a planetary signal. The use of incorrect temperature, model opacity wavelengths and model temperature-pressure profiles have important consequences for the least squares deconvolution procedure that we use to boost the signal-to-noise ratio (S/N) in our spectral time-series observations. We find that mismatches between the empirical and model planetary spectrum may weaken the significance of a detection by ∼30–60 per cent, thereby potentially impairing our ability to recover a planetary signal with high confidence.

1 INTRODUCTION

The field of exoplanet spectroscopy has advanced rapidly since the first observations of secondary eclipse events enabled the brightness temperatures of transiting planets to be determined (Charbonneau et al. 2005; Deming et al. 2005), thereby confirming heating due to stellar irradiation. Thanks largely to the Spitzer Space Telescope (Fazio et al. 2004; Houck et al. 2004; Rieke et al. 2004), the spectral energy distribution of a number systems has now been estimated (e.g. see Burrows, Budaj & Hubeny 2008). It has recently been suggested that close orbiting extrasolar giant planets (CEGPs) can be divided into two possible subgroups, one which exhibits absorption spectra and one which exhibits emission features due to the presence of a stratosphere (Burrows et al. 2008; Fortney et al. 2008). Individual planetary atmospheres are no doubt more complex with Burrows et al. (2008), for example, preferring to use parametrizations of the degree of stratospheric absorption and heat redistribution. Nevertheless, HD 209458b appears to fall broadly into the latter category (Knutson et al. 2008) while the spectral energy distribution of HD 189733b is found to be well fit by models where no stratosphere forms (Barman 2008).

Although planetary signatures are easier to detect at mid-infrared wavelengths, characterization of atmospheres is easier at shorter wavelengths where observational sensitivities are sufficient to probe the large amplitude signatures of species such as H2O and carbon bearing molecules. The first evidence for water and organic molecules in the atmospheres of CEGPs has come from space-based broad-band photometric observations and low resolution transit spectroscopy. While Barman (2007) and Tinetti et al. (2007) claimed independent detections of H2O at different wavelengths, Swain, Vasisht & Tinetti (2008) have now identified CO and CH4 in transit spectra which probe the terminator of the planet HD 189733b. Grillmair et al. (2008) reported a detection of water in the dayside spectrum through Spitzer Space Telescope observations of secondary eclipse events below 7.5 μm. Most recently, data taken with the Near Infrared Camera and Multi-Object Spectrometer on board the Hubble Space Telescope (NICMOS/HST) have been used to infer the additional presence of CO and CO2 in the dayside spectra of HD 189733b (Swain et al. 2009, hereafter S09).

Rather than attempting to identify molecules from their broad-band spectral signatures, we present a method which attempts to identify spectral structure through a statistical examination of the many thousand individual transitions found in a typical high resolution (_R_∼ 25 000–50 000) planetary spectrum. Specifically, in order to detect the faint planetary signal, spectral deconvolution is applied to individual spectra in a time series in order to derive mean absorption profiles with boosted S/N ratios. Modelling the phase dependent contrast ratio and radial velocity motion enables the maximum planet/star contrast ratio and velocity amplitude to be measured (from which the true mass of the planet may also be determined). Crucially, since the planet need not be transiting, this spectroscopic method increases the sample of planets which can potentially be studied with current instrumentation. With the aim of extending space-based mid-infrared measurements of planet/star contrast ratios into the near infrared, we search for the CEGP signature, manifested as H2O, CO and CO2 absorption in high resolution _K_-band spectroscopic time-series observations of HD 189733 (K1V−K2V).

2 DATA REDUCTION

2.1 Observations

_K_-band observations of HD 189733 were secured with NIRSPEC (McLean et al. 1998) at the Keck II Telescope on ut 2008 June 15 and 22. A total of 219 and 154 spectra were recorded, respectively, using a 1024 × 1024 InSb Aladdin-3 array. With the NIRSPEC-7 blocking filter, a wavelength span of 2.0311–2.3809 μm was achieved with a slit width of 0.432 arcsec, giving a resolution of _R_∼ 25 000. Our 60 s exposures comprised of 12 coadds, each of 5 s duration. The observations are summarized in Table 1. The seeing was mostly good at around 0.6–0.7 arcsec although observations were plagued by cloud for a period on June 15.

Table 1

Keck/NIRSPEC observations of HD 189733 for ut 2008 June 15 and 22.

ut Date ut start of first frame ut start of last frame Time per exposure (s) Number of co-adds per frame Number of observations
2008 June 15 08:36:29 15🔞33 5 12 219
2008 June 22 10:43:37 15:08:46 5 12 154
ut Date ut start of first frame ut start of last frame Time per exposure (s) Number of co-adds per frame Number of observations
2008 June 15 08:36:29 15🔞33 5 12 219
2008 June 22 10:43:37 15:08:46 5 12 154

Table 1

Keck/NIRSPEC observations of HD 189733 for ut 2008 June 15 and 22.

ut Date ut start of first frame ut start of last frame Time per exposure (s) Number of co-adds per frame Number of observations
2008 June 15 08:36:29 15🔞33 5 12 219
2008 June 22 10:43:37 15:08:46 5 12 154
ut Date ut start of first frame ut start of last frame Time per exposure (s) Number of co-adds per frame Number of observations
2008 June 15 08:36:29 15🔞33 5 12 219
2008 June 22 10:43:37 15:08:46 5 12 154

2.2 Data extraction

Pixel-to-pixel variations were corrected for each frame using flat-field exposures taken with an internal tungsten reference lamp. The worst cosmic ray events were removed at the pre-extraction stage using the Starlink figaro routine bclean (Shortridge 1993). Since we chose not to use an ABBA nodding sequence, in order to maximize stability in the K band, we carried out the same extraction procedure as detailed in Barnes et al. (2007b) for previous observations of HD 189733. The spectra were extracted using echomop's implementation of the optimal extraction algorithm developed by Horne (1986). echomop rejects all but the strongest sky lines Barnes et al. (2007b) and propagates error information based on photon statistics and readout noise throughout the extraction process.

In Barnes et al. (2007b), we reported observations of HD 189733 which were made in conditions of variable cloud and seeing. This resulted in a spectral time series in which the S/N ratio of individual spectra varied greatly. We found that rejecting the lowest S/N frames increased the sensitivity of our method. We thus carried out a similar rejection procedure, set at an arbitrary S/N = 350 to reject the lowest S/N ratio frames which formed a distinct distribution separate from those frames which were observed in the best seeing conditions. A total of 63 frames were rejected from observations made on June 15 (S/N ratios of the rejected frames were 245 ± 127) while none was rejected from observations made on June 22; a clear reflection of the more variable conditions prevailing on the first night. The S/N ratios of the residual time series were 497 ± 125 and 570 ± 77 on June 15 and 22, respectively. The phases of observations used during the subsequent analysis stages are represented visually in Fig. 1.

HD 189733b orbital phase diagram. Observations made on June 15 cover phases φ= 0.303–0.429 while observations made on June 22 cover φ= 0.498–0.581. Observations in the range φ= 0.498–0.517 were not used in the subsequent analysis since the planet is eclipsed during these phases. Only the phases used for analysis are shown.

Figure 1

HD 189733b orbital phase diagram. Observations made on June 15 cover phases φ= 0.303–0.429 while observations made on June 22 cover φ= 0.498–0.581. Observations in the range φ= 0.498–0.517 were not used in the subsequent analysis since the planet is eclipsed during these phases. Only the phases used for analysis are shown.

2.3 Residual spectra

Following the procedures outlined in Barnes et al. (2007b) and Barnes et al. (2008), we attempt to extract the planetary signature from time-series spectra by removing the dominant spectral contributions; namely the stellar spectrum and the telluric lines. A master spectrum is generated by combining all observed frames on a given night after nearest pixel alignment of each individual order in each spectrum to minimize blurring. The master template is subtracted from each spectrum in turn after shifting, scaling and blurring/sharpening (again on an order to order basis). This latter procedure is achieved by calculating the first to fourth order derivatives of the template spectrum. The master template is then scaled to each observed spectrum by using a spline with a fixed number of knots. Scale factors for the derivatives are similarly calculated. A model of each observed spectrum can thus be calculated by means of a Taylor expansion of the template spectrum. The procedure is described in detail in appendix A of Collier Cameron et al. (2002). Since we observe HD 189733b at phases close to phase φ= 0.5 when the planet spectrum is Doppler shifted on the steepest part of the radial velocity curve, the master template spectrum contains only a very weak, blurred out copy of the planetary spectrum. Hence subtraction of the template does not significantly attenuate the planetary signature.

2.4 Residual pattern noise removal

The resulting residual time series should thus ideally only contain noise and a copy of the planetary spectrum. However, residual pattern noise remains in the subtracted time series. Examination indicates that the dominant patterns are time variant residuals which are introduced when subtracting the scaled master template during the procedure outlined in Section 2.3. Although this procedure is effective at removing the stellar and telluric lines, variations between the strength of these lines throughout the night (due to changing airmass and hence telluric strength) are sufficient to preclude consistent results over time-scales of the variations.

In addition to pattern noise induced by telluric effects, noise may arise from time variable fringing effects, as found by Brown, Libbrecht & Charbonneau (2002) for _K_-band NIRSPEC observations for example. To remove these two effects, Brown et al. (2002) carried out a two stage procedure in the form of a regression and singular value decomposition filtering procedure on their spectra. Our implementation of the Taylor expansion algorithm, mentioned above and in detail in Collier Cameron et al. (2002), is analogous to the regression step. In a similar procedure to the singular value decomposition implemented by Brown et al. (2002), the remaining pattern nose can be removed by calculating the correlation matrix of the time variations in individual spectral bins for our set of residual spectra. The eigenvectors of the correlation matrix, which account for the largest fraction of the variance, are the principal components of the residual spectra. By subtracting only the principal components, any remaining pattern noise can be further removed (see appendix B of Collier Cameron et al. 2002 for further details). A balance must be struck between removing pattern noise and maintaining information in the residual spectra. We found that shifts of 1–2 pixels during a typical night of observations could be attributed to points during the night where the instrumental configuration was slightly modified to enable observations of a star for another project. We thus applied principal components analysis to continuous blocks (i.e. between observations of the other star, where the 1–2 pixel shifts occurred) of HD 189733 observations and found that 2–3 components were necessary to remove the remaining pattern noise in the data without attenuating the planetary signal (see Section 4).

3 DECONVOLUTION AND PLANETARY MODELS

In order to extract the planetary signature from the corrected residual spectra, we used least squares deconvolution (Donati et al. 1997) which requires the use of a model spectrum to describe the strengths (normalized profile depths) and wavelength positions of the strongest planetary opacities over the wavelength range of our observations. Our implementation of the algorithm (Barnes et al. 1998) propagates errors from the input spectra and has been used in reflected light searches in the optical by Collier Cameron et al. (1999, 2002) and Leigh et al. (2003a,b). For each residual spectrum a deconvolved profile is obtained, potentially containing a copy of the Doppler shifted planetary profile which can then be detected owing to the effective boost in S/N ratio. Typical boosts in S/N ratio of a few times to a few tens of times are achieved since several hundred to several thousand planetary absorption lines are used to deconvolve the planetary signature. The 2.04–2.06 μm and 2.36–2.38 wavelength ranges are omitted in all the following analyses due to the strong telluric features which dominate these regions of the spectra.

We have generated several models to represent the emergent spectrum of HD 189733b. Our standard HD 189733b model, an updated version of the model detailed in section 3.2 of Barnes et al. (2007b), is generated for an atmosphere with solar metallicity, a temperature, _T_= 1250 K and surface gravity, log _g_= 1.33 ms−1. For a detailed description of the model opacities and setup see Ferguson et al. (2005), Barman, Hauschildt & Allard (2001) and Barman, Hauschildt & Allard (2005) BHA05. The most recent models and their ability to fit HD 189733b observations made with the Spitzer Space Telescope are discussed at length in Barman (2008). A number of further models with adjusted temperature pressure profiles and relative chemical abundances will be discussed in the following sections. The S/N ratio after deconvolution with the standard model yields time series with mean profile S/N ratios of 7520 and 9450 on June 15 and 22, respectively. The mean deconvolved wavelength of the spectra depends on the wavelength dependent count rate and the relative strengths of the absorption lines. For the standard model, λmean= 2.19 μm. As reported in Section 4, the mean wavelength may change slightly when the line list used for deconvolution is modified.

4 RESULTS

A matched Gaussian filter approach is used to search for the planetary signal, where a model describing the radial velocity shift and planet/star contrast phase function is used to search for the deconvolved absorption profile of the planet in the spectral time series (Collier Cameron et al. 2002; Barnes et al. 2007a,b, 2008). Pairs of maximum contrast ratio (ε0) versus velocity amplitude (_K_p) are used in a 2D χ2 search to find the combination which yields the best improvement in χ2. The significance of candidate enhancements in χ2 are assessed by a bootstrap procedure which randomizes the order of the data within each night of observations (Collier Cameron et al. 2002). This process, carried out several thousand times, scrambles any planetary signature, but enables the data to retain the ability to give false χ2 enhancements due to systematics which may remain in the data above the photon noise. Reliable confidence levels can thus be plotted on the 2-parameter χ2 landscapes of log10(ε0) versus (_K_p) χ2.

Calibration of contrast ratios is achieved by injecting a ‘fake’ model planet spectrum into the time series and then recovering the signal using our matched filter method. The fake planetary spectrum is injected with known ε0 and _K_p after extraction of the spectra and before any of the subsequent steps described above are carried out. In this way, we are also able to assess our ability to correctly recover a planetary signature (Barnes et al. 2007b, 2008). We find that for a planet recovered with high significance, there may be a slight shift in the recovered _K_p velocity. This is most likely due to some removal of planetary signature when subtracting the template star and during principal components analysis. During these procedures, we chose parameters which strike a balance between removing residuals while not significantly affecting the planet signal. Essentially, an absorption signature located at phases close to φ= 0.25 will be attenuated most, since, in this region the planet shows the smallest radial velocity gradient with orbital phase. The resulting magnitude of the velocity amplitude uncertainty is typically <5 km s−1 for a planet simulated with greater than 99.9 per cent significance but may be as much as 10–20 km s−1 where a fake planetary signature is injected with ∼95.4 per cent significance. This is especially true if the data contain systematics above the photon noise level of the data.

4.1 Standard model

In Fig. 2 (top left), we present the phased deconvolved time series of the residual spectra (i.e. spectra with removed stellar spectrum and tellurics and containing a potential planetary signature) based on our standard model (see Section 3). For plotting purposes only, the time series has been normalized using the formal variances since some phases (particularly φ= 0.372–0.429 at the end of the first night of observations) are noisier than others. This enables the noise structure at all phases to be more clearly seen. The black and white grey-scale values are set at ±1.5σ in the plotted normalized time series. It should be stressed that the true formal variances are utilized when searching for the planetary signal in the un-normalized deconvolved time series. Hence spectra with lower S/N receive a lower weighting in our analysis. Since the planet undergoes eclipse during phases φ= 0.498–0.517, we do no use spectra taken during this interval in our analysis. These spectra are however plotted for completeness in Fig. 2.

Phased deconvolved time-series spectra of HD 189733b and corresponding 2D χ2 landscape plots for analysis using the standard model (top) and enhanced CO2 model (bottom). Left-hand column: the dashed line in the time-series plots represents the position of the planet as a function of phase. The black and white levels are set at the 1.5σ level of the time series. Right-hand columns: χ2 plot for matched filter combinations of maximum contrast ratio, log10(ε0) versus Kp. Observations covering phases φ= 0.498–0.517 are shown for completeness but were not used in the analysis since the planet undergoes eclipse at these phases. Black and white represent the best and worst improvements in χ2, respectively. The recovered enhancements in χ2 are likely caused mainly by systematics in the residual time-series spectra. No enhancement in χ2 is measured at the known velocity amplitude, Kp= 152.6 km s−1, marked by the dashed line. Upper confidence levels of 63.8, 95.4, 99 and 99.9 per cent (top to bottom solid lines) are plotted. The bold + symbols mark the mean value of the contrast ratio over the range of observations for the standard model (top right) and the corresponding mean value of the contrast ratio from the NICMOS/HST S09 observations (bottom) and are log10(ε0) =−3.163 and −3.286 (Fp/F*= 1/1460 and 1/1930), respectively.

Figure 2

Phased deconvolved time-series spectra of HD 189733b and corresponding 2D χ2 landscape plots for analysis using the standard model (top) and enhanced CO2 model (bottom). Left-hand column: the dashed line in the time-series plots represents the position of the planet as a function of phase. The black and white levels are set at the 1.5σ level of the time series. Right-hand columns: χ2 plot for matched filter combinations of maximum contrast ratio, log10(ε0) versus _K_p. Observations covering phases φ= 0.498–0.517 are shown for completeness but were not used in the analysis since the planet undergoes eclipse at these phases. Black and white represent the best and worst improvements in χ2, respectively. The recovered enhancements in χ2 are likely caused mainly by systematics in the residual time-series spectra. No enhancement in χ2 is measured at the known velocity amplitude, _K_p= 152.6 km s−1, marked by the dashed line. Upper confidence levels of 63.8, 95.4, 99 and 99.9 per cent (top to bottom solid lines) are plotted. The bold + symbols mark the mean value of the contrast ratio over the range of observations for the standard model (top right) and the corresponding mean value of the contrast ratio from the NICMOS/HST S09 observations (bottom) and are log10(ε0) =−3.163 and −3.286 (_F_p/F*= 1/1460 and 1/1930), respectively.

In Fig. 2 (top right) we present the χ2 landscape plot of log(ε0) versus _K_p based on our standard model. Dark features in the plot represent enhancements in χ2 whose significance can be measured relative to the plotted confidence levels. The large black feature representing the greatest enhancement in χ2 appears with low confidence (in the 68.3–95.4 per cent confidence region) at log (ε0) =−3.41 and _K_p= 85 km s−1. However, as can be seen in the phased time series, a number of low-level residual features are present. These appear as dark absorption areas, covering localized regions of velocity and phase. We believe that these features are responsible for the _K_p= 85 km s−1 signature and result from imperfect removal of telluric and stellar lines during our analysis, giving rise to false signals. As has been demonstrated previously (Barnes et al. 2008), a clear detection of the planetary signature would be expected to result in a more localized χ2 enhancement and greater significance than the _K_p= 85 km s−1χ2 enhancement. Nevertheless, _K_p is known for HD 189733b (since the system is eclipsing and the orbital inclination is known) and is indicated by the dashed lines in Fig. 2. A candidate planetary signature should thus appear at, or close to (see above), this velocity amplitude.

For HD 189733b, the standard model predicts a maximum contrast ratio of log10(ε0) =−3.163 or _F_p/F*= 1/1460 over the wavelength span of our observations. We are however unable to detect the planetary signature, at the mean deconvolved wavelength of 2.19 μm, with 68.3, 95.4, 99 and 99.9 per cent confidence levels of log10(ε0) =−4.065, −3.491, −3.366 and −3.193 or _F_p/F*= 1/11600, 1/3100, 1/2320 and 1/1560, respectively. In light of much better observing conditions, this is a significantly more sensitive result than our 2006 observations (Barnes et al. 2007b) permitted. Considerable care must be exercised if quoting sensitivities at contrasts ratios deeper than the 95.4 per cent level (this is investigated further in Section 4.3) owing to candidate signatures which arise from systematics at these levels. HD 189733b is therefore not detected at _K_p= 152.6 km s−1 at a contrast which is 2.1 and 1.1 times deeper (95.4 per and 99.9 per cent confidence, respectively) than the standard model predicts.

4.2 Enhanced CO2 model

Recent Spitzer/IRAC (Grillmair et al. 2008) and HST/NICMOS (S09) observations have enabled the dayside spectrum of HD 189733b to be measured using low resolution spectroscopy. In order to reliably fit the spectrum, a greater than expected abundance of CO2 is required with S09 reporting CO2 mixing ratios of _c_∼ 0.1–1 × 10−6. We have generated a model with augmented CO2 abundance which enables us to match the HST/NICMOS observations (Fig. 3). From here on, in Section 4, we only consider this model. We emphasize that our current model does not contain some of the hot CO2 bands which have been identified by Fourier Transform Spectroscopy carried out at the Jet Propulsion Laboratory and included in the latest edition of the High Resolution Transmission database (HITRAN) (see Rothman et al. 2009 and references therein). We therefore caution that in order to achieve the required level of absorption from CO2 in the 1.9–2.2 μm, the relative strengths of individual opacities in the output model are likely overestimated. Nevertheless, based on the S09 results, inclusion of such opacities would appear to give a more accurate representation of the expected opacities found in the spectrum of HD 189733b. Fig. 2 (bottom left and bottom right) presents the deconvolved time series and log(ε0) versus _K_pχ2 plot after deconvolution using our augmented CO2 model spectrum. At _K_p= 152.6 km s−1, the expected planetary signature is not detected with 68.3, 95.4, 99 and 99.9 per cent confidence levels of log10(ε0) =−4.074, −3.529, −3.400 and −3.283 or _F_p/F*= 1/11600, 1/3380, 1/2510 and 1/1920, respectively. The sensitivities are slightly greater than for the standard model, although we note that the enhanced CO2 sensitivities are quoted for a centroidal wavelength of 2.15 μm rather than 2.19 μm. This shift in mean wavelength results from greater normalized depths of the enhanced CO2 at the shorter wavelengths of our observations. In addition, although the mean planetary flux level is lower in the regions with enhanced CO2, the recorded count rate in the observed star + planet spectra is higher, leading to higher contrast confidence limits with this model. The mean 2.0–2.4 μm planet/star flux ratio reported by S09 is log10(ε0) =−3.286 or _F_p/F*= 1/1930 (i.e. almost identical to our 99.9 per cent confidence level of 1/1920), indicating that we are sensitive down to the 99.9 per cent level.

Comparison of the standard model (black) and augmented CO2 abundance model (grey) with the observed NICOMS/HST dayside spectrum of S09. The enhanced CO2 model shows a considerably improved fit in the 1.9–2.2 μm region relative to the standard model.

Figure 3

Comparison of the standard model (black) and augmented CO2 abundance model (grey) with the observed NICOMS/HST dayside spectrum of S09. The enhanced CO2 model shows a considerably improved fit in the 1.9–2.2 μm region relative to the standard model.

4.3 Wavelength splitting the data – a candidate signature?

While there is no clear candidate signature at the expected velocity amplitude of the planet in the enhanced CO2 deconvolved time series, we do detect a signal with relatively low confidence (95.4 per cent) at _K_p= 167.9 km s−1 and log10(ε0) =−3.63 (_F_p/F*= 1/4270). Although it is likely that this signature is again the result of alignment of systematic absorption features in the time series at this velocity amplitude, we have investigated splitting the time-series data into two spectral regions. Deconvolution was carried out on the first three orders (region 1: 2.03–2.18 μm) and on the second three orders (region 2: 2.21–2.36 μm) independently before carrying out a search for the planetary signature. With the data split in this manner, region 1 contains H2O and CO2 opacities while region 2 contains H2O and CO (bandhead at ∼2.29 μm) opacities. While region 1 did not reveal any candidate signature close to _K_p= 152.6 km s−1, region 2 has a more well defined candidate signature at _K_p= 147.8 km s−1 with log10(ε0) =−3.449 (1/2810) and 97.2 per cent confidence. The S09 contrast ratio for HD 189733b over the wavelength span of region 2 is (_F_p/F*∼ 1/1340). To assess the true nature of the region 2 signature, we refer the reader to Fig. 4 (left-hand column) which again indicates that there are a number of systematic features. Close examination reveals that the expected velocity position of the planet as a function of phase (indicated by the dashed line) appears to pass through, or near, a number of contiguous absorption regions. To investigate the contribution of these regions to the candidate signature, we have carried out three tests.

As for Fig. 2, but using only the wavelength range 2.21–2.36 μm. The corresponding NICMOS/HST S09 contrast ratio for this region is marked by a + symbol with Fp/F*∼ 1/1340. A candidate signature is detected with 97.2 per cent confidence close to the expected Kp velocity, although the contrast ratio of log10(ε0) =−3.449(Fp/F*= 1/2810) is lower than expected. See main text for details.

Figure 4

As for Fig. 2, but using only the wavelength range 2.21–2.36 μm. The corresponding NICMOS/HST S09 contrast ratio for this region is marked by a + symbol with _F_p/F*∼ 1/1340. A candidate signature is detected with 97.2 per cent confidence close to the expected _K_p velocity, although the contrast ratio of log10(ε0) =−3.449(_F_p/F*= 1/2810) is lower than expected. See main text for details.

  1. Analysis of the data observed on 15th June –φ= 0.303–0.429 alone. The result is a candidate signature with _K_p= 153.9 km s−1 and with log10(ε0) =−3.008 (1/1020) and 98.8 per cent confidence.
  2. Analysis of the data observed on 22nd June –φ= 0.517–0.581 alone. A candidate signature with _K_p= 141.8 km s−1 and with log10(ε0) =−3.528 (1/3372) and 94.8 per cent confidence is found.
  3. Contiguous dark regions omitted from the analysis by eye (the regions are not of sufficient amplitude to enable reliable sigma-clipping). No candidate signature within 26 km s−1 of the known _K_p= 152.6 km s−1 is apparent.

The varying velocity and contrast ratio of the candidate signals from tests 1 and 2 suggest that if any planetary signature contributes to the χ2 enhancements, it is biased by some other factor. Test 3, in which the contiguous regions are omitted from the analysis, has the effect of removing the candidate signature seen in Fig. 4 (right-hand column) completely. The confidence levels, at _K_p= 152.6 km s−1, with the omitted contiguous residual absorption regions are log10(ε0) =−3.964, −3.355, −3.214 and −3.098 or _F_p/F*= 1/9200, 1/2260, 1/1640 and 1/1250, respectively. In other words, the phases which do not show contiguous blocks of absorption residuals in the time series (66 per cent of the recorded spectra) along the radial velocity path of the planet do not possess the ability to recover a planetary signal.

We stress that the argument asserting that the contiguous residual absorption regions are wholly due to systematics and solely responsible for producing candidate signatures is however not strictly true. Any residual absorption features in the time series have the ability to modify the contrast ratio and velocity amplitude of a true planetary signal. Since the residual absorption features may be expected to vary in strength it is not unlikely that they would result in a planetary signature modified by differing degrees in tests 1 and 2. In our third test, removing 33 per cent of the data along the expected radial velocity curve of the planet leaves only regions which are consistent with the mean level or regions of contiguous ‘emission’ relative to the mean level. One might expect that this procedure would severely impair our ability to detect a planetary absorption signature. In the hypothesis that dark regions are artefacts of the data processing (i.e. imperfect telluric/stellar line removal), after their removal (test 3) we can rule out our ability to detect the planetary signal with 95.4 per cent confidence at a level of _F_p/F*= 1/2260. In light of this and our 99.9 per cent upper limit (Section 4.2) on the contrast ratio, we believe that since S09 detect the planet with _F_p/F*= 1/1930, further investigations of model dependency on our analysis are required.

5 MODEL DEPENDENCY EFFECTS

We are confident that the planetary signature is not severely attenuated (there is inevitably some attenuation as described in Section 4) during our analysis procedure since fake planetary signals which are injected before analysis are recovered. A cause of our inability to detect a planetary signature is likely to stem from a mismatch between the model planetary spectrum and observed planetary spectrum. The most likely direct causes of line strength mismatch and model line wavelength opacities were first highlighted in Barnes et al. (2007a). Line strength mismatches may arise from incorrect treatment of the model atmosphere, including uncertainty in the exact form of the temperature–pressure (T–P) profile. In addition, the precision of the calculated opacities is limited by the accuracy of the Einstein A coefficients. This latter effect may be true for important molecular species such as H2O (Barber et al. 2006) for instance. We investigate relative line depth, temperature and wavelength uncertainties below.

5.1 Relative line depths

Although the model planetary spectrum may show little variation as a result of T–P profile changes when observed at low resolution, the relative line depths may change significantly. In addition, the gradient of the T–P profile determines the absolute strength of the absorption lines. To investigate these effects, we have generated a series of ad hoc T–P profiles and resulting emergent spectra. Fig. 5 shows the models plotted for a short region of wavelength space. Steeper T–P gradients lead to the formation of deeper lines whereas the relative line strengths vary from model to model. These effects are important since mismatch of the model and observed spectra line depths will lead to a non-optimally deconvolved line profile, and therefore decrease in sensitivity. If all the lines are less deep, they will simply be harder to detect above the noise level.

Ad hoc T–P profiles and resulting spectra for 2.135–2.145 μm region. The black profile and spectrum represents the standard, standard model. While a steeper temperature profile at pressures where lines predominantly form (e.g. red/blue models) result in spectra with deeper absorption lines, note that the relative line strengths also differ from model to model.

Figure 5

Ad hoc T–P profiles and resulting spectra for 2.135–2.145 μm region. The black profile and spectrum represents the standard, standard model. While a steeper temperature profile at pressures where lines predominantly form (e.g. red/blue models) result in spectra with deeper absorption lines, note that the relative line strengths also differ from model to model.

To investigate the degree to which our ad hoc models affect the sensitivity of our procedure we used our standard model (black in Fig. 5) to inject a fake planet into the HD 189733 time series. The fake planet was then recovered and calibrated to match the contrast ratio at which it was injected. By deconvolving with line lists derived from each of the different models shown in Fig. 5, we find that we are able to recover the planetary signature in all cases. The planetary signature is however recovered with an incorrect contrast ratio and modified relative confidence. Fig. 6 shows the relative contrast ratio for a fake HD 189733b planet injected into the data with 99.9 per cent significance. Model 0 represents the standard model calibration to which the simulations are normalized. The contrast ratio is incorrectly recovered, with model 1 showing a 2.3 per cent overestimation of the contrast level and model 4 indicating a 66 per cent underestimation of the contrast ratio. In all cases however, the ad hoc models appear to recover the planet with increased significance. The effect is nevertheless relatively small, with models 2 and 3 showing the greatest increase in confidence. Model 2 indicates an increase in significance of 12.5 per cent relative to the 99 per cent confidence level. One might naturally expect a decrease in confidence to arise from mismatch of the line strengths during deconvolution rather than the counter-intuitive increase. We believe that the increase is most likely due to models 1–4 yielding strong lines which become even stronger and weak lines which become weaker relative to the standard model. If one of these models were a closer match to the empirical HD 189733b spectrum, we note that the relative significance of the standard model would decrease (with a maximum reduction in sensitivity of 14.8 per cent relative to the 99 per cent confidence level). In conclusion, the above ad hoc models alone are not able to explain the lack of true planetary signal (our observations are after all still sensitive enough to detect HD 189733b) through mismatch of line strengths. Relative to the confidence levels in Fig. 6, a true planetary signal could not change its confidence by more than half of the separation of the 99 and 99.9 per cent confidence levels.

Relative contrast ratio versus ad hoc deconvolution model. The standard model is represented by model 0. A fake planet was injected into the HD 189733b time series with 99.9 per cent confidence and recovered using our procedure. For each of models 1–4, the deconvolution step was carried out using the corresponding model. The solid line represents the recovered planetary signature relative to the standard model 0. The dashed lines are the 95.4, 99 and 99.9 per cent confidence levels (bottom to top), respectively. Note how the contrast ratio is incorrectly recovered with models 1–4, while the relative significance increases relative to the standard model.

Figure 6

Relative contrast ratio versus ad hoc deconvolution model. The standard model is represented by model 0. A fake planet was injected into the HD 189733b time series with 99.9 per cent confidence and recovered using our procedure. For each of models 1–4, the deconvolution step was carried out using the corresponding model. The solid line represents the recovered planetary signature relative to the standard model 0. The dashed lines are the 95.4, 99 and 99.9 per cent confidence levels (bottom to top), respectively. Note how the contrast ratio is incorrectly recovered with models 1–4, while the relative significance increases relative to the standard model.

5.2 Opacity wavelength uncertainties

Moderate wavelength uncertainties lead to an effective degradation of deconvolved resolution while model temperature uncertainties may lead to line strength mismatches with the observed spectra during deconvolution. We were able to use an unpublished improved version of the BT2 (Barber et al. 2006) water line list (with reduced wavelength uncertainties) to investigate these effects. Most of the strong lines which contribute to the deconvolution are transitions between states whose energies are experimentally known to very high accuracies. Consequently, by selecting only those lines which are greater in strength than 1/10 000 of the strongest line, we eliminate a large number of lines which are expected to have larger uncertainties in calculated positions and which in any case have negligible contribution to the deconvolved profile. Where the energies of the upper and lower states of a transition are both known experimentally, these are used for modified BT2 line frequencies, rather than the ab initio calculated values. At 1250 K, 75 per cent of the water lines in our trimmed list are transitions between experimentally known levels, and only 25 per cent employ BT2 ab initio frequencies/wavelengths (in all cases, however, BT2 Einstein A values are used in computing line strengths as these are generally more accurate than experimentally determined values).

Barber et al. (2006) state that a comparison of the BT2 ab initio frequencies with experimentally known transitions shows that the positions of ∼40 per cent of the lines tested are accurate to within 0.1 cm−1 and 91 per cent are within 0.3 cm−1. At 2.2 μm, this corresponds to resolutions of _R_= 45 500 and 15 150, respectively. Clearly, since water is the dominating opacity, these uncertainties will play a role in degrading the resolution of a deconvolved profile for data sets with resolutions of _R_≥ 15 000. These uncertainties are therefore applied to the 25 per cent of ab initio lines in our 1250 K list by using a Gaussian random uncertainty. This should represent a worst case scenario because we have removed those lines which are weaker than 1/10 000 of the strongest line and which are expected to exhibit the largest frequency/wavelength errors. We then carry out simulations by injecting a planetary signature into the HD 189733b spectra using our 1250 K spectrum and deconvolving firstly with the matching line list (case A), and the with an adjusted line list which models the wavelength uncertainties (case B). The mismatch (i.e. case B relative to case A) leads to a planet which is detected with a 6.5 per cent underestimation of contrast ratio and a 14.5 per cent decrease in significance relative to the 99.9 per cent confidence level.

5.3 Temperature uncertainty

The effect of using a model line list for deconvolution which varies from the observed spectrum in temperature alone is shown in Fig. 7. Here, a 1250 K planetary spectrum signature is recovered with 99.9 per cent confidence using a 1250 K modified BT2 line list. However both the recovered contrast ratio and significance change when deconvolved with line list temperatures which differ by ±250 and ±500 K from 1250 K. The effect is again relatively small for an underestimation of temperature (the increase in significance is likely due to over-weighting of strong lines and under-weighting of weak lines as described in Section 5.1) while slightly more significant for overestimation of temperature. In all instances a planetary signature is however recovered. For HD 189733b, the above effects alone are not sufficient to explain the lack of planetary signature (using our standard model) which is predicted at the log10(ε0) =−3.163 level. Combined wavelength and temperature mismatches should lead to a 99.9 per cent planetary signature appearing with 95.4 per cent confidence at worst. A ±250 K mismatch in model spectrum temperature results in a 20.5 per cent relative uncertainty in the confidence of a recovered signal.

Significance of recovered planet as a function of temperature. A planet with 99.9 per cent significance is simulated with a temperature of 1250 K and recovered with 99.9 per cent significance. By deconvolving the spectra with a range of temperatures, the recovered contrast ratio and significance of the planet are seen to vary.

Figure 7

Significance of recovered planet as a function of temperature. A planet with 99.9 per cent significance is simulated with a temperature of 1250 K and recovered with 99.9 per cent significance. By deconvolving the spectra with a range of temperatures, the recovered contrast ratio and significance of the planet are seen to vary.

5.4 Other sources of uncertainty

Additional model and observational uncertainties may contribute to an incorrect estimate of the planet/star contrast ratio or its relative significance. A possible source of error may arise from the planet ephemeris although this has been determined to high precision. Following our previous study of the HD 189733 system (Barnes et al. 2007b), we adopt the ephemeris of Winn et al. (2007), determined from Stromgen b and y passband observations [_Tt_= 2453988.80336(23) + 2.2185733(19)]. A more recent estimate of the ephemeris by Agol et al. (2009) using Spitzer 8 μm observations of planetary transit yields _Tt_= 2454279.436741(23) + 2.21857503(37). Since limb darkening and starspot effects are reduced at longer wavelengths, this has been claimed as the most precise measurement of the ephemeris to date. The predicted mid-transit time for our observations differs by 73 s when comparing the two ephemerides, a phase difference of 0.00038. The level of precision of ephemeris observations from HD 189733b is therefore now sufficient that new refinements have no measurable effect on the contrast or velocity amplitude of a planetary signal.

A more important consideration arises from global re-circulation patterns in the atmosphere of HD 189733b. For a tidally locked planet with a static atmosphere, one might expect the maximum planet/star contrast ratio to occur at orbital phase φ= 0.5 due to the highest effective irradiation of the planet at the substellar point. Knutson et al. (2007) found a difference in day and night side brightness temperatures of 238 K from 8 μm photometric Spitzer light curve variations. The 1212 ± 11 K dayside temperature was found to be displaced from the substellar point by 16°± 6°. This finding is in accordance with 3D circulation models Showman & Guillot (2002), Cooper & Showman (2005), Fortney et al. (2006) and Showman et al. (2008). More recently, Knutson et al. (2009) have re-analysed their 8 μm light curves and published 24 μm light curves of HD 189733b. Maximum brightness is found to occur at phase 0.396 ± 0.022 corresponding to a shift eastward of 20°–30° of the hottest region relative to the substellar point. We have carried out a simulation to estimate the effect of such a shift which is not accounted for in the preceding analysis. We created an artificial planetary signal which peaked 30° before secondary eclipse and recovered with a phase function which peaked at the same shift and also at secondary eclipse. The recovered planetary signature which did not account for the 30° shift was found to overestimate the contrast ratio by ∼10 per cent. The relative change in significance increases by 4.7 per cent since the contrast ratio must be increased to optimize the χ2 fit to the mis-aligned phase function. We note that this effect will be dependent on observational phase coverage and S/N ratio from night to night (i.e. shifting the phase function peak to a region of fewer observations or lower S/N ratio will reduce sensitivity).

We have assumed an effective v sin i for HD 189733b of 2.53 km s−1 which corresponds to a tidally locked planet. However there may be additional broadening as a result of the redistribution of heat. Showman et al. (2008) find that up to 3–4 km s−1 wind speeds are responsible for the advection of heat away from the substellar point. This shift is somewhat less than our resolution element of 11.99 km s−1. Although the wind speeds are effectively translational (an east-west flow) at the 100–1000 mbar levels from which the 2.2 μm spectrum is expected to predominantly arise (see fig. 4 of Showman et al. 2008), we have simulated an additional 4 km s−1 broadening of the spectral lines. Combined in quadrature with the rotational broadening, we simulate a planetary atmosphere which possesses lines broadened by 4.73 km s−1 rather than 2.53 km s−1 from rotational broadening alone as in the preceding sections. As expected this effect is also minor at a resolution of 25 000 with a 9.6 per cent drop in sensitivity.

5.5 A Semi-empirical approach – an L dwarf spectrum

In addition to model uncertainties, we have carried out a semi-empirical examination of our ability to recover the spectrum of a brown dwarf which closely matches the planetary temperature of 1250 K. _K_-band observations of an L3.5–L4.5 spectrum (Kirkpatrick et al. 2000; Knapp et al. 2004) were secured by with NIRSPEC at a spectral resolution of 22 000, covering redder wavelengths in each order than the HD 189733b observations. We were thus unable to use the L spectrum as a template which could be injected into our time series to mimic the signature of a fake ‘planet’. Instead, using the Taylor expansion scaling technique described in Section 2.3, we scaled the standard HD 189733b model spectrum to give the closest possible match to our observed L spectrum. Being an L dwarf, our spectrum exhibits significant rotation, with v sin _i_= 32 km s−1 (Zapatero Osorio et al. 2006; Reiners & Basri 2008). The same broadening was applied to our standard model prior to scaling it to the L spectrum. Deconvolution was then carried out using: (a) the standard model line list on the scaled standard model spectrum and, (b) the standard model line list on the L spectrum. Closely matched deconvolved profiles are recovered in both instances but with a smaller equivalent width for case b. Since v sin i is matched, the resulting profiles essentially differ in their depths only, with the case a profile being 55 per cent deeper than case b profile. It is difficult to assess wavelength mismatch effects given the broad nature of the profile; however, we can attribute the 55 per cent decrease in profile strength to line strength mismatches. Although there may be differences between a L spectrum and a planetary spectrum, this semi-empirical approach may be taken to represent an upper limit to our line depth sensitivity. The line depth, wavelength and temperature uncertainties in Sections 5.1, 5.2 and 5.3 yield a 28 per cent reduction in sensitivity when combined in quadrature. The semi-empirical analysis result may be equated with this combination of effects, and is almost twice the modelling estimate.

6 SUMMARY AND DISCUSSION

We have carried out a high-resolution search for the signature of the close orbiting extrasolar giant planet, HD 189733b. Our signal enhancement technique enables us to achieve the sensitivities required to detect the dayside spectrum of the planet that has already been observed at a mean contrast ratio of _F_p/F*= 1/1930 by S09 in the _K_-band region of our observations. Inclusion of augmented CO2 abundance is however not sufficient to detect the planet with a 99.9 per cent confidence level of _F_p/F*= 1/1920 (i.e. almost identical to the S09 result). A tentative candidate planetary signature is found at 15 km s−1 greater than the expected velocity amplitude of the planet at _F_p/F*= 1/4270. In light of the model uncertainties that have been investigated, finding a planetary signature with modified contrast ratio and velocity amplitude is reasonable. This prompted us to perform simulations in which planetary signals were injected at contrast ratio levels equivalent to those induced by contiguous absorption residuals. While these planetary signatures could be recovered, we found that the velocity amplitude may be uncertain by ±20 km s−1, further reflecting the difficulty of reliably extracting a real signal at the 95.4 per cent level. We note however that a planetary signature with 99.9 per cent confidence should easily be detected, as demonstrated in Barnes et al. (2008).

Splitting the data into two wavelength regions revealed that the 2.21–2.36 μm region (containing mainly H2O and CO opacities) yielded a candidate planetary signature with higher confidence. Analysing these subsetted data on a night by night basis however revealed that the signature was not stable in velocity amplitude or contrast ratio suggesting that it could result from a chance alignment of a number of systematic contiguous absorption residual features at the phase dependent velocity position of the planet. By removing these features, we found that the signature, close to the known _K_p= 152.6 km s−1, disappeared. The remaining 66 per cent of the data did not possess the ability to recover a planetary signature at the level determined by the results of S09 with between 95.4 and 99 per cent confidence. Since the remaining data contained contiguous regions with levels above the mean, this may not be surprising as we only search for absorption signatures. It is important to emphasize that the tests we have carried out do not rule out the possibility that a true planetary signal is contained within the spectra. The detected candidate features may be partially influenced by a true planetary signature, but at the 95.4 per cent levels, no confident claim for a detection can be made.

The effects of model opacity strength uncertainties, wavelength uncertainties, temperature mismatch, phase function mismatch and velocity field/broadening uncertainties contribute sensitivity uncertainties of 12.5, 14.5, 20.5, 4.7 and 9.6 per cent, respectively. Combining these effects in quadrature yields a total uncertainty in the significance of the result of 30 per cent. Further, if we take the semi-empirical 55 per cent uncertainty as an upper limit to our line depth, wavelength and temperature mismatches, the corresponding uncertainty is 56 per cent. Assuming that the L spectrum can provide a close match to that of HD 189733b, the semi-empirical result already shows that the model uncertainties may be significantly underestimated. Hence the 99.9 per cent confidence with which we reject a signal at the know _K_p could in fact be modified to a level with reduced significance, taking a candidate signal to contrast ratios that are plagued by systematic features.

While the current generation of models can adequately fit broad-band photometric and low resolution spectroscopic observations, it is clear that moving to higher resolution requires further model refinement. With the uncertainties investigated above, we cannot rule out the presence of the planet using our technique, especially if further model uncertainties remain unaccounted for. Only further observations which would bring about an increase in sensitivity, or more precise model atmospheres could increase our chances of detecting HD 189733b.

JRB was supported by a STFC funded research grant during the course of this work. TB acknowledges support from NASA's Origins of Solar System program and the NASA Advanced Supercomputing facility, and LP from NSF grant 04-44017. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. We would like to thank the referee for providing constructive suggestions for improving the manuscript.

REFERENCES

,

2009

, in

,

,

, eds,

Proc. IAU Symp. 253, Transiting Planets

.

Kluwer

, Dordrecht, p.

209

,

2007a

,

MNRAS

,

379

,

1097

et al.,

2005

,

ApJ

,

626

,

523

et al.,

2004

,

ApJS

,

154

,

10

et al.,

2008

,

Nat

,

456

,

767

et al.,

2004

,

ApJS

,

154

,

18

et al.,

2000

,

AJ

,

120

,

447

et al.,

2004

,

AJ

,

127

,

3553

et al.,

2007

,

Nat

,

447

,

183

et al.,

2009

,

ApJ

,

690

,

822

,

2003a

,

MNRAS

,

344

,

1271

et al.,

1998

,

Proc. SPIE

,

3354

,

566

et al.,

2004

,

ApJS

,

154

,

25

et al.,

2009

,

J. Quant. Spectrosc. Radiat. Transfer

,

110

,

533

,

1993

, in

,

,

, eds,

ASP Conf. Ser. Vol. 52, Astronomical Data Analysis Software and Systems II. Astron. Soc. Pac.

, San Francisco, p.

219

,

2009

,

ApJ

,

690

,

L114

(S09)

et al.,

2007

,

Nat

,

448

,

169

et al.,

2007

,

AJ

,

133

,

1828

© 2009 The Authors. Journal compilation © 2009 RAS

Citations

Views

Altmetric

Metrics

Total Views 402

241 Pageviews

161 PDF Downloads

Since 3/1/2017

Month: Total Views:
March 2017 1
April 2017 1
June 2017 3
July 2017 1
August 2017 1
September 2017 3
October 2017 1
November 2017 3
December 2017 2
January 2018 3
February 2018 5
March 2018 5
April 2018 3
May 2018 7
June 2018 9
July 2018 6
August 2018 7
September 2018 7
October 2018 3
November 2018 5
December 2018 3
January 2019 2
February 2019 6
March 2019 5
April 2019 10
May 2019 8
June 2019 3
July 2019 5
August 2019 6
September 2019 7
October 2019 8
November 2019 8
December 2019 6
January 2020 8
February 2020 4
March 2020 3
April 2020 2
May 2020 2
June 2020 6
July 2020 5
August 2020 2
September 2020 5
October 2020 3
November 2020 3
December 2020 5
January 2021 1
February 2021 4
March 2021 4
April 2021 3
May 2021 1
June 2021 4
July 2021 7
August 2021 4
September 2021 2
October 2021 4
November 2021 5
January 2022 3
February 2022 6
March 2022 1
April 2022 8
June 2022 2
July 2022 3
August 2022 3
September 2022 8
October 2022 16
November 2022 4
December 2022 2
January 2023 5
February 2023 3
April 2023 5
June 2023 1
August 2023 4
September 2023 2
November 2023 5
December 2023 6
January 2024 4
February 2024 8
March 2024 4
April 2024 8
May 2024 12
June 2024 10
July 2024 14
August 2024 6
September 2024 5
October 2024 4

Citations

19 Web of Science

×

Email alerts

Citing articles via

More from Oxford Academic