© ESO 2021

1 Introduction

Since the discovery of the first extra-solar planet by Mayor & Queloz (1995), the number of detected exoplanets increased enormously to over 4000 (The Extrasolar Planets Encyclopaedia¹; Schneider et al. 2011; from November 2020). The current generation of imaging instruments, for example the Gemini Planet Imager (GPI; Macintosh et al. 2008) or the SPectro-Polarimetric High-Contrast Exoplanet Research (SPHERE; Beuzit et al. 2019), potentially provide the opportunity to detect the infrared polarized signal from sufficiently bright (e.g., self-luminous) exoplanets. Furthermore, with modern polarimeters, such as the HIgh-Precision Polarimetric Instrument (HIPPI; Bailey et al. 2015) or the POlarimeter at Lick for Inclination Studies of Hot jupiters 2 (POLISH2; Wiktorowicz & Nofi 2015), it is now possible to measure the polarized flux at a parts-per-million level. In addition, new high-accuracy polarimeters are going into operation soon, such as the High-Precision Polarimetric Instrument-2 (HIPPI-2; Bailey et al. 2020), or are in preparation, such as POLLUX (Bouret et al. 2018).

Recently, Bott et al. (2018) reported observations of linear polarized radiation of the WASP-18 system, which harbors a massive planet (approximately 10 M_J) orbiting close to its star with an orbital period of p < 1 day. Although the measured polarization is dominated by the interstellar medium, the authors were able to set an upper limit of 40 ppm (99% confidence level) on the amplitude of a reflected polarized radiation planetary signal. Thus, they could rule out certain atmospheric models, such as optically thick atmospheres dominated by Rayleigh scattering clouds.

Since the emitted radiation of solar-type stars can be assumed to be unpolarized (Kemp et al. 1987), the reflected radiation of a planet is solely polarized due to scattering processes within its atmosphere. Different types and properties of atmospheric particles cause different characteristics in the polarization. Analyzing the reflected polarized radiation can therefore be a useful tool not only for detecting extra-solar planets, but also for characterizing their atmospheres. Hansen & Hovenier (1974) showed the benefits of this method by determining the cloud properties of the atmosphere of our neighboring planet Venus. With the previously mentioned instruments, detailed studies and characterization of exoplanetary atmospheres have come into reach.

Meanwhile, various theoretical studies of linear polarization resulting from Rayleigh scattering (Buenzli & Schmid 2009) and scattering by water clouds (Karalidi et al. 2012), circular polarization of cloudy exoplanets (Rossi & Stam 2018), the determination of the cloud coverage of exoplanets (Rossi et al. 2018), and the linear polarization of scattered radiation from self-luminous exoplanets (Stolker et al. 2017) have been performed. These studies have in common that they either focus on specific polarization mechanisms or assume certain simplifications. Examples of common simplifications are the assumption of a locally plane-parallel planetary atmosphere, the often missing cross-talk between latitudinal and longitudinal regions (i.e., locally horizontal radiation transport) in the scattering process (not in the case of Monte Carlo simulations), black planetary surfaces (i.e., absorbing surfaces), the often neglected state of circular polarization, and the assumption of either a point-like star or plane-parallel illumination on the planetary atmosphere.

Simulating the radiative transfer in a planetary atmosphere is a very complex problem since various parameters have an impact on the observed radiation. Furthermore, the planetary surface and its influence become relevant if the atmosphere is optically thin. The reflected (polarized) radiation caused by the surface depends on the material, such as land masses with or without biomass, or oceans. However, the reflected flux of an optically thick atmosphere is dominated by the wavelength-dependent scattering and absorption properties of the atmospheric particles.

Our goal is to develop a polarization radiative transfer tool that considers all the relevant contributions to the continuum polarization signal due to various polarization mechanisms while avoiding simplification inherent to previous approaches. Thus, we provide a tool for the comprehensive analysis of the polarized flux reflected by extra-solar planets. This development is based on the publicly available 3D Monte Carlo radiative transfer code POLARIS² (Reissl et al. 2016). The code is already optimized to handle the full spectrum of state-of-the-art polarization mechanisms (e.g., scattering and/or thermal re-emission, line and/or continuum polarization). It is well tested and has been applied to a broad range of astrophysical models: molecular clouds (Reissl et al. 2017; Pellegrini et al. 2020; Seifried et al. 2020), Bok globules (Brauer et al. 2016), and protoplanetary disks (Brauer et al. 2019; Heese et al. 2020; Brunngräber & Wolf 2020).

In Sect. 2 we define our radiative transfer model. The Stokes parameters that describe the intensity and state of polarization of the radiation are introduced in Sect. 2.1. In Sect. 2.2 the emission of photon packages by a spatially extended radiation source is briefly described, including our approach to reduce the run-time of the code by restricting the emissionof photon packages. In Sect. 2.3 the planetary model, its atmosphere, and its surface are briefly outlined. In Sect. 3 we verify the new routines and present the results for various atmospheric models including diffuse surface reflection. The results are compared to analytical solutions in the case of a simple diffuse reflecting sphere and already existing computations for various structured atmospheres. As an example case study (Sect. 4), we investigate the presence of an additional circumplanetary ring and its impact on the reflected (polarized) flux. Finally, our results are summarized in Sect. 5.

2 Radiative transfer model

We briefly introduce the Stokes parameters and the description of scattering that are used to model the radiation field and the radiativetransfer. In order to adapt the existing radiative transfer code POLARIS to the specific requirements for handling radiation scattering in planetary atmospheres, both the illuminating source and the atmosphere have to be considered in detail afterwards.

2.1 Stokes formalism

The radiation field is represented by photon packages that are defined by their wavelength-dependent Stokes parameters. The Stokes parameters are used to determine the intensity, state, and degree of polarization (see, e.g., Bohren & Huffman 1983); they are combined in a 4D vector $S={(I,Q,U,V)}^{\text{T}}$, where I is the total intensity, Q and U are the linear polarization, and V the circular polarization. Furthermore, the quantities $$P_{\text{l}}=\frac{\sqrt{Q^2+U^2}}{I},\tan (2\gamma )=\frac{U}{Q},P_{\text{c}}=\frac{V}{I},$$(1)

describe the degree of linear polarization P_l, the corresponding angle of linear polarization γ, and the degree of circular polarization P_c, respectively. If the photon package scatters, the change in polarization is obtained by multiplying the incoming Stokes vector S_in with a scattering matrix F(Θ, Φ), where Θ ∈ [0, π] and Φ ∈ [0, 2π] are the scattering angles. The resulting Stokes vector after scattering is then given by $$S_{\text{out}}\propto F(\Theta ,\Phi )\cdot L(\Phi )\cdot S_{\text{in}}.$$(2)

Here the rotation matrix L(Φ) rotates the Stokes vector into different frames (e.g., into the scattering plane or into the observers frame) if the photon package is detected. We assume that the atmospheric particles are spherical and distinguish between Rayleigh scattering (particle size r smaller compared to the wavelength λ) and the more general case of Mie scattering (arbitrary particle size). The general scattering matrix (or Müller matrix) for spherical particles has the following simplified structure (e.g., Bohren & Huffman 1983): $$F(\Theta )=\left(\begin{array}{cccc}{F}_{11}& F_{12}& 0& 0\\ F_{12}& F_{22}& 0& 0\\ 0& 0& F_{33}& F_{34}\\ 0& 0& -F_{34}& F_{44}\end{array}\right).$$(3)

In addition to their wavelength-dependent Stokes vector, the photon packages are characterized by their point of emission or last point of interaction as well as the direction of their propagation through the model space. Every photon package has its own reference frame that is transformed at every scattering event in order to describe the position and direction of propagation of the photon package in the 3D model space. A detailed description of the mathematically and geometrically random walk of the photon packages can be found in Fischer (1993), while an illustration of the rotation of the frame of the photon package due to scattering can be found in Fig. B.1. POLARIS tracks the position and path of every photon package through the model space until it is absorbed, or leaves the model space and is detected by the observer.

2.2 Illuminating source

The illuminating source is assumed to be a spherical, spatially extended radiation source with radius R_⋆. To simulate theemission of radiation, four independent angles are required to specify its random starting point on the stellar surface (photosphere) and direction of emission. The location of emission on the photosphere is defined by the angles θ₁ and ϕ₁. The polar angle θ₁ is defined with respect to the z-axis of the global model space and the azimuthal ϕ₁ is in the (x, y)-plane. The model space is illustrated in Fig. 1 . In this setup, the angle θ₁ has values in the range [0, π] and ϕ₁ in the range [0, 2π]. The angles θ₂ and ϕ₂ describe the direction of emission at that starting point. Here the polar angle θ₂ is the angle between the surface normal and the direction of propagation of the photon package, and ϕ₂ is the azimuthal angle in the surface plane. Thus, θ₂ has values in the range [0, π∕2] and ϕ₂ in the range [0, 2π], respectively. A further detailed description can be found in Cashwell & Everett (1959) or Niccolini et al. (2003), among others.

As the solid angle under which a planet is seen from the central star is very small (even in the case of hot Jupiters), we restrict the emission angles θ₁, θ₂, and ϕ₂ accordingly. Consequently, we sample only emission locations and directions to ensure that the photon package propagates towards the planet. To compensate the oversampling in that range, we have to weight the net energy of the photon package. Weighting the energy of the photon packages is a popular method for optimizing Monte Carlo simulations (see, e.g., Cashwell & Everett 1959; Yusef-Zadeh et al. 1984; Lucy 1999; Juvela 2005; Baes et al. 2016). These restrictions are necessary to avoid inefficient computations because only a small number of photon packages actually hit the planetary atmosphere if the four angles are calculated over their full range. For example, for an Earth-sized planet with radius R_p = 6.3781 × 10⁶ m at a distance of d_⋆ = 1 au to its central star, the solid angle under which the planet is seen from the star is approximately 5.7 × 10⁻⁹ sr. For a hot Jupiter with a radius of R_p = 7 × 10⁷ m and a distance of d_⋆ = 0.1 au, the solid angle is approximately 6.9 × 10⁻⁵ sr.

The allowed range for the emission angles has to be determined individually for every photon package. While the range of the angle θ₁ depends on the planetary radius R_p, stellar radius R_⋆, and the distance between the planet and the star d_⋆, the range of the angle θ₂ in addition depends on the position on the stellar photosphere that is defined by θ₁. The range of the angle ϕ₂ furthermore depends on the direction that is given by the previously sampled angle θ₂. A description of the calculation of the allowed range for the emission angles is given in Appendix A.

While a constant brightness distribution of the stellar photosphere as seen from the planet is certainly a valid approximation in most cases, the impact of limb darkening might become relevant for close-in planets. If the stellar brightness distribution is known, it can be considered directly during the sampling of the emission angles. For the case of restricted emission, we can simply adjust the weight of the photon package according to the stellar intensity profile.

2.3 Planetary model

After the emission the photon package travels towards the planet and experiences a random optical depth τ (e.g., in the atmosphere) before it interacts. The scattering and absorbing properties of our model atmosphere are defined by the optical properties of the gas particles (molecules, atoms) and cloud, aerosols, or dust particles. As shown in Fig. 1, the planet with radius R_p is in the center of the 3D model space.

2.3.1 Atmosphere

The atmospheric structure is described by a spherical grid with radial atmospheric boundaries and with polar and azimuthal boundaries. Therefore, both vertical and horizontal inhomogeneities of the atmosphere can be considered. Inside a spherical grid cell, the number density of particles is constant. A simple assumption for the radial atmospheric pressure profile or density profile, and thus the optical depth of the gas phase, is based on the condition of hydrostatic equilibrium and the equation of state for the ideal gas. Starting at the top of the atmosphere we can describe the increase in pressure towards the center as

Fig. 1 Illustration of the model space. The planet with radius Rp is located at the center of the model space, i.e., the origin of the coordinate system. The dashed lines represent the radial cell boundaries that divide the atmosphere above the surface into individual layers. The star, located at distance d* along the positive z-axis, has a radius of R*. The first two emission angles are sampled and define the starting position of a photon package on the stellar surface. The position is described by the polar angle θ1 with respect to the z-axis and the azimuthal ϕ1 in the (x, y)-plane. The second polar angle θ2 is defined with respect to the surface normal and the azimuthal angle ϕ2 is defined inthe surface plane at this location.

$$\frac{\text{d}p}{\text{d}h}=\rho g=\frac{p{M}_{\text{g}}g}{\mathcal{R}T}=\frac{p{m}_{\text{g}}g}{k_{\text{B}}T},$$(4)

where g is the gravitational acceleration, $\mathcal{R}$ the gas constant, k_B the Boltzmann constant, and T the temperature. The quantity M_g is the molar mass of the gas in kg mol⁻¹ and m_g is the mass of one gas particle (or molecule) in kg. By integrating this equation we can calculate the layer boundaries using $$h_{i+1}=h_i-\mathcal{H}(p_{i+1})\ln \left(\frac{p_{i+1}}{p_i}\right),\mathcal{H}=\frac{\mathcal{R}T(p)}{M_{\text{g}}g},$$(5)

where $\mathcal{H}$ is the scale height with a pressure-dependent temperature profile. Such a temperature profile could be given by the approximation of a gray atmosphere (e.g., Hansen 2008; Guillot 2010). The pressure bins are logarithmically distributed between p_min at the upper boundary (at h₀ = 0 m) of our model and p_max at the lower boundary at R_p. If we express the gas density following Eq. (4) as ρ = m_gn with a number density n of the gas particles, the optical depth τ_g of a given layer can be expressed as $$\frac{{\tau }_{\text{g}}}{C_{\text{ext}}}=n(h_{i+1}-h_i)=\frac{p_{i+1}-p_i}{m_{\text{g}}g},$$(6)

where C_ext is the extinction cross section of the gas particles (i.e., the scattering cross section plus the absorbing cross section).

For the gas particles (i.e., particles that are very small compared to the wavelength) the scattering properties can be described by the Rayleigh scattering theory including depolarization due to anisotropic molecules. The elements of the scattering matrix for Rayleigh scattering are given by Hansen & Travis (1974). The wavelength-dependent scattering cross section can be expressed by the Rayleigh cross section (Sneep & Ubachs 2005), $$C_{\text{sca}}(\lambda )=\frac{24{\pi }^3}{{\lambda }^4{n}_{\text{s}}^2}\frac{{(n^{\prime }{(\lambda )}^2-1)}^2}{{(n^{\prime }{(\lambda )}^2+2)}^2}\left(\frac{6+3{\rho }_{\text{d}}}{6-7{\rho }_{\text{d}}}\right),$$(7)

where n′(λ) is the wavelength-dependent real part of the refractive index and n_s is the number density at standard conditions; the density n_s is usually referred to as the Loschmidt constant. The real part of the refractive index can be approximated by $$n^{\prime }(\lambda )=A\left(1+\frac{B}{{(\lambda \times {10}^6)}^2}\right)+1,$$(8)

where A and B are real constants and depend on the type of the considered gas (Cox 2000). The wavelength-, pressure-, and temperature-dependent absorption cross sectionC_abs(λ, p, T) can be calculated by knowing the strength and width of the absorption line of the corresponding gas (see, e.g., Rothman et al. 1998; Tennyson et al. 2016).

In addition to the gas particles, we also consider clouds, aerosols, and dust particles in the atmosphere. For these larger particles the elements of the scattering matrix F as well as the particle cross sections C_ext, C_abs, and C_sca are calculated with MIEX (Wolf & Voshchinnikov 2004), which is based on the Mie scattering theory and takes into consideration arbitrarily large size parameters 2πr∕λ. The cloud particles have a size distribution described by Hansen (1971): $$N(r)\propto r^{(1-3v_{\text{eff}})/v_{\text{eff}}}{\text{e}}^{-r/(r_{\text{eff}}{v}_{\text{eff}})}.$$(9)

Here r is the radius of the cloud particle, r_eff the effective particle radius, and v_eff the effective variance (see Hansen 1971). Similar to the gas phase, the cloud particles, aerosols, or dust particles are constant in a grid cell defined by radial, azimuthal, and polar boundaries.

To account for the absorption of radiation, the energy of the interacting photon package is weighted by the single scattering albedo of the particle $$\omega =\frac{C_{\text{sca}}}{C_{\text{ext}}}=\frac{C_{\text{sca}}}{C_{\text{sca}}+C_{\text{abs}}}.$$(10)

2.3.2 Surface

If a photon package reaches the planetary surface, it is reflected back into the current layer or transmitted into the next layer, depending on the properties of the surface. Similar to the case of scattering by a particle, the Stokes vector of the photon package is transformed using a surface reflection (or transmission) matrix to describe the change in state of polarization.

For a diffuse (Lambertian) reflecting surface the reflection matrix element R₁₁ is equal to the planetary surface albedo, while the remaining elements are zero. Therefore, the reflected radiation does not depend on the incident angle and is fully depolarized. For a specular reflection, however, the angle of the reflected or transmitted photon package depends on the incoming angle and on the refractive indices of the optical media that are separated by the surface. In this case the reflection and transmission matrix is based on the Fresnel equations (see, e.g., Zhai et al. 2010; Garcia 2012).

3 Selected test cases

We start with various selected test cases to verify the correct implementation of the concepts outlined in Sect. 2. The test cases include a homogeneous Rayleigh scattering atmosphere, a diffuse reflecting sphere, and an inhomogeneous atmosphere with both gas and cloud particles. Our computations are compared to the results by Buenzli & Schmid (2009), the analytical solution by Russell (1916), and the results by Karalidi et al. (2012) for the Rayleigh scattering atmosphere, the diffuse reflecting sphere, and the inhomogeneous atmosphere, respectively. To make our results independent of the general planetary and stellar parameters, we normalize the detected flux by $$\pi B(\lambda ,T)\frac{R_{\star }^2{R}_{\text{p}}^2}{d_{\star }^2{d}_{\text{obs}}^2}.$$(11)

An overview of the general model parameters is given in Table 1.

3.1 Homogeneous cloud-free atmosphere

We considera simple homogeneous atmosphere consisting of gas particles with an optical depth τ_g above a diffuse reflecting surface with albedo ω_s, and compare our computations with the results by Buenzli & Schmid (2009). The grid is divided into one radial, one polar, and one azimuthal cell since the atmosphere is both vertically and horizontally homogeneous. As this atmospheric layer contains only gaseous particles, the scattering properties are described by the Rayleigh scattering theory. The depolarization factor is set to zero and we ignore the absorption by gas particles (C_abs = 0). Thus, the total radial optical depth of the gaseous layer with height H_atm can be calculated by $${\tau }_{\text{g}}=C_{\text{sca}}{n}_{\text{g}}{H}_{\text{atm}}.$$(12)

Here, C_sca is the scattering cross section of the particles, and n_g is the number density that can be scaled to simply adjust the optical thickness of the atmosphere. The radius of our planetary model is set to R_p = 7 × 10⁷ m with a height of the atmospheric layer of H_atm = 10⁵ m. The distance to the star is set to d_⋆ = 0.1 au. This model is referred to as model I.

Figure 2 shows the normalized reflected flux and degree of linear polarization as a function of the phase angle α for an optical depths τ_g = 5 with a surface albedo ω_s = 0 and ω_s = 1. We find a very good agreement for the reflected flux (relative difference |δI| < 0.008 at α = 7.5°) and the degree of linear polarization (relative difference |δP_l| < 0.002 at α = 92.5°) between our numerical calculations and the results by Buenzli & Schmid (2009).

An optical depth of τ_g = 0 corresponds to the case of a simple Lambertian reflecting sphere without an atmosphere above. Here, the computations are compared to the analytical solution derived by Russell (1916), $$I(\alpha )=\frac{2{\omega }_{\text{s}}}{3\pi }(\sin (\alpha )+(\pi -\alpha )\cos (\alpha )),$$(13)

where I(α) is the normalized reflected flux of the planet at a phase angle α. We find very good agreement between our numerical calculations and the analytical solution for the considered value of the surface albedo (relative difference |δI| < 0.004 at α = 0°).

Fig. 2 Normalized reflected flux (top) and degree of linear polarization (bottom) as a function of the phase angle α for various atmospheric models τg = 5, 0, 0.097; ωs = 1, 0.5, 0; fc = 0, 0.3, 0.7). Shown are the cloud-free homogeneous atmosphere (model I; orange and light blue lines), the Lambertian reflecting sphere without an atmosphere τg = 0; green line), and the cloudy inhomogeneous atmosphere (model II; yellow and dark blue lines) with various cloud fractions fc. Our numerical calculations (plus symbols) are compared with results by Buenzli & Schmid (2009) for the homogeneous atmosphere, with the analytical solution by Russell (1916) for the Lambertian reflecting sphere, and with the results by Karalidi et al. (2012) for the inhomogeneous atmosphere. The color-coding applies to both panels. See Sects. 3.1 and 3.2 for details.

3.2 Inhomogeneous cloudy atmosphere

In this section we investigate the impact of an inhomogeneous cloud cover in the atmosphere and compare our computations with the results by Karalidi et al. (2012). The atmosphere of the planetary model describes an Earth-like planet. Here the optical depth of the gas measured in radial direction is determined by given pressure and temperatureprofiles of a mid-latitude atmosphere tabulated by McClatchey (1972). The distance to the star is set to d_⋆ = 1 au and the planet has a radius of R_p = 6.3781 × 10⁶ m with a height of the atmosphere of 10⁵ m. The optical depth of each layer is calculated applying Eq. (6), the given pressure profile, and by using M_g = 2.8964 × 10⁻² kg mol⁻¹, which is typical for air (Cox 2000). The resulting total radial optical depth of the gas particles is thus τ_g ≈ 0.097. We use a depolarization factor of 0.028 (typical for air; Bates 1984) and ignore the absorption by the gas particles. The cloud particles are located between 3 and 4 km above the surface and have a total radial optical depth of τ_c = 2. They have an effective radius of r_eff = 2 μm and an effective variance of v_eff = 0.1 (Karalidi et al. 2012) with a refractive index of 1.335 + 10⁻⁵i (Karalidi et al. 2011) typical for water droplets. In contrast to the gas particles, the cloud layer has a horizontally inhomogeneous structure (i.e., a random patchy pattern) that is described by the cloud fraction f_c. The atmosphere is divided into 32 radial, 90 polar, and 180 azimuthal cells to define the 3D atmospheric model. The surface has an albedo of ω_s = 0, thus it absorbs all incident radiation. This model is referred to as model II.

The derived elements F₁₁ and F₁₂ of the Müller matrix of air and for the water cloud particles as a function of the scattering angle Θ are shown in Fig. 3. Since we use a size distribution for the cloud particles, the plotted values are averaged over all grain sizes ranging from r_min to r_max and weighted by the size distribution N(r) in order to represent the effective values: $$F_{ij}=\frac{{\displaystyle {\int }_{r_{\text{min}}}^{r_{\text{max}}}{F}_{ij}}(r)N(r)\text{d}r}{{\displaystyle {\int }_{r_{\text{min}}}^{r_{\text{max}}}N}(r)\text{d}r}.$$(14)

For incoming unpolarized radiation, the scattering matrix element F₁₁ indicates the fraction of radiation scattered towards the direction described by the phase angle α = π − Θ. The ratio F₁₂∕F₁₁ indicates the linear polarization degree after single scattering for initially unpolarized light. In the case of F₁₂ < 0 the polarization vectors are perpendicular to the scattering plane, while they are parallel for F₁₂ > 0. While a Rayleigh scattering profile is obtained in the case of molecular hydrogen, the scattering distribution of the cloud particles strongly increases towards smaller scattering angles, indicating strong forward scattering. In addition, the matrix element F₁₂ shows a maximum at a scattering angle of approximately 150° which is the characteristic rainbow feature of water droplets.

The normalized reflected flux and the degree of linear polarization as a function of the phase angle for various cloud fractions are shown in Fig. 2. Our numerical results are compared with the results by Karalidi et al. (2012). We find a good agreement for the reflected flux (relative difference of |δI| < 0.01 at α = 0°), but some deviation for the degree of linear polarization (relative difference of |δP_l | < 0.05 at α = 90°) between ournumerical calculations and the results by Karalidi et al. (2012). The deviation is larger for the smaller cloud fraction, so we assume that the deviation in the degree of polarization is either due to our random cloud coverage since equal cloud fractions can produce slightly different observing signals (Rossi & Stam 2017) and/or due to the additional horizontal transport of radiation at the terminator if the atmosphere is optically thin.

Fig. 3 Normalized matrix element F11 (solid lines) and the ratio − F12∕F11 (dashed lines) as a function of the scattering angle Θ = 0°: forward scattering) for molecular hydrogen (Rayleigh scattering with ρd = 0.02), air (Rayleigh scattering with ρd = 0.028), water cloud particles (Mie scattering, see Sect. 3.2 for details), and water-ice particles (Mie scattering, see Sect. 4 for details). Theelement F11 is normalized such that ∫F11 d Ω = 1. The solid orange line (H2) is not visible because it coincides with the solid light blue line (air). The plotted matrix elements of the cloud particles and of the water-ice particles are averaged over all grain sizes and are weighted by the particle size distribution using Eq. (14).

4 Impactof a circumplanetary ring

As a case study we investigate the impact of a circumplanetary ring on the net reflected flux and polarization. Circumplanetary rings are expected to be detectable through their influence on the light curve of the host star during a transit event of such a planet (e.g., Barnes & Fortney 2004). Ohta et al. (2009) derived specific predictions for corresponding photometric and spectroscopic signatures. While these studies are limited to the case of transiting planets, Sucerquia et al. (2020) developed an analytical model to estimate the photometric signatures of non-transiting exorings due to scattering. However, this model is restricted to the assumption of diffuse reflection by both the planetary atmosphere and the ring, so no preferential direction of scattered light is expected. However, the atmosphere and the circumplanetary ring consist of particles that have characteristic scattering properties (e.g., described by Mie scattering in the case of cloud or dust particles).

While the rings of Saturn mainly consist of water-ice (Nicholson et al. 2008), ices only exist if the separation to the star satisfies the relation $$d_{\star }\gtrsim {\left(\frac{L_{\star }}{16\pi \sigma T_{\text{sub}}^4}\right)}^{1/2}=2.7\text{au}{\left(\frac{L_{\star }}{L_{\odot }}\right)}^{1/2},$$(15)

where L_⋆ is the luminosity of the central star, σ the Stefan-Boltzmann constant, and T_sub = 170 K the sublimation temperature of water-ice (Gaudi et al. 2003). Therefore, rings of close-in planets (e.g., hot Jupiters) are expected to consist of dusty and rocky material. Following this line of reasoning, we choose a separation of d_⋆ = 3 au between the planet and the central star, and for the sake of simplicity water-ice particles as the only material in the circumplanetary ring with optical properties based on the complex refractive index taken from Warren & Brandt (2008). The ring ranges from 7.7 × 10⁷ to 1.5 × 10⁸ m, measured from the center of the planet, to mimic the width of Saturn’s main ring system (Colwell et al. 2009). With an opening angle of 0.2″ of the ring, we satisfy the upper limit of 200 m of Saturn’s ring thickness observed by Lane et al. (1982). The particles have a power-law size distribution N(r) ∝ r^−3.5, as expected to result from a collisional cascade (Dohnanyi 1969). The grain radii range from 0.1 μm to 1 mm. In this simple model the number density of water-ice particles is constant throughout the ring. The derived matrix elements F₁₁ and F₁₂ of the Müller matrix of water-ice are shown in Fig. 3. The particles show strong forward scattering. The high single scattering albedo of ω≈ 1 at the considered wavelength of 550 nm leads to a low absorption by the water-ice particles.

The optical depth of the ring, measured throughout the vertical coordinate of the ring at the outer edge, amounts to $${\tau }_{\text{r}}=C_{\text{ext}}{n}_{\text{r}}{H}_{\text{r}}.$$(16)

Here n_r is the number density of the water-ice particles and H_r the total height of the ring at the outer edge. We investigate an optically thin ring with τ_r = 0.1, a ring with an optical depth of τ_r = 1, and an optically thick ring with τ_r = 5. This corresponds to relatively low ring masses of M_r ≈ 2 × 10⁻⁶ M_S,r, M_r ≈ 2 × 10⁻⁵ M_S,r, and M_r ≈ 1 × 10⁻⁴ M_S,r, respectively, with M_S,r = 1.54 × 10¹⁹ kg being the mass of the rings of Saturn (Iess et al. 2019). However, the low mass ratio is due to the much smaller grains compared to those in the rings of Saturn (Zebker et al. 1985).

We assume a simple single layer homogeneous atmosphere consisting of molecular hydrogen (H₂) without any cloud particles for the planetary model. This model is referred to as model III. The scattering properties of molecular hydrogen are described by the Rayleigh scattering theory with a depolarization factor of 0.02 (Hansen & Travis 1974). We only consider scattering at gas particles and ignore absorption. The planetary radius is set to R_p = 7 × 10⁷ m with a height of the atmosphere of H_atm = 10⁵ m. The optical depth of the atmosphere is set to τ_g = 5 above a diffuse reflecting surface with an albedo of ω_s = 1.

Taking various inclinations of the ring i with respect to the orbital plane of the planet into account, we first consider an observer who is, as before, located in the plane of the planetary orbit (edge-on orbit). Subsequently, we choose an observer whose line of sight is perpendicular to the planetary orbit (face-on orbit). As half of the planetary disk is illuminated throughout the entire orbit of the planet, any variations in the (polarized) reflected flux are due to the circumplanetary ring. In addition, the polarization degree resulting from scattering in the atmosphere of the planet is largest since the scattering angle (in the case of single scattering) is approximately 90° towards the observer.

Fig. 4 Normalized reflected flux I and degree of linear polarization Pl of the planet–ring system for various inclination angles i (see Sect. 4.1 for details). The planetary orbit is seen edge-on by the observer. Following cases are presented: an optically thin ring with τr = 0.1 (top pair), a ring with an optical depth τr = 1 (middle pair) and an optically thick ring with τr = 5 (bottom pair). The solid black line pertains to a planet without a ring. The color coding applies to all figures.

Fig. 5 Spatially resolved maps of the reflected flux (top) and degree of linear polarization (bottom) for phase angles (from left to right). The direction and length of the white vectors in the lower figures represent the angle and degree of the linear polarization, respectively. The vectors are oriented perpendicular to the incoming radiation since the matrix element F12 is negative for both molecular hydrogen and water-ice particles at these scattering angles (see Fig. 3). The scale of the vector length is displayed on the right side between the figures. The ring has an inclination of i = 30° and an optical depth of τr = 1. The planetary orbit is seen edge-on by the observer. See Sect. 4.1 for details.

4.1 Edge-on orbit

Figure 4 shows the normalized reflected flux I and the degree of linear polarization P_l of the planet–ring system as a function of the phase angle α. As reference, the reflected flux and degree of linear polarization of the planet without a ring is represented as a solid black line.

For an inclination i = 0° the ring is in the plane of the orbit (i.e., seen edge-on by the observer). As the ring has a small vertical extension, the reflected flux and degree of linear polarization is similar to that of the planet alone. This is independent of the optical depth as well. However, if the ring is inclined, the illuminating and viewing conditions change. For all inclination angles, the strongest impact of the ring on the reflected flux is found at larger phase angles. Here, the reflected flux I strongly increases. While the reflected radiation from the planet is very small at large phase angles, the ring has a strong impact due to the forward scattering of water-ice particles. In addition, the degree of linear polarization decreases since the degree of polarization is small after single scattering at water-ice particles (see Fig. 3). If the ring faces the observer (i = 90°), the deviation of the reflected flux and degree of linear polarization compared to the case of a planet without a ring is largest because the observed surface area is largest at this inclination angle.

At small phase angles, the reflected radiation is dominated by radiation scattered in the planetary atmosphere. Thus, the reflected flux and degree of linear polarization is compareable to the case without a circumplanetary ring. While for an optically thin ring (τ_r = 0.1) most of the incoming radiation penetrates the ring, the amount of back-scattered radiation increases with increasing optical depth (τ_r = 1 and τ_r = 5). Thus, the reflected flux increases at smaller phase angles. In addition, the degree of linear polarization P_l decreases due to the increased likelihood of multiple scattering combined with a small polarization after single scattering at water-ice particles.

In addition to the properties of the particles, shadowing effects have an impact on the reflected (polarized) flux as well (see Fig. 5 where spatially resolved images for i = 30° and τ_r = 1 are shown). As the planet orbits the central star, the illuminated areas of planet and ring change, and thus the net polarization of the reflected flux changes as well. The planet always casts shadows on the circumplanetary ring if the inclination angle is small enough. In addition, the ring casts shadows on the planetary body (shown for α = 120° in Fig. 5 ), resulting in a lower reflected flux. Furthermore, the shadow of the ring also covers parts of the planetary body, and thus lowers the degree of linear polarization, for example at α = 90° where the polarization is at its maximum due to the Rayleigh scattering in the atmosphere.

4.2 Face-on orbit

In Fig. 6 the normalized reflected flux I and the degree of linear polarization P_l of the planet–ring system are shown. Similar to the previous case, we consider different inclination angles i of the ring with respect to the orbital plane of the planet and for different optical depths τ_r. As reference the reflected flux and degree of linear polarization of the planet without a ring is represented as a solid black line.

For an optically thin ring (τ_r = 0.1 in our model setup), the difference of the reflected flux I and degree of linear polarization P_l compared to a planet without a ring becomes negligible. The radiation reflected by the planet dominates the flux since most radiation penetrates the ring without scattering. Subsequently, the impact on the net flux is negligibly small for this geometrical setup combined with an optically thin ring.

For the case of τ_r = 1 and τ_r = 5, the reflected flux increases at smaller and larger phase angles. In contrast, the degree of polarization decreases. At α = 90°, the reflected flux and degree of linear polarization have its minimum and maximum, respectively, because the ring is illuminated from the side. Thus, the impact on the reflected (polarized) flux is lowest. This is also true if the ring is not inclined (i = 0°). Here, both I and P_l are constant throughout the entire orbit of the planet. In addition, the reflected flux is higher at small phase angles compared to larger phase angles. This is because in our geometrical setup the radiation has to penetrate the ring at larger phase angles, which lowers the observed flux with increasing optical depth.

Figure 7 shows spatially resolved images for i = 30° and τ_r = 1. As mentioned before, the reflected polarized radiation of the planet itself does not change while it orbits the central star. However, similar to the previous case, the planet can cast shadows on the circumplanetary ring, and vice versa. In addition, the total reflected flux and degree of polarization also decrease if the ring covers the planetary body.

4.3 Discussion

The polarization contrast C_pol, which is the ratio of the polarized flux to the total stellar flux (see, e.g., Hunziker et al. 2020), depends on the planetary radius R_p and the distance to the central star d_⋆: $$C_{\text{pol}}=P(\alpha )I(\alpha )\frac{R_{\text{p}}^2}{d_{\star }^2}.$$(17)

Here I(α) is the reflectivity (i.e., the normalized reflected intensity) and P(α) the degree of polarization. For our considered model (i.e., R_p = 7 × 10⁷ m and d_⋆ = 3 au) the ratio R_p ∕d_⋆ is approximately 2 × 10⁻⁸. Including polarization, the polarization contrast amounts to 10⁻⁹, which is out of reach for existing polarimetric devices (e.g., SPHERE; Beuzit et al. 2019). So far, only close-in planets (e.g., hot Jupiters) would produce a sufficiently high polarization contrast. In this case the observed polarized flux can reveal circumplanetary rings of extra-solar planets. However, rings around these close-in planets are expected to have relatively short lifetimes because multiple forces are at work (see, e.g., Goldreich & Tremaine 1982), such as the radiation pressure (Poynting-Robertson effect).

In addition to geometrical parameters characterizing the planet–ring and/or planet–star system investigated here, the reflected polarized flux of the circumplanetary ring also depends on the properties of the dust phase constituting the ring; however, this is outside the scope of the present study.

Fig. 6 Similar to Fig. 4, but the planetary orbit is seen face-on by the observer (see Sect. 4.2 for details). The color coding applies to all figures.

Fig. 7 Similar to Fig. 5, but the planetary orbit is seen face-on by the observer. See Sect. 4.2 for details.

5 Summary

In this study we presented a numerical tool that simulates the 3D polarized radiative transfer in planetary atmospheres based on pressure and temperature profiles; various light scattering species such as various types of gas, condensates, and dust; as well as the planetary surface that is absorbing or reflecting, and (de)polarizing the incoming radiation. The simulation software is based on the publicly available radiative transfer code POLARIS (Reissl et al. 2016) and has been equipped with dedicated numerical routines described in this paper. Our approach avoids the various simplifications commonly used in previous studies, for example the locally plane-parallel planetary atmosphere, the missing cross-talk between latitudinal and longitudinal regions, a point-like star, or plane-parallel illumination. In particular, we are able to consider a spatially extended radiation source and take an inhomogeneous brightness distribution into account, for example due to the effect of limb darkening or stellar spots. In addition to describing vertical and horizontal inhomogeneities in the structure of the atmosphere and on the surface, scattering material in the local planetary environment, such as planetary rings and exomoons can also be included. Furthermore, the code has been optimized to handle very small planetary cross sections.

We tested our numerical concepts and calculated the reflected and polarized flux for different types of planetary models. This includes an atmosphere-free planet with a Lambertian reflecting surface only, planets with a purely gaseous atmosphere considering various optical depths and surface albedos, and a planet with additional cloud particles in the atmosphere. The numerical results show very good agreement with already existing (analytically derived) computations for these atmospheric models.

Furthermore, we presented a study in which we investigated the influence of a circumplanetary ring. We considered the cases of an observer located in the plane of the planetary orbit and an observer located perpendicular to it. This study shows that an additional circumplanetary ring can have an essential impact on the observed polarized radiation.

For the transiting planet, the flux strongly increases at larger phase angles due to the forward scattering of water-ice particles. However, the degree of linear polarization decreases at these phase angles because the polarization degree after single scattering at water-ice particles is small. This is also true if the ring is not inclined (i = 0°). Here, the reflected flux and degree of linear polarization is similar to the case of a planet without a ring.

If the orbit of the planet is seen face-on by the observer, then a change in the reflected (polarized) flux is due only to the ring. The reflected (polarized) flux for a ring optical depth of τ_r = 0.1 differs slightly from the results of the same planetary model without a ring. For an inclined ring with an optical depth of τ_r = 1 and τ_r = 5, however, the reflected flux increases while the degree of linear polarization decreases. The characteristic profile is due to the additional scattering inside the ring and various shadowing effects on the planet and the ring as well. However, if the ring is not inclined, the observed polarized radiation remains constant as the planet orbits the star.

Acknowledgements

R.B. thanks the DFG for financial support under contract WO857/18-1. We thank the anonymous referee for very useful suggestions.

Appendix A Boundaries for restricted emission

The planet is located in the origin of the model space, while the radiation source is located on the positive z-axis. Thus, the first polar angle θ₁ is restricted from $${\theta }_{\text{1},\text{min}}=\pi -\arccos \left(\frac{R_{\star }-R_{\text{p}}}{d_{\star }}\right)$$(A.1)

to π. For the second polar angle θ₂, we determine the angle θ_c that is needed to rotate to the center of the planet, i.e, the center of the model space (see Appendix B for the corresponding calculation). The allowed range of values is then given by $${\theta }_{\text{2},\text{min}}={\theta }_{\text{c}}-\delta \theta ,{\theta }_{\text{2},\text{max}}={\theta }_{\text{c}}+\delta \theta ,$$(A.2)

where δθ is the deviation depending on R_p and d_⋆ with sin (δθ) = R_p∕d_⋆.

For the second azimuthal angle ϕ₂, we set a newcoordinate system $(\widehat{x},\widehat{y},\widehat{z})$ originating at the location of the emission point, i.e., the apex of the emission cone, with the ẑ -axis parallel to the surface normal of the stellarsurface and the $\widehat{x}$- and ŷ-axis perpendicular to the surface normal. In Fig. A.1 this coordinate space is indicated in orange. The ŷ -axis points into the drawing plane. In this new coordinate system the center of the planet is at $({\widehat{x}}_0,{\widehat{y}}_0,{\widehat{z}}_0)=(p\sin {\theta }_{\text{c}},0,p\cos {\theta }_{\text{c}})$, where p = |p| and p is the position vector of the origin of the new coordinate system in the global 3D model space. The height h_cone of the cone where the azimuthal angle ϕ₂ is largest, the radius r_cone at this height, and the length of the surface line l_cone at this height (see Fig. A.2) can be calculated by $$\begin{array}{l}{l}_{\text{cone}}=p\cos ({\theta }_{\text{c}}-{\theta }_2),\\ r_{\text{cone}}=l_{\text{cone}}\sin ({\theta }_2),\\ h_{\text{cone}}=l_{\text{cone}}\cos ({\theta }_2).\end{array}$$

The radius of the planet at this height is given by $$r_{\text{p}}^2=R_{\text{p}}^2-{(h_{\text{cone}}-{\widehat{z}}_0)}^2.$$(A.6)

The 3D problem is now reduced to a 2D problem where we have to consider an intersection of a circle with radius r_cone with a circle of radius r_p in the $(\widehat{x},\widehat{y})$-plane of the new coordinate space: $$\begin{array}{l}{\widehat{x}}^2+{\widehat{y}}^2=r_{\text{cone}}^2,\\ {(\widehat{x}-{\widehat{x}}_0)}^2+{\widehat{y}}^2={\widehat{x}}^2-2\widehat{x}{\widehat{x}}_0+{\widehat{x}}_0^2+{\widehat{y}}^2=r_{\text{p}}^2.\end{array}$$

Fig. A.1 Photon package position (defined by θ1 and ϕ1) on the stellar surface and distance p to the center of the planet. The emission cone (light blue line) with opening angle θ2 intersects with the planet in the range [ϕ2,min, ϕ2,max]. The new coordinate space at the location of the photon package is indicated in orange. See Appendix A for details.

Fig. A.2 Intersection of the emission cone with opening angle θ2 and the planet’s sphere. The maximum range of the azimuthal angle ϕ2 is at the height hcone of the new coordinate space where the planet has a radius rp and the conea radius rcone. See Appendix A for details.

Putting Eq. (A.7) into Eq. (A.8) leads to $$r_{\text{cone}}^2-r_{\text{p}}^2+{\widehat{x}}_0^2=2\widehat{x}{\widehat{x}}_0.$$(A.9)

By using polar coordinates in the new coordinate space with $\widehat{x}=r_{\text{cone}}\cos (\phi )$, the equation can be rewritten as $$\cos (\phi )=\frac{r_{\text{cone}}^2-r_{\text{p}}^2+{\widehat{x}}_0^2}{2r_{\text{cone}}{\widehat{x}}_0},$$(A.10)

where ± φ is the maximum deviation of the azimuthal angle ϕ₂. Finally, the boundaries for the azimuthal angle are $${\varphi }_{\text{2},\text{min}}={\varphi }_{\text{c}}-\phi ,{\varphi }_{\text{2},\text{max}}={\varphi }_{\text{c}}+\phi .$$(A.11)

Similar to the polar case, ϕ_c is the angle to rotate to the center of coordinate space (see Appendix B).

Appendix B: Determining the angles of rotation

If the direction of the incoming photon package and the direction of the outgoing photon package are given, we can calculate the angles that are needed to perform this rotation. This calculation is necessary, for example, to rotate to the center of coordinate space (ϕ_c and θ_c, see Appendix A).

Figure B.1 shows the rotation of the incoming (r_in, l_in, d_in) photon package frame into the outgoing (r_out, l_out, d_out) photon package frame. The angle θ_out can be calculated by using $$\cos {\theta }_{\text{out}}=d_{\text{in}}\cdot d_{\text{out}}.$$(B.1)

The azimuthal angle ϕ_out is determined by two equations due to the full range from 0 to 2π. With the relation $$r_{\text{out}}=-\frac{d_{\text{in}}\times d_{\text{out}}}{|d_{\text{in}}\times d_{\text{out}}|},$$(B.2)

the first equation is $$\sin {\varphi }_{\text{out}}=(r_{\text{in}}\times r_{\text{out}})\cdot d_{\text{in}}=-\frac{r_{\text{in}}\cdot d_{\text{out}}}{|d_{\text{in}}\times d_{\text{out}}|}.$$(B.3)

Fig. B.1 Rotation of the frame of the photon package before the scattering event (orange) into the frame after the scattering event (light blue). The planes represent the planes that are perpendicular to the r-axis (i.e., the scattering planes). See Appendix B for details.

The second equation is $$\cos {\varphi }_{\text{out}}=r_{\text{in}}\cdot r_{\text{out}}=\frac{l_{\text{in}}\cdot d_{\text{out}}}{|d_{\text{in}}\times d_{\text{out}}|}.$$(B.4)

References

All Tables

All Figures

Fig. 1 Illustration of the model space. The planet with radius Rp is located at the center of the model space, i.e., the origin of the coordinate system. The dashed lines represent the radial cell boundaries that divide the atmosphere above the surface into individual layers. The star, located at distance d* along the positive z-axis, has a radius of R*. The first two emission angles are sampled and define the starting position of a photon package on the stellar surface. The position is described by the polar angle θ1 with respect to the z-axis and the azimuthal ϕ1 in the (x, y)-plane. The second polar angle θ2 is defined with respect to the surface normal and the azimuthal angle ϕ2 is defined inthe surface plane at this location.

In the text

Fig. 2 Normalized reflected flux (top) and degree of linear polarization (bottom) as a function of the phase angle α for various atmospheric models τg = 5, 0, 0.097; ωs = 1, 0.5, 0; fc = 0, 0.3, 0.7). Shown are the cloud-free homogeneous atmosphere (model I; orange and light blue lines), the Lambertian reflecting sphere without an atmosphere τg = 0; green line), and the cloudy inhomogeneous atmosphere (model II; yellow and dark blue lines) with various cloud fractions fc. Our numerical calculations (plus symbols) are compared with results by Buenzli & Schmid (2009) for the homogeneous atmosphere, with the analytical solution by Russell (1916) for the Lambertian reflecting sphere, and with the results by Karalidi et al. (2012) for the inhomogeneous atmosphere. The color-coding applies to both panels. See Sects. 3.1 and 3.2 for details.

In the text

Fig. 3 Normalized matrix element F11 (solid lines) and the ratio − F12∕F11 (dashed lines) as a function of the scattering angle Θ = 0°: forward scattering) for molecular hydrogen (Rayleigh scattering with ρd = 0.02), air (Rayleigh scattering with ρd = 0.028), water cloud particles (Mie scattering, see Sect. 3.2 for details), and water-ice particles (Mie scattering, see Sect. 4 for details). Theelement F11 is normalized such that ∫F11 d Ω = 1. The solid orange line (H2) is not visible because it coincides with the solid light blue line (air). The plotted matrix elements of the cloud particles and of the water-ice particles are averaged over all grain sizes and are weighted by the particle size distribution using Eq. (14).

In the text

Fig. 4 Normalized reflected flux I and degree of linear polarization Pl of the planet–ring system for various inclination angles i (see Sect. 4.1 for details). The planetary orbit is seen edge-on by the observer. Following cases are presented: an optically thin ring with τr = 0.1 (top pair), a ring with an optical depth τr = 1 (middle pair) and an optically thick ring with τr = 5 (bottom pair). The solid black line pertains to a planet without a ring. The color coding applies to all figures.

In the text

Fig. 5 Spatially resolved maps of the reflected flux (top) and degree of linear polarization (bottom) for phase angles (from left to right). The direction and length of the white vectors in the lower figures represent the angle and degree of the linear polarization, respectively. The vectors are oriented perpendicular to the incoming radiation since the matrix element F12 is negative for both molecular hydrogen and water-ice particles at these scattering angles (see Fig. 3). The scale of the vector length is displayed on the right side between the figures. The ring has an inclination of i = 30° and an optical depth of τr = 1. The planetary orbit is seen edge-on by the observer. See Sect. 4.1 for details.

In the text

	Fig. 6 Similar to Fig. 4, but the planetary orbit is seen face-on by the observer (see Sect. 4.2 for details). The color coding applies to all figures.
In the text

	Fig. 7 Similar to Fig. 5, but the planetary orbit is seen face-on by the observer. See Sect. 4.2 for details.
In the text

	Fig. A.1 Photon package position (defined by θ1 and ϕ1) on the stellar surface and distance p to the center of the planet. The emission cone (light blue line) with opening angle θ2 intersects with the planet in the range [ϕ2,min, ϕ2,max]. The new coordinate space at the location of the photon package is indicated in orange. See Appendix A for details.
In the text

	Fig. A.2 Intersection of the emission cone with opening angle θ2 and the planet’s sphere. The maximum range of the azimuthal angle ϕ2 is at the height hcone of the new coordinate space where the planet has a radius rp and the conea radius rcone. See Appendix A for details.
In the text

	Fig. B.1 Rotation of the frame of the photon package before the scattering event (orange) into the frame after the scattering event (light blue). The planes represent the planes that are perpendicular to the r-axis (i.e., the scattering planes). See Appendix B for details.
In the text