Dynamic mineral clouds on HD 189733b II. Monte Carlo radiative transfer for 3D cloudy exoplanet atmospheres: combining scattering and emission spectra

3.7 MB·

by  Graham Lee

Astronomy&Astrophysicsmanuscriptno.RHD+MCRT (cid:13)cESO2017
January19,2017
Dynamic mineral clouds on HD 189733b
II. Monte Carlo radiative transfer for 3D cloudy exoplanet atmospheres:
combining scattering and emission spectra
G.Lee1(cid:63),K.Wood1,I.Dobbs-Dixon2,A.Rice1,andCh.Helling1
1 CentreforExoplanetScience,SUPA,SchoolofPhysicsandAstronomy,UniversityofStAndrews,NorthHaugh,StAndrews,Fife
KY169SS,UK
2 NYUAbuDhabi,POBox129188,AbuDhabi,UAE
Received:28Sep2016/Accepted:3Jan2017
7
1
0
ABSTRACT
2
n Context.Asthe3Dspatialpropertiesofexoplanetatmospheresarebeingobservedinincreasingdetailbycurrentandnewgenera-
tionsoftelescopes,themodellingofthe3Dscatteringeffectsofcloudformingatmosphereswithinhomogeneousopacitystructures
a
J becomesincreasinglyimportanttointerpretobservationaldata.
Aims.Wemodelthescatteringandemissionpropertiesofasimulatedcloudforming,inhomogeneousopacity,hotJupiteratmosphere
8
ofHD189733b.WecompareourresultstoavailableHubbleSpaceTelescope(HST)andSpitzerdataandquantifytheeffectsof3D
1
multiplescatteringonobservablepropertiesoftheatmosphere.WediscusspotentialobservationalpropertiesofHD189733bforthe

upcomingTransitingExoplanetSurveySatellite(TESS)andCHaracterisingExOPlanetSatellite(CHEOPS)missions.
]
P Methods.WedevelopedaMonteCarloradiativetransfercodeandappliedittopost-processoutputofour3Dradiative-hydrodynamic,
E cloudformationsimulationofHD189733b.Weemployedthreevariancereductiontechniques,i.e.nexteventestimation,survival
biasing,andcompositeemissionbiasing,toimprovesignaltonoiseoftheoutput.Forcloudparticlescatteringevents,weconstructed
.
h alog-normalareadistributionfromthe3Dcloudformationradiative-hydrodynamicresults,whichisstochasticallysampledinorder
p tomodeltheRayleighandMiescatteringbehaviourofamixtureofgrainsizes.
- Results.StellarphotonpacketsincidentontheeasterndaysidehemisphereshowpredominantlyRayleigh,single-scatteringbehaviour,
o
whilemultiplescatteringoccursonthewesternhemisphere.Combinedscatteredandthermalemittedlightpredictionsareconsistent
r
with published HST and Spitzer secondary transit observations. Our model predictions are also consistent with geometric albedo
t
s constraintsfromopticalwavelengthground-basedpolarimetryandHSTBBandmeasurements.Wepredictanapparentgeometric
a
albedoforHD189733bof0.205and0.229,intheTESSandCHEOPSphotometricbandsrespectively.
[ Conclusions.Modellingthe3Dgeometricscatteringeffectsofcloudsonobservablesofexoplanetatmospheresprovidesanimportant
2 contributiontotheattempttodeterminethecloudpropertiesoftheseobjects.ComparisonsbetweenTESSandCHEOPSphotometry
v mayprovidequalitativeinformationonthecloudpropertiesofnearbyhotJupiterexoplanets.
3
Key words. radiative transfer – planets and satellites: atmospheres – planets and satellites: individual: HD 189733b – methods:
8
numerical
9
0
0
1. Introduction ometric albedos (A > 0.1) in the Kepler bandpass (0.43-0.89
. g
1 µm),suggestingthepresenceofanopticaltonear-IRscattering
0The 3D spatial properties of exoplanet atmospheres are being
aerosol. B Band (0.29-0.45µm) and V Band (0.45-0.57µm) ge-
7observedinincreasingrateanddetail.Twokeyvariableforex-
ometric albedo measurements of HD 189733b by Evans et al.
1oplanet characterisation are the properties of their atmospheric
(2013) using HST STIS show a blueward slope in the optical,
:cloudcomplexandhowcloudparticlesize,compositions,global
v inferring the presence of wavelength-dependent, backscattering
idistributions and opacity affect the resultant observations. A cloudparticles.
XwestwardoffsetintheopticalphasecurvesofKepler-7bwasob-
Becauseofitsproximitytoitshoststar,HD189733bislikely
rservedbyDemoryetal.(2013)whichwasinterpretedascloud
a to be tidally locked, with one hemisphere in constant daylight.
particles that backscatter optical photons on the western hemi-
Infrared flux phase curve observations of HD 189733b using
sphere. But a reduced or absent cloud particle abundance on
Spitzer(Knutsonetal.2007;Charbonneauetal.2008;Knutson
the eastern hemisphere. More Kepler planets (e.g.. Kepler-12b,
etal.2009;Agoletal.2010;Knutsonetal.2012)inphotometry
Kepler-8b, Kepler-41b; Angerhausen et al. 2015; Esteves et al.
bandscentredat3.6µm,4.5µm,8.0µmand24µmshoweastward
2015;Shporer&Hu2015)werefoundtoexhibitsimilarphase
curve behaviour, revealing that offset optical phase curves may offsetsinpeakfluxfromthesub-stellarpoint.Globalcirculation
models and radiative-hydrodynamic models of this hot Jupiter
be a common feature of hot Jupiter atmospheres. Some Kepler
(e.g.Showmanetal.2009;Dobbs-Dixon&Agol2013;Kataria
fieldplanets(e.g..HAT-P-7b,Kepler-5b,Kepler-7b,Kepler-17b,
etal.2016)suggestthedevelopmentofanequatorialsupersonic
Kepler-41b, Keper-76b; Heng & Demory 2013; Angerhausen
jet that efficiently advects energy density from the dayside to
et al. 2015; Esteves et al. 2015) also have non-negligible ge-
nightside. Horizontally propagating gravity waves (Watkins &
(cid:63) E-mail:[email protected] Cho2010;Perez-Becker&Showman2013)havealsobeensug-
Articlenumber,page1of17
A&Aproofs:manuscriptno.RHD+MCRT
gestedtocontributetothedayside-nightsideenergydensityad- sphere(Hellingetal.2016b),importanttotheionchemistryin
vection.Theseprocessesresultinaneast-westhemisphereinho- exoplanetatmospheres(Rimmer&Helling2016).
mogeneity in the temperature profiles for HD 189733b, in line Inthisstudy,wepresenta3DMonteCarloradiativetransfer
withtheSpitzermaximumfluxoffsets. simulation for use in post-processing GCM/RHD atmospheric
These observational and modelling investigations suggest results with a strong scattering component such as clouds. Our
that, for cloud forming hot Jupiters, optical and near-IR wave- paperoutlineisasfollows;Section2explainsourMCRTmodel
length instruments such as HST STIS and Kepler are more framework and our approach to modelling the HD 189733b at-
sensitive to reflected/scattered light from cloud particles in the mosphereusingMCRT.Section3detailsourmodelset-up,input
atmosphere; while longer wavelength IR instruments such as quantities and equations for relating the MCRT output to ob-
HST WF3, NICOMS, and Spitzer are more sensitive to the servable quantities. In Sect. 4 we present our model results for
thermal emission from the planet itself. Future optical to near- scatteredandemittedlightaswellascombiningthetwocompo-
IR wavelength photometry observational programmes such as nents.WealsoproduceKepler,TESS,andCHEOPSphotomet-
TESS (0.6-1.0µm) and CHEOPS (0.4-1.1µm), will likely char- ric band predictions for HD 189733b based on our simulation
acteriseacombinationofbothscatteredlightandthermalemis- output.Section5containsthediscussionandSect.6summarises
sionfornearbyhotJupiters.Examiningdatafrombothofthese keyresultsandconclusions.
instrumentscouldthereforeprovideinformationaboutthecloud
propertiesoftheiratmospheres,similartostudiesoftheKepler
fieldplanets. 2. Approach
In recent years, theoretical and modelling efforts have ex-
We present a MCRT code based on the Hood et al. (2008)
plored into exploring how the properties of cloud particles in
model framework, with significant additions and alterations to
exoplanet atmospheres affect their observed phase curves and
theschemetocaptureamorerealistic3DhotJupiteratmosphere.
geometric/spherical albedos. Marley et al. (1999) used 1D at-
InLeeetal.(2016)(PaperIinthisseries)wesimulatedtheat-
mospheric models, combined with a parameterised equilibrium
mosphere of HD 189733b using a 3D radiative-hydrodynamic
cloudmodeltopredictthealbedospectraofexoplanets.Cahoy
[RHD] model coupled with a self-consistent cloud formation
etal.(2010)updatedandextendedtheMarleyetal.(1999)ap-
module.Inthepresentpaper,wepost-processourcloudforming
proach and presented model scattering albedos for Jupiter and
RHDresultsusingtheMCRTmodeltocompareoursimulation
Neptune-likeexoplanetsatdifferentorbitaldistancesfromtheir
to available observational data and make observational predic-
host stars. Webber et al. (2015) employed a similar method to
tionsfortheTESSandCHEOPSphotometricbands.
Cahoyetal.(2010)tomodelawest-eastpatchycloudscenario
forKepler-7btofittheoffsetopticalphasecurve(Demoryetal.
2013).Morleyetal.(2015)investigatedthescatteringandther- 2.1. MCRTmodellingframework
mal emission spectra for photochemical haze and cloud form-
MonteCarloradiativetransferisamicrophysical,stochasticap-
ingSuper-Earthatmospheres,comparingtheirmodelstotheflat
proach to solving radiative transport problems. MCRT simu-
transmissionspectraofGJ1214b(Kreidbergetal.2014).Kop-
lates the random walk of a luminosity packet (henceforth "L-
parlaetal.(2016)appliedamultiplescatteringvectorradiative
packet"or"packet")throughamediumanddescribestheinter-
transferschemetointerpretthepolarimetryobservationsofHD
actions(scattering,absorption)ofthepacketwiththesurround-
189733bbyWiktorowiczetal.(2015).
ingmediumbyprobabilisticsampling.InMCRT,eachinitialised
Recently, studies using 3D global circulation models
packet carries a luminosity fraction (cid:15) / ∆t [erg s−1 ], propor-
[GCMs]haveexploredmodellingtheatmosphericreflectedlight 0
tional to a source (e.g.. parent star or planetary atmosphere) of
phase curves of Kepler field objects. Oreshenko et al. (2016)
luminosityL [ergs−1].IfNL-packetsareinitialisedbythe
performed a GCM experiment using the orbital parameters of source
source of luminosity, the luminosity carried per packet is then
Kepler-7b, and applied an equilibrium cloud scheme to model
(e.g.Lucy1999)
the optical and infrared phase curves offsets of the planet. Par-
mentieretal.(2016)performedseveralGCMmodelexperiments (cid:15) L
over a grid of equilibrium temperatures (Teq = 1000−2200 K), ∆0t = soNurce, (1)
correspondingtotheKeplerhot-Jupiterplanetparameterspace.
Theypost-processedtheirmodeloutputandappliedanequilib- where ∆t [s] is the time over which the MCRT experiment is
riumcloudschemewithatwo-streamradiativetransfermethod performed;usually∆t=1sfortimeindependentMCRTsuchas
tomodelthereflectivityofcloudsacrosstheKeplerbandpass. the current study. By tracking the proportion of the luminosity
Modelling the effects of cloud particles on observable ofeachL-packetthatescapesduringthesimulationinacertain
properties of exoplanets using Monte Carlo radiative transfer direction,thetotalluminosityescapingtowardsaparticularob-
[MCRT] techniques began with Seager et al. (2000), which servationdirectioncanbefoundbysummingupthecontribution
mainly focused on producing phase and polarised light curves ofeachNL-packettothetotalescapingluminosity.Ifthetotal
from a model atmosphere. Hood et al. (2008) applied a pa- luminositybudgetinthesimulationisnormalisedtothesource
rameterised cloud model to a 3D idealised atmosphere of HD luminosity (i.e. L = 1), the total contribution of all pack-
source
209458b and modelled scattered light curves and geometric ets is then the fraction (i.e. f ∈ (0, 1)) of source luminosity
tot
albedos of varying cloud heights and scattering properties. de thatescapedthesimulationboundaries.Toretrievetheluminos-
Kok & Stam (2012) used a 3D Monte Carlo scattering code ity that escaped, the source luminosity is calculated and multi-
to investigate the influence of forward scattering particles on pliedbythefractionofluminositythatescapedcarriedbytheL-
transit spectroscopy. Garcia Munoz & Isaak (2015) used a packets.Thisisdifferentfromnumericallysolvingtheradiative
Pre-conditioned Backward Monte Carlo method to constrain transfer equation via, for example, the two-stream approxima-
the cloud particle scattering properties of Kepler-7b. Recently, tion(e.g.Toonetal.1989;Marley&McKay1999;Fortneyetal.
Monte Carlo transport methods have also been used to model 2006) or ray-tracing method (e.g. Rijkhorst et al. 2006). Since
the path and decay of cosmic rays through a Jupiter-like atmo- MCRT tracks the random walk of a packet and its interactions,
Articlenumber,page2of17
G.Leeetal.:DynamicmineralcloudsonHD189733b
radiative transfer through complicated 3D geometries, inhomo- 2.2. L-packetinteractionswithdustandgas
geneousopacities,andhighlymultiplescatteringregionscanbe
Oncethecoordinatesspecifiedbytheopticaldepthsamplingand
modelled.MCRTmethodshavebeenextensivelyusedtoinves-
directionoftravelhavebeenreachedbythepacket,ascattering
tigateradiativetransferindustymedia,forexample,protoplan-
orabsorptioneventisdeterminedstochastically.Incloudform-
etary disks (e.g. Whitney et al. 2003; Harries et al. 2004; Pinte
inghotJupiteratmospheres,thepacketinteractswithtwocom-
etal.2006;Minetal.2009),theISM(e.g.Robitaille2011)and
ponents: the gas and cloud particles. An interacting packet ex-
dustygalacticscaleproblems(e.g.Wood&Loeb2000).MCRT
hibits different scattering and absorption behaviours depending
has also been applied to the study and retrieval of Earth based
what component it is interacting with. The local probability of
cloud properties (e.g. Mayer 2009; Stap et al. 2016). Whitney
thepacketinteractingwiththegasphaseP isgivenby
(2011) and Steinacker et al. (2013) provide a review of MCRT gas
methodsandmodelsusedforastrophysical-basedproblemsand
κ
Mayer (2009) gives details on Earth based applications. The P = ext,gas , (6)
gas κ +κ
heavily irradiated, inhomogeneous opacity, and thermal struc- ext,cloud ext,gas
turesofcloudforminghotJupiteratmospheresisanidealappli-
whereκ [cm2g−1]isthetotalopacityofthegasandκ
cationforaMCRTapproach. ext,gas ext,cloud
ThebeginningoftheMCRTL-packettrackingprocessstarts [cm2g−1](definedbyEq.28)thetotalopacityofthecloudcom-
bystochasticallyassigninginitialstartingcoordinatesanddirec- ponent. Should ζ1 < Pgas, (where the indexed random number
tion of a given packet, depending on the location of the lumi- (e.g.. ζ1,ζ2) denotes a sequence of random, independent, num-
nosity source that is emitting the packet. For example, a stellar bers for a particular scheme) the packet is assumed to interact
L-packetisinitialisedrandomlyatthetopoftheatmosphereon withthegasphase.Thetypeofinteraction,scatteringorabsorp-
thedaysideoftheHD189733bcomputationalgridinourmodel. tionisdeterminedbysamplingthesinglescatteringalbedoωgas
The wavelength-dependent luminosity carried by these stellar of the gas phase, which is the probability of the packet under-
packets is the total monochromatic luminosity incident on the goingascatteringinteraction.SinceweonlyconsideraH2scat-
atmosphere,dividedbythenumberofsimulatedpackets((cid:15) / tering gas component n the MCRT scheme (Sect. 3), the local
0,inc
∆t = Linc,λ / Ninc). The L-packets emitted from the atmosphere singlescatteringalbedoofthegasωgasisgivenby
itselfareassignedarandomstartingpositionwithinacomputa-
κ
tional cell volume, where the luminosity is proportional to the ω = sca,H2, (7)
total monochromatic luminosity of the atmosphere itself ((cid:15)0,atm gas κgas,ext
/ ∆t = L / N ). Initialised stellar packets are assumed to
atm,λ atm
travel in a plane-parallel direction towards the planet, while at- whereκsca,H2 [cm2g−1]isthescatteringopacityofH2andκgas,ext
mospheric,thermallyemittedpacketsaregivenanisotropicini- =κ +κ [cm2 g−1]thetotalextinctionforthegascom-
gas,abs sca,H2
tial direction. The random walk of the L-packet is then deter- ponent.Shouldζ <ω thepacketisRayleighscatteredbyH ,
2 gas 2
mined by stochastically sampling probability distributions that otherwiseitisabsorbedbythegasphase.
governthebehaviourofthepacket(distancetravelled,scattering Iftheinteractioniswiththecloudcomponent(i.e.ζ >P ),
1 gas
directionsetc.). thetypeofinteractionisgivenbysamplingthelocalcloudpar-
From the Beer-Lambert law in radiative transfer theory, ticlesinglescatteringalbedoω givenby
cloud
the probability of a photon packet passing, without interaction,
κ
throughamediumoftotalopticaldepthτλisgivenby ω = scat,cloud, (8)
cloud κ
ext,cloud
P(τλ)=e−τλ, (2)
where κ [cm2 g−1] is the scattering opacity of the cloud
scat,cloud
wheretheopticaldepthτλisdefinedas particles and κext,cloud [cm2 g−1] the cloud particle total extinc-
tion. Should ζ < ω the packet is assumed to be scattered
τ =(cid:90) lmaxρ (l)κ (λ,l)dl, (3) bythecloudco3mponcelonutdtowardsanewdirection,otherwiseitis
λ gas ext,total
absorbedbythecloudparticles.
0
After each scattering event, the next interaction location is
where ρgas [g cm−3] is the local gas density, κext,total [cm2 g−1] determinedbyEq.(5)andthesamplingprocessisrepeated.Itis
(Eq. 26) the local total extinction opacity, including absorption importanttonotethatpacketscanundergomultiplescatteringin-
and scattering, and l [cm] the path length. The probability of a teractionsbeforegettingabsorbed,especiallyinhigh(ω
gas/cloud
packetinteractingwiththesurroundingmediumisthen >0.9)singlescatteringalbedoregions.
P(τλ)=1−e−τλ. (4)
2.3. ScatteringofL-packetsbydustandgas
The MCRT method stochastically samples this probability dis-
tributionwiththeuseofa(pseudo)randomnumber,P(τ )=ζ ∈ Duringascatteringevent,anewdirectionoftravelisgiventothe
λ
packet,whichisdependentonthegeometricandcompositional
[0,1]. An optical depth, determined stochastically, encountered
properties of the interacting material as well as the wavelength
bythepacketbeforeaninteractionisgivenby
oftheinteractingpacket.Akeyquantitytoconsiderwhenmod-
τ =−ln(1−ζ). (5) elling a scattering event is the size parameter of the scattering
λ
particle given by x = 2πa/λ, where a [cm] is the particle size
Using the definition of τ (Eq. 3), the distance that the packet and λ [cm] the wavelength of the L-packet. The scattering be-
λ
travelled through the simulation can then be computed should haviourofaninteractingpacketcanbesplitintotwocategories:
the density and extinction opacity in the path of the packet be
known. After the packet has travelled the distance given by the 1. x(cid:28)1-Rayleighscattering
sampledτ ,aninteractionwiththesurroundingmediumoccurs. 2. x(cid:38)1-Miescattering
λ
Articlenumber,page3of17
A&Aproofs:manuscriptno.RHD+MCRT
For Raleigh scattering, the particle size is magnitudes smaller whereαisthemainlyforwardscatteringcomponentwithg >0,
a
than the wavelength of light interacting with it, resulting in the andβamainlybackscatteringcomponentwithg <0.Forallpa-
b
well-known double lobed scattering phase function where the rameters,therelationα+β=1hastobesatisfied.Weapplythe
probabilityofscatteringintheforwardandbackwarddirections parametersinCahoyetal.(2010)whereg =g ,g =-g /2andα
a λ b λ
areequal.IntheMieregime,theparticlesizeandwavelengthof =1-g2,β=g2.TheTTHGformiscommonlyusedinvariousat-
b b
light are of similar magnitudes; the light is typically more for- mosphericradiativetransferscatteringcodes(e.g..Marleyetal.
wardscatteringwithasmallbackscatteredcomponentatoptical 1999;Cahoyetal.2010;Barstowetal.2014;GarciaMunoz&
wavelengths. For a scattered packet in MCRT, the new travel Isaak2015;Dyudinaetal.2016)andmorequalitativelycaptures
directionisgeneratedstochasticallyfromthescatteringproper- theMiebackscatteringlobeatopticalwavelengths.
ties of the interacting material. This is determined by sampling TodeterminethescatteringangleofaTTHGevent,theprob-
the(normalised)scatteringphasefunctionΦ (θ,φ)oftheinter- abilityofforwardscatteringissimplygivenbyζ<α(Pfeiffer&
scat
acting scattering particle, which is the probability of scattering Chapman2008)afterwhichEq.(13)withg =g canbeused,
λ a
towardsangle(θ,φ). otherwisethescatteringisabackscatteringeventinwhichg =
λ
For small cloud particles with small size parameters x (cid:28) 1 g isapplied.Althoughαandβareusedtodetermineaforward
b
andH gasparticlescattering,weapplytheRayleighscattering or backward scattering event, due to the probabilistic sampling
2
phasefunctiongivenby of the HG function scattering angle, sampling the α term does
notnecessarilyalwaysresultinaforwardscatteringevent,norβ
3
Φ (θ,φ)= (1+cos2θ). (9) abackwardscatteredevent.
RS 16π Since there is a distribution of cloud particle sizes in each
Theθscatteringdirectionissampledrandomlyfromthisdistri- cell,thesizeofanindividualcloudparticleinteractingwiththe
bution,wherearejectionmethodisappliedfollowingtheadvice packetmustbestochasticallydeterminedfromthesizedistribu-
inWhitney(2011).Theφscatteringangleangleisgivenby tionproperties.Thisisrequiredasthesizedistributionmaycon-
tain a combination of cloud particles in the Rayleigh and Mie
φ=2πζ, (10) sizeparameterregimes,resultinginamixtureofpacketscatter-
ingbehaviour.Tocapturethisscatteringbehaviour,weperform
which is a uniform sampling of the φ direction across the 2π
astochasticsamplingfromthecumulativedistributionfunction
radian circular circumference. Rayleigh scattering events are
[CDF] of the size distribution, which is the fractional contribu-
wavelengthindependentandassumedtobeelastic.
tionofeachparticlesizetothetotalareaoftheparticleensemble.
Whenapacketscattersoffacloudparticlewithsizeparam-
Theareadistributionisused,ratherthanthesizedistribution,as
eterx(cid:38)1,weapplyaHenyey-Greenstein[HG]scatteringphase
this better captures the contribution of each cloud particle ra-
function(Henyey&Greenstein1941).TheHGphasefunctionis
dius to the scattering cross-sectional opacity of the size distri-
ananalyticapproximationtotheMiescatteringphasefunction,
bution (since κ ∝ f(a)a2, Eq. 28). The cloud particle mean
givenby sca
grainsize,(cid:104)a(cid:105)[cm],andmeangrainarea,(cid:104)A(cid:105)[cm2],ofthecloud
1 1−g2 particlesineachRHDcellcanbefoundfromthedustmoment
ΦHG(gλ,θ,φ)= 4π[1+g2−2g λcosθ]3/2, (11) values, Lj (Woitke&Helling2003,2004;Leeetal.2015).As-
λ λ suming an arithmetic log-normal distribution of cloud particle
whereg isthewavelength-dependentscatteringasymmetrypa- sizes(e.g.Ackerman&Marley2001),themeanµandstandard
λ
rameter in the range −1 to 1, given by the results of the Mie deviation σ (more specifically their natural logarithms) of the
theory.Thisisdefinedasthemeanscatteringcosineanglefrom ensemble of particles in each cell is then calculated (e.g. Stark
therelation etal.2015).Thesurfacearealog-normaldistributionisgivenby
(cid:90) (e.g.Heintzenberg1994)
gλ =(cid:104)cosθ(cid:105)= ΩΦHG(gλ,θ)cosθdΩ. (12) A(a)= aσA√02πexp−(ln(a2)σ−2µA)2 (15)
Avalueofg<0indicatesapreferenceforpacketbackscattering, A A
g = 0 an equal backward/forward scattering and g > 0 forward which is the log-normal size distribution multiplied by the sur-
facearea4πa2,withthetotalareaA =4πn exp(2µ+2σ2)(e.g.
scattering. The scattering angle cosθ is sampled stochastically 0 d
Zender2015),wherea[cm]isthesampledgrainsize,n [cm−3]
fromthisdistributionusingthefollowingform d
the total cloud particle number density, µ the arithmetic mean
cosθ= 21g 1+g2λ−1−1g−+g2λ2g ζ2. (13) odfistthriebuarteioand.iTsthreibruetliaotnioanndbeσtwAetehnesthtaenmdaeradAndaenvdiasttiaonndoafrdthdeeavrieaa-
λ λ λ tionofthesizedistribution(µ,σ)andareadistribution(µ ,σ )
A A
aregivenby(e.g.Heintzenberg1994)
The sampled φ direction is given by Eq. (10). We assume elas-
ticscatteringforMiescatteringevents.TheHGprobabilitydis- µ =µ+2σ2,
A
tribution has been shown to be a reasonable approximation to σ =σ. (16)
A
Miescatteringwithintheopticalandnear-IRwavelengthregime
The cloud particle area cumulative distribution function of the
(Draine2003).However,asmallbutnotinsignificantprobabil-
log-normal size distribution A ∈ (0,1) is then constructed in
ityofbackscatteringisnotcompletelycapturedbythisapproxi- CDF
eachcellfrom(e.g.Zender2015)
mationatopticalwavelengths(e.g.Kattawar1975;Draine2003;
(cid:34) (cid:35)
Hoodetal.2008;Dyudinaetal.2016).Toaddressthis,weapply 1 1 ln(a)−µ
a Two-term Henyey-Greenstein [TTHG] function (e.g.. Pfeiffer ACDF(a)= 2 + 2erf √2σ A , (17)
&Chapman2008;Cahoyetal.2010;Dyudinaetal.2016)given A
whereerf istheerrorfunction.Weconstruct A for100log-
by CDF
spacedcloudparticlesizeabinsbetweenaminimumseedparti-
ΦTTHG(ga,gb,θ)=αΦHG(ga,θ)+βΦHG(gb,θ), (14) clesizeaseed∼0.001µmandamaximumofamax=10·aeff µm,
Articlenumber,page4of17
G.Leeetal.:DynamicmineralcloudsonHD189733b
whereaeff istheeffectivecloudparticleradius(Eq.29).During stellar(left)andatmosphericemitted(right)luminosityfractions
thesimulation,bysamplingarandomnumberζ ∈(0,1),thepar- at1.0µmforaviewingangleofφ=60◦,θ=0◦.
ticle size from the distribution interacting with the packet can Foratmosphericemission,assumedtobeisotropic,thefrac-
thereforebestochasticallydeterminedbysampling A (i.e.ζ tionofthepacketluminositypeeledofftowardsobservationdi-
CDF
→a). rectionαisgivenby
We assume two scattering regime limits in our model for
1
a certain cloud particle size parameter (x = 2πa/λ). For x < W (α)= e−τ2, (18)
0.1weassumeaRayleighscatteringevent,whileifx≥0.1we po 4π
assume a TTHG scattering event. We assume TTHG scattering where1/4πisthenormalisationfactorforisotropicemissioninto
eventstooccuratthepropertiesoftheeffectivemeanradiusval- 4πsteradiansandτ (λ)(whereindexedτterms(e.g.τ ,τ etc.)
ues (e.g. gλ(aeff), Sect. 3.2). The effect of this scheme is that if indicate optical dep2ths calculated using the ray-tracing1m2ethod
themeanandvarianceskewtheareadistributiontowardslarger to the simulation boundaries, and non-indexed τ denotes the
cloud particle sizes, then the fractional contribution of particle stochasticallysampledopticaldepthforphotonpackets)theop-
areas,wherex≥0.1increases,andsotheprobabilityofaTTHG ticaldepththroughthecomputationalzoneatwavelengthλ,to-
eventisincreased.However,ifthedistributionisskewedtowards wardstheobservationaldirection.
smallercloudparticlesizesthenagreaterfractionofcloudparti- For non-isotropic scattering events such as Rayleigh or
cleareasatisfiesx<0.1,increasingthelikelihoodofaRayleigh Henyey-Greenstein scattering, the fraction of peeled off lumi-
event. nosity is related to the probability of the packet scattering to-
wards the observational direction. For Rayleigh scattering, the
weightofpeeledoffluminosityisgivenby
2.4. MCRTvariancereductiontechniques
W (α)=Φ (α)e−τ2, (19)
In a basic MCRT scheme, to produce observable quantities for po RS
a particular observation direction, the MCRT code would track
andforTwo-termHenyey-Greenstein,
how many packets of a particular luminosity entered the simu-
lation and how many escaped. The ratio of these values would Wpo(α)=ΦTTHG(ga,gb,α)e−τ2. (20)
thengivethefractionofluminosityescapingthesimulationdo-
main towards a particular direction. However, the probability Phase curve information can be calculated by adding the con-
of an individual packet escaping by chance towards a particu- tribution of each event to an image in the desired observation
lardirectioncanbeincrediblysmall,dependingonthegeomet- direction (e.g.. 360◦ in longitude). Each weight is added into a
ricpropertiesofthesimulation,andcanproducenoisyresultsif square image pixel array of (x, y) = (201, 201) pixels and then
notenoughMonteCarlopacketsaresimulated(e.g.seeDupree normalisedbythenumberofinitialisedpacketsafterthesimula-
& Fraley 2002). Therefore, variance reduction techniques were tioniscomplete.Theluminosityscattered/emittedintheobser-
developed for MCRT, which allows for a much improved sig- vationaldirectioncanthenbefoundbymultiplyingthefraction
naltonoiseandgreatercomputationalefficiencythanthebasic of escaping luminosity by the total luminosity of the emitting
scheme. We employ three techniques in this study as follows: source.
nexteventestimation,survivalbiasing,andcompositeemission
biasing.Thevarianceandconvergencepropertiesourmodelout-
2.4.2. Survivalbiasing
putispresentedinSect.4.5.
Fordiffuseemissionfromwithinthesimulationboundaries,the
packetisforcedtoscatterateachinteractionandnotterminated
2.4.1. Nexteventestimation upon absorption (see Sect. 3). In order to conserve energy, the
weightoftheluminositycarriedbythepacketisreducedateach
Next event estimation (e.g. Yusef-Zadeh et al. 1984; Wood &
interactionbythesinglescatteringalbedo
Reynolds1999)isavariancereductionmethodthatusesaray-
tracingtechniqueinconjunctionwiththeMonteCarloscheme. W =ωW, (21)
new
Ateveryscatteringoremissionevent,afractionofthepacket’s
luminosity, W (cid:15) , (Eq. 1) is peeled off towards an observation where,inaninteractionwiththegasphase,ω=ω (Eq.7)and
po 0 gas
directionusingaray-tracingmethod.Thefractionofluminosity withthecloudparticlesω=ω (Eq.8).Thefractionoflumi-
cloud
is given by the normalised phase function of the type of inter- nositylostbythephotonpacketduetoabsorptionisthentaken
action, weighted by the probability of the packet escaping the intoaccountinthepeeledoffimagesasthepacketinteractswith
simulationtowardsthisobservationdirection.Thepeeledofflu- thesurroundingmedium.Sothatapacketdoesnotscatterindef-
minosityiswavelengthdependentandtheluminositycarriedby initelywitheverdecreasingweight,aRussianRoulettescheme
the packet is not altered in the peeling off scheme. The peeled (e.g. Dupree & Fraley 2002) is applied to stochastically termi-
off luminosity is then stored in a 2D pixel image array at this nate packets below a predefined weight cut-off. The packet is
observationdirection,whichcanbeusedtoderiveobservational thengivena1in10chanceofsurviving,withthenewweightof
quantities. thesurvivingpacketsW =10·W.Inourmodel,packetswith
new
Multiplepixelimagescanbeconstructedformanyviewing weightW<10−3areenteredintotheRussianRoulettescheme.
angles and wavelengths, which allows a single photon packet
interaction to contribute to multiple output images at different
2.4.3. Compositeemissionbiasing
planetaryphases.Sinceeverypacketnowcontributestotheout-
put images, this improves the signal to noise of the output as ForhotJupiteratmospheres,thetotaldaysideemissionluminos-
images can be produced without relying on packets emerging ity can be orders of magnitude greater than the emission from
bychancefromthesimulationboundariestowardstheobserver. thenightsideoftheplanet.Innon-biasedMCRT,thenumberof
Figure 1 illustrates the resultant next event estimator scattered packets N emitted from cell i is given by the fraction of the
i,λ
Articlenumber,page5of17
A&Aproofs:manuscriptno.RHD+MCRT
200 200 4.0e-06
180 1.2e-05 180
3.5e-06
160 160
1.1e-05
3.0e-06
140 140
9.0e-06
120 120 2.5e-06
7.5e-06 α) α)
y100 λ, y100 2.0e-06 λ,
6.0e-06 (y (y
80 x, 80 x,
f 1.5e-06 f
4.5e-06
60 60
1.0e-06
40 3.0e-06 40
20 1.5e-06 20 5.0e-07
0 0.0e+00 0 0.0e+00
0 20 40 60 80 100120140160180200 0 20 40 60 80 100120140160180200
x x
Fig.1.Incidentscattered(left)andatmosphericemitted(right)lightimagearrayoftheL-packetluminosityfractionsf (λ,α)∈(0,1)escaping
x,y
fromtheHD189733bsimulation.Theseimagesareproducedbythenexteventestimatormethodat1.0µmforaviewingangle/phaseofα=
60◦(easterndayside).Thesharptransitionfrombluetoblackatpixelnumbersx=120−150istheeasternterminatorinthescatteredlightimage,
includinga"twilight"effectwherepacketsscatterpasttheterminatorline.Theeffectofthegridstructureontheimagesismoreapparentinthe
emittedlight,sincepacketsareinitialisedfromwithinthecellvolumes.
luminosityofthecelltothetotalluminosityoftheatmospheric 3. Inputandoutputquantities
cellsatthiswavelength(e.g.Pinteetal.2006),
In thissection, we providedetails of ourinputquantities to the
MCRT scheme and the relationships between the MCRT out-
L
Ni,λ = (cid:80) iL,λ Natm, (22) put, which is in normalised luminosity fraction units, and ob-
i i,λ servational quantities (geometric albedos, flux ratios, etc.). We
perform all scattering and emission calculations for 250 wave-
(cid:80)
where Li,λ is the luminosity of the cell, iLi,λ the total lumi- lengths in the range 0.3 and 5.0 µm, linear in wave-number
nosity of the considered cells, and Natm the total number of at- space.Table1showstheHD189733starandplanetparameters
mospheric L-packets emitted at wavelength λ. Low luminosity usedinthesimulation.
(cid:80)
regions (e.g.. thenightside), where L (cid:28) L can therefore
i,λ i i,λ
be under sampled in the emission scheme, which increases the Table1.HD189733systemparametersadoptedforthisstudy.Weuse
varianceofresultsderivedfromtheseregions. theparameterspresentedinTorresetal.(2008).
To address this, we implement the multiple component,
Parameter Value
emission composite biasing method from Baes et al. (2016),
R 0.76R
where a linear combination of the unbiased emission (Eq. 22) (cid:63) (cid:12)
probabilitydistributionfunctionandauniform(incellnumber) Teff,(cid:63) 5040K
a 0.031au
distribution function is applied. This results in the number of
R 1.138R
emittedpacketspercellgivenby p Jup
(cid:34) (cid:35)
L η
Ni,λ = Natm (1−η)(cid:80) iL,λ + (cid:80) i , (23)
i i,λ i λ
3.1. Summaryof3DRHDresults
whereη∈(0,1),isaparametergoverningthelinearcombination
(cid:80) We briefly provide a summary of the results of our 3D cloud
and i the total number emitting cells at wavelength λ. To
i λ forming RHD simulation of HD 189733b (Paper I), which are
conservethetotalluminosity,eachemittedpacketisgivenacell
used as input to the MCRT post-processing simulation. The
dependentstartingweightof
mean cloud particle sizes in the upper atmosphere (0.1 - 100
mbar) range from seed particle sizes (cid:104)a(cid:105) ∼ 0.001 µm at the
1
W = , (24) hottest regions of the dayside to (cid:104)a(cid:105) ∼ 0.1 µm at higher lati-
em (1−η)+η(cid:104)Li,λ(cid:105)/Li,λ tudesandthecoolernightside.Cloudparticlestotheeastofthe
sub-stellarpointaregenerallysmallerthanthosetothewestof
where thesub-stellarpoint,with∼1orderofmagnitudedifferencesbe-
tweencloudparticlesizesattheeastandwestterminatorregions.
(cid:80)
L
(cid:104)Li,λ(cid:105)= (cid:80)i ii,λ, (25) Ctuldoeu;dTpiOart[isc]lericcohmgpraoisnistiodnomalisnoatveatrhieeshwotittehrldaatiytusiddeeaenqdualtoonrgiai-l
i λ 2
regions, while silicate materials (SiO[s], SiO [s], MgSiO [s],
2 3
istheaveragecellluminosityatwavelengthλ.Weapplyafixed Mg SiO [s])dominatethecoolerhigherlatitudesandnightside
2 4
η=0.5throughoutoursimulations.Baesetal.(2016)providean regions.TiO [s]richgrainsarealsomoreprevalentontheeast-
2
indepthdescriptionofthiscompositeemissionbiasingmethod. ern terminator region, while silicates dominate the western re-
Articlenumber,page6of17
G.Leeetal.:DynamicmineralcloudsonHD189733b
gions.Thismaterialcompositiondisparity,duetodifferencesin N
the thermal stability between TiO [s] and silicates, is the main
2
reasonwhyeachterminatorregioncontainsdifferentcloudpar-
ticlesizes.
FortheMCRTpost-processingsimulation,weusethesame
spherical RHD simulation dimensions of radial, longitude and
latitude sizes (R, φ, θ) = (100, 160, 64) respectively. Latitude
θ cells defined near the pole edges in the RHD simulation are
assumed to extend to −π and π radians respectfully. This is re-
W
quired so that no regions remain undefined across the whole -80 Sub-stellar 80 E
sphereintheMCRTsimulation.Cellthermodynamicquantities -60 point 60
(Tgas,ρgas)andcloudmomentproperties(Lj)aretakenfromthe -40 40
snapshot(∼60Earthdays)simulationofPaperI.Thecellther- -20 0 20
modynamic states are kept constant at the values given by the
α [ ]
◦
RHDmodelthroughouttheMCRTsimulations.
3.2. Cloudandgasopacity
S
The input opacities are critical to the accuracy of the MCRT
scheme as they govern the distance travelled by the stochasti-
Fig.2.Diagramvisualisingthedefinitionofequatorialviewingangle
cally sampled τ , the probability of packets scattering, and the
λ α [◦] used in this study. Viewing angles east of the sub-stellar point
type of scattering events. The total extinction opacity κ
ext,total are positively defined, while regions west of the sub-stellar point are
[cm2g−1]ineachcellisgivenby
negativelydefined.
κ =κ +κ +κ , (26)
ext,total abs,gas sca,H2 ext,cloud
ratio of the third (L ) and second (L ) cloud moments (Hansen
where κ [cm2 g−1] is the absorption due to gaseous atoms 3 2
abs,gas &Travis1974)
andmolecules,κ [cm2 g−1]isthecontributionofRayleigh
scatteringbymolsecac,uHl2arH ,andκ =κ +κ is (cid:114)
tchleest.otalextinction(absor2ption+esxct,calottuedring)abds,ucleoutdocloscuadt,clpoaudrti- aeff = 3 43πLL23, (29)
Forgasabsorptionopacitiesκ ,weadoptthesolarmetal-
abs,gas
licity,chemicalequilibriumSharp&Burrows(2007)opacityta- whichistheareaweightedaveragecloudparticlesizeacrossthe
bles,withoutTiOandViOopacities,linearlyinterpolatedtothe grain size distribution. The single scattering albedo ωcloud (Eq.
RHDcellthermodynamicproperties.Weperformk-distribution 8)andscatteringasymmetryparametergλ (Eq.12)arealsocal-
(e.g. Goody et al. 1989; Lacis & Oinas 1991; Grimm & Heng culatedattheeffectivecloudparticlesize.Effective(n,k)coef-
2015) band averages from the high-resolution opacity tables ficients of the cloud particles required as input for Mie theory
across 250 bins, linear in wave-number between (wavelengths) arecalculatedformixedmaterialcloudparticlesusingeffective
0.3and5µm.EightGaussianquadraturepointswereusedwith mediumtheoryinthesamewayasLeeetal.(2015).
asplitat95%ofthecumulativeopacity,inthesamemanneras
Showmanetal.(2009).
3.3. Scatteredlightphasecurvesandalbedospectra
FortheH scatteringopacityweusethecross-section,wave-
2
lengthrelationfromDalgarno&Williams(1962)givenby To produce scattered light phase curves, we track the random
walks of stellar L-packets incident on the HD 189733b simu-
σ (λ)= 8.14×10−13 + 1.28×10−6 + 1.61, (27) latedatmosphere.Thestellarilluminationisassumedtobeplane
sca,H2 λ4 λ6 λ8 parallel from the direction of the host star, with the initial lati-
tudeandlongitude(φ,θ)topofatmospherepositionsdetermined
wophaecrietyλbiystihneÅrealnatdioσnscκas,cHa2,Hi2n(λc)m=2.nTHh2σisscisa,Hc2o(nλv)e/rρtegdasi.nTtoheaambuanss- sctuolcahraasntincuallulysforofmthaerdeajeycstiidoenfmaceethoodf,tchaerrpieladnoeut.tIancrtohsissthmeocdier-,
dance of molecular H2 (nH2 [cm−3]) is calculated at each cell shouldthepacketbeabsorbedbythegasorcloudparticles,itis
assumingchemicalequilibriumusingtheroutinesfromHelling terminatedandnolongercontributestothesimulation.Thesur-
etal.(2016a). vivalbiasingvariancereductiontechniqueisnotrequiredinthis
For cloud particle absorption and scattering opacities we modesincethemajorityofincidentstellarpacketsarescattered
calculate a mono-disperse absorption κabs,cloud and scattering orabsorbedintheatmosphere.Therefore,thefractionalweight
κsca,cloudopacityofthecloudparticlesineachcellgivenby ofthepacketluminosityremainsunchangedinthismode(W=
1).Wechoseaviewingangleofθ=0◦ latitude(equator)while
κabs,cloud(λ)=πa2Qabs(λ,a)nd/ρgas, alteringthelongitudeviewingangleαtocaptureplanetaryphase
κ (λ)=πa2Q (λ,a)n /ρ , (28) information. We define α = 0◦ as the viewing angle at the sub-
sca,cloud sca d gas stellar point (i.e.. full phase) with +◦α eastward and −◦α west-
respectively, where a [cm] is cloud particle size, Q and Q wardfromthesub-stellarpoint.Figure2showsaschematicvi-
abs sca
aretheabsorptionandscatteringefficiencyfactorsgivenbythe sualisationofthisαdefinition.Asaconsequence,weignorepos-
MieX code of Wolf & Voshchinnikov (2004), and n [cm−3] is sibleviewingeffectsofnon-negligibletransitimpactparameters,
d
the total cloud particle number density. We calculate the cloud similartootherpost-processingRHD/GCMmethods(e.g.Fort-
opacities at the effective grain size, aeff, which is given by the neyetal.2006;Showmanetal.2009;Amundsenetal.2016).
Articlenumber,page7of17
A&Aproofs:manuscriptno.RHD+MCRT
Thetotalfractionofluminosity,carriedbythephotonpack- In this scheme the packets are not terminated upon absorp-
ets,escapingtowardsviewingangleαatwavelengthλisgiven tion but are forced to scatter with reduction in their weight at
by summing the image pixel array produced by the peeling off each event. Survival biasing and Russian Roulette techniques
method(Sect.2.4.1,Fig.1) are therefore applied (Sect. 2.4.2). Composite emission bias-
(cid:88) ingisalsoappliedtodecreasethenoiseerrorofresultsderived
f (λ,α)= f (λ,α), (30) from nightside regions (Sect. 2.4.3). We employ a wavelength-
total x,y
x,y dependent dark zone (Pinte et al. 2006) of longitude- and
latitude-dependent, radially integrated optical depth τ (θ,φ) =
where f is the total fraction of luminosity escaping towards r,λ
total 30(Eq.3).Nopacketsareemittedfromtheseregionsandpack-
viewing angle α at wavelength λ and f is the fraction of lu-
x,y etsareterminatediftheycrossintoregionsdeeperintheatmo-
minositycontainedina(x,y)pixelimagearray.Themonochro-
sphere than the cell radial index given by this condition. This
maticphasedependentalbedoA (α)isdefinedas
λ preventstrackingpacketsthroughlowluminosityweightregions
f (λ,α)L (λ) which contribute negligibly to the photospheric emission. The
Aλ(α)= totalL (λ)inc = ftotal(λ,α), (31) contribution of cells below this depth to the total cell emission
inc luminosityisalsoignored(i.e.L =0).Theemittedmonochro-
i,λ
where Linc [erg s−1] is the incident luminosity onto the day- maticluminosityLp,em,λ [ergs−1 cm−1]fromtheplanetatview-
side face of the exoplanet atmosphere from the host star. The ingangleαisthengivenby
monochromatic apparent geometric albedo A is defined as
g,λ
the fraction of scattered light at zero phase angle (α = 0◦) to
(cid:88)
an equivalent spherical Lambertain surface (e.g.. Seager 2010; L (α)= f (λ,α) L . (37)
p,em,λ total i,λ
Madhusudhan & Burrows 2012). In the MCRT scheme A is
g,λ i
derivedfromtheluminosityfractionsbytheequation
f (λ,α=0◦)L (λ) f (λ,α=0◦) As in the scattered light case, we assume the same number of
Ag,λ = totalf L (λi)nc = total 2/3 , (32) packets (Natm = 108) for each (pseudo) monochromatic wave-
Lambert inc length.
where f = 2/3 is the scattering fraction of the theoretical
Lambert
Lambertianplanetatzerophase.
The wavelength-dependent classical phase function (e.g.
4. Results
Seager 2010; Madhusudhan & Burrows 2012) can be con-
structedfromthescatteringfractions In this section, we present the scattered and emitted light re-
sults of our MCRT post-processing method. We chose 36 dif-
f (λ,α)
Φ (α)= total , (33) ferentviewinganglestocaptureorbitalphaseinformationofthe
λ ftotal(λ,α=0◦) scattered and emitted light, which gives a phase resolution of
∆α = 10◦. Sections 4.1 and 4.2 present our incident scattered
whichistheratioofreflectedluminosityatviewingangleαcom-
pared to the fraction reflected at α = 0◦. The monochromatic light apparent geometric albedo Ag, and albedo spectra/phase
curves,respectively.Section4.3presentsouratmosphericemit-
planetary luminosity due to reflected starlight L [erg s−1
p,scat,λ tedluminosityspectraandphasecurves.InSect.4.4wecombine
cm−1]asafunctionofviewingangleisgivenby
our scattered and emitted light results and compare these com-
Lp,scat,λ(α)= ftotal(λ,α)L(cid:63),λ4Ra2p2, (34) bTiEnSedSraensdulCtsHwEiOthPcSuprrheonttoombesterricvaotbiosnesrv;awtieonaslsfoorprHedDic1t8K97ep3l3ebr,.
Finally, in Sect. 4.5 we discuss the variance and convergence
propertiesofourMCRTmodelresults.
where R is the radius of the planet and a the semi-major axis.
p
Themonochromaticluminosityofthestarisgivenby
4.1. Geometricalbedo
L(cid:63),λ =4π2R2(cid:63)Bλ(T(cid:63),eff,λ), (35)
Figure 3 shows the resultant geometric albedo of our post-
where R is the radius of the star and B the Planck function,
(cid:63) λ processed cloud forming HD 189733b 3D RHD model, which
whichisdependentonthestellareffectivetemperatureT(cid:63),eff.To wecomparetotheHSTSTISobservationsofEvansetal.(2013).
ensuregoodsignaltonoise,weemitthesamenumberofpackets
OurresultsintheBbandareconsistentwiththeindividualSTIS
(N =107)foreach(pseudo)monochromaticwavelength.
inc data points, but underestimate the convolved B Band geomet-
ricalbedo.OurresultsarealsoconsistentwiththeWiktorowicz
3.4. Emittedlightphasecurvesandspectra etal.(2015)groundbasedB+Vbandpolarimetrymeasurement
of an upper limit of A < 0.4. Our V Band result compares
g,λ
WefollowasimilarschemetotheprotoplanetarydiskandISM
poorertotheobservations,howeveritisclearthatNaabsorption
MCRTbasedstudiesofPinteetal.(2006)andRobitaille(2011)
isresponsibleforthedownwardtrendnear0.6µm.Thisisper-
wherethemonochromaticluminosityL [ergs−1cm−1]ofcom-
i,λ hapsbecauseofthepresenceofanunknownabsorbingmaterial
putationalcelliis
intheVband(Sect.5).
L =4πρ V(κ +κ )B (T), (36) Figure 3 also shows the individual contributions from H2
i,λ i,gas i i,λ,abs,cloud i,λ,abs,gas λ i
Rayleigh scattering, cloud particle Rayleigh scattering, and
where V is the volume of the cell, ρ the gas density, cloudparticleTTHGscatteringtothetotalfractions.TheTTHG
i i,gas
κ the cloud absorption opacity, κ the gas ab- scatteringcontributesthegreatestfractiontothetotalgeometric
i,λ,abs,cloud i,λ,abs,gas
sorption opacity, and B the Planck function, which is depen- albedo across all wavelengths. A small H Rayleigh scattering
i,λ 2
dentonthecelltemperatureT. componentisalsopresentatshorteropticalwavelengths.
i
Articlenumber,page8of17
G.Leeetal.:DynamicmineralcloudsonHD189733b
0.6 0.6
Evansetal.(2013) Evansetal.(2013)
0.5 HSTSTIS 0.5 HSTSTIS:B-band
Total HSTSTIS:V-band
H2-Rayleigh Total
0.4 0.4
Cloud-Rayleigh H2-Rayleigh
Cloud-TTHG Cloud-Rayleigh
λ0.3 λ0.3 Cloud-TTHG
,g ,g
A A
0.2 0.2
0.1 0.1
0.0 0.0
0.3 0.4 0.5 0.6 1 2 3 4 5
λ [µm] λ [µm]
Fig.3.ScatteredlightapparentgeometricalbedoA (Eq.32)ofourHD189733bsimulationoutputcomparedtoEvansetal.(2013)’sindividual
g,λ
HSTSTISmeasurements.Thetotalcontributionisindicatedbythesolidblueline;TTHGcloudparticlescatteringbythedashedorangeline;
Rayleighcloudparticlescatteringbythetealdash-dottedline;H Rayleighscatteringbythebrownsolidline.Left:(0.29-0.6µm)TheRayleigh&
2
TTHGmodelisgenerallyconsistentwiththeBBandHSTSTISmeasurements,butdoesnotpredictthelowalbedooftheVband.Right:(0.29-5
µm)Extendedspectrafrom0.3-5µm,showingtheconvolvedBBandandVbandmeasurementsfromEvansetal.(2013).A becomesnegligible
g,λ
beyond5µm.
4.2. Scatteredlightphasecurves phase curves around the sub-stellar point, while longer 1-5µm
can show asymmetry in their scattering properties. Longer IR
During a single orbit of HD 189733b around its parent star, wavelengthscanalsoshowwestwardoffsetsof10-20◦,depend-
different atmospheric regions come in and out of view of the
ingonwavelength.Thisisduetolargercloudparticlesresiding
detector, dependent on the planetary phase being observed. If
on the western hemisphere of the dayside (Paper I), increasing
eastwardandwestwardhemispheresofthedaysideatmospheric
theopacityandscatteringprobabilityofIRpacketstravellingin
properties are different, an asymmetric signal as a function of
theseregions.
phase (α) is expected to be observed. In order to extract this
We compare our classical phase function to the Lambertian
phasebehaviourfromour3DRHDresults,weproduceobserv-
phase function (Fig. 5, dashed black line), which is the phase
ables at various viewing angles (longitude from the sub-stellar
function of a theoretical, perfectly isotropic scattering sphere
point)αduringtheMonteCarlosimulation.
(e.g.Seager2010;Madhusudhan&Burrows2012).Thisisuse-
Figure4showsthealbedoAλ(α)spectra(Eq.31)asafunc- ful for interpreting the type of scattering behaviour at east and
tionofphase.Slightdifferencesinthefractionofscatteredlight westdaysidehemispheres(Madhusudhan&Burrows2012).The
between the regions east and west of the sub-stellar point are eastern hemisphere shows a classical phase function behaviour
present,duetothedifferencesintheeast-west3Dcloudproper- typical of Rayleigh-like scattering across all wavelengths, sug-
tiesinourHD189733bRHDmodel.Thegreatestdifferencein gestingthatthemajorityofopticalandnear-IRpacketsundergo
Aλ (∆ Aλ = 0.005) in the B and V bands between the east and a single scattering event in these regions. This region contains
west hemispheres occur between viewing angles of ∼ 80◦and thesmallestcloudparticlesizesintheRHDsimulationandlow-
100◦. This corresponds to the east and west day-night termina- est cloud opacity, increasing the likelihood of the packet es-
tor regions, respectively, which have the largest differences in caping the boundaries without further interaction. The western
cloud properties. For example, the eastern terminator generally hemisphereshowsamixofisotropicandRayleigh-likescatter-
has mean particles sizes (cid:104)a(cid:105) ∼ 0.01 µm, composed of a mix of ing classical phase function behaviour for optical, near-IR and
Si-OmaterialsandTiO2.Whilethewesternterminator(cid:104)a(cid:105)∼0.1 IR wavelengths. This suggests that packets experience multi-
µmwithlargervolumefractionsofMgSiO3[s]andMg2SiO4[s] ple scattering events, which results in packets escaping in an
(Paper I). Wavelengths > 5 µm show little or no scattering be- isotropicmanner.Theseregionscontainsomeofthelargercloud
haviour, as the single scattering albedo for both gas and cloud particlesonthedaysideandanincreasedcloudscatteringprob-
componentsbecomenegligible. ability.
Figure 5 (left) presents the albedo Aλ(α) (Eq. 31) phase Figure5(right)alsoshowsIR(3-5µm)classicalphasefunc-
curvesat25differentwavelengthsintherangeconsideredinthe tion asymmetry that arises from the differences in the opac-
MCRT post-processing. Optical and some near-IR wavelengths itystructurebetweenthetwoday-nightterminatorregions.The
showscatteringfractionsAλ(α)>0.1across−90◦...90◦viewing westward limb (α = −90◦) contains larger cloud particle sizes,
angles.Thissuggeststhatasignificantpercentageoftheplanets opacity,andscatteringprobabilityintheIRcomparedtotheeast-
opticalscatteredluminosityremainsobservable,whileamajor- ernlimb(α=90◦).ThebehaviouroftheIRwavelengthclassi-
ity of the dayside hemisphere remains in view. Scattered light calphasefunctionssuggeststhatpacketsdonotinteractwiththe
fractionsdropbelow10%forλ>3µm. easternterminatorcloudstructuresandpassthroughwithoutin-
Figure 5 (right) shows the classical scattered light phase teraction, are absorbed, or interact with cloud particles on the
function Φ (α) (Eq. 33). This adds emphasis to any asymmet- nightside hemisphere since there is a return to the Lambertian
λ
ric scattering behaviour as a function of phase. Optical and functionatnightsideviewingangles.Forthewestwardlimb,the
near-IR wavelengths 0.3-1µm show very symmetric scattering IR cloud opacity is higher, packets travel shallower into the at-
Articlenumber,page9of17
A&Aproofs:manuscriptno.RHD+MCRT
mospherebeforeinteractingandaremorelikelytobescattered. The total monochromatic luminosity from the planet is the
Thesedifferencesincloudstructuresbetweentheterminatorre- sumofthereflectedlightandemittedlight
gionsresultinaskewingoftheIRclassicalphasefunction.The
L = L +L . (38)
westernsideoftheplanetisthereforetypicallybrighterintheIR p,tot,λ p,scat,λ p,em,λ
by10-20%inscatteredlightthantheeasternsideatcomparable The apparent geometric albedo (which includes scattered and
phases. emittedlight)isthen
L (α=0◦)
A = p,tot,λ , (39)
4.3. Emittedlightspectraandphasecurves g,λ L 2/3
inc,λ
Atsecondarytransit,observationsarealsosensitivetothether- in which the numerator and denominator are integrated across
malemissionfromtheplanetaryatmosphereitself.Becausethe the Kepler, TESS, and CHEOPS bandpasses. For our compar-
star-planetcontrastislargeratIRwavelengthsthanintheopti- isonstoobservationaldataandobservationalpredictions,instead
cal,theplanet’satmosphericemissionfeaturesareeasiertode- of assuming the host star radiates as a black body in previous
tect. sections, the stellar monochromatic luminosity of HD 189733
Figure 6 shows the spectral energy distribution (SED) λ (L ) is taken from the Kurucz 1 stellar atmosphere model for
(cid:63),λ
L [erg s−1] (Eq. 37) as a function of wavelength for the HD189733.
p,emλ
atmospheric emitted light. We find that the SED follows simi- Figure 8 shows the spectral luminosity of the planet as a
lar trends across all viewing angles with ∼ 1−2 magnitude dif- functionofphaseincludingboththescatteringandemittedcom-
ferences between dayside and nightside emission luminosities. ponents.Onthedaysideoftheplanet,theluminosityfromscat-
The peak of emission occurs at ∼ 2 µm, where our model HD teredlightisatasimilarmagnitudetotheIRemission.Fornight-
189733b is brightest in emitted light. The least emission lumi- sideprofiles,thescatteredlightfractiondropsoffandtheemitted
nosityintheIRoccursataviewingangleofα=−140◦,corre- luminositydominates.
spondingtothecoldestregionsofthenightside. Figure 9 presents the flux ratio L /L of our com-
p,tot,λ (cid:63),λ
Figure 7 (left) presents the planetary atmosphere SED λ bined scattered and emitted light, compared to the HST sec-
L (Eq. 37) as a function of viewing angle α. We find that ondaryeclipsedatafromSwainetal.(2009);Evansetal.(2013);
p,em,λ
emissionisdominatedbythelongerIRwavelengths2-5µmfor Crouzetetal.(2014);Barstowetal.(2014)andSpitzerdatafrom
all phases. Emission at optical wavelengths is confined to the Knutsonetal.(2012).SimilartoSect.4.1,oursimulationiscon-
hottest dayside regions of the planet, however, this emission is sistentwiththeBBandHSTSTISdatafromEvansetal.(2013)
> 4 magnitudes in luminosity smaller than IR wavelengths. A but overestimates the planet-star flux ratio for the V Band. The
dip in luminosity is seen at ∼ −135◦for all wavelengths, corre- resultsareconsistentwiththeHSTWFC3(Crouzetetal.2014)
sponding to the coldest nightside regions of the RHD simula- and NICMOS (Swain et al. 2009; Barstow et al. 2014) spectral
tion. Figure 7 (right) shows the normalised emitted light phase trends,withoffsetsintheamplitudeofthesefeatures.Severalre-
curves, which emphasises the shape of the phase curve. Most trievalmodelsonthisandsimilarobservationaldata(e.g.Mad-
wavelengths show eastern offsets from the sub-stellar point in husudhan&Seager2009;Leeetal.2012;Lineetal.2012,2014)
peakemissionat10-25◦withsomeIRwavelengthsshowingoff- suggest a possible sub-solar H2O abundance on the dayside of
sets∼45◦fromthesub-stellarpoint. HD189733b,whichwouldlowertheamplitudeoftheH2Ofea-
We find that the effects of atmospheric scattering by emit- tures in our results. Our results are also consistent with the 3.6
tedpacketsonthephasecurvestobenegligible.Thisisbecause µmand4.5µmSpitzerIRACphotometryobservationsofKnut-
of the low gas phase single scattering albedos (ω < 10−3) at sonetal.(2012).Overall,ourcurrentmodelreproducesthesec-
gas
IRwavelengths,whichlowerstheweightsofthepacketstoneg- ondarytransitopticaltonear-IRobservationaltrendswell.
ligible proportions after one or two scattering interactions. The Figure 10 (left) shows the dayside flux ratio predictions for
majorityofpacketsarethenlikelytobeterminatedbytheRus- Kepler, TESS, and CHEOPS bandpasses. We predict a ∼ 10 %
sian Roulette scheme, and those that survive will contain low difference in peak flux ratio between the CHEOPS and TESS
luminosityweightsinfutureinteractions. bandpasses.ThisisdirectlyduetothesensitivityoftheCHEOPS
bandpasstotheopticalscatteringcomponentofourRHDmod-
elledcloudparticles,whiletheTESSbandpassisunaffectedby
4.4. Totalluminositiesandobservationalpredictions thiscomponentandmoresensitivetothenear-IRthermalemis-
sionoftheHD189733bmodel.
The observed luminosity from an exoplanet atmosphere is a
Figure10(right)presentsthefluxratiophasecurvesforthe
combinationofthereflectedstarlightandtheemissionfromthe
Kepler,TESS,andCHEOPSbandpasses.Ourmodellingresults
planetitself.Thecomponentthatdominatestheobservationwill
showthatweexpectHD189733btohaveazeroorsmallwest-
be wavelength- and phase-dependent owing to the 3D inhomo- wardoffsetofnolessthan-5◦ fromthesub-stellarpointforthe
geneous scattering opacities (from cloud particles) and temper-
Kepler,TESS,andCHEOPSbands.Thissuggeststhatthecloud
atureprofiles(fromatmosphericcirculation)oftheplanet.Such particledifferences(size,composition,etc.)attheeastandwest-
inhomogeneityispresentintheresultsoftheHD189733bRHD ernhemispheresinourRHDsimulationarenotradicallydiffer-
cloudformingsimulation.Inthissection,wecombineourscat- entenoughtoproducethelarger(<−10◦)westwardoffsetsseen
teringandthermalemissionresultsabovetoproduceobservable
for some Kepler planets. The TESS and CHEOPS photometry
quantitiesintheTESS(Rickeretal.2014)andCHEOPS(Broeg
become comparable at greater longitudes as the cloud particle
et al. 2013) bandpasses for our HD 189733b model. We also
scattering component becomes less dominant compared to the
comparetoreadilyavailableobservationaldatafromcurrent/past
thermalemission(Fig.8).
missions.WealsoproduceKeplerbandpredictions,despiteK2
If the extra absorption component from the Evans et al.
notscheduledtoobserveHD189733btoallowacomparisonto
(2013) HST STIS measurements at ∼ 0.4-0.5 µm is taken into
otherstudiesfocusedonmodellingKeplerobjects(e.g.Parmen-
tieretal.2016). 1 http://kurucz.harvard.edu/stars/hd189733/
Articlenumber,page10of17