Polarization of Rotationally Oblate Self-luminous Exoplanets with Anisotropic Atmospheres

The polarimetric technique has been gaining momentum over the past few years in the field of exoplanet characterization, especially for young directly imaged planets. The use of adaptive optics coronagraphic systems has enabled us to conduct direct photometric, spectroscopic, and polarimetric observations of substellar companions (e.g., Bryan et al. 2018; Miles-Páez et al. 2019; Jensen-Clem et al. 2020). Polarimetric observations of brown dwarfs and the directly imaged planets using highly sensitive instruments have already been reported by Millar-Blanchaer et al. (2020), Jensen-Clem et al. (2020), and van Holstein et al. (2021), among others. The increasing number of reports of polarization observations of substellar mass objects calls for a better understanding of the atmospheric processes that give rise to polarization. The thermal radiation of the objects becomes linearly polarized due to scattering by atmospheric molecules and cloud particles. The detected net nonzero disk-integrated polarization, on the other hand, is attributed to the asymmetry of the visible disk, for various reasons such as oblateness, inhomogeneous or patchy cloud coverage in the atmosphere, and gravitational darkening, etc. (e.g., Sengupta & Kwok 2005; Sengupta & Marley 2009, 2010; Marley & Sengupta 2011; de Kok et al. 2011; Stolker et al. 2017; Sanghavi & Shporer 2018). Polarimetric observations of the red dwarf stars and brown dwarfs (Miles-Páez et al. 2015; Millar-Blanchaer et al. 2020) strongly suggest that the observed polarimetric variations correlate with the rotation-induced oblateness of those objects. Hence, in this paper, we focus on the estimation of the polarization caused solely by the rotation-induced oblateness of young self-luminous gas giant planets.

Theoretical computation of the flux and polarization observable from a directly imaged exoplanet plays a pivotal role in guiding future polarimetric missions. These forward models can be used to relate the observed flux and polarization from the planets to their different physical and atmospheric properties. The high-resolution spectroscopic studies (e.g., Snellen et al. 2014; Bryan et al. 2018; Xuan et al. 2020) provide information about the line-of-sight (LOS) component of the equatorial rotation velocity (${v}_{e}\sin i$) of the brown dwarfs and the directly imaged giant exoplanets. To date, however, the values of the inclination angles (i, not to be confused with the orbital inclination angle) of the rotation axis with respect to the observer cannot be found from such studies. Polarimetric observations can complement such spectroscopic as well photometric observations, and break the degeneracies among the estimated properties, such as the rotation speed (v_e ), the inclination angle of the rotation axis (i), and the surface gravity (g), among others.

The net nonzero disk-averaged polarization of a substellar object can be attributed to its rotation-induced oblateness even for a cloud-free atmosphere, as in the case of T dwarfs (e.g., Sengupta & Marley 2009). However, polarization from a cloud-free atmosphere is predominantly caused by the Rayleigh scattering of the thermal emission by the atoms and molecules of the gases present in the atmosphere. Hence, such a polarization signal is detectable only in the visible wavelength region. Given the presence of condensate cloud or haze particulates in the atmospheres, a significant amount of polarization arises in the infrared wavelength region, where the brightness of the objects peaks (e.g., Sengupta & Marley 2010; de Kok et al. 2011; Marley & Sengupta 2011; Sanghavi & Shporer 2018). Photometric and spectroscopic observations of brown dwarfs and the directly imaged gas giants (e.g., Burgasser et al. 2002; Marois et al. 2008; Zhou et al. 2016) indicate the presence of clouds in their atmospheres. In the present investigation, therefore, we consider the cloudy atmospheres of the young gas giants in order to estimate the observable polarization.

Sengupta & Marley (2009, 2010) and Marley & Sengupta (2011) presented theoretical models based on the spherical harmonic expansion (SHE) technique in order to estimate the detectable polarization of a fast-rotating substellar object with a certain oblateness. In order to estimate the disk-integrated polarization, these authors followed the technique prescribed by Simmons (1982). They have presented polarization models for different properties of the substellar objects using self-consistent radiative–convective equilibrium models for the atmospheres. The model presented by Jensen-Clem et al. (2020) calculates polarization from the rotation-induced oblateness using the same technique, but also includes evolution models to self-consistently calculate the moment of inertia, size, and surface gravity of a rotating substellar object at a given mass and age. de Kok et al. (2011) presented numerical models for the polarization observable from the directly imaged planets due to various sources of asymmetry, such as oblateness, banded clouds, and hot spots, etc. They followed a technique similar to that of Sengupta & Marley (2009, 2010) and Marley & Sengupta (2011), of correlating the rotation-induced oblateness with the observable polarization. Again, Sanghavi & Shporer (2018) and Sanghavi & West (2019) included factors such as the gravitational darkening and flattening of the planetary disks in order to calculate the polarization from the fast-rotating brown dwarfs. On the other hand, Stolker et al. (2017) presented a Monte Carlo-based technique of calculating the polarization from a self-luminous exoplanet arising from its oblateness and due to the presence of banded and patchy clouds. However, none of these works consider the effect of the anisotropy in the atmosphere of the ellipsoidal planet across the disk, except for the recent work by Mukherjee et al. (2021) who adopted a Monte Carlo-based method to model the polarization observable from the binary brown dwarfs Luhman 16 A and B by considering inhomogeneous atmospheres.

In this paper, we present a technique for calculating the polarization detectable from a substellar object with an oblate spheroid shape by incorporating the effects of both the vertical (layerwise) inhomogeneity and horizontal (across the disk) anisotropy of the atmospheres of the planets. The atmosphere of an oblate planet is usually inherently anisotropic, even if the gas molecules or the cloud particles across the disk are distributed homogeneously. This anisotropy in the atmosphere arises due to the fact that the depth of the atmospheric shell itself varies across the disk, depending on its latitude and the inclination angle i (Simmons 1982). We calculate the disk-resolved and disk-integrated emergent radiation field by using the same numerical recipes described in Chakrabarty & Sengupta (2021). However, the atmospheric models and the radiative transfer equations that we solve for the present cases are different, since we focus on the thermal emission from the young wide-orbit giant exoplanets, and not the reflected flux. In order to calculate the layerwise and wavelength-dependent properties of the atmospheres of such exoplanets, we use the state-of-the-art Sonora models (Marley et al. 2018,2021a, 2021b), which are available as grids of effective temperature (T_eff), surface gravity (g), and metallicity. These calculations are based on the radiative–convective equilibrium and evolutionary models (Saumon & Marley 2008; Marley & Sengupta 2011; Jensen-Clem et al. 2020) of the substellar objects. We have also incorporated a generalized Henyey–Greenstein–Rayleigh (HGR) phase matrix for a better representation of the scattering due to the cloud particulates present in the atmospheres. This new technique will allow us to model the flux and polarization observable from substellar objects with inhomogeneous or patchy atmospheres in our follow-up investigation.

We have chosen the self-luminous directly imaged exoplanets β Pic b and ROXs 42B b, whose projected rotation velocities ${v}_{e}\sin i$ are reported to be 25 ± 3 km s⁻¹ (Snellen et al. 2014) and ${9.5}_{-2.3}^{+2.1}$ km s⁻¹ (Bryan et al. 2018), respectively. Also, Jensen-Clem et al. (2020) reported the J-band polarimetric observations of these planets. They could not measure the polarization of these planets precisely, but could only set upper limits on the observed polarization that can be attributed to these planetary sources. However, a number of future observations and dedicated missions are being planned for the precise polarimetric characterization of the exoplanets. Hence, in the present study, we have developed models for the polarization detectable from these planets in four different wavelength bands of the SPHERE-IRDIS instrument, to present the maximum amount of polarization that can be expected from these planets in those wavelength bands, mainly due to their oblate structures.

Section 2 shows the relation between the rotation rate and the induced oblateness of a fast-rotating planet, as well as the effect of this oblateness on the shape of the atmosphere. Section 3, which is further divided into four subsections, provides a detailed description of our overall approaches, such as the atmospheric model adopted, the numerical technique used to calculate the local Stokes vectors, the detailed calculation of the HGR phase matrix in regard to the scattering by clouds, and the numerical recipe developed to calculate the disk-integrated flux and polarization from an oblate spheroid planet. We elaborate the specific band-averaged models that we have developed for the planets β Pic b and ROXs 42B b, for different values of ${v}_{e}\sin i$ and i, in Section 4. The results are discussed in detail in Section 5, and the investigation is concluded in Section 6.

The oblateness of a fast-rotating young giant planet can be expressed as f = 1 − R_p /R_e , where R_p and R_e are the polar and the equatorial radii of the planet, respectively. The oblateness (f) depends on the rotation rate (Ω), the surface gravity (g), and the mass of the planet (M_P), as dictated by the Darwin–Radau relationship (Barnes & Fortney 2003), given by

$$\begin{eqnarray}&&f=\displaystyle \frac{{{\rm{\Omega }}}^{2}{R}_{e}}{g}{\left[\displaystyle \frac{5}{2}{\left(1-\displaystyle \frac{3K}{2}\right)}^{2}+\displaystyle \frac{2}{5}\right]}^{-1},\end{eqnarray} \tag{ 1 }$$

where K = I/(M ${}_{{\rm{P}}}{R}_{e}^{2}$), I being the moment of inertia of the planet. As I does not depend on the oblateness to the first order (Barnes & Fortney 2003), I can be calculated by assuming the planet to be spherical. The moments of inertia of the self-luminous and cloud-free substellar objects can be obtained from the Sonora Bobcat tables (Marley et al. 2021b) available online. ⁴ However, for most of our calculations for the cloudy self-luminous giant planets, we have assumed the interior of the planet to be a stable polytropic gas of index n = 1, by following Marley & Sengupta (2011; hereafter, MS11), and hence set K = 0.261 in Equation (1) (Chandrasekhar 1933). From the oblateness–rotation speed relationship (see Figure 15) of Jensen-Clem et al. (2020), the value of K for the planets β Pic b and ROXs 42B b can be found to be 0.276 and 0.297, respectively. Clearly, our assumption of K = 0.261 provides a valid representation of the moments of inertia of the self-luminous giant planets. However, while calculating the models for those planets (see Section 4), we have assumed the specific values of K mentioned above.

The atmosphere of a fast-rotating planet can be considered to be an oblate spheroidal shell (see Figure 2(b) of Simmons 1982). The outermost radius of the planet at a colatitude Θ_col across the disk can be expressed as R(Θ_col) = r(Θ_col)R_e . The factor r(Θ_col) can be defined as (Simmons 1982; Sengupta & Marley 2009)

$$\begin{eqnarray}&&r({{\rm{\Theta }}}_{\mathrm{col}})=\displaystyle \frac{1}{{\left[1+({A}^{2}-1){\cos }^{2}{{\rm{\Theta }}}_{\mathrm{col}}\right]}^{1/2}},\end{eqnarray} \tag{ 2 }$$

where A = 1/(1 − f) = R_e /R_p . Clearly, the shell has a varying thickness across the disk, and the thickness at a colatitude Θ_col is equal to r(Θ_col) times the thickness along the equatorial plane.

Observations (e.g., Snellen et al. 2014) suggest that the young gas giants and the brown dwarfs do not have any systematic differences in their rotation rates. However, the relatively low surface gravity of the young gas giants allows them to attain higher rotation-induced oblateness, as high as ∼0.2–0.3, compared to the brown dwarfs. In all of our calculations, we consider the oblateness of the planets to be ≲0.44, as beyond this limit (this upper limit is for n = 1; for higher values of n, the upper limit is even lower) the atmospheres of the gas giant planets or the brown dwarfs are likely to become unstable (e.g., James 1964; Marley & Sengupta 2011; Sanghavi & Shporer 2018).

3.1. Model Atmosphere of Self-luminous Exoplanets

In the present work, we use the Sonora grid of pressure–temperature (P–T) profiles (e.g., Ackerman & Marley 2001; Marley et al. 2002; Saumon & Marley 2008; Sengupta & Marley 2009, 2010; Marley & Sengupta 2011; Jensen-Clem et al. 2020; Marley et al. 2021a), which are calculated for the atmospheres of the substellar mass objects by using the radiative–convective equilibrium, the chemical equilibrium, and evolutionary models self-consistently. We also use the precalculated Sonora grids for the layerwise properties of the atmospheres, such as mass extinction coefficients, mass absorption coefficients, and the scattering asymmetry parameters. These models are based on the opacity calculations of Freedman et al. (2008, 2014) and the molecular abundance calculations of Lodders (2020, 2010). In all of our derivations, we assume solar metallicity and the solar system value for the carbon-to-oxygen ratio. Also, we assume the pressure and temperature to be uniform along each of the stratified oblate spheroid layers, calculated over a range of effective temperature (T_eff) and surface gravity (g) that describe the atmospheres of the young gas giants. Figure 2 shows the P–T profiles of the young gas giants for different values of T_eff and g.

Zoom In Zoom Out Reset image size

To start with, we consider the cloud-free atmosphere of a Jupiter-sized T brown dwarf or methane dwarf with T_eff = 1000 K and g = 1000 ms⁻² to benchmark our calculations based on Rayleigh scattering. We use the corresponding P–T profile (see Figure 2), the mass extinction coefficients, and the mass absorption coefficients from the Sonora Bobcat model for cloud-free atmospheres (Sengupta & Marley 2009; Marley et al. 2018, 2021a, 2021b). We consider an equatorial speed (v_e ) of 90 kms⁻¹ that corresponds to a rotation period (P_rot) of 1.386 hr and i = 90°. We then compare our results with the results presented by Sengupta & Marley (2009).

Next, we use the cloudy Sonora models to calculate the polarization detectable from the young self-luminous giant planets (Marley & Sengupta 2011; Jensen-Clem et al. 2020). We consider young planets with T_eff ranging between 1000 K and 2068 K and g between 10 ms⁻² and 140 ms⁻². For the cloud models, we assume a sedimentation efficiency, f_sed = 2 (Ackerman & Marley 2001). In order to solve the vector radiative transfer equations, we apply the δ-Eddington approximation (Joseph et al. 1976; Batalha et al. 2019) on the optical depth, single-scattering albedo, and asymmetry parameters, before calculating the HGR phase matrix for the cloud particles.

3.2. Calculations of the Stokes Vectors at Local Points of the Atmosphere

Similar to the case of brown dwarfs (Sengupta & Marley 2009, 2010), MS11 use the SHE technique to calculate the local intensity vectors (or Stokes vectors) across the disk of the self-luminous gas giant exoplanets. In order to account for the effect of the oblateness on the polarization of these planets, they follow the formalism presented by Simmons (1982), which has also been explained in detail by Sengupta & Marley (2009). While the formalism prescribed by Simmons (1982) is a first-order approximation and uses less computational effort, a better and more accurate approach is needed in order to calculate the disk-integrated polarization arising purely due to the oblateness of the planets. We have developed a new technique to incorporate the effect of the anisotropy of the medium over the disk by using a three-dimensional approach.

As explained in Chakrabarty & Sengupta (2021), we divide the observable planetary disk into a grid of M and Φ, where M denotes the cosine of the polar angular position and Φ denotes the azimuthal angular position on the disk defined with respect to the disk center M = 1. We consider Φ = 0 along the Y-axis. The rotation axis of the planet is assumed to lie along the YZ-plane (see Figure 1), depending on the inclination angle of the rotation axis, i, with respect to the observer. When i = 90°, the poles lie on the Y-axis, and for i = 0°, the poles lie on the Z-axis, i.e., on the LOS of the observer. Clearly, when i = 0° (polar view), the disk exhibits circular symmetry (uniform along Φ), and hence the net disk-integrated polarization becomes zero in this case. On the other hand, the effect of the oblateness is maximum for i = 90°, i.e., for the equatorial view. We transform the disk-centric coordinates (M, Φ) to the corresponding colatitudes (Θ_col) of the planet to calculate the effect of the oblateness of the planet, as shown in Equations (1)–(2). For an inclination angle, i, we calculate the colatitudes by

$$\begin{eqnarray}&&\cos {{\rm{\Theta }}}_{\mathrm{col}}=M\cos i+\sqrt{1-{M}^{2}}\sin i\cos {\rm{\Phi }}.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Unlike the case presented by Chakrabarty & Sengupta (2021), we here assume the planetary atmospheric shell to be stratified and homogeneous along the concentric oblate spheroid layers. We divide the direction of the propagation of light at each location on the disk into a grid of μ and ϕ, where μ denotes the cosine of the polar angular direction and ϕ denotes the azimuthal angular direction defined with respect to the local meridian. While solving the radiative transfer equations for the case of thermal emission, we consider the radiation field to be circularly symmetric about the local normal at each location on the disk. However, in the present study, we assume circular symmetry about the radial line joining the center of the planet and the point on the disk, in order to avoid further complications, leveraging the fact that ω < 0.44. As a result, we can assume the incoming and outgoing radiation to be uniform along the ϕ-direction. Hence, we adopt the axisymmetric form of the vector radiative transfer equation, including multiple scattering of the internal radiations at each location (M, Φ) and at each layer, with an optical depth, τ(M, Φ), defined from the top of atmosphere (TOA). This is given by

$$\begin{eqnarray}&&\begin{array}{l}\mu \displaystyle \frac{d{{\boldsymbol{I}}}_{{\boldsymbol{mer}}}(\tau ,M,{\rm{\Phi }},\mu )}{d\tau }={{\boldsymbol{I}}}_{{\boldsymbol{mer}}}(\tau ,M,{\rm{\Phi }},\mu )-\displaystyle \frac{\omega (\tau )}{2}\\ \times \ {\int }_{-1}^{1}{{\boldsymbol{P}}}_{{\boldsymbol{mer}}}(\mu ;\mu ^{\prime} ){{\boldsymbol{I}}}_{{\boldsymbol{mer}}}(\tau ,M,{\rm{\Phi }},\mu ^{\prime} )\,d\mu ^{\prime} -(1-\omega ){\boldsymbol{B}}(\tau ).\end{array}\end{eqnarray} \tag{ 4 }$$

The wavelength-dependent intensity vector (or Stokes vector) I _mer = [I_mer, Q_mer, U_mer, V_mer] (see, e.g., Chandrasekhar 1960; Chakrabarty & Sengupta 2021) is a function of the optical depth (τ), the disk location (M, Φ), and the direction (μ). The subscript "mer" signifies that this vector is defined with respect to the local meridian plane at (M, Φ). B (τ) denotes the unpolarized internal radiation source function from a layer with optical depth, τ, from the TOA. Therefore, B (τ) = [B(τ), 0, 0, 0]. ω denotes the single-scattering albedo of the atmosphere at τ.

For the linear polarization due to scattering, V_mer = 0. ${P}_{\mathrm{mer}}(\mu ;\mu ^{\prime} )$ denotes the ϕ-averaged phase matrix from the direction $\mu ^{\prime}$ to the direction μ. We consider the Rayleigh phase matrix (see Equation (5) of Chakrabarty & Sengupta 2021) for the cloud-free atmosphere, and for a cloudy atmosphere we add the HGR phase matrix (Liu & Weng 2006) to the Rayleigh phase matrix in the ratio of their single-scattering albedos, ω_cld and ω_mol, respectively (MS11; Chakrabarty & Sengupta 2021). In our axisymmetric calculations, both the Rayleigh and HGR phase matrices do not contribute to the U-component, and hence we can set U_mer = 0. Further details of the calculations of the HGR phase matrix are described in Section 3.3.

We solve Equation (4) using Discrete Space Theory (e.g., Sengupta & Marley 2010, 2009; Marley & Sengupta 2011; Sengupta et al. 2020; Chakrabarty & Sengupta 2021, 2020), and, finally, calculate the intensity vector at the TOA toward the observer by

$$\begin{eqnarray}&&{{\boldsymbol{I}}}_{{\boldsymbol{mer}},{\boldsymbol{obs}}}({\boldsymbol{M}},{\boldsymbol{\Phi }})={{\boldsymbol{I}}}_{{\boldsymbol{mer}}}(\tau =0,M,{\rm{\Phi }},\mu =M).\end{eqnarray} \tag{ 5 }$$

For a validity check, we first adopt the SHE formalism used by Sengupta & Marley (2009), and compare the results with those calculated using the new technique developed and described in this paper. For this purpose, we considered the case for a cloud-free atmosphere of a T dwarf with T_eff = 1000 K and g = 1000 ms⁻², and the case for a cloudy atmosphere of a young giant planet with T_eff = 1000 K and g = 30 ms⁻². We first consider the object to be perfectly spherical while calculating the intensity vectors ( I _{mer , obs} (M, Φ)). The flux and polarization calculated at all the local points of the planetary disk are then integrated over an oblate spheroid by using the SHE formalism. This is done by including the factor r(Θ_col), as explained in Sengupta & Marley (2009).

Next, in order to account for the effect of the anisotropy of the atmosphere due to the rotation-induced nonsphericity of the object, we introduce the optical depth (d τ) at each stratified layer as a function of the colatitude (Θ_col), and hence a function of M and Φ (see Equation (3)), which is given by

$$\begin{eqnarray}&&d\tau ({{\rm{\Theta }}}_{\mathrm{col}}(M,{\rm{\Phi }}))=r({{\rm{\Theta }}}_{\mathrm{col}}(M,{\rm{\Phi }}))\,d{\tau }_{e},\end{eqnarray} \tag{ 6 }$$

where d τ_e = d τ(Θ_col = 0), which is the optical depth along the equatorial plane. By substituting Equation (6) into Equation (4), we calculate I _{mer , obs} ( M , Φ) at each location on the disk. Thus, in our present formalism, the optical depth is considered to be anisotropic.

Since I _{mer , obs} ( M , Φ) is defined with respect to the local meridian at each location on the disk, we need to transform I _{mer , obs} ( M , Φ) to a common plane of reference all over the disk before integrating them. We choose the YZ-plane, as this is a common plane of reference. Accordingly, we apply the rotation matrix, L (Φ) (see Equation (2) of Chakrabarty & Sengupta 2021), to calculate the intensity vector I _obs ( M , Φ) defined with respect to the YZ-plane as

$$\begin{eqnarray}&&{{\boldsymbol{I}}}_{{\boldsymbol{obs}}}({\boldsymbol{M}},{\boldsymbol{\Phi }})={\boldsymbol{L}}({\rm{\Phi }})\,{{\boldsymbol{I}}}_{{\boldsymbol{mer}},{\boldsymbol{obs}}}({\boldsymbol{M}},{\boldsymbol{\Phi }}).\end{eqnarray} \tag{ 7 }$$

From Equation (7), we calculate the components of I _obs ( M , Φ) as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\mathrm{obs}}(M,{\rm{\Phi }}) & = & {I}_{\mathrm{mer},\mathrm{obs}}(M,{\rm{\Phi }})\\ {Q}_{\mathrm{obs}}(M,{\rm{\Phi }}) & = & {Q}_{\mathrm{mer},\mathrm{obs}}(M,{\rm{\Phi }})\cos 2{\rm{\Phi }}\\ {U}_{\mathrm{obs}}(M,{\rm{\Phi }}) & = & -{Q}_{\mathrm{mer},\mathrm{obs}}(M,{\rm{\Phi }})\sin 2{\rm{\Phi }}.\end{array}\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

3.3. The Generalized Henyey–Greenstein–Rayleigh Phase Matrix

The phase matrix P_mer, which determines the angular distribution of photons before and after scattering, strongly dictates the amount of polarization. In the case of cloudy self-luminous exoplanets, MS11 (see also Sengupta & Marley 2010 and Jensen-Clem et al. 2020) have used an HGR phase matrix (Liu & Weng 2006), where they have assumed a Rayleigh-dominated symmetric form for the phase function (See Figure 4). This is valid for the axisymmetric radiation field. However, this approach tends to overestimate the polarization for higher values of the asymmetry parameter g_as, because the amount of light scattered in the forward and backward directions differs greatly at higher values of g_as, as shown in Figure 4. For this reason, we have used the azimuth-dependent expression of the HGR phase matrix as given in Equation (9) of Liu & Weng (2006). Therefore, in order to derive the azimuth-dependent phase matrix ${{\boldsymbol{P}}}_{{\boldsymbol{HGR}}}(\mu ,\phi ;\mu ^{\prime} ,\phi ^{\prime} ;{g}_{\mathrm{as}})$, we replace the term $\cos {\rm{\Theta }}$ by

$$\begin{eqnarray}&&\cos {\rm{\Theta }}=\mu \mu ^{\prime} +\sqrt{(1-{\mu }^{2})(1-{\mu }^{{\prime} 2})}\cos \,(\phi ^{\prime} -\phi ),\end{eqnarray} \tag{ 9 }$$

where $\phi ^{\prime}$ and ϕ are the azimuthal angles that determine the directions of the incoming (before scattering) and the outgoing (after scattering) photons, respectively. However, we have adjusted for the change in the reference planes of the intensity vectors during the scattering process from the direction ($\mu ^{\prime} ,\phi ^{\prime}$) to (μ, ϕ), following Chakrabarty & Sengupta (2021):

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{P}}}_{{\boldsymbol{mer}},{\boldsymbol{HGR}}}(\mu ,\phi ;\mu ^{\prime} ,\phi ^{\prime} ;{g}_{\mathrm{as}})\\ \ \ =\ {\boldsymbol{L}}(-(\pi -{i}_{2})){{\boldsymbol{P}}}_{{\boldsymbol{HGR}}}(\mu ,\phi ;\mu ^{\prime} ,\phi ^{\prime} ;{g}_{\mathrm{as}}){\boldsymbol{L}}({i}_{1}),\end{array}\end{eqnarray} \tag{ 10 }$$

where i₁ denotes the angle between the local meridian plane along $(\mu ^{\prime} ,\phi ^{\prime} )$ and the plane of scattering and i₂ denotes the angle between the plane of scattering and the local meridian plane along (μ, ϕ). P _{mer , HGR} denotes the HGR phase matrix defined with respect to the local meridian at (M, Φ). The axisymmetric form of the phase matrix, i.e., ${{\boldsymbol{P}}}_{{\boldsymbol{mer}},{\boldsymbol{HGR}}}(\mu ;\mu ^{\prime} ;{g}_{\mathrm{as}})$, is calculated by averaging ${{\boldsymbol{P}}}_{{\boldsymbol{mer}},{\boldsymbol{HGR}}}(\mu ,\phi ;\mu ^{\prime} ,\phi ^{\prime} ;{g}_{\mathrm{as}})$ over $\phi ^{\prime}$ numerically.

Zoom In Zoom Out Reset image size

Figure 5 shows the spectra, at a local point, of the intensity I_obs and the degree of polarization Q_obs/I_obs in the direction toward the observer, which are calculated by assuming the planet to be a perfectly spherical body. It demonstrates the sole effect of the modified HGR phase matrix on the local intensity vector. The intensity and polarization are shown at two different polar angular points, viz. M = 0.24 and M = 0.9, by adopting the HGR phase matrix of MS11 and the generalized phase matrix.

Zoom In Zoom Out Reset image size

3.4. Integration of the Local Stokes Vectors over the Oblate Planetary Disk

By integrating the intensity vector I _obs over the solid angles of the disk projected to the observer, we get the disk-integrated flux vector F = [F, F_Q, F_U, F_V], defined with respect to the YZ-plane. Here, F denotes the total disk-integrated flux, while F_Q and F_U denote the polarized disk-integrated flux (F_V] = 0 for linear polarization). F_U is also zero because of axisymmetry (see Equation (8)) with respect to the rotation axis. The disk-integrated degree of polarization is then given by

$$\begin{eqnarray}&&P=\displaystyle \frac{{F}_{{\rm{Q}}}}{F}.\end{eqnarray} \tag{ 11 }$$

When we use the SHE method, we employ Equations (46) and (47) of Sengupta & Marley (2009) to calculate F and F_Q. On the other hand, when we apply the present technique, we calculate F and F_Q by using the following expression:

$$\begin{eqnarray}&&{\boldsymbol{F}}=\displaystyle \frac{{R}_{e}^{2}}{{D}^{2}}{\int }_{0}^{2\pi }{\int }_{-1}^{1}{{\boldsymbol{I}}}_{{\boldsymbol{obs}}}(M,{\rm{\Phi }})\,{MdMd}{\rm{\Phi }},\end{eqnarray} \tag{ 12 }$$

where D is the distance between the planet and the observer. We consider $\tfrac{{R}_{e}^{2}}{{D}^{2}}=1$ in Equation (12), so that the emergent flux and its state of polarization are estimated at the TOA of the planet.

Figure 3 shows the disk-integrated flux F and polarization P estimated by using both the SHE method (Sengupta & Marley 2009) and the technique described here (i.e., Equation (12)). On the other hand, Figure 6 shows a comparison between F and P derived by using the formalism of MS11 and by using the present technique with and without the modified HGR phase matrix (see Section 3.3).

Zoom In Zoom Out Reset image size

All further derivations are done by solving Equation (12) and employing the generalized HGR phase matrix. The corresponding color maps of the intensity and the polarization [I_obs(M, Φ) and Q_obs(M, Φ)] at all points across the disk of a self-luminous giant planet with T_eff = 1000 K and g = 30 ms⁻² for two different rotation periods P_rot = 5 hr and P_rot = 10 hr are shown in Figures 7 and 8, respectively. These figures help us to visualize the asymmetry in the values of Q_obs(M, Φ) over the Φ-direction that gives rise to the net nonzero disk-integrated polarization, the amount of which depends on i and on the oblateness, and hence on the rotation rate. Similarly, Figure 9 presents the detectable F and P of a self-luminous Jupiter-sized gas giant for i = 90° and for different values of T_eff and g. The corresponding atmospheric models are explained in Section 3.1. Figures 10–11 show P at different wavelength bands as a function of P_rot for planets with different values of T_eff and g and different values of i, respectively. In order to calculate the flux and polarization at different wavelength bands, we have chosen the response functions of the broadband filters Y, J, H, and K of the SPHERE-IRDIS instrument placed at the backend of the Very Large Telescope.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

High-resolution spectroscopic observations have been used to estimate the LOS components of the equatorial rotation velocities (${v}_{e}\sin i$) of young directly imaged planets such as β Pic b (Snellen et al. 2014) and ROXs 42B b (Bryan et al. 2018), among others. However, these studies cannot provide us with any information about the inclination angle of the rotation axis with respect to the observer. Jensen-Clem et al. (2020) have reported the polarimetric observations of these two planets in the J-band, along with those of five brown dwarfs. The polarimetric observations of β Pic b and ROXs 42B b have only been able to set upper limits on the degree of linear polarization.

In Figures 12 and 13, we present updated models for the polarization of the planets β Pic b and ROXs 42B b, respectively, as functions of ${v}_{e}\sin i$ for different values of i. The oblatenesses for different equatorial rotation speeds have been calculated following Jensen-Clem et al. (2020; see Figure 15). The amounts of polarization have been constrained by their observed values of ${v}_{e}\sin i$. In order to estimate the degree of polarization at different wavelength bands, we consider the response functions of the broadband filters corresponding to the Y-, J-, H-, and K-bands ⁵ of SPHERE-IRDIS. We discuss the results in the next section.

We have presented new models for the scattering polarization detectable from a self-luminous directly imaged exoplanet arising due to the asymmetry caused by its rotation-induced oblateness. We have introduced a new formalism to account for the effect of the anisotropic distribution of the scatterer in the atmosphere of an oblate planet. It is shown that the emergent flux of an oblate substellar mass object does not differ much from that of a spherical object. In fact, the flux can be calculated with an insignificant error by assuming the object to be spherical. On the other hand, polarization, being a measure of the anisotropy in the radiation field, strongly depends on the distribution of the scatterers and on the shape of the visible disk of the object. This is clearly demonstrated in Figure 3, which shows that the predicted polarization of a T dwarf can increase significantly when the visible disk is nonspherical, while the flux does not alter much. The T dwarf has been chosen to focus our calculations for a cloud-free atmosphere in order to highlight only the effect of our new three-dimensional approach. The polarization of a cloud-free object arises by atomic and molecular scattering, and the angular distribution of the scattered photons is described by the Rayleigh phase matrix.

In the case of cloudy self-luminous giant planets, we have adopted the atmospheric models of MS11. However, we have introduced a generalized HGR phase matrix in order to treat the scattering due to the cloud particulates more accurately. Figure 4 demonstrates that the generalized scattering phase matrix described in this paper provides a better representation of the asymmetric nature of the angular distribution of the photons before and after the scattering of condensate cloud particles. Figure 5 shows that although the total intensity does not alter significantly, the polarized intensity is slightly overestimated if the Rayleigh-dominated symmetric approximation of the HGR phase matrix is used. This in turn alters the disk-integrated polarization, as shown in Figure 6. Figure 6 also shows the flux and polarization estimated by using the generalized HGR phase matrix, as well as by using the new algorithm of integrating the polarization over the rotation-induced oblate disk.

In Figures 7 and 8, we present color maps of the intensity vectors normalized to that at the disk center for a self-luminous Jupiter-sized planet with T_eff = 1000 K and g = 30 ms⁻² rotating with periods P_rot = 5 hr and P_rot = 10 hr, respectively, at different values of i. These figures help us to visualize the patterns of the total intensity (I_obs) and the polarized intensity (Q_obs) across the planetary disk, and to understand how they contribute to the disk-integrated flux and polarization. The variation of the intensities in the M-direction exhibits limb darkening and an increase in polarized intensity from the disk center to the limb, which agrees with Figure 5.

Zoom In Zoom Out Reset image size

The variation in the polarization at the local points across the disk along the Φ-direction determines the amount of the net nonzero disk-integrated polarization. At an inclination of i = 0°, both Figures 7 and 8 show that I_obs is uniform along the Φ-direction and that Q_obs exhibits an exact cosinusoidal pattern (see Equation (8)) along the Φ-direction. The amount of positive polarization along the Y-axis is exactly nullified by the amount of negative polarization along the X-axis, causing the net disk-averaged polarization to be zero. This occurs because of the perfect symmetry about the disk center when the inclination angle is i = 0°. However, for a given value of P_rot, the asymmetry in the Φ-direction increases with an increasing value of i as the poles shift from M = 0. This effect is more pronounced as the rotation speed increases (a lower value of P_rot), and hence the oblateness increases. Figure 7, for P_rot = 5 hr, represents such a case of highly oblate object. In this case, the asymmetry is maximum at i = 90°, and hence the cancellation is minimum because the positive polarization (along the Y-axis) is more than the negative polarization (along the X-axis). This causes the disk-integrated polarization to be maximum at i = 90° for a value of the rotation period. This effect, however, reduces for the case of P_rot = 10 hr, since the rotation-induced oblateness is significantly small (f ∼ 0.05) at this low rotation rate. Figure 8 does not show significant variation in the patterns of I_obs and Q_obs over the disk for different values of i. Although in this case, the oblateness is close to that of Jupiter (∼0.064), detecting the disk-averaged degree of polarization ( ∼ 0.005% for i = 30° and ∼ 0.02% for i = 90°; see Figure 11) of such slowly rotating directly imaged exoplanets may be quite challenging due to the extremely low signal-to-noise ratio.

Zoom In Zoom Out Reset image size

Figure 9 shows the spectra of the flux and polarization observable from young self-luminous giant planets with different effective temperature and surface gravity. It shows that the polarization significantly reduces with the increase in surface gravity for the same values of T_eff and P_rot, as pointed out by MS11. This is also demonstrated in Figure 10, where the detectable disk-averaged polarizations at different wavelength bands for planets with different effective temperature T_eff and surface gravity g are presented as functions of the spin rotation period P_rot. On the other hand, Figure 11 presents disk-integrated polarization at different wavelength bands that may be detectable from a self-luminous Jupiter-sized planet with T_eff = 1000 K and g = 30 ms⁻². In this figure, the amounts of polarization of the planet with different rotation periods and different inclination angles (i) of the rotation axis with respect to the observer are presented. Clearly, with an increase in i, the polarization increases for a given value of P_rot. The polarization predicted is maximum for i = 90°, i.e., for an equatorial view when the disk asymmetry is maximum, and close to zero for i = 0°, i.e., for a polar view (see Figure 1) when disk asymmetry is almost absent.

Figures 12 and 13 demonstrate that the polarization due to the oblateness of the planets β Pic b and ROXs 42B b can be constrained with the observational reports of their ${v}_{e}\sin i$, viz. 25 ± 3 km s⁻¹ (Snellen et al. 2014) and ${9.5}_{-2.3}^{+2.1}$ km s⁻¹ (Bryan et al. 2018), respectively. For a given value of ${v}_{e}\sin i$, a lower value of i implies a higher value of the equatorial rotation speed (v_e ), which causes the planet to be more oblate, resulting in a higher value of the detectable degree of polarization. Conversely, for a given value of v_e , a decreasing value of i implies declining asymmetry of the visible planetary disk, leading to a lower degree of polarization, as evident from Figure 11. However, the former effect outweighs the latter, causing an overall drop in the detectable polarization with an increase in i. The minimum value of i is determined by the stability limit of the spin rotation speed of a planet for a given observed value of ${v}_{e}\sin i$. These lower limits of i for the two planets are found to be ∼20° and ∼14°, respectively. Thus, we have been able to determine both the upper and lower limits of the detectable disk-integrated polarization of these planets by using the reported values of ${v}_{e}\sin i$ from observations.

Zoom In Zoom Out Reset image size

Clearly, the degrees of polarization of the two planets predicted by our present model are much below the 1σ nondetection upper limits, viz. 0.18% and 0.19%, reported by Jensen-Clem et al. (2020). These predictions are subject to the particular atmospheric models and cloud models that we adopt for a particular planet. With a change in the atmospheric model, e.g., if we assume chemical inequilibrium (Madhusudhan et al. 2020), or if we choose a different value of the sedimentation parameter f_sed for the cloud model, the predicted degree of polarization can alter. So, if future observations confirm the amount of polarization to be in the order of the nondetection upper limits set by Jensen-Clem et al. (2020), we need to adopt a different atmospheric or cloud model or to identify other sources of polarization, such as cloud band or surface inhomogeneity. Our new technique is capable of performing such analysis by including the inhomogeneity of atmospheres, which we will explore in our future work. However, the present study only focuses on the new three-dimensional technique of accounting for the anisotropy in the atmosphere of the planet introduced by the departure from sphericity due to spin rotation and on the generalized form of the HGR phase matrix on the detectable polarization.

It is foreseeable that imaging polarimetry and, in the future, spectropolarimetry in synergy with photometric and spectroscopic techniques will be able to open the door to the unexplored regimes of exoplanets. At present, polarimetry is the only prospective tool that can convey information about the axial tilt of a directly imaged planet. This technique is also useful for excavating information about the deeper layers of atmospheres, especially about the condensate clouds that cannot be fully probed with the spectroscopic technique. However, the correct interpretation of the polarimetric observations demands self-consistent and extensive work on the development of a forward model that can describe the atmospheric processes correctly and explain the polarization arising from those processes.

All of our calculations are centered on the atmospheric models and cloud models (fixing f_sed at 2) of Marley & Sengupta (2011) and Jensen-Clem et al. (2020). With changes in the models adopted, the predicted polarizations would also change. Consequently, polarimetric observations can be utilized to distinguish these models. In the present study, we have explored only a few cases in order to demonstrate our new approach of accounting for the effect of the oblateness of fast-rotating directly imaged planets on their observable polarization. This study also presents our generic methodology for calculating the phase matrices for a more accurate representation of the scattering process that predominantly dictates the polarization observable from a cloudy exoplanet. Our work provides a three-dimensional view of the atmosphere of an oblate planet, which allows us to calculate the polarization arising from the asymmetry caused by rotation-induced nonsphericity more accurately. This technique will be further applied in our follow-up work: calculating the polarization owing to the inhomogeneous atmospheres of substellar objects, for example, due to the presence of banded or patchy clouds.

To bring this paper into proper shape, we have intensively used the high-performance computing facility (Delphinus) of the Indian Institute of Astrophysics, Bangalore. We are thankful to the computer division of the Indian Institute of Astrophysics for the help and cooperation extended for the present project. We are thankful to the reviewer for the critical reading of the manuscript and for providing useful suggestions.