Claudio Quarti
email: claudio.quarti@umons.ac.be

Narrow and broadband light emission in layered organic lead halide perovskites: interplay between weak electron-lattice interactions and defect-related effects

Keywords: 2-dimensional hybrid lead halide perovskites, light-emission, electron-phonon coupling, point-defects

2-Dimensional metal Halide Perovskites (2D-HP) are at the limelight for their potential exploitation in light-emission related applications. In particular, the most-investigated <001>-terminated 2D HP family shows dominant narrow light emission, with reduced Stokes shift, of great interest for display applications. In parallel, these systems often show additional largely Stokes-shifted emission, with reduced spectral resolution, of interest for the lighting application, e.g. development of white light-emitting diodes. Clarifying the emission mechanisms in 2D-HP and explaining the coexistence of these two contrasting emission regimes is greatly coveted, for the further exploitation of this class of semiconductors. Here, Density Functional Theory (DFT) simulations estimate total electron-phonon interaction in 2D-HP in the order of few tens of meV, in agreement with ps-resolved UV-vis measurements, consistent with the reported narrow emission. On the other hand, such small coupling significantly contrasts with the assignment of broadband emission to some form of intrinsic (defect-free) self-trapped exciton. Additional DFT simulations rather assign broadband emission to extrinsic self-trapping, associated to point defects, halide interstitials in particular. This result is in line with similar findings for the parental 3D halide perovskite and highlights, on the one hand, the role of halide defects in lead-halide perovskitoids frames, and explains, on the other hand, the apparent contrasting nature of narrow and broadband emission, as due to two different emission mechanisms. Defect engineering protocol are therefore suggested, to optimize broadband emission in 2D-HP materials.

INTRODUCTION

With first reports tracking back to late nineteenth century, 1,2 layered 2-Dimensional metal Halide Perovskites (2D-HP) are currently facing a second golden era, due to their promising exploitation in opto-electronics. This class of materials has attracted renewed attention following the recent success of their 3-Dimensional (3D) lead halide perovskites analogues in the field of photovoltaic, characterised by an unprecedented and continuous improvement in the corresponding light-to-current conversion efficiency, recently reaching 25.2%. 3 In this frame, 2D-HP are currently under intense investigations for various applications, not solely related to photovoltaics, [START_REF] Tsai | High-efficiency two-dimensional Ruddlesden-Popper perovskite solar cells[END_REF] but also displaying, [START_REF] Dou | Atomically thin two-dimensional organic-inorganic hybrid perovskites[END_REF][START_REF] Tsai | Stable Light-Emitting Diodes Using Phase-Pure Ruddlesden-Popper Layered Perovskites[END_REF] lasing, [START_REF] Leyden | Amplified spontaneous emission in phenylethylammonium methylammonium lead iodide quasi-2D perovskites[END_REF][START_REF] Liang | Lasing from Mechanically Exfoliated 2D Homologous Ruddlesden-Popper Perovskite Engineered by Inorganic Layer Thickness[END_REF] sensors, [START_REF] Deng | Aligned Single-Crystalline Perovskite Microwire Arrays for High-Performance Flexible Image Sensors with Long-Term Stability[END_REF] etc. 3D perovskites possess a AMX 3 composition and a specific crystalline structure, corresponding to an extended network of MX 6 octahedral connected in corner-shared fashion, M being a metal, X a halide and A a small organic cation or alkali atom, filling the spaces between octahedra. 2D perovskites have composition A' 2 A n-1 M n X 3n+1 and can be ideally obtained by separating extended corner-shared XM 6 planes, as cutting the parental 3D perovskite structure along specific crystallographic directions, as shown in Figure 1a, by partially or fully substituting the A alkali atom or small cation with bulky A' cation (common examples are phenylethylammonium or alkylammonium). [START_REF] Smith | Tuning the Luminescence of Layered Halide Perovskites[END_REF][START_REF] Katan | Quantum and Dielectric Confinement Effects in Lower-Dimensional Hybrid Perovskite Semiconductors[END_REF] Altogether, the spatial confinement of the semiconducting MX 6 structure down to a length-scale close to the De-Broglie wavelength imparts these systems with unique opto-electronic features. [START_REF] Katan | Quantum and Dielectric Confinement Effects in Lower-Dimensional Hybrid Perovskite Semiconductors[END_REF] This, along with their improved stability to moisture, [START_REF] Tsai | High-efficiency two-dimensional Ruddlesden-Popper perovskite solar cells[END_REF] as induced by the hydrophobic character of the bulky A' organic spacer, makes 2D-HP extremely appealing for wide range of applications. In particular, spatial and dielectric confinement, the latter arising from the contrast in the dielectric susceptibility of the inorganic layer (e=5-6) with that of the organic (e =2.5-3.5), result in enhanced electron-hole interaction, [START_REF] Katan | Quantum and Dielectric Confinement Effects in Lower-Dimensional Hybrid Perovskite Semiconductors[END_REF] and consequent formation of stable excitons with sizable (350 -450 meV [START_REF] Blancon | Scaling law for excitons in 2D perovskite quantum wells[END_REF] ) and tunable electron-hole binding energy. [START_REF] Smith | Decreasing the electronic confinement in layered perovskites through intercalation[END_REF] Resulting light emission is characterised by considerable photoluminescence quantum-yield, [START_REF] Kawano | Effects of Organic Moieties on Luminescence Properties of Organic-Inorganic Layered Perovskite-Type Compounds[END_REF] narrow emission (full-width-half-maximum between 20 and 40 meV) [START_REF] Smith | Tuning the Luminescence of Layered Halide Perovskites[END_REF] and reduced Stokes shift (few tens of meV) [START_REF] Hong | Photoconductivity and Electroluminescence in Lead Iodide Based Natural Quantum Well Structures[END_REF] compared to absorption, as shown in Figure 1b, which makes 2D-HP ideal for display applications. Narrow emission is generally found for inorganic emitters, [START_REF] Saito | A narrow photoluminescence linewidth of 21 meV at 1.35 µm from strain-reduced InAs quantum dots covered by In 0.2 Ga 0.8 As grown on GaAs substrates[END_REF][START_REF] Kim | Metal halide perovskite light emitters[END_REF] where the electron/exciton phonon interaction is weak [START_REF] Yu | Fundamentals of Semiconductors[END_REF] and defect concentration severely controlled, via the use of expensive high-temperature production processes. The neat exciton emission observed in 2D-HP comes therefore with surprise, considering that these materials are produced from solution-based techniques. In this frame, it is worth to mention that several recent investigation pointed out towards additional broadband emission for 2D-HP materials, characterised by very large Stokes shift (600 meV) and reduced spectral resolution (full width half maximum, FWHM~1 eV), [START_REF] Dohner | Self-Assembly of Broadband White-Light Emitters[END_REF][START_REF] Yangui | Optical Investigation of Broadband White-Light Emission in Self-Assembled Organic-Inorganic Perovskite (C6H11NH3)2PbBr4[END_REF][START_REF] Cortecchia | Polaron self-localization in white-light emitting hybrid perovskites[END_REF][START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] as shown in Figure 1c. Originally observed in <110> terminated 2D-HP, [START_REF] Dohner | Self-Assembly of Broadband White-Light Emitters[END_REF] this broadband emission has been observed later-on also in <001>-terminated 2D-HP, although under specific temperature conditions, typically limited to low (<100 K) temperature, [START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] hence demonstrating this feature is not exclusive of a precise octahedral connection pattern but is a general property of this class of emitters.

In this perspective, the presence of both narrow exciton light emission and broadband emission is surprising, as it points out towards the coexistence of two different regimes in the optical behaviour of these organic-inorganic semiconductors, specifically a weak exciton-lattice-defect interaction regime, associated to narrow emission, along with a strong exciton-lattice-defect interaction regime, associated to broadband emission. Overall, the deep understanding of the optical behaviour of layered 2D-HP requires a theoretical frame which takes into account not only the electronic structure and exciton features of these systems but which considers also for the perturbations induced by the lattice dynamics (phonons) and defects, as well as the back-perturbation of the excited states on the lattice and defects themselves. Along the direction for the achievement of such ambitious task, I will show here some recent steps forward in the general understanding of the exciton properties of layered halide perovskites. Here, I will discuss the abovementioned two regimes of 2D-HP. The emerging picture is that 2D HP are generally weakly perturbed to phonons, unlike organic semiconductors, [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] or recently synthesised halide double perovkites. [START_REF] Steele | Giant Electron-Phonon Coupling and Deep Conduction Band Resonance in Metal Halide Double Perovskite[END_REF][START_REF] Zelewski | Revealing the nature of photoluminescence emission in the metal-halide double perovskite Cs2AgBiBr6[END_REF] However, 2D HP can be sensitive to lattice defects, mainly halide interstitials, that induce deep trap states and that can be at the origin of complex photochemical processes, [START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] leading to broadband emission. 1b reproduced with permission from Ref. [START_REF] Dohner | Self-Assembly of Broadband White-Light Emitters[END_REF] © 2014 American Chemical Society. Figure 1b reproduced with permission from Ref. [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] © 2017 American Chemical Society.

RESULTS AND DISCUSSIONS

Narrow emission in 2D HP and weak electron-phonon coupling

Narrow emission with reduced Stokes shift is associated to materials with reduced electron-/exciton-phonon coupling, along with low defect density, as typically the case of conventional inorganics semiconductors. [START_REF] Yu | Fundamentals of Semiconductors[END_REF] In this sense, recent efforts have been directed in clarifying the nature and the strength of the electron-and exciton-phonon coupling in 2D-HP. [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF][START_REF] Straus | Direct Observation of Electron-Phonon Coupling and Slow Vibrational Relaxation in Organic-Inorganic Hybrid Perovskites[END_REF][START_REF] Ni | Real-Time Observation of Exciton-Phonon Coupling Dynamics in Self-Assembled Hybrid Perovskite Quantum Wells[END_REF] In particular, Ref. [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] provides an elegant demonstration of such interaction, along with an estimate of its strength on the order of few tens of meV, using joint experimental/computational approach. Pico-second resolved differential UV-vis transmission measurements have been performed at 5 K for two of the most studied 2D-HP, namely buthylammonium (NBT) and phenylethylammonium (PEA) lead iodide perovskites. In the upper panel of Figure 2a, the steady-state UV-vis absorption of PEA 2 PbI 4 2D-HP is showed close to the exciton resonance. Corresponding differential transmission spectrum, as measured after pumping at 3.06 eV, is also reported in the Figure 2a, showing the typical negative feature across the A band and a positive feature at the B band. This line-shape is characteristic of a superposition of excitation-induced shift and broadening of the exciton transition by many-body interactions along with the bleach of higher lying electronic states. In the following pico-seconds, we observe the usual exponential decay of the optical signals, as related to the de-excitation of the system, but also a periodic modulation of the signal, which we identify as the signatures of coherent phonons generated via Resonant Impulsive Stimulated Raman Scattering (RISRS) induced by the ultrashort pump pulse. As matter of fact, the present measurement is the time-domain analogue of a Resonant Raman measurement, probing the phonons which most strongly couple with the given electronic transition. Resonant Raman scattering is often interpreted on the basis of popular displaced harmonic oscillators model in Figure 2b. [START_REF] Coropceanu | Charge Transport in Organic Semiconductors[END_REF] The potential energy surface of the ground and excited state are treated as two harmonic oscillators with same curvature along the coordinate Q (typically corresponding to the normal modes, in molecules) and displaced by a quantity D. According to the Franck-Condon principle, the ground-to-excited state takes place at the geometry of equilibrium of the ground state, which in the most general case does not correspond to that of the excited state (D¹0), followed by relaxation towards the equilibrium of the excited state, in the fs-ps time-domain. Corresponding relaxation energy l rel depends on the strength of the coupling along the coordinate Q and on the displacement D of the two oscillators:

! "#$ = &' ( (* + ) &* - 4 & -' /0/ (* + ) &* - (Eq.1)
Eq. 1 requires to evaluate the derivative of the ground to excited transition energy (E g ) with respect to the normal mode Q, as estimated at the ground state equilibrium structure Q 0 (enumerator) and the curvature of the potential energy surface of the ground/excited state with respect to the normal mode (denominator). In the case of systems having inversion symmetry, selection rules predict zero Raman activity (either resonant or non-resonant) for those normal modes with ungerade (odd) character, which reflects the fact that along these vibrations the two oscillators are not displaced (D=0).

Consistently with the above-mentioned theoretical frame, we performed periodic Density Functional Theory calculations, within the Linear Combination of Atomic Orbitals formalism, as implemented in the CRYSTAL17 package, [START_REF] Dovesi | Quantum-mechanical condensed matter simulations with CRYSTAL[END_REF] and we computed the quantities needed in Eq. 1, namely, normal modes of vibration at the electronic ground state, corresponding curvatures, and the derivative of the ground-to-excited state electronic transition, with respect to the displacement patter for each normal mode. In the present case, the ground-to-excited state electronic transition was roughly approximated with the lowest energy band-to-band transition. Our theoretical predictions for NBT 2 PbI 4 2D-HP in Figure 2c nicely parallel the measured RISRS signals from differential transmission spectrum, as probed at an energy close to the band-to-band transition, at 2.67 eV, supporting the general theoretical frame. We found two sets of phonons which couple to the band-to-band electronic transition, around 3-4 meV (~25-32 cm -1 ), and around 5 meV (~40 cm -1 ), respectively. Displacement patter associated to the most interacting normal mode at 3 meV is showed in Figure 2d. Computed frequencies for the low energy phonons in Figure 2c are slightly underestimates, compared to experiment, but it is worth to stress that the accurate prediction of normal modes at such low frequency range is inherently complicated by the shallow form of the potential energy surface and by anharmonic effects. As the main interest in the optical properties of 2D HP is related to their excitonic properties, the RISRS signal was also measured at probe energy corresponding to the A and B signals in Figure 2a. In particular, tuning the RISRS signal for probe at 2.36 eV, corresponding to the lowest energy excitonic A state, the highest phononic component at 5 meV was shown to be dominant, while the component at 3-4 meV was much less important. [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] From a theoretical and computational viewpoint, estimating the deformation-of-potential electron-phonon interaction via Eq. 1 can be conceptually performed, inherently considering electron-hole interaction, via for instance Time Dependent-Density Functional Theory (TD-DFT) or the solution of the Bethe-Salpeter Equation (BSE). However, TDDFT calculations were shown to not accurately account for the electron-hole screening, with final exciton binding energy being highly sensitive with respect to the amount of exact exchange used in the definition of the exchange-correlation kernel. [START_REF] Quarti | Tuning the Optoelectronic Properties of Two-Dimensional Hybrid Perovskite Semiconductors with Alkyl Chain Spacers[END_REF] In contrast, ab-initio solution of the Bethe-Salpeter Equations (BSE) performed on electronic structure from GW approximation was shown to provide exciton binding energy nicely matching the experimental datum, at least for the lowest energy 1s-component of the Rydberg series. [START_REF] Giorgi | Nature of the Electronic and Optical Excitations of Ruddlesden-Popper Hybrid Organic-Inorganic Perovskites: The Role of the Many-Body Interactions[END_REF] On the other hand, the computational cost of this methodology makes it unsuitable for similar study as the one reported in Ref. [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] , as the latter requires several BSE calculations, via distorting the crystalline structure along the different normal modes. In this sense, the availability of simplified methods for the simulation of excited states would be greatly needed, to further clarify these problems. [START_REF] Blancon | Scaling law for excitons in 2D perovskite quantum wells[END_REF] Altogether, total relaxation energy can be estimated from the sum of the relaxation energies computed for each phonon, which sums up to ca. 46 meV, that is ca. 2 times kT, in fair agreement with the measured Stokes shift. The general picture coming out is that, although not being completely transparent to lattice dynamics, excitons in 2D HP are just weakly perturbed by phonons. This seems not the case, for instance, for recently reported halide double perovskites. [START_REF] Steele | Giant Electron-Phonon Coupling and Deep Conduction Band Resonance in Metal Halide Double Perovskite[END_REF][START_REF] Zelewski | Revealing the nature of photoluminescence emission in the metal-halide double perovskite Cs2AgBiBr6[END_REF] Cortecchia and co-workers provide slightly contrasting results from DFT calculations, [START_REF] Cortecchia | Polaron self-localization in white-light emitting hybrid perovskites[END_REF] but the nature of the structural perturbation is slightly more arbitrary than in the case of Figure 2d, where the displacement is along convenient set of coordinates (e.g. phonons) and the amount of distortion is in line with curvature of the potential energy surface. In other words, the perturbation of the electronic structure in Ref. [START_REF] Cortecchia | Polaron self-localization in white-light emitting hybrid perovskites[END_REF] can be overestimated, as result of an unrealistic distortion of the lattice.

Broadband emission in 2D HP as due to defect induced electron-lattice interaction

In light of the previous results, presence of broad-band emission in 2D-HP comes with even more surprise, as this openly contrasts with the RISRS signals measured by Thouin et al. in Ref. [START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] and the corresponding weak electron phonon coupling computed from deformation-of-potential. Of course, other electron-phonon coupling mechanisms, as Frölich interactions, can be invoked, which, in light of different selection rules, may not be probed via RISRS. [START_REF] Yu | Fundamentals of Semiconductors[END_REF] On the other hand, it is worth to stress that the assumption of large electron-phonon interaction is in open contrasts with the observed narrow emission, which at least for <001>-terminated 2D HP is ubiquitous with the temperature. In other words, if the lattice were going to sizably rearrange in response to the formation of photogenerated species (either free-charges or electron-hole exciton pairs), one would expect such relaxation to be general and independent from external factors, for instance the temperature, and should therefore result in the suppression of the narrow emission. In this sense, in their review Karunadasa and co-workers, nicely distinguished between intrinsic and extrinsic self-trapping, as responsible for the broadband emission in 2D-HP. [START_REF] Smith | Tuning the Luminescence of Layered Halide Perovskites[END_REF] Intrinsic self-trapping does not requires imperfections in the material structure and is therefore characteristic of the lattice itself, that means, it implicitly assumes sizable electron-phonon coupling, as discussed above. Extrinsic self-trapping instead requires the presence of defects or imperfections in the lattice to kick-off the mechanism of self-trapping for the excitons, hence potentially opening new channels for the exciton decay in 2D-HP (either radiative or non-radiative), alternative to the narrow emission, as schematically depicted in Figure 3a. In this sense, the presence of different types of defects and/or different defect configurations in the material potentially explains the emission broadening, as due to not a single and well defined decay mechanism but to a manifold of decay mechanisms, with slightly different energetics. Aiming to investigate the possible role of point defects on the optical properties of 2D HP, I performed periodic DFT calculations for supercell models of hexylammonium lead iodide perovskite, containing different types of point defects, namely metal/halide vacancies/interstitials. [START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] Resulting atomic projected Density Of State for the non-defective model, as computed using standard PBE exchange correlation functional, is shown in Figure 3b-d. The DFT band gap computed for the non-defective model (Figure 3b) is obviously underestimated, because of the unbalanced error induced by the lack of spin-orbit coupling and inaccurate treatment of electronic many-body effects (as typical of the employed PBE exchange correlation functional). [START_REF] Quarti | Tuning the Optoelectronic Properties of Two-Dimensional Hybrid Perovskite Semiconductors with Alkyl Chain Spacers[END_REF] Still, the Valence/Conduction Band Edge of the material are related to outer-shell p-states of the iodine and to 6s-states of lead, as widely reported in the literature. [START_REF] Quarti | Tuning the Optoelectronic Properties of Two-Dimensional Hybrid Perovskite Semiconductors with Alkyl Chain Spacers[END_REF][START_REF] Even | Electronic model for self-assembled hybrid organic/perovskite semiconductors: Reverse band edge electronic states ordering and spin-orbit coupling[END_REF][START_REF] Umebayashi | Electronic structures of lead iodide based low-dimensional crystals[END_REF] The introduction of charged (+1) halide vacancy in the model (V I +1 ) only weakly perturbs the electronic structure, mainly in correspondence of the conduction band, as shown in Figure 3c, while the introduction of charged (-1) iodine interstitial defect (I I -1 ) interestingly results in the formation of intragap deep trap states, associated to the interstitial halide itself, as shown in Figure 3d. Further tests were performed on the iodine-related Frenkel defect as composed by the vacancy-interstitial iodide pair, to verify that the DFT results are not biased by the used of background charges. Additionally, effect of Spin-Orbit Coupling (SOC) and of Hubbard correction on the band gap were checked, always finding deep trap states associated to the interstitial iodine species. [START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] Interstitial lead defects have been also found to significantly perturb the electronic structure of the system, via the formation of PbI 2 molecular species and consequent introduction of deep trap states into the electronic structure of 2D-HP. [START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] Interestingly, the results in Figure 3 go on pair with similar results for their 3D HP analogues, highlighting the potential role of interstitial halides as electron acceptors. [START_REF] Meggiolaro | Iodine chemistry determines the defect tolerance of lead-halide perovskites[END_REF] Meggiolaro and coworkers in fact computed the defect formation energy for various types of point defects (metal/halide interstitials/vacancies), finding that interstitial iodine defects form similar deep trap states for widely investigated methylammonium lead iodide 3D perovskite, as in Figure 3d. Therefore, these authors proposed possible photochemical-mechanisms related to charge trapping and subband gap light emission in 3D HP, with reference to experimental UV-vis measurements. [START_REF] Meggiolaro | Iodine chemistry determines the defect tolerance of lead-halide perovskites[END_REF][START_REF] Hoke | Reversible photo-induced trap formation in mixed-halide hybrid perovskites for photovoltaics[END_REF] In this frame, more investigations are required to unveil the detailed mechanism behind the broadband emission, which can potentially involve the trapping of the photogenerated holes into the interstitial halide deep trap states or impose more complicated pathways, but this common role of iodine-related point defects in 2D in Ref. [START_REF] Booker | Formation of Long-Lived Color Centers for Broadband Visible Light Emission in Low-Dimensional Layered Perovskites[END_REF] and 3D HP in Ref. [START_REF] Meggiolaro | Iodine chemistry determines the defect tolerance of lead-halide perovskites[END_REF] is clearly of interest to understand the photophysics of halide perovskite systems. Namely, the connection between defects and broadband emission can have important technological consequences, as it paves the way for the development of a "defectengineering" protocol, for the optimization of this type of emission for light illumination applications.

CONCLUSIONS

The optoelectronic response of semiconductors is often the result of a complex interplay between direct/indirect electronic structure, symmetry related properties, and effects related to the perturbation of the lattice associated to the lattice dynamics (i.e. phonons). In the case of solution-based processable materials, further effects are often related to the almost unavoidable presence of structural defects. 2-Dimensional, metal Halide Perovskites (2D-HP) perfectly epitomize this situation, showing, at least for the family of <001>-terminated 2D-HP (see Figure 1a), two regimes of light emissions, one characterized by small Stokes-shifted, narrow emission and the other by large Stokes-shifted, broadband emission. In the present contribution, I try to provide a global perspective about the mechanisms at the atomistic scale dictating these two features in the photophysics of 2D-HP. Narrow emission is consistent with weak exciton/electronphonon coupling in the range of few tens of meV, as obtained from Density Functional Theory (DFT) calculations and Resonant Impulsive Stimulated Raman Scattering (RISRS) measurements. Such weak exciton/electron-phonon coupling regime resembles the typical behavior of conventional inorganic semiconductors and strongly criticizes the assignment of the broadband, Stokes-shifted emission to intrinsic self-trapped excitons. Furthermore, these results inherently call for the effect of some defect or perturbation in the lattice, as starter for the self-trapping mechanism. In this frame, our Density Functional Theory simulations, as performed on defective crystalline supercell models, highlight the presence of deep-trap states associated to point defects, interstitial halides in particular, with energetics consistent with the measured Stokes shift. Further works will be needed to clarify the mechanism behind the broadband emission in 2D HP. Still, the present contribution smoothly connects the (apparently contrasting) narrow emission of 2D-HP with observed broadband emission. Furthermore, it also links the role of halide defects in 2D-HP and in their 3D analogues, as highlighted by Meggiolaro and coworkers. [START_REF] Meggiolaro | Iodine chemistry determines the defect tolerance of lead-halide perovskites[END_REF] The proposed connection of broadband emission of point-defects can be at the basis of the optimization of broadband emission via "defect-engineering" protocol, aimed to stabilize the defects with the temperature and minimize non-radiative recombination.

Figure 1 .

 1 Figure 1. a) derivation of 2D layered metal halide perovskites from the parental 3D perovskite structure; b) Absorbance (blue) and photoluminescence (red) spectra of <001>-terminated methylpropane-1,3-diammonium lead bromine perovskite; c) Temperature resolved PL of (C 12 H 25 NH 3 ) 2 PbI 4 perovskite.Figure1breproduced with permission from Ref.[START_REF] Dohner | Self-Assembly of Broadband White-Light Emitters[END_REF] © 2014 American Chemical Society. Figure1breproduced with permission from Ref.[START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] © 2017 American Chemical Society.

Figure

  Figure 1. a) derivation of 2D layered metal halide perovskites from the parental 3D perovskite structure; b) Absorbance (blue) and photoluminescence (red) spectra of <001>-terminated methylpropane-1,3-diammonium lead bromine perovskite; c) Temperature resolved PL of (C 12 H 25 NH 3 ) 2 PbI 4 perovskite.Figure1breproduced with permission from Ref.[START_REF] Dohner | Self-Assembly of Broadband White-Light Emitters[END_REF] © 2014 American Chemical Society. Figure1breproduced with permission from Ref.[START_REF] Thouin | Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites[END_REF] © 2017 American Chemical Society.

Figure 2 .

 2 Figure 2. a) steady-state (upper panel) UV-vis spectrum of phenylethylammonium lead iodide 2D-HP, measured at 5K (upper panel) and time-resolved differential transmission spectrum, with pump at 3.06 eV; b) double-well model for potential energy surfaces of ground and excited electronic states, with respect to a general deformation Q. D is the displacement between the minima of the two potentials and l rel is the relaxation energy from the vertical transition to the minimum of the excited state potential surface; c) comparison of the relaxation energies and Franck-Condon factors from DFT simulations and Fourier-transform spectrum of the oscillatory component of the TA spectrum of buthylammonium lead iodide, at probe energy 2.67 eV; d) Normal mode displacement patter associated to the signal measured at 3 meV. Figure 2a-c reproduced with permission from Ref. 22 . © 2019 Springer Nature.

Figure 3

 3 Figure 3. a) energy level scheme associated to the narrow emission (ħw emi ) and broad band emission (ħw broad ); b-d) crystalline model and corresponding projected Density Of States (DOS) for 2D (C 6 H 13 NH 3 ) 2 PbI 4 HP without defects (b), containing one iodide vacancy with +1 charge (c) and containing an interstitial iodide with -1 charge (interstitial, circled in red).Figure 3b-d reproduced with permission from Ref. 21 © 2017 Americal Chemical Society.

  Figure 3. a) energy level scheme associated to the narrow emission (ħw emi ) and broad band emission (ħw broad ); b-d) crystalline model and corresponding projected Density Of States (DOS) for 2D (C 6 H 13 NH 3 ) 2 PbI 4 HP without defects (b), containing one iodide vacancy with +1 charge (c) and containing an interstitial iodide with -1 charge (interstitial, circled in red).Figure 3b-d reproduced with permission from Ref. 21 © 2017 Americal Chemical Society.

ACKNOWLEDGEMENTS

The author acknowledges support from Agence Nationale pour la Recherche (MORELESS project).


Sean P Matt
email: sean.matt@cea.fr

Olivier Do Cao

Benjamin P Brown

Allan Sacha Brun

Convection and Differential Rotation Properties of G & K Stars Computed with the ASH Code

Keywords:

Key words convection -hydrodynamics -Stars:

Introduction

The sun and sun-like stars with convection zones in their outer envelopes, have long been known to exhibit emission line and X-ray activity, associated with hot gas in chromospheres, transition regions, and coronae (e.g., Pizzolato et al., 2003;Strassmeier et al., 1990;Wright et al., 2011). Understanding this activity, and understanding how the solar case relates to the activity observed accross the HR diagram is a long-standing puzzle. It is clear that the existence of hot gas above the photosphere is related to magnetic processes associated with the convection zone itself. The magnetic field of low-mass main sequence stars is generally believed to be generated by dynamo processes, which derive their properties from convective motions, differential rotation, and meridional circulations in the convection zones, as well as from the interaction between the convection zone and the radiative interior.

Recent observations, by either spectropolarimetry (e.g., Donati et al., 2003), doppler imaging (e.g., Barnes et al., 2005), or monitoring of various activity indicators (e.g., Baliunas et al., 1995;Donahue et al., 1996;Lovis et al., 2011;Oláh et al., 2009;Saar & Brandenburg, 1999), show that solar-like stars possess activity cycles and differential rotation, analogous to the sun. Solar analogues are even starting to be discovered (Petit et al., 2008).

In order to gain a better theoretical understanding of how convective properties depend upon stellar parameters, we carry out 3D numerical dynamical simulations of convective envelopes of solar-like stars. As a first step, we model here the convective regions only, neglecting the effects of the interface region with a radiative envelope, and restrict ourselves to relatively slow (solar) rotation rates. Specifically, we simulate the convection dynamics for 4 main sequence stars with masses of 0.5, 0.7, 0.9, and 1.1 solar masses, spanning spectral types G0 to K7. This mass interval exhibits a large range of the physical properties of convective envelopes (such as the depth, physical size, mass, and density), as well as in the overall stellar luminosity transported by convection. These types of stars are also targets for asteroseismic studies (e.g., Verner et al., 2011), which have the potential to give precise measurements of stellar properties for large numbers of stars. The goal here is to determine how the convection, differential rotation, and meridional circulation is influenced by stellar mass, and to see if general trends or scaling laws can be extracted that will guide a deeper understanding of the inner hydrodynamics of these stars. The present study lays the groundwork for later studies to consider (e.g.) faster rotation rates, convection-radiation zone interface dynamics, and the dynamo generation of magnetic fields in stars in this mass range.

Section 2 contains a description of our simulation method and presents a comparison of the overall structures of each star. Section 3 describes the main results of our 3D simula-tions, focusing on both the convective turbulence properties, as well as the differential rotation and meridional circulation flows. A summary and brief discussion is contained in section 4.

Simulation method

We use the anelastic spherical harmonic (ASH) code (Clune et al., 1999) to compute the 3-dimensional and turbulent flows in convectively unstable stellar envelopes. This code has been extensively tested and used for computing several aspects of the solar interior (e.g., Browning et al., 2006;Brun et al., 2004;Brun & Toomre, 2002;DeRosa et al., 2002;Miesch et al., 2006), rapidly rotating young stars (Ballot et al., 2007;Brown, 2009;Brown et al., 2008Brown et al., , 2011)), the convective cores of massive stars (Browning et al., 2004a,b;Featherstone et al., 2009), fully convective low mass stars (Browning, 2008), red giant stars (Brun & Palacios, 2009), and pre-main-sequence stars (Bessolaz & Brun, in this volume;Bessolaz & Brun, 2011). We briefly describe the basic aspects of the code here, but the reader can find further details of the code in those previous works (see especially, Brun et al., 2004;Clune et al., 1999).

The code solves the fluid equations, under the anelastic approximation, in a computational domain consisting of a spherical shell and in a rotating reference frame. Under the anelastic approximation, sound waves are filtered out and assumed to have a negligible effect on the dynamics, in order to allow for a larger computational timestep. This approximation is appropriate (e.g.) in the interiors of stars because typical motions are highly sub-sonic. The thermodynamic variables are linearized with respect to a spherically symmetric background state with a density ρ, pressure P , temperature T , and specific entropy S and fluctuations about the background state of ρ, P , T , and S. The timedependent equations describe the conservation of mass, momentum, and entropy expressed as

∇ • (ρv) = 0, (1)
ρ ∂v ∂t + (v • ∇)v + 2Ω o × v = -∇P + ρg -[∇ P -ρg] -∇ • D, (2)
ρ T ∂S ∂t + ρ T v • ∇( S + S) = + ∇ • [κ r ρc p ∇( T + T ) + ρ T (κ∇S + κ 0 ∇ S)] + 2ρν e ij e ij -1/3(∇ • v) 2 , (3)
where v = (v r , v θ , v φ ) is the velocity in the rotating frame in spherical coordinates, Ω 0 = Ω 0 êz is the angular rotation rate of the reference frame, g is the acceleration due to gravity, κ r is the radiative diffusivity, and c p is the specific heat at constant pressure. The term D is the viscous stress tensor, with the components

D ij = -2ρν[e ij -1/3(∇ • v)δ ij ], (4)
where e ij is the strain rate tensor, and ν, κ, and κ 0 are effective eddy diffusivities. The code also uses a linerized equation of state, ρ ρ = P P -

T T = P γ P - S c p , (5)
and the ideal gas law,

P = Rρ T , ( 6
)
where γ is the ratio of specific heats (we use γ = 5/3), and R is the gas constant. For all 4 stars in our study, the luminosity has reached a value of more than 99.9% of the total stellar luminosity by the base of the convection zone, so we do not need to include any energy generation by nuclear burning within our computational domain.

We set up four different models in which the domain boundaries and stratification coincides with 1D models of the convection zones of main sequence stars with different masses. Section 2.1 describes the global properties of these stars, and sections 2.2 and 2.3 describe the initial and boundary conditions used in our ASH models, as well as the method for evolving the simulations to a fully-convective, statistical steady-state.

1D stellar structure

In order to define the background structure for our 3D models, we use the 1D stellar evolution code CESAM (Morel, 1997). With CESAM, we computed the evolution of four stars with masses of 0.5, 0.7, 0.9, and 1.1 M until the age of 4.6 Gyr. This age is approximately equal to that of the sun, so our cases can be compared with the many previous results of solar studies. Also, this age is appropriate for the slow (solar) rotation rates considered here. For all four stars, we assumed the same initial metallicity of (X, Y, Z) ≈ (0.71, 0.27, 0.02) and a mixing length parameter of 1.77, chosen to best represent the solar case. Table 1 lists the key global properties of our four stars, at the age of 4.6 Gyr, which are also graphically represented in Figures 123.

Figure 1 shows the total luminosity (solid line), as well as the effective temperature (dashed line) as a function of mass, for stars in the mass range of our models. The "X" symbols mark the values for the 4 stars modeled here, which are also listed in table 1. Also indicated on the plot are the approximate spectral types, corresponding to the effective temperature. The Figure demonstrates the steep rise in lumonosity with mass expected for main sequence stars (the stars here approximately follow L * ∝ M 4.6 * ), such that the 1.1 M star is 40 times more luminous than the 0.5 M star.

Figure 2 shows the photospheric radius R * (upper solid line) and the radial location of the base of the convection zone R cz (lower solid line) as a function of mass. The stellar radius is strong function of stellar mass in the plotted range, with the 1.1 M star being 2.8 times bigger than the 0.5 M star. The radial extent of the convection zones for these stars is represented by the region between the two Table 1 Global properties of the 4 stars used in our ASH models, computed with the CESAM stellar evolution code and at an age of 4.6 Gyr. We adopt M = 1.989 × 10 33 g, R = 6.9599 × 10 10 cm, and L = 3.846 × 10 33 erg s -1 .

Mass Radius

L * T eff SpT Mcz Rcz T (Rcz) ρ(Rcz) (M ) (R ) (L ) (K) (M , M * ) (R , R * ) (K) (g cm -3 ) 0.
T eff SpType = K7 K4/K5 G8 G0
Fig. 1 Stellar luminosity in solar units (solid line, left scale) and effective temperature (dashed line, right scale) as a function of stellar mass, for the mass range spanned by our models, and as computed using the CESAM stellar evolution code. The "X" symbols indicate the values for the 4 stars used in our 3D ASH simulations. The approximate spectral types are also indicated next to each of these.

solid lines, and the region below contains a convectively stable radiation zone. It is clear that, in the mass range shown, the thickness of the convection zone increases slightly with increasing stellar mass. However, since the stellar radius increases more rapidly, the fractional size of the convection zone decreases with increasing mass. Thus, the convection zone thickness ranges from 44% to 25% of the stellar radius, for the 0.5 and 1.1 M stars, respectively (for the sun, this value is approximately 30%). The "X" symbols indicate the upper and lower boundaries of our 3D simulation domains for the 4 stars simulated (discussed below).

Figure 3 shows the extent of convection zones, expressed in mass coordinates and normalized to the mass of each star, as a function of stellar mass. The lines show the total mass enclosed by the stellar surface (upper line) and enclosed by the location of the base of the convective envelope (lower line). The amount of mass contained in the convection zone, M cz , is the difference between these two lines and is listed in Table 1. It is clear that the convection zone mass is a

0.4 0.6 0.8 1.0 1.2 M * / M sun 0.0 0.2 0.4 0.6 0.8 1.0 1.2 radius (R sun ) R cz R * Radiation Zone C o n v e c t i o n Z o n e
Fig. 2 Stellar radius (upper line) and radial location of the base of the convection zone (lower line) in solar units and as a function of stellar mass, as computed by the CESAM stellar evolution code. The convection zone exists between these two lines, while a convectively stable radiation zone lies below. The "X" symbols indicate the upper and lower radial boundaries of the computational domain used in our ASH simulations, for the 4 cases considered.

strong function of stellar mass in the range considered. Stars with slightly lower mass will be fully convective, while significantly more massive stars will not have convective envelopes. The convection zone mass varies from 36% to 1% of the stellar mass, for the 0.5 and 1.1 M stars, respectively (the solar value is approximately 3%). The "X" symbols indicate the simulation domain boundaries, expressed in mass coordinates, for the 4 stars simulated (discussed below).

The last two columns of table 1 show the temperature and the mass density at the base of the convection zone for each case. Although it is true that the temperature and density increases with mass at the stellar center (R = 0), these quantities decrease with increasing mass at the location of R cz . In other words, for increasing stellar mass, the convection zone generally occupies a more tenuous and cooler outer layer of the star.

Table 2 Simulation parameters for each case. The nuber of radial, latitudinal, and longitudinal gridpoints are N r , N θ , and N φ . The outermost radius of the simulated spherical shell is R out , and H ρ (R out ) is the density scale height there. The radial size of the domain is L = R out -R cz . The viscosity at R out is ν top and it varies in the domain with the inverse square root of the background density. All stars have a Prandtl number P r ≡ ν/κ = 0.25 and rotate at the solar rate, Ω 0 = 2.6 × 10 -6 rad/s, corresponding to a rotation period of 2π/Ω 0 ≈ 28 days. Also listed are the Rayleigh number R a ≡ (-∂ρ/∂S)∆SgL 3 /(ρνκ), the Taylor number T a ≡ 4Ω 2 L 4 /ν 2 , and the convective Rossby number R oc ≡ (R a /T a P r ) 1/2 , all evaluated at the midlevel of the domain. Radiation Zone

Mass

Nr, N θ , N φ Rout Hρ(Rout) L νtop Ra Ta Roc (M ) (R , R * ) (Mm) (Mm) (cm 2 s -1 ) 0.

C o n v e c ti o n Z o n e

Fig. 3 Mass contained beneath the stellar surface (upper line) and location of the base of the convection zone (lower line) as a function of stellar mass, as computed by the CE-SAM stellar evolution code. The convection zone exists between these two lines, while a convectively stable radiation zone lies below. The "X" symbols indicate the upper and lower radial boundaries of the computational domain used in our ASH simulations, for the 4 cases considered.

Initial and boundary conditions for 3D simulations

We use the 1D stellar structure models described above for the spherically symmetric, initial conditions of the the 3D ASH models. Table 2 lists the key properties of each of the 4 simulations. In the ASH models presented here, the domain consists of the stellar convection zone only. Thus, the lower radial boundary of our simulated spherical shells coincides with the base of the convection zone, R cz , for each star (see table 1 and the lower row of "X" symbols in Figs. 2 and3). We chose the location of the outer boundary such that the mass density there is a factor of 100 smaller than at the base of the convection zone. This places the outer boundary somewhat below the photosphere of the star. The location of the outer boundary R out is listed in table 2 and shown as the upper row of "X" symbols in and Figures 2 and3.

In each of our 4 cases, the gravity (g) is taken directly from the corresponding CESAM model, and the entropy gradient (∂ S/∂r) is initialized with a constant value equal the (superadiabatic) entropy gradient in the middle of the convection zone in the CESAM model. The background thermodynamic variables (ρ, P , T , and S) are set according to this entropy gradient and to be in hydrostatic balance with the stellar gravity, while the fluctuating thermodynamic variables (ρ, P , T , and S) are initially zero. The initial entropy gradient and choice of diffusivities (see below) ensures that the initial state is convectively unstable and has a Rayleigh number near or above the critical Rayleigh number necessary for convection. The velocity in the rotating frame (v) is given an initial random perturbation, with negligible kinetic energy, in order to initiate convective motions.

For both the inner and outer boundaries, we use stressfree and impenetrable boundary conditions on the velocity and hold the entropy gradient fixed at the initial value. The radiative diffusivity κ r is chosen so that the ASH model has the same radiative flux at all radii as the corresponding CE-SAM model. The initial stellar structure and boundary conditions ensure that the energy flux into the domain at the bottom of the convection zone (given entirely by the radiative flux) is constant in time. At the top of the domain, the density scale height becomes small (see table 2), and since this approximately determines the convection cell size scale, it becomes numerically challenging to resolve the convective motions there. Furthermore, the impenetrable boundary condition precludes any convective enthalpy flux from escaping the top of the domain. To address both of these issues we introduce a diffusive energy flux (the term proportional to κ 0 in eq. [3]) that is assumed to represent an unresolved and spherically symmetric enthalpy flux carried by small-scale convection near the top of the domain (hereafter "unresolved eddy flux"). In all 4 cases, we chose κ 0 so that this flux is negligible in the bulk of the convection zone, but increases as a smooth function near the outer boundary, such that the flux leaving the domain equals the stellar energy flux and is constant in time. In each case, we chose the viscosity at the outer part of the domain, ν top , in order to ensure that the resulting convective turbulence has a significant Reynolds number. As shown in section 3, a higher stellar luminosity leads to larger convective velocities, and thus a larger viscosity is needed to keep the Reynolds numbers comparable (explaining the trend of ν top apparent in table 2). Within the domain, the viscosity in each case depends only upon radius, and it varries as the inverse square root of the background density. Scaling the viscosity in this way ensures a significant turbulence level at all radii and is used in most previous ASH models in the literature. The thermal diffusivity κ equals 4ν everywhere in the domain, giving a constant Prandtl number of 0.25 for all cases. The rotation rate of the reference frame (Ω 0 ), which equals the rotation rate of the star, is set to the solar rate for all 4 stars.

The Rayleigh numbers (calculated using the steady-state value of the entropy in the simulations), Taylor numbers, and convective Rossby numbers for each case are also given in table 2. The Rayleigh and Taylor numbers, show a strong decrease with increasing stellar mass, which is primarily explained by the trend in the diffusivities. If the diffusivities were constant, these numbers would increase with mass for these stars. The convective Rossby number generally increases with mass, meaning that for the same rotation rate, the coriolis force has less influence on the dynamics.

Reaching a statistical steady-state and energy flux balance

At the start of the simulation, the star is in a quiescent state. Thus, there is no significant enthalpy flux, and the system is not initially in radial energy flux balance. In this situation, the evolution of the system, according to equations ( 1)-( 3), ensures that the entropy gradient evolves toward an energy flux balanced state. Furthermore, since the gas in the computational domain is unstably stratified (∂ S/∂r < 0), with a large Rayleigh number (see table 2), significant convective motions begin rapidly after the start of the simulations. Once vigorous convection begins, an energy flux balance is achieved within a few convective turnover times (months, typically) and maintained for the duration of the simulations. In this state, the net energy transport across the domain is constant at all radii, when averaged over several convective turnover times.

Figure 4 shows the spherically-and time-averaged luminosity, as a function of radius in the whole domain, for the highest and lowest mass stars in this study. These are shown after the simulated stars have evolved for several years from the initial state. The boundary conditions ensure a fixed radiative energy flux into the domain at the inner boundary and a fixed unresolved eddy flux out of the domain at the outer boundary, both with luminosities equal to that of the modeled stellar luminosity. Within the domain, the radiative energy flux is nonzero, but it decreases in importance with increasing radius in the convection zone. The enthalpy flux, associated with convective motions, carries the bulk of the remaining energy flux from the lower boundary, across the convection zone, to the outer boundary. Due to the asymmetry between broad, slow upflows and narrow, fast downflows, there is generally a net negative (downward) flux of kinetic energy in the convection zone (e.g., as in the right panel of Fig. 4; Cattaneo et al., 1991;Hurlburt et al., 1986). As a result, the steady-state enthalpy luminosity (dash-triple-dotted line) peaks at a value significantly larger than L * , in order that the net energy flow across the convection zone is constant and equal to the stellar luminosity. As demonstrated by the two extreme cases in Figure 4, the presence of a super-luminal enthalpy flux is increasingly important with mass, in our 4 simulations. It is not clear whether this is due to an intrinsic property of the stars (such as the luminosity), or whether it only depends upon (e.g.) the viscosity, which is systematically different in each of our cases (table 2). However, it does seem that the faster flows found in the more massive stars (discussed below) yields a larger negative kinetic energy flux, since the difference between the broad-slow upflows and narrow-fast downflows is accentuated.

Although the energy flux balance is established relatively rapidly in the convection zone, the statistical properties of the turbulence can take longer to settle down, as they evolve in response to (e.g.) Reynolds stresses and the dissipation of energy on small scales. In addition, these systems typically exhibit global, axisymmetric flows such as differential rotation and meridional circulation. These flows are also influenced by the Reynolds stresses, and the meridional circulation typically has velocities smaller than the convective velocities. Thus, the global flows settle down on timescales of the order of 10 times the convective turnover time. In order to ensure the simulations have time to establish these global flows and are near their turbulent statistical steady-state, we have run each model for a timescale comparable to the global viscous dissipation timescale (L 2 /ν).

As an example, Figure 5 shows the globally integrated kinetic energies in the 0.5 M case, as a function of time from the start of the simulation. The evolution during the first year is characterized by an initially linear (exponential) growth of the convective instability, followed by a nonlinear saturation phase, then by relaxation oscillations until it settles down into a statistically steady state. The longer term evolution exhibits a growth in the differential rotation, which also reaches a steady state. It is clear from Figure 5 that there is little change during the last 5000 days of evolution. The results in Figure 4 and presented below are shown during this mature state of the systems.

3D simulations

The ending time of each of the 4 simulations (t f ) is listed in table 3. In this section, we present the convective properties and the global flows existing in each case.

Properties of convection

Figure 6 shows the radial velocity on slices of constant radius at slightly below the top of the domain and at the middle of the domain, for the 0.5 (top), 0.7 (middle), and 1.1 (bottom) M cases. Also shown are the slices of temperature at the middle of the domain (right panels). The temperature slices generally appear less structured than the radial velocity slices because the Prandtl number is less than unity (0.25), so that the thermal difussivity (κ) is larger than the kinematic viscosity (ν).

The convective patterns and characteristic eddy sizes are influenced strongly by the density scale height, as well as by the contrast between horizontal and vertical velocity (see, e.g., Bessolaz & Brun, 2011). Near the surface, where the density scale height is smallest, the shell slices exhibit typical convective patterns characterized by hot, broad/patchy upflows surrounded by cool, narrow downflow lanes (e.g., Cattaneo et al., 1991;Miesch et al., 2008). Deeper in the convection zone, the downflow lanes often merge and generally lose their connectivity with respect to the surface patterns, and the convection has a less patchy appearance. The near-surface shell slices (left panels of Fig. 6) are shown at a location where the unresolved eddy flux dominates the other fluxes (see Fig. 4). The unresolved eddy flux is a spherically symmetric quantity, and it does not affect the shape of the convective patterns, but it does influence the level of (resolved) turbulence, since the unresolved eddy flux reduces the convective driving. However, at the location of the nearsurface shell slices shown, all cases are still turbulent, with a significant convective (enthalpy) flux and Reynolds numbers in the range of 10 -30.

In all of our cases, the convective patterns evolve on typical timescales of a few weeks to a few months, with cells merging, splitting, and disappearing. The evolution of the convective motions in the 1.1 M case generally evolve on timescales of a factor of a few times shorter than in the 0.5 M case. The convective patterns are partly advected by the local shear and, depending on whether the flow is prograde of retrograde, are tilted to the right or the left near the equator. At higher latitude, beyond the tangent cylinder (an imaginary circle crossing the upper surface with a cylindrical radius of R cz ), the convective cells are less aligned with the rotation axis and exhibit a more patchy behavior.

Table 3 Simulation results for each case. The amount of evolution time since the start of the simulation is denoted t f . Temporal averages of the rms components of velocity ṽr , ṽφ , ṽθ , speed ṽ, and fluctuating velocities ṽ φ and ṽ are evaluated at the midlevel of the domain and given in units of m s -1 . The temporally averaged rms temperature fluctuation at midlevel is T . Also listed are the rms Reynolds number Re ≡ ṽ L/ν, Rossby number Ro ≡ ṽ /(2Ω 0 L), and Péclet number Pe ≡ ṽ L/κ, evaluated at the midlevel, and the convective turnover timescale τ to ≡ L/ṽ r . The ∆Ω is the difference in angular rotation rate at the outer boundary, between latitudes of 0 • and 60 • . The ∆T is the difference in the temporally and azimuthally averaged temperature at the base of the convection zone (R cz ), between latitudes of 60 • and 0 • . Finally, ṽmc is the rms meridional circulation speed, calculated from temporally and azimuthally averaged poloidal velocity and by taking the rms of all values at a constant midlevel radius. As evident in the Figure, the 0.5 M case exhibits convective patterns near the equatorial region that are inhomogeneous, consisting of vigorous convection over some range of longitudes and relatively quiescent regions at other longitudes. Such "active nest" convective patterns have been discussed by Ballot et al. (2007) and Brown et al. (2008) and are likely related to our choice of low Prandtl number. The remaining 3 cases, exhibit vigorous convection within the entire volume of the convection zone, as seen for the 0.7 and 1.1 M cases shown in the bottom two rows of Figure 6.

Mass t f ṽr ṽθ ṽφ ṽ φ ṽ ṽ T Re Ro Pe τto ∆Ω/Ω0 ∆T ṽmc (M ) (years) (m/s) (m/s) (m/s) (m/s) (m/s) (m/s) (K) (days) (%) (K) (m/
Table 3 lists the rms velocities at the middle of the computational domain for each case. There is a clear trend of increasing convective velocities with increasing stellar mass. As evident in the right panels of Figure 6 and in the values of the rms temperature ( T ) listed in table 3, there is also a strong trend of increasing temperature contrast between up-and downflows, with increasing mass. This trend is steeper if one considers the temperature fluctuation relative to the mean temperature, since the mean background temperature in the 1.1 M convection zone midlevel is ∼3 times colder than for the 0.5 M star (see §2.1). Furthermore, the average mass density in the convective envelope decreases with increasing stellar mass (the midlevel density in the 1.1 M case is 280 times lower than in the 0.5 M case), which decreases the efficiency of convective energy transport. In spite of the lower density, the convection in higher mass stars is able to carry a larger luminosity than the lower mass stars, due to a combination of larger stellar radii, higher convective velocities, and higher temperature contrasts between up-and downflows.

Table 3 also lists the Reynolds, Rossby, and Péclet numbers for each model. The Reynolds (and Péclet) have a range of a factor of a few, with no trend in mass. These values pri-marily reflect our choice of ν top for each case (tab. 2) rather than any intrinsic property of the stars. For our choices of ν top , the 0.5 M case has the highest Reynolds number. On the other hand, the Rossby numbers do not depend strongly on our choice of diffusivities, and they show a strong dependence on mass, primarily due to the trend in convective velocity with mass, at a constant rotation period. Finally, table 3 lists a characteristic convective turnover timescale for each case. This also shows a strong trend with mass, primarily due to the trend in convective velocities.

Differential rotation and meridional circulation

Figure 7 shows the differential rotation in all 4 cases. The 2D plots in the Figure exhibit isorotation contours that are nearly aligned on cylinders in all cases. This behavior is typical for simulations such as these that have no enhanced latitudinal entropy gradient present at the lower boundary of the convection zone (Ballot et al., 2007;Miesch et al., 2006).

The angular rotation rate in the 0.5 M case has a minimum value at mid latitudes, a local maximum in the outer equatorial region, and a global maximum near the pole. A similarly "banded" differential rotation pattern is evident in a number of previously published ASH simulations (e.g. Bessolaz & Brun, 2011;Browning, 2008). The differential rotation in the 0.7 and 0.9 M cases are the most solar-like, with the fastest angular rotation rate at the equator and slowest at higher latitudes. However, the 1.1 M case exhibits anti-solar rotation, where the slowest angular velocity is at the equator. A reversal of the sense of differential rotation was also observed in the simulations of Bessolaz & Brun (2011) to be an effect of the thickness of the convection zone. In that study, the stellar structure and luminosity was held fixed, while only the thickness of the convection zone was varried. The thickness of the convection zone in our 1.1 M case (as a fraction of R * ) is comparable to the models of Bessolaz & Brun (2011) showing anti-solar rotation. This would suggest that the switch from solar to anti-solar may have more to do with convection zone thickness than with stellar structural properties. At the same time, our choice of parameters may also have influenced this outcome. In particular, the 1.1 M case has the largest Rossby numbers (see tabs. 2 and 3), which means that the convection is less influenced by rotation than in the other cases.

Each case exhibits a time-averaged, axisymmetric, latitudinal temperature gradient. For the 3 highest mass cases, this is characterized by a monotonic change from equator to pole. The 1.1 M case has a hotter equator, while the 0.7 and 0.9 M cases have hotter poles (as in the sun). The 0.5 M case has a more complex temperature pattern with a maximum at mid latitudes, while the latitudinal entropy gradient in all cases follows a more monotonic behavior from equator to pole. The second-to-last column of table 3 shows the temporally and azimuthally averaged temperature difference between the equator and a latitude of 60 • , at the base of the convection zone. It is evident that the temperature gradients are steaper for higher mass stars. The more complex temperature structure of the lowest mass case is consistent with the banded differential rotation pattern, and the reversed polarity of the temperature variation in the highest mass case is consistent with anti-solar differential rotation.

The fractional differential rotation of each case, measured at the top of the domain between latitudes 0 and 60 degrees, is listed in the final column of table 3. These values of ∆Ω/Ω 0 in the table, show that the magnitude of the differential rotation ranges from 30 to 42% for all cases.

Figure 8 shows the meridional circulation for all 4 cases. The plots are shown to scale, according to the physical size of each of the stars. For the 2 highest mass stars, the circulation pattern is dominated by a global, unicellular flow pattern from equator to pole near the surface and from pole to equator near the base of the convection zone. However, as evident in the Figure, there is a trend for this flow to break up into a multi-cellular network of circulation patterns roughly aligned on cylinders, for the lower mass stars. The tendancy for the meridional flow to breakup in this way for increasing convection zone depth was also seen in the simulations of Bessolaz & Brun (2011), suggesting that this effect is mostly due to the size of the convection zone, rather than with other stellar properties.

The last column of table 3 shows the characteristic meridional circulation flow speed evaluated at the midlevel of the domain for each case. As with the convective velocities, there is a clear trend of increasing meridional circulation flow speed with increasing stellar mass. A similar trend is also seen in 2D models (e.g., Küker & Rüdiger, 2008). Since the physical size of the circulating patterns (evident in Fig. 8) also increases strongly with mass, this suggests that the characteristic timescale for these circulations is not strongly dependent on stellar mass. The values listed in table 3 also suggest that the ratio of the meridional circulation speed to the convection speed (e.g., ṽmc /ṽ) range from 2% for the 3 lowest mass stars to 9% for the 1.1 M case.

Summary & discussion

We have presented results from 3D dynamical simulations of the convection zones of 4 sun-like stars, spanning a mass range from 0.5 to 1.1 M and rotating at the solar rate. As presented in section 2.1, these stars cover a range in luminosities of a factor of 40 and have convection zone masses ranging from 18% to 1.1% of a solar mass. The physical sizes of the convection zones increase slightly with mass (ranging from 19% to 31% of a solar radius), although the fractional size decreases with mass (ranging from 44% to 24% of a stellar radius).

We found that the convective velocities and temperature contrast between up-and downflows is a strong function of mass. This is due to the strong luminosity dependence on mass, coupled with the fact that higher mass stars have less dense convective envelopes and thus require more vigorous convection. As a consequence of the trend in convective velocities, the Rossby number (for a constant Ω 0 ) increases and the convective turnover time decreases significantly with increasing mass.

The convective pattern in the 0.5 M case exhibits "active nests" on a finite range of longitudes in the equatorial region. The presence of these active nests appears to be related to both the Rossby and Prandtl numbers, in a sense that they have been seen in simulations with low R o number and with P r = 0.25 (Ballot et al., 2007;Brown et al., 2008) but not in cases with P r = 1.0 (Ballot et al., 2007;Bessolaz & Brun, 2011). This may be due to the relative amplitude of the local shear, which tends to be stronger for lower P r . Stronger shear disrupts convective motions and can locally inhibit radial motions.

The angular rotation rate has a peak value at the surface and equator for the 0.5, 0.7, and 0.9 M cases. However, our simulations exhibit a differential rotation t ¯hat is "banded," with a minimum at mid latitudes and another maximum at the poles, for the 0.5 M case. The 0.7 and 0.9 M cases exhibit the most solar-like differential rotation, with a maximum at the equator and minimum near the poles, but the 1.1 M case exhibits anti-solar differential rotation. The magnitude of the differential rotation is comparable in all cases.

The 4 simulations presented here exhibit isorotation contours that are nearly aligned on cylinders, approximately following the Taylor-Proudman constraint. From simulations of the solar convection zone (e.g., Miesch et al., 2006), we know that the Taylor-Proudman constraint can be broken by a thermal wind driven by very small latitudinal gradients of temperature or entropy, with differences of a few Kelvin sufficing to move cylindrical profiles to radially-aligned rotation profiles. Such gradients and thermal winds do occur in the simulations presented here, but the presence of the domain boundary at the base of the convection zone somewhat inhibits this effect. By moving the lower domain boundary deep into the radiation zone, Brun et al. (2011) demonstrated that the presence of a dynamically self-consistent tachocline alleviates the effects of the lower boundary on the convection zone dynamics and results in a stronger thermal wind, isorotation contours more aligned to the radial direction, and more realistic meridional circulations and convective penetration. In future work, similar models will be developed for the stars in the present study.

The meridional circulation is dominated by a single, global circulation pattern in each hemisphere for the most massive cases, but breaks up into smaller scale patterns, aligned with the cylindrical z-direction, for lower mass stars.

Future work with these models should explore cases that include, for example, (a) a part of the radiation zone in the computational domain, in order to determine how the presence of a tachocline affects the convection zone properties, (b) the effects of different rotation rates, in particular faster rotation, to understand how the properties of these stars differ at younger ages, and (c) magnetic fields, in order to address the magnetic dynamo properties of these stars.

Fig. 4

 4 Fig.4Energy flow in 0.5 (left) and 1.1 M (right) cases as a function of radius in the computational domain. Shown are the spherically-and time-averaged luminosities, after the simulations have reached a statistical steady-state.

Fig. 5

 5 Fig. 5 Integrated kinetic energies in the domain, divided by the total volume, as a function of time, for the 0.5 M case. Shown are the kinetic energies associated with the differential rotation (DRKE), meridional circulation (MCKE), convective motions (CKE), and the total (solid line).

Fig. 6

 6 Fig.6Slices at constant radius (Mollweide view) for the 0.5 M case (top), 0.7 M case (middle), and 1.1 M case (bottom), showing the radial velocity near the surface (left panels) and in the midlevel (middle panels) of the domain and the temperature at midlevel (right panels). In the temperature slices, the azimuthal average has been subtracted, to emphasize the temperature fluctuations. The dotted line corresponds to the stellar surface and the dashed line to the equator.

Fig. 7

 7 Fig.7Azimuthally and temporally averaged Rotation frequency in a meridional slice through the domain (contour/color plots; rotation axis is vertical) and shown along radial lines at different latitudes (line plots, latitude given in degrees). Shown are the 0.5 (top left), 0.7 (top right), 0.9 (bottom left), and 1.1 (bottom right) M cases.

Fig. 8

 8 Fig. 8 Contours of the stream function of the meridional flow for the 0.5, 0.7, 0.9, and 1.1 M cases (left to right), where the physical sizes of the stars are shown to scale. The contours follow streamlines in the meridional plane from temporally and azimuthally averaged data. The colors indicate clockwise (blue) and counter-clockwise (red) circulation, and the stellar rotation axis is vertical.

www.an-journal.orgc 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

c 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.an-journal.org

c 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Acknowledgements. SPM, ODC, and ASB were supported by the ERC through grant 207430 STARS2. BPB is supported in part by NSF Astronomy and Astrophysics postdoctoral fellowship AST 09-02004. CMSO is supported by NSF grant PHY 08-21899.