WATER FORMATION DURING THE EPOCH OF FIRST METAL ENRICHMENT

Water is an essential ingredient for life as we know it (Kasting 2010). In the interstellar medium (ISM) of the Milky Way and also in external galaxies, water has been observed in the gas phase and as grain–surface ice in a wide variety of environments. These environments include diffuse and dense molecular clouds, photon-dominated regions (PDRs), shocked gas, protostellar envelopes, and disks (see the review by van Dishoeck et al. 2013).

In diffuse and translucent clouds, H₂O is formed mainly in gas-phase reactions via ion–molecule sequences (Herbst & Klemperer 1973). The ion–molecule reaction network is driven by cosmic-ray or X-ray ionization of H and H₂, which leads to the formation of H⁺ and ${\rm H}_{3}^{+}$ ions. These interact with atomic oxygen and form OH⁺. A series of abstractions then lead to the formation of H₃O⁺, which makes OH and H₂O through dissociative recombination. This formation mechanism is generally not very efficient, and only a small fraction of the oxygen is converted into water; the rest remains in atomic form or freezes out as water ice (Hollenbach et al. 2009).

Sonnentrucker et al. (2010) showed that the abundance of water vapor within diffuse clouds in the Galaxy is remarkably constant, with ${{x}_{{{{\rm H}}_{2}}{\rm O}}}\sim {{10}^{-8}}$, that is, ∼0.1% of the available oxygen. Here, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ is the H₂O number density relative to the total hydrogen nuclei number density. Toward the Galactic center, this value can be enhanced by up to a factor of ∼3 (Monje et al. 2011; Sonnentrucker et al. 2013).

At temperatures $\gtrsim 300$ K, H₂O may form directly via the neutral–neutral reactions, O + H₂ $\to$ OH + H, followed by OH + H₂ $\to$ H₂O + H. This formation route is particularly important in shocks where the gas heats up to high temperatures and can drive most of the oxygen into H₂O (Draine et al. 1983; Kaufman & Neufeld 1996).

Temperatures of a few hundred kelvin are also expected in very low metallicity gas environments, with elemental oxygen and carbon abundances of $\lesssim {{10}^{-3}}$ solar (Bromm & Loeb 2003; Omukai et al. 2005; Glover & Clark 2014), associated with the epochs of the first enrichment of the ISM with heavy elements, in the first generation of galaxies at high redshifts (Loeb & Furlanetto 2013). At such low metallicities, cooling by fine structure transitions of metal species, such as the [Cii] 158 μm line, and by rotational transitions of heavy molecules, such as CO, becomes inefficient, and the gas remains warm.

Could the enhanced rate of H₂O formation via the neutral–neutral sequence in such warm gas compensate for the low oxygen abundance at low metallicities?

Omukai et al. (2005) studied the thermal and chemical evolution of collapsing gas clumps at low metallicities. They found that for models with gas metallicities of ${{10}^{-3}}-{{10}^{-4}}$ solar, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ may reach 10⁻⁸, but only if the density, n, of the gas approaches ${{10}^{8}}$ cm⁻³. Photodissociation of molecules by far-ultraviolet (FUV) radiation was not included in their study. While at solar metallicity, dust grains shield the interior of gas clouds from the FUV radiation; at low metallicities, photodissociation by FUV becomes a major removal process for H₂O. H₂O photodissociation produces OH, which is then itself photodissociated into atomic oxygen.

Hollenbach et al. (2012) developed a theoretical model to study the abundances of various molecules, including H₂O, in PDRs. Their model included many important physical processes, such as freezeout of gas species, grain-surface chemistry, and also photodissociation by FUV photons. However, they focused on solar metallicity. Intriguingly, Bayet et al. (2009) report a water abundance close to 10⁻⁸ in the optically thick core of their single PDR model for a low metallicity of 10⁻² (with $n={{10}^{3}}$ cm⁻³). However, Bayet et al. did not investigate the effects of temperature and UV intensity variations on the water abundance in the low-metallicity regime.

Recently, comprehensive studies of molecular abundances for the bulk ISM gas as functions of the metallicity were conducted by Penteado et al. (2014) and Bialy & Sternberg (2015, hereafter BS15). These models, however, focused on the "low-temperature" ion–molecule formation chemistry.

In this Letter, we present results for the H₂O abundance in low-metallicity gas environments for varying temperatures, FUV intensities, and gas densities. We demonstrate the importance of the onset of the neutral–neutral formation sequence and explore the role of the FUV field. We find that for temperatures T in the range of 250–350 K, H₂O may be abundant, comparable to or higher than that found in diffuse Galactic clouds, provided that the FUV intensity to density ratio is smaller than a critical value. In Section 2, we discuss our physical and chemical model. We present our results in Section 3. Finally, we summarize and discuss our conclusions in Section 4.

We calculate the abundance of gas-phase H₂O for low-metallicity gas parcels that are exposed to external FUV radiation and cosmic-ray and/or X-ray fluxes. Given our chemical network, we solve the steady state rate equations using our dedicated Newton-based solver, and we obtain ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ as function of T and the FUV intensity to density ratio.

We adopt a 10⁵ K diluted blackbody spectrum, representative of radiation produced by massive Pop-III stars. The photon density in the 6–13.6 eV interval is ${{n}_{\gamma }}\equiv {{n}_{\gamma ,0}}\;{{I}_{{\rm UV}}}$, where ${{n}_{\gamma ,0}}=6.5\times {{10}^{-3}}$ photons cm⁻³ is the photon density in the interstellar radiation field (ISRF; Draine 2011) and ${{I}_{{\rm UV}}}$ is the "normalized intensity." Thus, ${{I}_{{\rm UV}}}=1$ corresponds to the FUV intensity in the Draine ISRF.

Cosmic-ray and/or X-ray ionization drive the ion–molecule chemical network. We assume an ionization rate per hydrogen nucleon ζ (s⁻¹). In the Galaxy, Dalgarno (2006) and Indriolo & McCall (2012) showed that ζ lies within the relatively narrow range ${{10}^{-15}}-{{10}^{-16}}$ s⁻¹. We therefore introduce the "normalized ionization rate" ${{\zeta }_{-16}}\equiv (\zeta /{{10}^{-16}}\;{{{\rm s}}^{-1}})$. The ionization rate and the FUV intensity are likely correlated, as both should scale with the formation rate of massive OB stars. We thus set ${{\zeta }_{-16}}={{I}_{{\rm UV}}}$ as our fiducial case but also consider the cases ${{\zeta }_{-16}}={{10}^{-1}}\;{{I}_{{\rm UV}}}$ and ${{\zeta }_{-16}}=10\;{{I}_{{\rm UV}}}$.

Dust shielding against FUV radiation becomes ineffective at low metallicities. However, self-absorption in the H₂ lines may significantly attenuate the destructive Lyman Werner (11.2–13.6 eV) radiation (Draine & Bertoldi 1996; Sternberg et al. 2014), and high abundances of H₂ may be maintained even at low metallicity (BS15). In the models presented here, we assume a H₂ shielding column of at least $\sim {\mkern 1mu} {{10}^{22}}$ cm⁻². (For such conditions, CO is also shielded by the H₂.) The LW flux is then attenuated by a self-shielding factor of ${{f}_{{\rm shield}}}\sim {{10}^{-8}}$, and the H₂ photodissociation rate is only $5.8\times {{10}^{-19}}\;{{I}_{{\rm UV}}}$ s⁻¹. With this assumption, H₂ photodissociation by LW photons is negligible compared to cosmic-ray and/or X-ray ionization as long as ${{I}_{{\rm UV}}}\lt 85\;{{\zeta }_{-16}}$.

However, even when the Lyman Werner band is fully blocked, OH and H₂O are photodissociated because their energy thresholds for photodissociation are 6.4 and 6 eV, respectively. For the low metallicities that we consider, photodissociation is generally the dominant removal mechanism for H₂O and OH. We adopt the OH and H₂O photodissociation rates calculated by BS15.

We assume thermal and chemical steady states. In the Milky Way, the bulk of the ISM gas is considered to be at approximate thermal equilibrium, set by cooling and heating processes. We discuss the relevant chemical and thermal timescales in Section 2.2.

Given the aforementioned assumptions, the steady state solutions for the species abundances depend on only two parameters, the temperature T and the intensity to density ratio ${{I}_{{\rm UV}}}/{{n}_{4}}$. Here, ${{n}_{4}}\equiv (n/{{10}^{4}}\;{\rm c}{{{\rm m}}^{-3}})$ is the total number density of hydrogen nuclei normalized to typical molecular cloud densities. T and ${{I}_{{\rm UV}}}/{{n}_{4}}$ form our basic parameter space in our study.

2.1. Chemical Model

We consider a chemical network of 503 two-body reactions between 56 atomic, molecular, and ionic species of H, He, C, O, S, and Si. We assume a cosmological elemental helium abundance of 0.1 relative to hydrogen (by number). For the metal elemental abundances, we adopt the Asplund et al. (2009) photospheric solar abundances, multiplied by a metallicity factor Z' (i.e., $Z^{\prime} =1$ is solar metallicity). In our fiducial model, we assume $Z^{\prime} ={{10}^{-3}}$, but we also explore cases with $Z^{\prime} ={{10}^{-2}}$ and $Z^{\prime} ={{10}^{-4}}$. Since our focus here is on the very low metallicity regime, where dust grains play a lesser role, we neglect any depletion on dust grains and dust-grain chemistry (except for H₂, as discussed further below). Direct and induced ionizations and dissociations by the cosmic-ray/X-ray field ($\propto \zeta$) are included. For the gas-phase reactions, we adopt the rate coefficients given by the UMIST 2012 database (McElroy et al. 2013).

The formation of heavy molecules relies on molecular hydrogen. We consider two scenarios for H₂ formation: (a) pure gas-phase formation and (b) gas-phase formation plus formation on dust grains. In gas phase, radiative-attachment

$$\begin{eqnarray}&&{\rm H}+e\to {{{\rm H}}^{-}}+\nu \end{eqnarray} \tag{ 1 }$$

followed by associative-detachment

$$\begin{eqnarray}&&{{{\rm H}}^{-}}+{\rm H}\to {{{\rm H}}_{2}}+e\end{eqnarray} \tag{ 2 }$$

is the dominant H₂ formation route.

H₂ formation on the surface of dust grains is considered to be the dominant formation channel in the Milky way. We adopt a standard rate coefficient

$$\begin{eqnarray}&&R\simeq 3\times {{10}^{-17}}\ T_{2}^{1/2}\ Z^{\prime} \ {\rm c}{{{\rm m}}^{3}}\ {{{\rm s}}^{-1}}\end{eqnarray} \tag{ 3 }$$

(Hollenbach et al. 1971; Jura 1974; Cazaux & Tielens 2002), where ${{T}_{2}}\equiv (T/100\;{\rm K})$. In this expression, we assume that the dust-to-gas ratio is linearly proportional to the metallicity Z'. Thus, in scenario (b), H₂ formation on dust grains dominates even for $Z^{\prime} ={{10}^{-3}}$. Scenario (a) is the limit where the gas-phase channel dominates, as is appropriate for dust-free environments or for a superlinear dependence of the dust-to-gas ratio on Z' (see also BS15).

2.2. Timescales

The timescale for the system to achieve chemical steady state is dictated by the relatively long H₂ formation timescale. In the gas phase (scenario (a)), it is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{t}_{{{{\rm H}}_{2}}}} & = & \frac{1}{{{k}_{2}}\;n\;{{x}_{e}}}\;\approx \;8\times {{10}^{8}}\ \zeta _{-16}^{-1/2} \\ {} & {} & \times n_{4}^{-1/2}\;T_{2}^{-1}\ {\rm yr}, \\ \end{array}\end{eqnarray} \tag{ 4 }$$

where ${{x}_{e}}=2.4\times {{10}^{-5}}{{({{\zeta }_{-16}}/{{n}_{4}})}^{1/2}}T_{2}^{0.38}$ is the electron fraction as set by ionization recombination equilibrium and ${{k}_{2}}\approx 1.5\times {{10}^{-16}}T_{2}^{0.67}$ cm³ s⁻¹ is the rate coefficient for Reaction (1).

For formation on dust grains ($\propto Z^{\prime}$), the timescale is generally shorter, with

$$\begin{eqnarray}&&{{t}_{{{{\rm H}}_{2}}}}\;=\;\frac{1}{R\;n}\;\approx \;{{10}^{8}}\ T_{2}^{-1/2}\ n_{4}^{-1}\;\left( \frac{{{10}^{-3}}}{Z^{\prime} } \right)\ {\rm yr}.\end{eqnarray} \tag{ 5 }$$

Gas clouds with lifetimes of $t\gg {{t}_{{{{\rm H}}_{2}}}}$ will reach chemical steady state.

The relevant timescale for thermal equilibrium is the cooling timescale. For low-metallicity gas with $Z^{\prime} ={{10}^{-3}}$, the cooling proceeds mainly via H₂ rotational excitations (Glover & Clark 2014). If the cooling rate per H₂ molecule (in erg s⁻¹) is $W(n,T)$, then the cooling timescale is given by

$$\begin{eqnarray}&&{{t}_{{\rm cool}}}\ =\ \frac{{{k}_{{\rm B}}}\ T}{W(n,T)}.\end{eqnarray} \tag{ 6 }$$

Here, k_B is the Boltzmann constant. For $n={{10}^{4}}$ cm⁻³ and T = 300 K, $W\approx 5\times {{10}^{-25}}({{x}_{{{{\rm H}}_{2}}}}/0.1)$ erg s⁻¹ (Le Bourlot et al. 1999), and the cooling time is very short, $\approx 2\times {{10}^{3}}(0.1/{{x}_{{{{\rm H}}_{2}}}})$ yr. For densities much smaller than the critical density for H₂ cooling, $W\propto n$ and ${{t}_{{\rm cool}}}\propto 1/n$. In the opposite limit, W saturates and ${{t}_{{\rm cool}}}$ becomes independent of density. We see that generally ${{t}_{{\rm cool}}}\ll {{t}_{{{{\rm H}}_{2}}}}$.

Because the free fall time,

$$\begin{eqnarray}&&{{t}_{{\rm ff}}}\ =\ {{\left( \frac{3\pi }{32\;G\rho } \right)}^{1/2}}\ =\ 5\times {{10}^{5}}\ n_{4}^{-1/2}\ {\rm yr},\end{eqnarray} \tag{ 7 }$$

is generally much shorter than ${{t}_{{{{\rm H}}_{2}}}}$, chemical steady state may be achieved only in stable, non-collapsing clouds, with lifetimes $\gg {{t}_{{\rm ff}}}$. Obviously, both ${{t}_{{{{\rm H}}_{2}}}}$ and ${{t}_{{\rm cool}}}$ also must be shorter than the Hubble time at the redshift of interest.

Next, we present and discuss our results for the steady state, gas-phase H₂O fraction ${{x}_{{{{\rm H}}_{2}}{\rm O}}}\equiv {{n}_{{{{\rm H}}_{2}}{\rm O}}}/n$ as function of temperature T, and the FUV intensity to density ratio ${{I}_{{\rm UV}}}/{{n}_{4}}$.

3.1.  ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ as a Function of T and ${{I}_{{\rm UV}}}/{{n}_{4}}$

Figure 1 shows ${{{\rm log} }_{10}}({{x}_{{{{\rm H}}_{2}}{\rm O}}})$ contours for the two scenarios described in Section 2.1. In one, H₂ forms in the pure gas phase (scenario (a)—left panel), and in the other, H₂ also forms on dust grains (scenario (b)—right panel). Our fiducial parameters are $Z^{\prime} ={{10}^{-3}}$ and ${{\zeta }_{-16}}={{I}_{{\rm UV}}}$. In the upper left region of the parameter space, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ is generally low $\lesssim {{10}^{-9}}$. In this regime, H₂O forms through the ion–molecule sequence, that is operative at low temperatures. In the lower right corner, the neutral–neutral reactions become effective and ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ rises.

Zoom In Zoom Out Reset image size

In both panels, the solid line highlights the ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ = 10⁻⁸ contour, which resembles the H₂O gas-phase abundance in diffuse and translucent Galactic clouds. This line also delineates the borderline between the regimes where H₂O forms via the "cold" ion–molecule sequence and the "warm" neutral–neutral sequence. The temperature range at which the neutral–neutral sequence kicks in is relatively narrow, ∼250–350 K, because the neutral–neutral reactions are limited by energy barriers that introduce an exponential dependence on temperature.

The dependence on ${{I}_{{\rm UV}}}/{{n}_{4}}$ is introduced because the FUV photons photodissociate OH and H₂O molecules and therefore increase the removal rate of H₂O and at the same time suppress formation via the OH + H₂ reaction. This gives rise to a critical value for ${{I}_{{\rm UV}}}/{{n}_{4}}$ below which H₂O may become abundant.

For pure gas-phase H₂ formation (scenario (a)—left panel), the gas remains predominantly atomic, and H₂O formation is less efficient. In this case, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ $\gtrsim \;{{10}^{-8}}$ only when ${{I}_{{\rm UV}}}/{{n}_{4}}$ is smaller than a critical value of

$$\begin{eqnarray}&&\left( {{I}_{{\rm UV}}}/{{n}_{4}} \right)_{{\rm crit}}^{({\rm a})}\ =\ 2\times {{10}^{-2}}.\end{eqnarray} \tag{ 8 }$$

However, when H₂ formation on dust is included (scenario (b)—right panel), the hydrogen becomes fully molecular, and H₂O formation is then more efficient. In this case, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ may reach 10⁻⁸ for ${{I}_{{\rm UV}}}/{{n}_{4}}$ smaller than

$$\begin{eqnarray}&&\left( {{I}_{{\rm UV}}}/{{n}_{4}} \right)_{{\rm crit}}^{({\rm b})}\ =\ 3\times {{10}^{-1}},\end{eqnarray} \tag{ 9 }$$

an order of magnitude larger than for the pure gas-phase formation scenario.

3.2. Variations in Z' and ${{\zeta }_{-16}}/{{I}_{{\rm UV}}}$

In Figure 2, we investigate the effects of variations in the values of Z' and ${{\zeta }_{-16}}/{{I}_{{\rm UV}}}$. The figure shows the ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ = 10⁻⁸ contours for scenarios (a) (left panels) and (b) (right panels). As discussed above, H₂O is generally more abundant in scenario (b) because the hydrogen is fully molecular in this case, and therefore the 10⁻⁸ contours are located at higher ${{I}_{{\rm UV}}}/{{n}_{4}}$ values in the right panels.

Zoom In Zoom Out Reset image size

The upper panels show the effect of variations in the metallicity value Z' for our fiducial normalization ${{\zeta }_{-16}}={{I}_{{\rm UV}}}$. In both panels, the oxygen abundance rises, and ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ increases with increasing Z'. Thus, at higher Z', the 10⁻⁸ contours shift to lower T and higher ${{I}_{{\rm UV}}}/{{n}_{4}}$ and vice versa. An exception is the behavior of the $Z^{\prime} ={{10}^{-2}}$ curve, for which the metallicity is already high enough so that reactions with metal species dominate H₂O removal for ${{I}_{{\rm UV}}}/{{n}_{4}}\lesssim {{10}^{-2}}$. The increase in metallicity then results in a decrease of the H₂O abundance, and the 10⁻⁸ contour shifts to the right. For ${{I}_{{\rm UV}}}/{{n}_{4}}\gtrsim {{10}^{-2}}$, removal by FUV dominates, and the behavior is similar to that of the $Z^{\prime} ={{10}^{-3}}$ and $Z^{\prime} ={{10}^{-4}}$ cases.

The lower panels show the effects of the variations in the ionization rate normalization ${{\zeta }_{-16}}/{{I}_{{\rm UV}}}$ for our fiducial metallicity value of $Z^{\prime} ={{10}^{-3}}$. First, we consider the pure gas-phase formation case (scenario (a)—lower left panel). For the two cases ${{\zeta }_{-16}}/{{I}_{{\rm UV}}}=1$ and 10⁻¹, the H₂O removal is dominated by FUV photodissociation, and therefore is independent of ζ. As shown by BS15, the H₂O formation rate is also independent of ζ when the H₂ forms in the gas phase. Therefore, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ is essentially independent of ζ, and the contours overlap. For the high ionization rate ${{\zeta }_{-16}}/{{I}_{{\rm UV}}}=10$, the proton abundance becomes high, and H₂O reactions with H⁺ dominate H₂O removal. In this limit, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ decreases with ζ, and the 10⁻⁸ contour moves down.

When H₂ forms on dust (scenario (b)—lower right panel), the H₂O formation rate via the ion–molecule sequence is proportional to the H⁺ and H$_{3}^{+}$ abundances, which rise with ζ. Since the gas is molecular, the proton fraction is low, and the removal is always dominated by FUV photodissociations (independent of ζ). Therefore, in this case, ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ increases with ${{\zeta }_{-16}}/{{I}_{{\rm UV}}}$, and the curves shift up and to the left toward lower T and higher ${{I}_{{\rm UV}}}/{{n}_{4}}$.

In this Letter, we have demonstrated that the H₂O gas-phase abundance may remain high even at very low metallicities of $Z^{\prime} \sim {{10}^{-3}}$. The onset of the efficient neutral–neutral formation sequence at $T\sim 300$ K may compensate for the low metallicity and form H₂O in an abundance similar to that found in diffuse clouds in the Milky Way.

We have considered two scenarios for H₂ formation, representing two limiting cases: one in which H₂ is formed in the pure gas phase (scenario (a)), and one in which H₂ forms both in the gas phase and on dust grains, assuming that the dust abundance scales linearly with Z' (scenario (b)). Recent studies by Galametz et al. (2011), Herrera-Camus et al. (2012), and Fisher et al. (2014) suggest that the dust abundance might decrease faster than linearly with decreasing Z'. As shown by BS15, for $Z^{\prime} ={{10}^{-3}}$ and dust abundance that scales as $Z{{^{\prime} }^{\beta }}$ with $\beta \geqslant 2$, H₂ formation is dominated by the gas-phase formation channel. Therefore, our scenario (a) is also applicable for models in which dust grains are present, with an abundance that scales superlinearly with Z'. For both scenarios (a) and (b), we have found that the neutral–neutral formation channel yields ${{x}_{{{{\rm H}}_{2}}{\rm O}}}$ $\gtrsim \;{{10}^{-8}}$, provided that ${{I}_{{\rm UV}}}/{{n}_{4}}$ is smaller than a critical value. For the first scenario, we have found that this critical value is ${{({{I}_{{\rm UV}}}/{{n}_{4}})}_{{\rm crit}}}=2\times {{10}^{-2}}$. For the second scenario, ${{({{I}_{{\rm UV}}}/{{n}_{4}})}_{{\rm crit}}}=3\times {{10}^{-1}}$.

In our analysis, we have assumed that the system had reached a chemical steady state. For initially atomic (or ionized) gas, this assumption offers the best conditions for the formation of molecules. However, chemical steady state might not always be achieved within a cloud lifetime or even within the Hubble time. The timescale to achieve chemical steady state (from an initially dissociated state) is dictated by the H₂ formation process and is generally long at low metallicities. For $Z^{\prime} ={{10}^{-3}}$ and our fiducial parameters, the timescales for both scenarios are of the order of a few 10⁸ years (see, e.g., Bell et al. 2006) and are comparable to the age of the universe at redshifts of ∼10. The generically high water abundances we find for warm conditions and low metallicities will be maintained in dynamically evolving systems so long as they remain H₂ shielded.

Our results might have interesting implications for the question of how early could have life originated in the universe (Loeb 2014). Our study addresses the first step of H₂O formation in early enriched, molecular gas clouds. If such a cloud is to collapse and form a protoplanetary disk, some of the H₂O may make its way to the surface of forming planets (van Dishoeck et al. 2014, p. 835). However, the question of to what extent the H₂O molecules that were formed in the initial molecular clouds are preserved all the way through the process of planet formation is beyond the scope of this paper.

S.B. acknowledges support from the Raymond and Beverly Sackler Tel Aviv University–Harvard/ITC Astronomy Program. A.L. thanks the Sackler professorship at Tel Aviv University for generous support, as well as the NSF grant AST-1312034. This work was also supported in part by the DFG via German–Israeli Project Cooperation grant STE1869/1-1/GE625/15-1 and by a PBC Israel Science Foundation I-CORE Program grant 1829/12.