Neutral Metals in the Atmosphere of HD 149026b

We analyzed the archival data for transmission spectroscopy of HD 149026, which were taken from the SMOKA system (Baba et al. 2002). The data were obtained on the night of 2009 May 11 with the High Dispersion Spectrograph (HDS; Noguchi et al. 2002) on the Subaru telescope. A total of 41 frames were observed using the standard Ra (StdRa) mode with no iodine cell. The exposure time per frame ranged from 480 to 720 s, and the mean was ∼550 s. The total exposure was 22,680 s, and the exposure during the transit was 10,320 s. The slit width was 0.2 mm, which yields a wavelength resolution of R ∼ 90,000. The spectra were taken by two (blue and red) CCDs, which contained 26 and 17 orders and covered the wavelength ranges 4923–6227 and 6340–7660 Å, respectively.

We conducted a standard reduction of the data with IRAF tools and hdsql, provided by the HDS team. ¹⁰ First, we performed an overscan correction, removed the bias, and converted the analog digital unit to electron numbers. We made mask frames of the bad columns for each CCD. These processes were performed by hdsql. Scattered light was then subtracted by the IRAF task apscatter. Using 50 flat frames, we generated a median flat frame with a nonlinearity correction (Tajitsu et al. 2010). To remove the fringe pattern on the spectrum, we normalized the flat frame by apnormalize and divided the object frames by the normalized flat frame. Then we extracted a 1D spectrum by apall. The grid of the wavelength was derived by the IRAF task ecidentify using thorium–argon arc lamp spectra, which were taken at the beginning and end of the night. We assigned the wavelength grid to each object frame by refspectra and dispcor. We conducted these reductions to the median flat frame, as well as the object frames, to estimate the blaze function of each order. We divided the object spectra by the estimated blaze function using sarith. The spectrum contains several regions affected by bad columns that were not completely removed by the reduction with the CL script. Because these bad regions can be problematic for analysis, we identified them by visual inspection and masked them.

As described in previous research, the blaze function of HDS varies during an observation (e.g., Winn et al. 2004; Narita et al. 2005). To correct for this variation, we conducted the following procedure.

First, we divided the spectrum of each frame by that of a reference frame. We used the 21st frame as the reference frame because its signal-to-noise ratio (S/N) was the highest. We took the raw ratio of each spectrum and the reference spectrum and its change during the observation (Figure 1). Then we performed a 3σ clipping and smoothing of the ratio to remove the effects of residual absorption lines. After that, we fitted the raw ratio with the Chebyshev 15th-order function. The fitted ratio determined in this way represents the variation of the blaze function for the frame compared with the reference frame. We multiplied the reference spectrum by the raw ratio and divided it by the fitted ratio, which was a spectrum with the same blaze function as the reference spectrum. Finally, we derived the continuum spectrum of the reference spectrum and normalized it with the IRAF task continuum. Because the spectra of all of the frames have the same blaze function, we normalized the spectra of the other frames by this continuum spectrum.

Zoom In Zoom Out Reset image size

As described in Section 2.3, we obtained the normalized raw spectra S_raw for 41 frames. For further analysis, we converted S_raw to the transmission spectrum matrix in both the Earth and stellar rest frames with common wavelength grids (denoted by F_Earth and F_stellar, respectively).

Instrumental instability during the observation causes a time-varying shift of the spectrum on the detector. We used the telluric absorption lines in S_raw to measure the instrumental shift relative to the reference frame. The shift was estimated from the peak of the cross-correlation between S_raw and the model telluric transmission spectrum ¹¹ for each echelle order ¹² . The cross-correlation function (CCF) was computed from the velocity difference of −20 to 20 km s⁻¹ with a 0.1 km s⁻¹ step. The peak of the CCF was derived by fitting a Gaussian to the CCF. The instrumental shift derived by the above procedure is shown in Figure 2. The shift correlates with the dome temperature. We shifted S_raw with spline interpolation. We conducted a 5σ clipping to remove bad pixels and cosmic rays and obtained F_Earth.

Zoom In Zoom Out Reset image size

Then, we generated F_stellar by a cross-correlation analysis of F_Earth with the theoretical stellar spectrum. The detailed procedure is as follows. We measured the apparent radial velocity of the star in each frame by the CCF with the stellar model spectrum (Coelho et al. 2005). We fitted a Gaussian to the CCF and used the center position of the best-fit Gaussian as the stellar radial velocity in each frame. The apparent radial velocity variation is ∼600 m s⁻¹ between the first and last frames, which is consistent with the barycentric radial velocity variation of the Earth toward HD 149026. Then we shifted the spectra according to the stellar radial velocity of each frame. We note that this correction includes both the barycentric velocity variation and the stellar Rossiter–McLaughlin (RM) effect.

We also obtained the systemic radial velocity in this reduction. The systemic radial velocity used throughout this paper is −17.92 ± 0.04 km s⁻¹.

We performed a correction of the stellar line profile and detrending of the quasi-static signals with the SYSREM algorithm as follows. In the transmission spectrum matrix, the dominant signal is the stellar lines. The line profile for the stellar spectrum is time-varying. We calculated a CCF matrix C of the transmission spectrum matrix. Here C is a function of the transmission spectrum matrix F_xxx and the stellar model spectrum M and is defined as

$$\begin{eqnarray}&&{C}_{\mathrm{xxx}}^{\mathrm{st}}[i,\delta v]={F}_{\mathrm{xxx}}[i]\,\star \,M(\delta v),\end{eqnarray} \tag{ 1 }$$

where i denotes the frame number, ⋆ is the cross-correlation operator, and M(δ v) is the model spectrum shifted by δ v. We used a grid with steps of 1 km s⁻¹ as δ v. We denote F_stellar normalized by the mean of the out-transit as F_norm. This quantity ideally becomes the transmission spectrum matrix without the stellar signals. However, we found significant residuals originating from stellar signals in F_stellar. Figure 3 shows the CCF matrices, ${C}_{\mathrm{stellar}}^{\mathrm{st}}$ and ${C}_{\mathrm{norm}}^{\mathrm{st}}$. Residuals remain in ${C}_{\mathrm{norm}}^{\mathrm{st}}$ at the stellar radial velocity. The strong residuals in ${C}_{\mathrm{norm}}^{\mathrm{st}}$ were correlated with the mean count and the resolution of each frame (Figure 4). Therefore, we assume that these residuals originate from an imperfect correction of the nonlinearity of the CCD and instrumental profile (IP) variations. These residuals are problematic for detecting planetary signals because the stellar radial velocity is close to that of the planet in the planetary radial velocity–systemic velocity (V_sys) plane.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a first step to correct these residuals, we performed a deconvolution of the IP in S_raw. We estimated the IP of the observation night by a flat frame with an iodine cell, obtained just before the observation of the targets. We shifted each deconvolved spectrum and aligned them to the stellar rest frame, as shown in Section 2.4. Then, we generated a template stellar spectrum by averaging the shifted spectra. We fitted

$$\begin{eqnarray}&&a\,[T\,\ast \,\mathrm{IP}](\lambda )\end{eqnarray} \tag{ 2 }$$

to each raw value of F_stellar, where a is a scaling factor, T is the template stellar spectrum, and ∗ indicates the convolution operator. The scaling factor a (for the nonlinearity) and the IP are free parameters in this fitting. The IP was modeled by nine satellite Gaussian functions (Valenti et al. 1995). The positions and widths of the Gaussians were fixed. The only free parameters were their amplitudes. We normalized each spectrum by dividing each frame by the fitted template spectrum. The iodine cell generates numerous deep and sharp absorption lines in the 5000–5800 Å region in the 19 echelle orders (100–118) we used. We made this correction to these 19 orders only. For the other orders, we normalized them by dividing by the mean spectrum. We denote the transmission spectrum matrix generated by this procedure by F_IP. Figure 3 also shows ${C}_{\mathrm{IP}}^{\mathrm{st}}$, which is a CCF matrix generated from F_IP and the model stellar spectrum. As shown, the correction suppressed the noise at the stellar radial velocity.

Next, we used the SYSREM algorithm to remove the residuals from the stellar and telluric absorptions (Tamuz et al. 2005; Mazeh et al. 2007). SYSREM is an algorithm based on principal component subtraction and was developed to remove systematic variations in large numbers of transit light curves. Previous studies have applied this technique to remove telluric and stellar absorption lines in high-resolution spectra and succeeded in detecting signals of planetary atmosphere (e.g., Birkby et al. 2013, 2017; Nugroho et al. 2017; Hawker et al. 2018; Alonso-Floriano et al. 2019; Sánchez-López et al. 2019). In our analysis, the wavelength bins of the spectrum matrix were treated as numerous light curves. The systematic variation was removed by minimizing

$$\begin{eqnarray}&&S=\displaystyle \sum _{{ij}}\displaystyle \frac{{\left({r}_{{ij}}-{c}_{i}{a}_{j}\right)}^{2}}{{\sigma }_{{ij}}^{2}},\end{eqnarray} \tag{ 3 }$$

where r_ij is the average-subtracted stellar magnitude of the wavelength bin i, c_i is the effective extinction coefficient of the wavelength bin i, a_j is the effective air mass of the jth frame, and σ_ij is the uncertainty of the wavelength bin i of the jth frame. Each data point was weighted by its uncertainty, and the most significant systematic component was fitted and removed. Therefore, c_i and a_j are not necessarily the real extinction coefficient and air mass. The SYSREM algorithm can remove not only the telluric and stellar absorption lines but also other systematic effects (e.g., variations of air mass or water vapor column density, variations in instrumental conditions) by iterating the SYSREM reduction to the residual. We applied SYSREM to F_IP with 15 iterations.

Thus, we obtained 3D F_SYS[s, i, λ] for all of the echelle orders, with dimensions of 16 (denoted by s and ranging from zero to 15 SYSREM iterations) × 41 (frame number, denoted by i) × the number of wavelength bins.

The center position of the planetary absorption lines changes during observations according to the orbital motion, while the dominant noises from the stellar and telluric absorption lines are quasi-static. SYSREM removes this static noise, leaving the planetary signals. However, the planetary signals are also eliminated as the iterations proceed too much. When the S/N reaches the highest value, we regard the iteration number as the optimized value. Because the noise level may be different for each order, the optimal SYSREM iteration number might be different for each order. The injection of artificial planetary signals was often employed (e.g., Birkby et al. 2017) to find the optimal iteration number for each order, although this method may result in some bias (Cabot et al. 2019). Thus, recent studies use a common iteration number for all orders with the highest S/N (e.g., Nugroho et al. 2020b; Turner et al. 2020). Therefore, we determined a common iteration number for all orders with the highest S/N.

We generated the model transmission spectra for the cross-correlation template of Sc i, Ti i, V i, Cr i, Mn i, Fe i, and Co i for atomic species and TiO for molecular species. These atoms/molecules have many absorption lines in the wavelength range of our data. Each model spectrum included absorption by a single atom or molecule plus continuum absorption. We included a continuum absorption cross section of H⁻ and Rayleigh scattering by H and H₂. Here we describe the procedure to generate them.

We assumed a cloud-free, isothermal atmosphere at a temperature of 2000 K. We calculated the cross section using Python for Computational ATmospheric Spectroscopy (Py4CAtS; Schreier et al. 2019) for atomic species. We used partition functions taken from Barklem & Collet (2016) with a spline interpolation to obtain the partition function at the atmospheric temperature. We used the line list given by Kurucz (2018). For TiO, we used HELIOS-K (Grimm & Heng 2015) to compute the cross section. We used two different line lists for TiO, provided by Plez (1998; updated and corrected in 2012, hereafter TiO-Plez) and McKemmish et al. (2019; hereafter TiO-Toto). We considered the five most abundant isotopes, ⁴⁶Ti¹⁶O, ⁴⁷Ti¹⁶O, ⁴⁸Ti¹⁶O, ⁴⁹Ti¹⁶O, and ⁵⁰Ti¹⁶O, and combined them with the natural isotope ratio.

The abundance of each atom and molecule in each layer was calculated using FastChem (Stock et al. 2018), which calculates the abundances of atoms and molecules in the chemical equilibrium. We assumed a metallicity of 0.36 (Sato et al. 2005). We computed the line profile based on the Voigt profile with natural and thermal broadening because the transmission spectroscopy mainly probes a low-pressure region of the planetary atmosphere, and thus the contribution of pressure broadening is much smaller than that for other broadening mechanisms. We assumed a 1D plane-parallel atmosphere with 100 layers, which were evenly spaced in the log pressure scale from 10 to 10⁻¹⁵ bars.

The planetary radius at a wavelength of λ, R_p (λ), and the model transmission spectrum, $\mathrm{Tr}(\lambda )$, were calculated by

$$\begin{eqnarray}\begin{array}{rcl}{R}_{p}^{2}(\lambda ) & = & {R}_{{p}_{0}}^{2}+{{dR}}_{p}^{2}(\lambda )\\ & = & \ {R}_{{p}_{0}}^{2}+2{\displaystyle \int }_{{R}_{{p}_{0}}}^{{R}_{{p}_{0}}+{R}_{\max }}(1-\exp (-{\tilde{\tau }}_{\lambda ,\tilde{r}}))\tilde{r}d\tilde{r},\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{Tr}(\lambda )=1\ -\ {\left(\displaystyle \frac{{R}_{p}(\lambda )}{{R}_{s}}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

where ${R}_{{p}_{0}}$ is the planetary radius with white light, ${R}_{\max }$ is the maximum height of the model atmosphere, τ_λ,r is the optical depth at λ and the coordinates r, and R_s is the stellar radius. The symbols with tildes indicate physical quantities integrated over the transmission chord direction.

We generated model spectra with a radial velocity shift from −262.5 to 224.5 km s⁻¹ with a step size of 1 km s⁻¹, corresponding to a systemic velocity (V_sys) ranging from −120 to 80 km s⁻¹ and a semiamplitude of the planetary radial velocity (K_p ) from −100 to 300 km s⁻¹ (see Equation (7)).

We computed C^atom[s, i, δ v] as a function of F_SYS and a shifted model spectrum of each atom (denoted by M_atom). We subtracted the continuum of the model spectra to zero. We applied two high-pass filters to the spectrum of each frame to remove the effects of low-order variations before computing the cross-correlation. The cross-correlation was computed by

$$\begin{eqnarray}&&{C}^{\mathrm{atom}}[s,i,\delta v]=\displaystyle \sum _{k}\displaystyle \frac{{F}_{\mathrm{SYS}}[s,i,k]{M}_{\mathrm{atom}}(\delta v)[k]}{{\sigma }_{k}^{2}},\end{eqnarray} \tag{ 6 }$$

where σ_k is the error value at the kth wavelength bin. We assumed that σ_k is the square root of the sum of the variance of each wavelength bin and frame of F_norm. Then, we summed C^atom for the orders in which the planetary signal is strong. The criterion for this selection process is that the strongest line depth in the order is deeper than 2e⁻⁵ to the continuum, which corresponds to ${{dR}}_{p}\sim \tfrac{H}{2}$ in Equation (4), where H is the scale height of the planetary atmosphere. In this way, we obtained C^atom, whose dimensions are 41 (frame number) × 488 (radial velocity from −262.5 to 224.5 km s⁻¹) for each SYSREM iteration number. We computed the CCF map in the K_p –V_sys plane as follows. We integrated C^atom[s] (in the frame–radial velocity plane) along the path of expected planetary radial velocity RV_p ,

$$\begin{eqnarray}&&{\mathrm{RV}}_{p}={K}_{p}\ \sin (2\pi \phi (t))+{V}_{\mathrm{sys}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\phi (t)=\displaystyle \frac{t-{T}_{o}}{P},\end{eqnarray} \tag{ 8 }$$

where ϕ(t) is the orbital phase of the planet, and t is the Barycentric Julian Date (BJD) of each exposure. We only used the in-transit frames for this calculation. The orbital phase of each frame was calculated by the system parameters (Table 2) using the batman python software (Kreidberg 2015). We also calculated the expected amount of starlight blocked by the planet for each frame with batman and used it as a weighting coefficient with a K_p –V_sys integration. To make the S/N map, we divide the CCF map by the standard deviation of the map after excluding the expected position of the planetary signal, as well as the region with ∣K_p ∣ < 50 km s⁻¹, since the signal with low ∣K_p ∣ approaches zero as the SYSREM iteration proceeds.

Table 2. Basic Parameters of HD 149026 and HD 149026b

Parameter	Value
HD 149026
Radius (R⊙)	$1.41{\pm }_{0.03}^{0.03}$ a
Mass (M⊙)	$1.42{\pm }_{0.33}^{0.33}$ a
Teff (K)	$6179{\pm }_{15}^{15}$ a
log g	$4.37{\pm }_{0.04}^{0.04}$ a
Spectral type	G0 IV b
$v\sin i$ (km s−1)	$6.0{\pm }_{0.5}^{0.5}$ b
[Fe/H]	$0.36{\pm }_{0.05}^{0.05}$ b
HD 149026b
Radius (RJup)	$0.74{\pm }_{0.02}^{0.02}$ a
Mass (MJup)	$0.38{\pm }_{0.012}^{0.014}$ a
T0 (BJD)	$2,454,597.70713{\pm }_{0.00016}^{0.00016}$ c
Period (days)	$2.8758916{\pm }_{0.000002}^{0.000002}$ c
Inclination (deg)	$84.55{\pm }_{0.58}^{0.58}$ a
Kp (km s−1)	$169{\pm }_{20}^{20}$

Notes. K_p means the radial velocity semiamplitude.

^a Stassun et al. (2017). ^b Sato et al. (2005). ^c Bonomo et al. (2017).

Download table as:  ASCIITypeset image

Table 3. Positions, S/N, and Significance of Detected Peaks in the K_p –V_sys Plane

Species	Kp	ΔV	S/N	σ
Fe i	$\ =163.4{\pm }_{7.8}^{7.8}$	$\ =-0.7{\pm }_{0.6}^{0.6}$	2.8	2.8
Ti i	$\ =146.9{\pm }_{5.1}^{5.1}$	$\ =-3.2{\pm }_{0.4}^{0.4}$	5.0	4.4

Note. Here ΔV is the difference between the V_sys for the detected peak and the stellar radial velocity (−17.92 km s⁻¹). The detection significance was estimated by the phase-shuffling method (Esteves et al. 2017).

Download table as:  ASCIITypeset image

To estimate the errors in K_p and V_sys for the detected peaks, we calculated the likelihood maps (Brogi & Line 2019; Gibson et al. 2020; Nugroho et al. 2020b) in the K_p –V_sys plane. The log likelihood is defined as

$$\begin{eqnarray}&&\mathrm{log}L=-\displaystyle \frac{N}{2}\mathrm{ln}\left[\displaystyle \frac{1}{N}\displaystyle \sum _{i}\left(\displaystyle \frac{{f}_{i}^{2}}{{\sigma }_{i}^{2}}+\displaystyle \frac{{m}_{i}^{2}}{{\sigma }_{i}^{2}}-2\displaystyle \frac{{f}_{i}{m}_{i}}{{\sigma }_{i}^{2}}\right)\right],\end{eqnarray} \tag{ 9 }$$

where N is the number of data points. The third term can be represented by the CCF maps in the K_p –V_sys plane. The likelihood map is calculated with the exponential of $\mathrm{log}L$ after subtracting its maximum value. We fitted a Gaussian function to the slice of the likelihood map at the maximum position. Finally, the errors in K_p and V_sys were estimated from the standard deviation of the best-fit Gaussian function. We also calculated the detection significance with the phase-shuffling method (e.g., Esteves et al. 2017). We randomly shuffled the in-transit CCF map in the phase axis and integrated them with the K_p of the detected peak. We repeated this calculation 10,000 times. The noise at each V_sys was estimated by the standard deviation of the integrated CCF value.

We made the first detection of neutral titanium (Ti i) in a planetary atmosphere. We also detected a marginal signal of neutral iron (Fe i). The S/N for the detected signals was 5.0 for Ti i and 2.8 for Fe i. The detection significances calculated by the phase-shuffling method were 4.4 and 2.8. Table 3 surmmarizes the information of the detected signals.

Figure 5 shows the S/N map in the K_p –V_sys plane. The noise level was estimated from the standard deviation of the CCF values, excluding the region around the expected positions of the planetary signal.

Zoom In Zoom Out Reset image size

The peak values are ${K}_{p}=146.9{\pm }_{5.1}^{5.1}$, ${\rm{\Delta }}V=-3.2{\pm }_{0.4}^{0.4}$ km s⁻¹ for Ti i and ${K}_{p}=163.4{\pm }_{7.8}^{7.8}$, ${\rm{\Delta }}V=-0.7{\pm }_{0.6}^{0.6}$ km s⁻¹ for Fe i. The values of ΔV mean the difference between the V_sys of the detected peak and the stellar radial velocity (−17.92 km s⁻¹). We note that we only used echelle orders 100–118 for Fe i because the stellar noise was significant in orders without the line profile correction (see Section 2.5). We excluded frames 22 and 23 (the two closest frames to the mid-transit) for both atomic species in order to avoid the effects of the residual of the stellar signal because the planetary signal was expected to be overlapped with the center of the stellar absorption lines in these frames. Figure 6 shows the S/N for the detected peak with SYSREM iteration numbers from zero to 15. The SYSREM iteration number for the highest S/N is 9 for Fe i and 11 for Ti i. Figure 6 also shows the CCFs for the model transmission template spectrum and F_stellar. Figure 6 shows that the absorption lines for Fe i in the stellar spectrum were stronger than those for Ti i. The marginal detection of Fe i might be due to the strong contamination from the stellar spectrum.

Zoom In Zoom Out Reset image size

We could not find planetary signals for any other atomic or molecular species. For the nondetected species, we conducted an injection test to investigate the detection limit in our analysis (Figures 7 and 8). We injected a model template spectrum of each of the species into the spectrum matrix before SYSREM detrending. We convoluted a model spectrum with a rotation kernel to simulate line broadening by planetary tidally locked rotation by ∼1.3 km s⁻¹ and then injected the broadened model spectrum assuming K_p = −150 and V_sys = −20 km s⁻¹. To avoid the effect of unseen planetary signals from the planetary atmosphere, K_p was inverted (Merritt et al. 2020). The detection significance for the injected signal was calculated by the phase-shuffling method. In this paper, we consider species to be undetectable with our analysis if the significance of the injected signal is ≲3σ. For Sc i, Cr i, Mn i, and Co i, the expected signal was ≲3σ, comparable to the noise, and thus we could not confirm the existence of these atomic/molecular species. Our analysis was able to detect V i with ≳5σ, showing that V i is depleted in the HD 149026b atmosphere with ≳5σ, with the assumption of a cloud-free, 2000 K isothermal atmosphere and chemical equilibrium.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. Titanium Oxide

Assuming a solar-like abundance and chemical equilibrium, titanium mainly exists in the form of titanium oxide (TiO) at a temperature of ∼2000 K. Figure 9 shows the abundance of Ti i and TiO in chemical equilibrium calculated using Fastchem. We note that the TiO-Toto line list was confirmed to be accurate by cross-correlation with an M dwarf spectrum for λ > 4400 Å (McKemmish et al. 2019). To consider the implication of the nondetection of TiO in our analysis for both line lists, TiO-Plez and TiO-Toto, we performed injection tests for TiO.

Zoom In Zoom Out Reset image size

As shown in Figure 8, if the planet had a cloud-free, 2000 K isothermal atmosphere and chemical equilibrium, our analysis could detect TiO signals with >6σ for both line lists, assuming the line lists are completely accurate. This is indicative because we detected neutral Ti but not TiO. We discuss the implications in the next section in more detail.

Both Fe and Ti can condense and take the form of clouds below the condensation temperature. The condensation temperatures for typical condensates at a pressure of 1 mbar are 1650 K for FeO and 1582 K for CaTiO₃ (Wakeford & Sing 2015). Our detection of Ti i (and the hint of Fe i) in a gaseous phase indicates that the temperature of HD 149026b is higher than the condensation temperatures for Ti- and Fe-bearing condensates. We note that the equilibrium temperature for HD 149026b is ∼1680 K, consistent with our results, if we assume the received energy is evenly redistributed to the whole planet and the Bond albedo = 0.1. This temperature is higher than the condensation temperature for most chemical species that can form clouds (Wakeford & Sing 2015). Given this high temperature and the fact that we detected atomic species in the atmosphere, the effect of clouds would not be significant for HD 149026b. However, we cannot rule out the cloud hypothesis because several other chemical species have condensation temperatures higher than ≥2000 K (e.g., see Figure 2 in Mbarek & Kempton 2016).

A possible explanation for the nondetection of TiO is a supersolar C/O ratio in the planetary atmosphere, as reported for several other hot Jupiters (Brewer et al. 2017), as the abundance of TiO strongly depends on the C/O ratio. Other scenarios that have been proposed to explain the absence of TiO in the planetary atmospheres are less likely in this case, as described below. Our detection shows that titanium is not depleted in the atmosphere of HD 149026b; thus, titanium depletion by the cold-trap effect (Spiegel et al. 2009; Parmentier et al. 2013) would not be significant in the atmosphere of HD 149026b. Strong chromospheric activity of the host star might destroy the TiO in planetary atmospheres, though the activity of HD 149026 is low based on the analysis of Ca ii H & K lines (Knutson et al. 2010). Figure 9 also shows the volume mixing ratio of Ti i and TiO with C/O values of 0.55 and 1.0. As the C/O ratio increases, the abundance of TiO significantly decreases. Figure 8 shows the result of an injection test assuming a C/O ratio of 0.55 and 1.0. In the case of a C/O of 1.0, our analysis will not detect a TiO signal with ≲3σ; therefore, the scenario of C/O ≳ 1 is consistent with our results. We also estimated the expected significance of the Ti i signal in the same way (Figure 10). The predicted significance of Ti i with a C/O value of 0.55 and 1.0 was ∼3.3 and 4.5, respectively. The significance of the detected signal was 4.4; thus, our result is consistent with the scenario that the C/O ratio is higher than the solar value. However, we note that the predicted significance strongly depends on atmospheric properties such as thermal profile or scattering properties such as clouds, because the depth of the absorption lines in the model transmission spectrum is affected by these properties.

Zoom In Zoom Out Reset image size

The C/O ratio for the planetary atmosphere provides information on its formation history. Some hot Jupiters have been reported to have supersolar C/O ratios (e.g., HD 209458b; Madhusudhan & Seager 2009). A supersolar C/O ratio suggests that the planet swept up a large portion of its atmosphere from the gas in the protoplanetary disk outside the water snow line (gaseous oxygen is reduced); that is, the planet formed outside the water snow line and migrated to the current orbit (e.g., Öberg et al. 2011). Therefore, our detection of Ti i and nondetection of TiO suggest that HD 149026b formed outside the disk and migrated to the current orbit. Object HD 149026b is considered to have a massive core (∼67 M_⊕) like WASP-59b (Hébrard et al. 2013). The formation history of hot Jupiters with massive cores is still in debate. It is thought that such hot Jupiters are formed in an environment with a high dust grain density (Ikoma et al. 2006). Kanagawa et al. (2018) showed that a planetary-induced gap in the protoplanetary disk can make a "dust ring" at the outer edge of the disk, and a gas giant with a massive core can be formed in the dust ring. In this scenario, HD 149026b would have migrated by dynamical scattering. Our detection of neutral metals and the implication of a high C/O ratio provide clues on the formation history of hot Jupiters with massive cores. Line et al. (2014) reported that the atmospheric C/O ratio for HD 149026b is 0.45–1.0 (68% confidence), determined by Spitzer photometry of a secondary eclipse. The C/O ratio of ∼1 hypothesis does not contradict their results. However, the detection limit for atom/molecule abundance discussed in this paper depends on the properties of the planetary atmosphere, such as scattering and vertical mixing. In our model, we assumed Rayleigh scattering by H and H₂, continuum absorption by H⁻, and a cloud-free, isothermal atmosphere under chemical equilibrium, and the detection limits we set should be interpreted in this context only.

The detected signals were slightly blueshifted compared with the central star. A typical wind speed from the dayside to nightside in a tidally locked hot Jupiter's atmosphere is a few kilometers per second (e.g., Kataria et al. 2016), which is consistent with the blueshift of the detected signals. However, if we assume the planet is tidally locked, the rotational velocity at the equator is ∼1.3 km s⁻¹. The asymmetric distribution of atomic species in the planetary atmosphere also shifts the planetary signals (Ehrenreich et al. 2020). Since the blueshift of the detected signals is also consistent with the planetary rotation, we could not confirm the cause of the blueshift of the detected signals.

In this paper, we presented the results of high-resolution transmission spectroscopy of HD 149026b at visible wavelengths using Subaru/HDS. We detected the signal of neutral titanium (4.4σ) and marginal signal of neutral iron (2.8σ) in the atmosphere of HD 149026b. Our study is the first detection of neutral titanium in the planetary atmosphere.

The detection of Ti i and nondetection of TiO suggests a supersolar C/O ratio in the atmosphere of HD 149026b.

The formation scenario for hot Jupiters is of great interest in exoplanetary science. Object HD 149026b is a hot Jupiter with a massive core, and its formation history is still controversial. Our detection of neutral metals and the implication of a supersolar C/O ratio add important clues for understanding the unique properties of this hot Jupiter.