Dust Polarization Maps from TIGRESS: E/B power asymmetry and TE correlation

TL;DR

This paper presents the first large-scale synthetic dust polarization maps from a self-consistent MHD simulation, revealing E/B asymmetry and TE correlation consistent with Planck observations, with detailed analysis of fluctuations and observer dependence.

Contribution

It introduces a comprehensive set of synthetic dust polarization maps derived from TIGRESS MHD simulations, capturing turbulence, magnetic fields, and star formation effects.

Findings

Prevalent E/B power asymmetry with EE>BB

Positive TE correlation consistent with observations

Significant fluctuations depending on observer position and star formation bursts

Abstract

We present the first large set of all-sky synthetic dust polarization maps derived directly from a self-consistent magnetohydrodynamics simulation using the TIGRESS framework. Turbulence in this simulation is predominantly driven by supernova explosions, with rates that are self-consistently regulated by feedback loops. The simulation covers both the outer scale and inertial range of turbulence with uniformly high resolution. The shearing-box utilized in the simulation in concert with resolved supernova-driven turbulence enables to capture the generation, growth, and saturation of both turbulent and mean magnetic fields. We construct polarization maps at 353 GHz as seen by observers inside a model of the multiphase, turbulent, magnetized interstellar medium (ISM). To fully sample the simulated ISM state, we use 350 snapshots spanning over ~350Myr (more than six feedback loops) and nine representative observers. The synthetic skies show a prevalent E/B power asymmetry (EE>BB) and positive TE correlation in broad agreement with observations by the Planck satellite. However, the ranges of EE/BB~1.4-1.7 and TE/(TT*EE)^{1/2}~0.2-0.3 are generally lower. We find large fluctuations of E/B asymmetry and TE correlation depending on the observer's position, and temporal fluctuations of interstellar medium properties due to bursts of star formation. The synthetic maps are made publicly available to provide novel models of the microwave sky.

Figures

Equations

$$v_{\mathrm{turb}}\equiv {\left(\frac{{\sum }_{2p}\rho |\delta \mathbf{v}{|}^{2}dV}{{\sum }_{2p}\rho dV}\right)}^{1/2},$$

$$v_A\equiv {\left(\frac{{\sum }_{2p}|\mathbf{B}{|}^2/(4\pi )dV}{{\sum }_{2p}\rho dV}\right)}^{1/2},$$

$$c_s\equiv {\left(\frac{{\sum }_{2p}PdV}{{\sum }_{2p}\rho dV}\right)}^{1/2}.$$

$${\mathcal{M}}_s\equiv v_{\mathrm{turb}}/c_s,$$

$${\mathcal{M}}_A\equiv v_{\mathrm{turb}}/v_A,$$

$$\beta \equiv \frac{{\sum }_{2p}PdV}{{\sum }_{2p}|\mathbf{B}{|}^2/(8\pi )dV}=\frac{2c_s^2}{v_A^2}.$$

$$\begin{array}[]{llllllll}\hat{\mathbf{X}}=&\cos\theta\cos\phi&\hat{\mathbf{x}}&+&\cos\theta\sin\phi&\hat{\mathbf{y}}&-&\sin\theta\hat{\mathbf{z}}\\ \hat{\mathbf{Y}}=&-\sin\phi&\hat{\mathbf{x}}&+&\cos\phi&\hat{\mathbf{y}}&&\\ \hat{\mathbf{Z}}=&\sin\theta\cos\phi&\hat{\mathbf{x}}&+&\sin\theta\sin\phi&\hat{\mathbf{y}}&+&\cos\theta\hat{\mathbf{z}}.\end{array}$$

$$I_{\nu }$$

$$Q_{\nu }$$

$$U_{\nu }$$

$$s_{\min }=\frac{\Delta x}{\pi }{\ell }_{\max }=163\mathrm{pc}\left(\frac{\Delta x}{4\mathrm{pc}}\right){\left(\frac{{\ell }_{\max }}{128}\right)}^{-1}.$$

$$D_{\ell }^{XY}=A_{80}^{XY}(80/\ell {)}^{{\alpha }_{XY}+2},$$

Full text

galaxies: ISM – ISM: magnetic fields – turbulence – polarization – magnetohydrodynamics (MHD) – methods: numerical

††journal: Astrophysical Journal††software: Athena (Stone et al., 2008; Stone & Gardiner, 2009), yt (Turk et al., 2011), astropy (Astropy Collaboration et al., 2013), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), IPython (Perez & Granger, 2007), pandas (McKinney, 2010), magnetar (Soler et al., 2013), healpy (Górski et al., 2005).

\turnoffedit

1 Introduction

Interstellar magnetic fields are believed to play a crucial role in controlling star formation either by directly providing additional support to the self-gravitating entities on large and small scales (Shu et al., 1987) or by altering the characteristics of turbulence (Elmegreen & Scalo, 2004; McKee & Ostriker, 2007). As a consequence, understanding interstellar magnetic fields is of paramount importance in the field of star formation and the study of the interstellar medium (ISM).

Measurements of the polarization of the microwave sky by the Planck satellite have recently delivered a wealth of information about interstellar magnetic fields. For example, because non-spherical interstellar dust grains are expected to align tightly with local magnetic fields (e.g. Hoang & Lazarian, 2008; Andersson et al., 2015), the polarized thermal dust emission directly traces the interstellar magnetic fields shaped by a variety of nonlinear physical processes, including large scale energy injection from star formation feedback, turbulence cascade, and galactic dynamo (Haverkorn, 2015, for a review).

Planck Collaboration et al. (2016, PXXX hereafter) characterized the $E$-mode and $B$-mode angular power spectra of the 353 GHz polarization maps over the multipole range $40<\ell<600$. $E$ and $B$ modes form a basis of linear polarization, and are expected to be equal for randomly oriented polarization (e.g., Kamionkowski et al., 1997; Zaldarriaga & Seljak, 1997; Zaldarriaga, 2001). However, the observations show that the $E$-mode (gradient) power is about a factor two larger than the $B$-mode (curl) power (see also Planck Collaboration et al., 2018, for recent updates). In the same analysis, $E$-mode polarization shows a positive correlation with temperature anisotropy ($TE$ correlation). $EB$ correlations are statistically insignificant but there are recent claims of a weak positive $TB$ correlation (Planck Collaboration et al., 2018; Bracco et al., 2019).

These Planck observations provide new constraints that models of the local ISM must satisfy. An asymmetry in the $E$- and $B$-mode powers may emerge naturally from the perspective of magnetohydrodynamic (MHD) turbulence in the ISM because turbulence is not random motion, but follows a specific scaling relation (e.g., Kolmogorov, 1941; Burgers, 1972; Beresnyak, 2011) with scale-dependent anisotropy in the presence of strong magnetic fields (e.g., Goldreich & Sridhar, 1995; Cho & Vishniac, 2000; Cho & Lazarian, 2003; Boldyrev, 2006).

Using an analytic description of the statistical properties of MHD turbulence, Caldwell et al. (2017) claimed that the observed $E/B$ power asymmetry and positive $TE$ correlation are not easily reproduced by simple MHD turbulence prescriptions decomposed into linear MHD waves. The model adopted an arbitrary functional form for an anisotropy, and the authors questioned a narrow parameter range for an anisotropy parameter (as well as the plasma beta defined by the ratio between thermal and magnetic pressures) required by the observed constraints. Kandel et al. (2017) performed a similar analysis and reached a different conclusion. Kandel et al. (2017) replaced the anisotropy parameter with Alfvén Mach number using an anisotropy function derived from the Goldreich-Sridhar theory (Goldreich & Sridhar, 1995; Cho & Lazarian, 2003; Lazarian & Pogosyan, 2012). In this model, the observed $E/B$ power asymmetry can be reconciled with sub-Alfvénic turbulence (magnetic energy is larger than turbulence kinetic energy), while the positive $TE$ correlation additionally requires uncorrelated magnetic field and density (Kandel et al., 2018). The desired parameter range for turbulence is still narrow (especially if as suggested in Kandel et al. (2017) slow and Alfvén waves dominate), but these conditions appear plausible for the high-latitude, low-density ISM analyzed in PXXX. Kandel et al. (2018) further suggested that an asymmetry would exist in synchrotron polarization as well in sub-Alfvénic turbulence, consistent with both the WMAP (Page et al., 2007) and Planck data.

These studies give some hints how MHD turbulence could be responsible for the observed results. However, a direct application of this analysis to the real ISM is not entirely straightforward in part because it is difficult to develop an analytic theory of MHD turbulence for the multiphase ISM. The real ISM is subject to a variety of cooling and heating processes, resulting in distinct thermal phases with typical densities and temperatures that span many orders of magnitude (e.g., McKee & Ostriker, 1977; Wolfire et al., 2003; Kim & Ostriker, 2017).

A complementary approach to understand the observed dust polarization is to use direct numerical simulations of MHD turbulence for the multiphase ISM. Kritsuk et al. (2018) investigated MHD simulations of driven multiphase turbulence in a periodic box. In these simulations, turbulence is driven by large-scale random solenoidal forcing, and the simulations include major cooling and heating processes for the cold ($T\sim 10^{2}\;{\rm K}$) and warm ($T\sim 10^{4}\;{\rm K}$) medium. The synthetic dust maps based on these simulations of the two-phase ISM show good agreement with the observed constraints. The results depend on the “turbulence parameters” (e.g., Alfvén Mach number, sonic Mach number), which are given as initial conditions in this type of simulations. For strong magnetic fields, the simulations show a weak dependence on the choice of projection with respect to the mean field direction. More targeted simulations informed by observations and/or large-scale simulations are required to refine the comparison with data.

Polarized microwave emission also plays a crucial role for next-generation cosmic microwave background (CMB) experiments that aim to detect the unique B-mode signal induced by primordial gravitational waves because the primordial B-mode signal may be weaker by many orders of magnitude than the galactic foreground signals from dust and synchrotron (Choi & Page, 2015; Kamionkowski & Kovetz, 2016). The importance of galactic foregrounds was highlighted by the claim of a detection of the primodial gravitational waves by the BICEP2 collaboration (BICEP2 Collaboration et al., 2014) and retraction after incorporating the polarized emission from galactic dust (BICEP2/Keck Collaboration and Planck Collaboration et al., 2015), which had been underestimated (e.g., Flauger et al., 2014; Mortonson & Seljak, 2014). A better understanding of the ISM can thus aid in the interpretation of current and future microwave data.

In addition, simulations of the microwave sky play an important role in the design of future CMB experiments. The most commonly used models in design studies are based on templates of synchrotron and dust emission derived from the WMAP and Planck data. Because future experiments are significantly more sensitive than WMAP and Planck, the data is smoothed to remove instrumental noise. This also removes the small scale galactic emission which is then reintroduced in the form of Gaussian random fluctuations (e.g., Delabrouille et al., 2013; Hervías-Caimapo et al., 2016; Thorne et al., 2017; Martínez-Solaeche et al., 2018). The ISM is highly non-Gaussian, and models derived from realistic MHD simulations that naturally capture the statistical properties of the ISM are likely to play a key role in future design studies and ultimately the analysis and interpretation of data (e.g., Delabrouille et al., 2018; Simons Observatory Collaboration et al., 2018).

Analytic theory and periodic box simulations allow to study the statistical properties applicable to the inertial range of turbulence. In addition, the scales of interest for future CMB missions for B-mode observations may include the energy injection scale (the reionization bump at \ell\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10; Kamionkowski & Kovetz 2016). On large scales, the turbulence that shapes the distribution of the density and magnetic fields in the ISM is not the result of self-similar cascades, but is primarily driven by energy injection mechanisms and determined by interaction with large-scale physical structure. The predominant energy source of ISM turbulence is star formation feedback (mainly supernova explosions; Mac Low & Klessen 2004; McKee & Ostriker 2007). The star formation rates (SFRs) are highly time-dependent and spatially and temporally correlated. Therefore, comprehensive direct numerical simulations are necessary to produce realistic polarized signals from the multiphase, turbulent, magnetized ISM in galactic disks.

In this paper, we utilize a numerical simulation framework of the star-forming ISM called TIGRESS (Three-phase ISM in Galaxies Evolving with Resolved Star formation and Supernova feedback; Kim & Ostriker 2017). In the TIGRESS simulation framework, we model local patches of vertically-stratified galactic disks, including self-consistent treatment of star formation and feedback, interstellar cooling and heating, magnetic fields, and galactic differential rotation. This local shearing box approach enables us to achieve uniformly high resolution to resolve key physics in the diffuse ISM. As a result, star formation rates, turbulence, magnetic fields, and thermal phases are self-consistently regulated by star formation feedback (Kim & Ostriker, 2017). Most importantly, in the context of dust polarization, the small-scale turbulent and large-scale mean magnetic fields are generated, grown, and saturated by supernova driven turbulence and galactic differential rotation (e.g., Gent et al., 2013; Kim & Ostriker, 2015a).

In what follows, we construct a large number of synthetic dust polarization maps using a fiducial solar neighborhood model of the TIGRESS simulation suite (Kim & Ostriker, 2017, 2018), analyze angular power spectra, and compare dust polarization characteristics observed in the simulations with the ISM properties observed by Planck. In Section 2, we provide an overview of the numerical simulation (§ 2.1) and a method to synthesize dust polarization maps (§ 2.2). We then calculate power spectra as described in Section 3.1, show example spectra in Section 3.2, and overall statistics in Section 3.3. Section 4 presents the full time evolution of the $EE/BB$ ratio and the $TE$ correlation seen in the synthetic maps. Sections 4.1, 4.2, and 4.3 respectively delineate the time evolution of dust polarization characteristics focusing on instantaneous, quasi-periodic, and secular variations in the context of the temporal fluctuations of the ISM properties. In Section 5, we summarize the results and discuss astrophysical and cosmological implications. Our synthetic maps are made publicly available111http://tigress-web.princeton.edu/~changgoo/TIGRESS-dustpol-maps/index.html.

2 Synthetic Dust Polarization Maps

2.1 TIGRESS MHD simulation

Numerical MHD simulations that serve for polarized foreground studies for CMB experiments should provide realistic magnetized ISM structure as generated by nonlinear interactions between turbulence, sheared galactic rotation, self-consistent star-formation and supernova feedback. All of these elements are incorporated in the newly-developed TIGRESS framework (Kim & Ostriker, 2017), in which we model local patches of vertically-stratified galactic disks using the Athena MHD code (Stone et al., 2008).

The advantage of this shearing tall-box setup is twofold. First, by modeling only a portion of galactic disks, we can achieve a uniformly high resolution to resolve key physics in the diffuse ISM. Our fiducial simulation adopts uniform spatial resolution $\Delta x=4\;{\rm pc}$ using a number of zones of $(N_{x},N_{y},N_{z})=(256,256,1792)$ with a box size of $(L_{x},L_{y},L_{z})=(1024,1024,7168)\;{\rm pc}$. Adaptive mesh refinement or (semi-)Lagrangian techniques like moving-mesh and smoothed particle hydrodynamics usually follow mass elements so that they under-resolve the diffuse ISM, which is crucial for the present purpose. Second, since it is global in the vertical direction, the realistic vertical large-scale structure of the ISM is well captured and allows a straightforward line-of-sight integration to construct reliable maps, especially toward the high-latitude sky. See Figure 1 for a sample snapshot showing three dimensional magnetic field line structures along with slices of density and temperature.

The current TIGRESS framework includes: (1) star formation in self-gravitating gas with sink particles, (2) numerically well-resolved supernovae (Kim & Ostriker, 2015b) from star clusters and runaway OB stars, (3) photoelectric heating from a time-dependent FUV field, and (4) galactic differential rotation. Kim & Ostriker (2017) have shown that a fiducial TIGRESS simulation with parameters similar to those of the solar neighborhood reaches a fully self-consistent quasi-steady state, with more than ten self-regulation cycles of star formation and feedback. The simulation shows realistic gas properties, including mass and volume fractions of cold, warm, and hot phases, turbulence velocities, magnetic field strengths (both mean and turbulent components), and disk scale heights. We refer the reader to Kim & Ostriker (2017) for more complete descriptions of the TIGRESS framework. Since the ISM properties and SFRs in TIGRESS simulations are self-regulated for given galactic conditions, no specific parameter choice for turbulence and magnetic field strength and geometry is necessary.

To overview the key ISM properties in the simulation, Figure 2 presents the time evolution of (a) surface density of star formation and supernova rates, (b) mass-weighted velocity dispersions, and (c) sonic and Alfvén Mach numbers and plasma $\beta$. The gray area is time range used to construct synthetic observations and analyze power spectra (with an interval of $\sim 1\;{\rm Myr}$). The SFR surface density is calculated from the mass of young star cluster particles with the age of star cluster particles younger than $t_{\rm bin}=10\;{\rm Myr}$ (roughly equivalent to SFRs traced by H$\alpha$ emissions; Kennicutt & Evans 2012). We count the total number of supernovae within the last $10\;{\rm Myr}$ to calculate the areal supernova rate $\Sigma_{\rm SNR}$ and plot in Figure 2(a) by multiplying the mass of young stars formed per SN, $m_{*}=95.5\;M_{\odot}$ (Kim & Ostriker, 2017). The recent star formation and supernova rates are self-regulated to maintain an equilibrium state, $\Sigma_{\rm SFR,10}=2.4^{+1.8}_{-0.9}\times 10^{-3}\;M_{\odot}\;{\rm kpc}^{-2}\;{\rm yr}^{-1}$ (ranges indicate the first and third quartiles), that is consistent with the observed $\Sigma_{\rm SFR}$ in the solar neighborhood (Fuchs et al., 2009), with large temporal fluctuations.222Since our model does not have inflows, the gas surface density (and total gas mass) decreases gradually as the gas is gradually converted into the stars and escapes as galactic winds (Kim & Ostriker, 2018). The surface density during the time range of interest varies from $\Sigma=10.5\;M_{\odot}\;{\rm pc^{-2}}$ to $8\;M_{\odot}\;{\rm pc^{-2}}$, covering the solar neighborhood condition.

The mass-weighted turbulent333For the turbulent velocity dispersion, we subtract the background shear arising from galactic differential rotation, $\mathbf{v}_{\rm shear}=-q\Omega x\hat{\mathbf{y}}$), and use the perturbed velocity of gas. $\delta\mathbf{v}\equiv\mathbf{v}-\mathbf{v}_{\rm shear}$. Here, we adopt a flat rotation curve with the shear parameter $q\equiv|d\ln\Omega/d\ln R|=1$ and galactic rotational speed $\Omega=28\;{\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$., Alfvén, and thermal velocity dispersions are respectively defined by

[TABLE]

[TABLE]

[TABLE]

Here, we further decompose Alfvén velocity into mean and turbulent components using the horizontally-averaged mean field $\overline{\mathbf{B}}$ and turbulent field $\delta\mathbf{B}\equiv\mathbf{B}-\overline{\mathbf{B}}$. To obtain representative values of velocity dispersions where most dust emission arises, we exclude the hot gas using a temperature criterion $T<2\times 10^{4}\;{\rm K}$. We use the subscript “2p” in the summation to denote this “two-phase” medium (or the warm and cold medium). It is also useful to calculate dimensionless ratios between velocity dispersions, defining sonic and Alfvénic Mach numbers and plasma $\beta$ as

[TABLE]

[TABLE]

[TABLE]

We initially impose a uniform azimuthal magnetic field with a moderate strength ($\beta=10$). This gives an initial magnetic field at the midplane with $B_{y,0}(z=0)=2.3~{}\mu G$. As soon as turbulence is driven, the turbulent component of magnetic fields is rapidly created by tangling the initial fields. The turbulent Alfvén velocity dispersion is well correlated with the turbulent velocity dispersion with a slight offset as stronger turbulence makes magnetic fields more turbulent, while the mean Alfvén velocity smoothly grows through dynamo action (e.g., Gent et al., 2013; Kim & Ostriker, 2015a). As the mean field slowly grows with the time scale comparable to the orbital period ($2\pi/\Omega=220\;{\rm Myr}$), the turbulent magnetic fields grows as well, and the total magnetic energy saturates at a level similar to the kinetic energy.

In our simulation, saturation of the turbulent fields occurs after $300\;{\rm Myr}$, while the mean field growth saturates after $500\;{\rm Myr}$ (see Figure 2(b)). We also find that the saturation level is independent of the initial field strength, but it takes longer to saturate with weaker initial field strength (e.g., Kim & Ostriker, 2015a). At saturation, the strengths of the mean and turbulent magnetic fields become $|\overline{B}|=3.6~{}\mu G$ and $|\delta B|=3.9~{}\mu G$, respectively. In the period of interest (the gray area in Figure 2), the secular evolution of magnetic fields enables us to cover the global gas properties $\mathcal{M}_{A}\sim 1.6-1$ and $\beta\sim 0.7-0.2$, while the turbulent velocity fluctuations results in $\mathcal{M}_{s}\sim 1.8-3.2$.

2.2 Synthetic Dust Polarization Maps

Our simulation is carried out in local Cartesian coordinates. To synthesize all-sky dust polarization maps, we place observers in our simulation box. For a given observer’s position, we tile the sky using HEALPix (Górski et al., 2005) with $N_{\rm side}=128$. The HEALPix angles $(\phi,\theta)$ are the angles measured in the $x-y$ plane from $\hat{\mathbf{x}}$-axis in the counterclockwise direction and from the $\hat{\mathbf{z}}$-axis downward, respectively. At a given line-of-sight (LOS), the unit vectors along the LOS ($\hat{\mathbf{Z}}$) and the plane of the sky ($\hat{\mathbf{X}}$ and $\hat{\mathbf{Y}}$) are defined by rotations of the unit vectors of the Cartesian coordinates $(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}})$ by $\phi$ about $\hat{\mathbf{z}}$-axis and $\theta$ about $\hat{\mathbf{y}}$-axis (e.g., the HEALPix Primer444https://healpix.sourceforge.io/pdf/intro.pdf Fig. 5). The relationship between the two unit vectors are then given by

[TABLE]

Along each LOS, we sample the physical quantities from the observer’s position to the maximum path length $s_{\rm max}$ with an interval of $\Delta s\equiv\Delta x=4\;{\rm pc}$. We apply the shearing-periodic boundary conditions (e.g., Hawley et al., 1995) whenever we march a ray outside the simulation domain horizontally. At every position along the LOS, we extract the gas density and the magnetic field vector from the original data cube via a trilinear interpolation using the eight nearest grid zones. At every snapshot, we construct nine realizations corresponding to nine observers at $x=-256,0,256\;{\rm pc}$ and $y=-256,0,256\;{\rm pc}$ in the $z=0$ plane.

Following the HEALPix convention, the synthetic Stokes I, Q, and U parameters at frequency $\nu$ are (e.g., Planck Collaboration et al., 2015)

[TABLE]

Here, $n_{H}=\rho/(1.4271m_{p})$ is the hydrogen number density, $B_{X}$ and $B_{Y}$ are the plane of the sky components of the magnetic field vectors projected by Equation (7), and the magnitude of the magnetic field is $B=(B_{X}^{2}+B_{Y}^{2}+B_{Z}^{2})^{1/2}$. Following Planck Collaboration et al. (2015), we adopt an intrinsic polarization fraction $p_{0}=0.2$, the blackbody source function $B_{\nu}$ with an uniform dust temperature $T_{d}=18\;{\rm K}$, and the dust opacity at 353 GHz $\sigma_{d,353}=1.2\times 10^{-26}\;{\rm cm}^{-2}$.

Since our simulation utilizes uniform cubic grid zones, the coordinate conversion and integration along the radial directions cause non-trivial resolution effects. For example, due to the large covering area of nearby cells, angular resolution can become poor with an existence of nearby dense cells. If one aims to resolve modes up to a multipole moment $\ell_{\rm max}$, the minimum distance to the cell is

[TABLE]

In practice, many observers are located in a real bubble in the simulation created by SN(e), reducing such contamination from nearby cells. For the sake of uniformity, however, we take $s_{\rm min}=160\;{\rm pc}$, implying that we assume an artificial local bubble with radius of $160\;{\rm pc}$ surrounding an observer, similar to the Local Interstellar Bubble (e.g., Lallement et al., 2018). We use $s_{\rm max}=3.5\;{\rm kpc}$ to cover a majority of the gas in the simulation. However, most of the signal comes from the gas within the disk scale height ($H\sim 400\;{\rm pc}$) so that choosing any $s_{\rm max}>1\;{\rm kpc}$ results in similar synthetic maps.

The right panel of Figure 1 shows example Stokes I, Q, and U maps along with the drapery patterns calculated by the line integral convolution (Cabral & Leedom, 1993) representing the apparent magnetic field orientation in the plane of the sky.

3 Angular Power Spectra

3.1 Method

We compute the angular power spectra of our synthetic maps using the nominal curved sky pseudo-$C_{\ell}$ method, which accounts for masks imposing sky cuts (e.g., Hivon et al., 2002; Kogut et al., 2003; Brown et al., 2005). This power spectrum estimator code has been thoroughly tested and was also used in analyzing the WMAP and Planck maps to study polarized galactic synchrotron and dust emission (Choi & Page, 2015). We do not use the pure $B$-mode estimator in the following discussion as the synthetic maps have comparable power in $E$- and $B$-modes, but we have verified that the conclusions remain unchanged if a pure estimator is used.

We use $\ell$ bins of width 20 centered at 29.5, 49.5, 69.5, etc. The bin-to-bin correlation depends on the width of the $\ell$ bins, the sky fraction of the mask, and the smoothness of the mask apodization. To estimate the level of the bin-to-bin correlation, we generate 2000 Gaussian simulations with power law $EE$ and $BB$ spectra. For the mask with the smallest sky fraction described below, we find the average bin-to-bin correlation to be $<5\%$. We compute $D_{\ell}$, where $D_{\ell}\equiv[\ell(\ell+1)/2\pi]C_{\ell}$.

3.2 An Example

Figure 3 shows example $TE$ , $EE$, and $BB$ power spectra derived from synthetic maps presented in Figure 1. We construct masks for high dust intensity regions with four threshold values $I_{\rm th}=0.05$, 0.1, 0.2, and 0.3 MJy/sr, which cover a range of the dust intensity adopted in PXXX. We also apply a latitude mask, masking $|b|<5.7^{\circ}$, where long path length with periodic boundary conditions would sample the same gas artificially multiple times. For this particular realization, from the higher to lower threshold value, the visible sky fractions are $f_{\rm sky}=0.23$, 0.4, 0.56, and 0.66, respectively.

Following PXXX, we conduct least-square fits to the power spectra with a power-law function

[TABLE]

where $XY$ can be $TE$ , $EE$, and $BB$. The fits use multipole ranges of $(\ell_{\rm min},\ell_{\rm max})$, where $\ell_{\rm min}\equiv 8/f_{\rm sky}$ and $\ell_{\rm max}=128$ (see Equation 11).

We define the $EE/BB$ ratios at a given multipole and get the mean $EE/BB$ ratio (or the $E/B$ ratio, in short) $\mathcal{R}_{EB}\equiv\left\langle{D_{\ell}^{EE}/D_{\ell}^{BB}}\right\rangle$, where $\left\langle{}\right\rangle$ stands for the mean over the multipole range used in the power law fit. We also calculate the $E/B$ ratio using the amplitudes of the power-law fits, i.e., $A_{80}^{EE}/A_{80}^{BB}$, as in PXXX. The two values are consistent with each other, implying that the $EE/BB$ ratios are not strongly scale dependent; in general, $E$-mode power spectra have slightly shallower slopes than that of $B$-mode power spectra, although both power-law slopes are steeper in our simulations than those in Planck observations due to the projection effect (see Section 3.3).

We also calculate the mean value of the dimensionless $TE$ correlation coefficient (or the $TE$ correlation, in short) at a given multipole $r_{TE}\equiv\left\langle{D_{\ell}^{TE}/(D_{\ell}^{TT}D_{\ell}^{EE})^{1/2}}\right\rangle$. For the example realization shown in Figures 1 and 3, the $E/B$ ratios and the $TE$ correlations are close to the Planck measurements, but, as we shall show later, this is a rare temporal case in the suite of synthetic maps.

3.3 Overall Statistics and Convergence

We construct total 3,150 synthetic maps using 350 snapshots during $t\in(300,650)\;{\rm Myr}$ with $\Delta t=1\;{\rm Myr}$ for 9 observers in the midplane for each snapshot.555All synthetic maps are available at We then calculate power spectra using masks with intensity thresholds ($I>0.05$, 0.1, 0.2, and 0.3 MJy/sr) and latitude cuts ($|b|<5.7^{\circ}$) as described in the previous section. We further select the cases with $f_{\rm sky}>0.16$ to ensure that the sky coverage is not too small to determine meaningful power spectra within the $\ell$ range of interest. We finally get 9,175 sets of power spectra, and at least 16 sets of power spectra at a given time. We refer to this result as fiducial in what follows.

For comparison, we also synthesize 900 dust polarization maps for the first 100 snapshots using the warm medium only by applying a temperature cut, $5\times 10^{3}\;{\rm K}<T<2\times 10^{4}\;{\rm K}$, while for the calculation of power spectra, we apply the same masks used in the fiducial maps. We refer to this result as warm-only.

In addition, we create 900 maps using a lower resolution simulation ($\Delta x=8pc$) for 100 snapshots during $t\in(300,400)$. We refer to this result as low-res. Note that overall evolution and ISM properties are statistically converged at this resolution (Kim & Ostriker, 2017), but this does not imply exactly the same ISM realization with lower resolution because the evolution of the star-forming ISM is chaotic in nature. Therefore, the resolution comparison should be understood in a statistical sense.

Figure 4 shows the overall statistics of the $E/B$ ratios and $TE$ correlations as well as power law slopes of $EE$ and $BB$ power spectra. Note that we only show results from the fiducial maps from the first 100 snapshots to enable direct comparison with the warm-only and the low-res maps.

As noted in Figure 3, the slopes of the power spectra are generally steeper than those in the Planck observation. This is largely caused by the limited coverage of the inertial range of turbulence in the simulation in combination with the sky projection. and a dissipation range with a steep drop of power (e.g., Vestuto et al., 2003, for an inviscid, isothermal, compressible turbulence simulation). In inviscid numerical simulations (i.e., without explicit viscosity and resistivity), the dissipation range is set by numerical schemes and resolution, typically happening at scales of $\sim 10-20$ resolution elements depending on the numerics.

For projections along the Cartesian axes, there is a one-to-one correspondence between wavenumber and multipole moment, and the inertial and dissipation ranges are well separated. As a consequence, one can infer the slope of power spectra even for a relatively narrow inertial range (e.g., Kritsuk et al., 2018). For the radial projection of the simulations onto the sky used here, there is no correspondence between wavenumber and angular scale because even small nearby structures appear large in the sky. This artificially steepens the power-law slopes significantly compared to their values in the inertial range.

To investigate the effect more closely, we have also considered projections along the Cartesian axes. We take Equations (8)-(10), but consider the vertical axis ($\hat{\mathbf{z}}$) instead of the radial axis ($\hat{\mathbf{Z}}$). We integrate toward the upper boundary $z_{\rm max}=L_{z}/2$ by varying $z_{\rm min}=0$, 128, $\cdots$, 1024 pc. The $E/B$ ratio and slopes in the $EE$ and $BB$ power spectra are essentially unchanged if $z_{\rm min}>128\;{\rm pc}$ (equivalent to the intensity masks to exclude the cold, dense medium near the midplane).

We find shallower slopes ranging between $-2.9<\alpha<-2.6$ significantly closer to the slopes observed in Kritsuk et al. (2018) and the Planck data. In addition, since Kritsuk et al. (2018) observe an $E/B$ ratio that approaches unity on small scales where numerical dissipation is significant, one may also expect a bias on the $E/B$ ratio and $TE$ correlation in the radial projections. However, we find that the observed $E/B$ ratio and $TE$ correlation coefficients for the projections along the Cartesian coordinates are similar to those obtained for the radial projections.

This is further supported by comparisons between fiducial (blue) and low-res (green) results shown in Figure 4 which demonstrate that the systematic biases in the measurements of $\mathcal{R}_{EB}$ and $r_{TE}$ are less problematic, while the steepening of power spectra at lower resolution is evident in both $EE$ and $BB$ power spectra.

The warm-only maps also show consistent values of $\mathcal{R}_{EB}$ and $r_{TE}$ (slightly higher than the fiducial maps), while the slopes are the same. As we already masked out low-latitude, high-intensity skies, our fiducial maps are mostly dominated by the warm medium, which is better resolved than the cold medium. Overall, we conclude that our synthetic maps are converged in the sense that the measured $E/B$ ratio and $TE$ correlation coefficient are not sensitive to the projection and resolution, allowing us to investigate these two polarization characteristics in more detail.

The effects of the radial projection on the spectral slope can be mitigated by increasing the resolution near the observer, but we leave this for future work.

4 Time Evolution of Polarization Characteristics

Our numerical simulation provides a long-term evolution of the ISM for conditions similar to those in the solar neighborhood, with self-consistent changes of turbulence energy injection due to self-regulation of SFRs as well as corresponding changes of the ISM state (see Figure 2). Therefore, it is interesting to explore time evolution of the dust polarization characteristics rather than analyzing time averaged “stacked” power spectra.

Figure 5 shows $E/B$ ratios $\mathcal{R}_{EB}$ and $TE$ correlations $r_{TE}$ for all power spectra as a function of time. Note that we have a maximum of 36 (9 observers and 4 masks) power spectra at a given time, but the number can be as small as 16 after discarding cases with low sky fraction ($f_{\rm sky}<0.16$). The median and 10th and 90th percentiles smoothed by the rolling mean with the time window of $10\;{\rm Myr}$ are shown as the solid line and shaded region, respectively. The results from the warm-only maps are shown as red in the same figure.

First of all, we emphasize that we find ubiquitous $E/B$ power asymmetry and positive $TE$ correlation in the synthetic maps as seen by PXXX. However, the median values of $\mathcal{R}_{EB}\sim 1.5$ and $r_{TE}\sim 0.23$ are smaller than those measured by Planck. Figure 6 presents a heat map between $\mathcal{R}_{EB}$ and $r_{TE}$, showing that $\mathcal{R}_{EB}$ and $r_{TE}$ are generally correlated.

Given the generally higher $E$-mode power than $B$-mode power and positive $TE$ correlation, in what follows, we focus on temporal fluctuations shown in the synthetic maps. There are three main features in the time evolution of $\mathcal{R}_{EB}$ and $r_{TE}$ that we separately address in the following subsections: (1) large instantaneous variations, (2) quasi-periodic fluctuations, and (3) secular evolution.

4.1 Large Instantaneous Variations

We observe large variations in both measured characteristics $\mathcal{R}_{EB}$ and $r_{TE}$ from the same snapshot depending on observer’s positions. The mean levels of fractional variations within the same snapshot are $\delta\mathcal{R}_{EB}=0.18$ and $\delta r_{TE}=0.45$, where $\delta q\equiv(q_{90}-q_{10})/q_{50}$ is defined using the $n$-th percentile of a quantity $q$, $q_{n}$. Note that the levels of fractional variations arising from different masks (intensity thresholds) are $\delta\mathcal{R}_{EB}=0.085$ and $\delta r_{TE}=0.24$, about a half of the total fractional variation levels.

It is natural to expect that synthetic observations from the same snapshot would share similar dust polarization characteristics since they are looking through the same ISM state. However, the dust polarization signals mainly come from the nearby medium (density-weighted). In other words, the scale length of signal should be comparable to the (projected) scale height of the warm/cold medium ($H\sim 400\;{\rm pc}$ in the simulation; Kim & Ostriker 2017). The energy injection scale of turbulence is also similar to the scale height of the disk. This means that the synthetic skies seen by different observers separated by a few hundreds of parsec are effectively independent realizations.

4.2 Quasi-periodic Fluctuations

In addition to large instantaneous variations, Figure 5 also shows quasi-periodic fluctuations of mean values of both $\mathcal{R}_{EB}$ and $r_{TE}$, similarly with ISM properties presented in Figure 2. The period of these quasi-periodic fluctuations is roughly $\sim 40-50\;{\rm Myr}$, corresponding to the period of the duty cycle of star formation self-regulation; when SFRs are enhanced, SN feedback expels gas and destroys dense clouds near the midplane, and star formation is quenched until the expelled gas falls back and is recollected by gravity and large scale flows into dense clouds. Here, the governing time scale of this process is the vertical dynamical time scale, $t_{\rm ver}\sim 2H/v_{z}\sim 40-50\;{\rm Myr}$ (Kim et al., 2013; Kim & Ostriker, 2015a, 2017).666 The star formation duty cycle is also related to the star cluster lifetime (or duration of feedback); since SN rates from one star cluster are slowly decreasing as a function of time, one burst can persist longer than the vertical dynamical time. Governing time scale for the star formation duty cycle is longer one between the vertical dynamical time and feedback lifetime. For our adopted population synthesis model (Starburst99; Leitherer et al. 1999), the half-life of SN rates is about $20\;{\rm Myr}$ (see Kim & Ostriker, 2017).

The duty cycle of self-regulation also results in offsets between SFRs and thermal, turbulent, and magnetic components of velocity dispersion (or pressure). To show this clearer, Figure 7 plots SFR surface density and velocity dispersions shown in Figure 2 as a function of time along with the median values of $\mathcal{R}_{EB}$ shown in Figure 5, with an emphasis on fluctuations by adjusting plotting ranges. From top to bottom, we compare $\mathcal{R}_{EB}$ with SFR surface density, sound speed, turbulent velocity, and turbulent Alfvén velocity. One can clearly see that the patterns in velocity dispersions are shifted progressively from $c_{s}$ to $v_{\rm turb}$ to $v_{\rm A,turb}$ as each component is controlled by different physical processes of star formation feedback (see Ostriker et al., 2010; Kim et al., 2011; Ostriker & Shetty, 2011; Kim et al., 2013; Kim & Ostriker, 2015a).

If SFRs are higher than their equilibrium values, the warm medium gas fraction ($f_{w}$) increases almost immediately since FUV radiation provided by massive stars younger than $\sim 10\;{\rm Myr}$ enhances the photoelectric heating rate from small grains, which is a dominant heating process in the diffuse ISM (e.g., Wolfire et al., 1995). The mass-weighted mean sound speed $c_{s}\approx f_{w}c_{w}$ (e.g., Ostriker et al., 2010) thus simply follows the recent SFRs. After $\sim 4-5\;{\rm Myr}$ from the birth of star clusters, massive stars begin to die and explode as SNe, driving turbulence. The turbulent velocity dispersion peaks after $\sim 10\;{\rm Myr}$ from the peak of $\Sigma_{\rm SFR}$ and $c_{s}$. Then, turbulent magnetic fields are generated by tangling of existing magnetic fields due to turbulence. The peak of the turbulent Alfvén velocity is expected to be further delayed after reaching the peak of turbulence.

At a glance, $\mathcal{R}_{EB}$ shows a positive correlation with $\Sigma_{\rm SFR}$ and $c_{s}$, especially for $t<500\;{\rm Myr}$. Kritsuk et al. (2018) pointed out that dense gas, which is masked out with a number density threshold $n_{\rm th}=70\;{\rm cm}^{-3}$ in their synthetic dust polarization maps before projection, can randomize polarization to significantly reduce $\mathcal{R}_{EB}$. When $c_{s}$ is lower and the cold medium fraction is higher, the intervening cold medium is more likely to reduce the overall $\mathcal{R}_{EB}$.

A comparison between $\mathcal{R}_{EB}$ from the fiducial and warm-only maps can provide evidence for the general agreement with this expectation (see the second panel of Figure 7). $\mathcal{R}_{EB}$ in the warm-only maps are larger than that in the fiducial maps (Figure 4 also shows slight increase of $\mathcal{R}_{EB}$ for the warm medium only case), especially when $\mathcal{R}_{EB}$ is low (in other words, when the cold medium dominates). Therefore, the dips in $\mathcal{R}_{EB}$ are less prominent in the warm-only maps, in which less randomization effects are expected by the cold medium. However, at the peak of $\mathcal{R}_{EB}$, there is almost no difference between fiducial and warm-only maps as the volume is already dominated by the warm medium.

While reduced, the warm-only maps still show similar temporal fluctuations. Simple considerations for the variation of the warm medium fraction alone cannot fully explain the quasi-periodic oscillations in $\mathcal{R}_{EB}$. Although there is another apparent (anti-)correlation between $\mathcal{R}_{EB}$ and the turbulent Alfvén velocity, we do not push this analysis further with just simple, volume-integrated turbulence properties. A more careful analysis on characteristics of MHD turbulence in compressible, multiphase ISM driven by SN feedback like our simulation is needed to show any causal relations between SFRs, ISM properties, and dust polarization characteristics more robustly.

4.3 Secular Evolution

In the simulation, the total gas mass is decreasing due to star formation and mass loss, and the mean magnetic fields are growing due to galactic dynamo (Figure 2). They drive the secular evolution of SFRs and the ISM state, which potentially drive the secular evolution of the dust characteristics as well. Although local and temporal variations are much more prominent, there is a weak secular decreasing trend in $\mathcal{R}_{EB}$ but not in $r_{TE}$ (Figure 5).

To gauge the relative importance between thermal, turbulent, and magnetic states of the ISM for the dust characteristics, Figures 8 and 9 respectively plot $\mathcal{R}_{EB}$ and $r_{TE}$ measurements as a function of the dimensionless parameters such as Alfvén Mach number $\mathcal{M}_{A}$ (Eq.(5); top), sonic Mach number $\mathcal{M}_{s}$ (Eq.(4); middle), and plasma $\beta$ (Eq.(6); bottom). Overall, $\mathcal{R}_{EB}$ and $r_{TE}$ are not strongly correlated with these global ISM parameters.

While less distinctive, $\mathcal{R}_{EB}$ increases with both $\mathcal{M}_{A}$ and $\beta$, implying that $\mathcal{R}_{EB}$ is generally higher with weaker magnetic fields. This is consistent with the apparent anti-correlation between $\mathcal{R}_{EB}$ and turbulent Alvén velocity seen in Figure 7. It is also possible that $\mathcal{R}_{EB}$ is just higher with stronger turbulence (may relate to higher SFRs).

Intriguingly, $r_{TE}$ does not show such a correlation due to the secular evolution although there is a general correlation between $\mathcal{R}_{EB}$ and $r_{TE}$ as shown in Figure 6. Since similar temporal fluctuations are shown in Figure 5 for both $\mathcal{R}_{EB}$ and $r_{TE}$, the correlation between $\mathcal{R}_{EB}$ and $r_{TE}$ should arise from short-term variation of the ISM state that is responsible for quasi-periodic fluctuations rather than secular variation of turbulence properties.

In addition, comparing the state of the ISM at early and late times in the simulation (Figure 2), we see that turbulence gets weaker (both kinetic and magnetic components) as the mean field grows, but the sound speed is unchanged. The lack of secular variation of $r_{TE}$ with still substantial quasi-periodic fluctuations thus implies that $r_{TE}$ may be a more sensitive probe of the thermal state of the ISM than of the turbulence state, while both ISM properties can be responsible for $\mathcal{R}_{EB}$.

5 Summary & Discussion

We have presented the first set of full-sky synthetic dust polarization maps constructed from a self-consistent MHD simulation. Our simulation utilizes the TIGRESS framework (Kim & Ostriker, 2017), in which the magnetized ISM in a local patch of galactic disks evolves self-consistently under the influence of interstellar cooling and heating, galactic differential rotation, self and external gravity, and turbulence driven by SN feedback. We have constructed a large number of synthetic maps seen by nine representative observers in the galactic plane for time snapshots covering $\sim 350\;{\rm Myr}$ ($\sim 1.5$ orbital time; 6-7 feedback cycles). Our analysis of these maps demonstrates that the recent Planck observations of $E/B$ power asymmetry and positive $TE$ correlation emerge naturally in MHD simulations of the star-forming ISM. Our main results are summarized in more detail as follows.

From our large number of synthetic skies (3,150), we find that the $EE/BB$ ratios ($\mathcal{R}_{EB}$) are larger than unity and that the dimensionless $TE$ correlation coefficients ($r_{TE}$) are positive, in broad agreement with the Planck observations. However, the typical values seen in our synthetic skies are $\mathcal{R}_{EB}\sim 1.4-1.7$ and $r_{TE}\sim 0.2-0.3$, somewhat lower than the values reported by Planck, except for short periods during which the Planck values are reached. 2. 2.

Numerical simulations are inherently limited in the range of scales that are resolved, and numerical dissipation occurs at scales near the grid resolution, where we expect a steep drop of turbulence power. The radial projection onto the sky necessary to generate the full-sky maps mixes up the inertial and dissipation ranges of turbulence so that the angular power spectra for our synthetic skies are affected by numerical dissipation. This results in steeper slopes (see Figure 3) than in Kritsuk et al. (2018) and PXXX. To quantify this effect, we also consider a projection along the Cartesian axes and find that the slopes measured in the inertial range are significantly closer to those in both Kritsuk et al. (2018) and PXXX. However, the values for the $E/B$ ratio and $TE$ correlation we obtain for the projections along the Cartesian axes are consistent with those obtained for the radial projection and are somewhat lower than those observed by Planck. This gives us confidence in the measurements of these two characteristics from our simulations, although we cannot fully exclude the possibility that numerical dissipation and imperfect separation of the dissipation range introduces some bias. To eliminate the mixing of scales as a possible bias, the simulations should be further developed to increase the resolution near the observers. Higher resolution simulations are generally also required to better resolve the diffuse cold neutral medium and to better understand whether values of the $E/B$ ratio and the $TE$ correlations seen by Planck are indeed rare. We leave this for future work. 3. 3.

The dynamic ISM model, which includes self-consistent time evolution due to self-regulation of star formation rates, shows a distinctive temporal evolution of both $\mathcal{R}_{EB}$ and $r_{TE}$. $\mathcal{R}_{EB}$ shows a quasi-periodic temporal fluctuations with period similar to the star formation duty cycle and a secular decreasing trend following the growth of the mean fields and/or decreasing SFRs. $\mathcal{R}_{EB}$ and $r_{TE}$ are generally correlated, but $r_{TE}$ does not show the secular decrease. 4. 4.

Caldwell et al. (2017) and Kandel et al. (2017) have provided methods to calculate polarization power spectra for dust emission from MHD turbulence theory (see also Kandel et al. 2018 for synchrotron polarization). Although the interpretations differed, the results imply that the Planck constraint on the $E/B$ power asymmetry suggests a dominance of either isotropic fast waves or sufficiently anisotropic slow and Alfvén waves in a low-$\beta$ plasma. Caldwell et al. (2017) raised the questions whether such conditions are realistic, but Kandel et al. (2017) argued that the Goldreich & Sridhar (1995) theory which appear to be common in the high-latitude sky probed by PXXX.

It is potentially interesting to attempt to use the $E/B$ power asymmetry and $TE$ correlation coefficient as a probe of turbulence parameters, e.g., Alfvén Mach number and plasma $\beta$. For example, Kandel et al. (2017) claimed that $\mathcal{R}_{EB}$ would be higher as the magnetic fields gets stronger in sub-Alfvénic, low-$\beta$ turbulence. However, a clear and sensitive dependence of dust polarization statistics to global turbulence parameters are not evident in our synthetic maps. Turbulence in our simulations is driven by realistic driving mechanisms (supernovae) and rates (self-regulated by self-consistent feedback loops). The spatially and temporally localized driving interacts with large-scale structures of the ISM developed by gravitational collapse, vertical stratification, thermal instability, and galactic differential rotation. This is drastically different from “turbulence-in-a-box” simulations with prescribed driving schemes and isothermal equation of state (e.g., Cho & Lazarian, 2002, 2003; Kowal & Lazarian, 2010; Federrath et al., 2010). Therefore, the characterization of local turbulence parameters in our simulation is itself another challenge. A more careful analysis is an essential prerequisite to successfully connect dust polarization statistics with local and/or global turbulence properties in the multiphase, stratified ISM. 5. 5.

MHD simulations are a promising tool to generate realistic maps for galactic foregrounds, especially for polarized dust emission. Further improvements of foreground maps based on MHD simulations are possible and are currently underway. We are, for example, incorporating the adaptive ray tracing module available in the Athena code (Kim et al., 2017) to calculate interstellar radiation fields at every grid point. This will be combined with more sophisticated dust models (Hensley and Draine, in preparation) to obtain a spatially-varying spectral energy distribution and to produce multi-frequency maps. These simulations will allow to assess the degree of frequency decorrelation arising from variation of dust properties in inhomogeneous structures along the line-of-sight, which is a key parameter for future CMB experiments.

We are grateful to the referee for helpful report. CGK acknowledges Eve Ostriker for her continuous encouragement and Juan Soler, Susan Clark, Brandon Hensley, Colin Hill, and François Boulanger for useful discussions. Resources supporting this work were provided in part by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and in part by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center. This work took part under the program Milky-Way-Gaia of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. The Flatiron Institute is supported by the Simons Foundation. CGK and RF are supported in part by the NASA TCAN grant No. NNH17ZDA001N-TCAN. CGK also acknowledges support from the NASA ATP grant No. NNX17AG26G and from the Simons Foundation Award No. 528307. RF is grateful for support from the Alfred P. Sloan Foundation, the Department of Energy under grant de-sc0009919, and the NASA ATP Grant No. 80NSSC18K0561. SKC acknowledges support from the Cornell Presidential Postdoctoral Fellowship.