A new scale in the bias expansion

TL;DR

This paper introduces a new scale in galaxy bias expansion caused by ionizing radiation's mean free path, affecting large-scale galaxy clustering and potentially impacting cosmological measurements.

Contribution

It proposes a non-perturbative model to incorporate radiative-transfer effects into galaxy bias, enhancing the robustness of large-scale structure analyses.

Findings

Radiative transfer introduces a significant new scale in galaxy clustering.

Higher-derivative bias terms may be non-negligible due to ionizing radiation effects.

A non-perturbative model can preserve the predictive power of galaxy clustering measurements.

Abstract

The fact that the spatial nonlocality of galaxy formation is controlled by some short length scale like the Lagrangian radius is the cornerstone of the bias expansion for large-scale-structure tracers. However, the first sources of ionizing radiation between $z\approx 15$ and $z\approx 6$ are expected to have significant effects on the formation of galaxies we observe at lower redshift, at least on low-mass galaxies. These radiative-transfer effects introduce a new scale in the clustering of galaxies, i.e. the finite distance which ionizing radiation travels until it reaches a given galaxy. This mean free path can be very large, of order $100\,h^{-1}\,{\rm Mpc}$. Consequently, higher-derivative terms in the bias expansion could turn out to be non-negligible even on these scales: treating them perturbatively would lead to a massive loss in predictivity and, for example, could spoil the determination of the BAO feature or constraints on the neutrino mass. Here, we investigate under what assumptions an explicit non-perturbative model of radiative-transfer effects can maintain the robustness of large-scale galaxy clustering as a cosmological probe.

Tables

Table 1. Table 1 : This table shows which assumptions are necessary to predict the scale dependence of the bias from RT effects, assuming no inhomogeneities in the optical depth. Δ ​ η em Δ subscript 𝜂 em \Delta\eta_{\rm em} and Δ ​ η G Δ subscript 𝜂 𝐺 \Delta\eta_{G} are the interval over which ε em subscript 𝜀 em \varepsilon_{\rm em} is non-vanishing and the typical extent in time of the galaxy response, respectively. By ✓ we mean that the higher-derivative terms can be resummed into well-defined functions of k 2 ​ λ eff 2 superscript 𝑘 2 superscript subscript 𝜆 eff 2 k^{2}\lambda_{\rm eff}^{2} , each multiplied by an RT bias coefficient and an increasing power of ℋ ​ λ eff ℋ subscript 𝜆 eff {\cal H}\lambda_{\rm eff} . The case of an instantaneous galaxy response is discussed in Appendix A.4 .

$\delta \tau =0$	$\Delta {\eta }_G\ll {\mathscr{H}}^{-1}$	$\Delta {\eta }_G\sim {\mathscr{H}}^{-1}$
$\Delta {\eta }_{\mathrm{em}}\ll {\mathscr{H}}^{-1}$	✓	✓
$\Delta {\eta }_{\mathrm{em}}\sim {\mathscr{H}}^{-1}$	✓	✗

Table 2. Table 2 : Same as Tab.  1 , but taking into account the inhomogeneities in the optical depth. Note that having an inhomogeneous τ 𝜏 \tau is not a problem if the galaxy response is instantaneous: we refer to Appendix A.4 for more details.

$\delta \tau \ne 0$	$\Delta {\eta }_G\ll {\mathscr{H}}^{-1}$	$\Delta {\eta }_G\sim {\mathscr{H}}^{-1}$
$\Delta {\eta }_{\mathrm{em}}\ll {\mathscr{H}}^{-1}$	✓	✗
$\Delta {\eta }_{\mathrm{em}}\sim {\mathscr{H}}^{-1}$	✓	✗

Equations

$${\delta }_h(\eta ,\mathbf{x})=\int \mathrm{d}{\eta }^{\prime }{\mathrm{d}}^{3}y{F}_h(\eta ,{\eta }^{\prime },|\mathbf{y}|)\delta ({\eta }^{\prime },\mathbf{x}+\mathbf{y}),$$

$${\delta }_h(\eta ,\mathbf{x})=\int \mathrm{d}{\eta }^{\prime }{\mathrm{d}}^{3}y{F}_h(\eta ,{\eta }^{\prime },|\mathbf{y}|)\delta ({\eta }^{\prime },\mathbf{x})+\frac{1}{6}\int \mathrm{d}{\eta }^{\prime }{\mathrm{d}}^{3}y{F}_h(\eta ,{\eta }^{\prime },|\mathbf{y}|)|\mathbf{y}{|}^2{\nabla }^2\delta ({\eta }^{\prime },\mathbf{x})+\dots .$$

$${\delta }_h(\eta ,\mathbf{x})=b_1(\eta )\delta (\eta ,\mathbf{x})+b_{{\nabla }^2\delta }(\eta ){\nabla }^2\delta (\eta ,\mathbf{x})+\dots ,$$

$$\frac{\mathrm{D}\mathcal{I}}{\mathrm{d}\lambda }=\frac{3\mathcal{I}}{E}\frac{\mathrm{D}E}{\mathrm{d}\lambda }+\mathcal{I}{\sigma }_{\mathrm{ab}}{j}_{\mathrm{ab}}^{\mu }{P}_{\mu }-\frac{{\rho }_{\mathrm{em}}{\epsilon }_{\mathrm{em}}{U}_{\mathrm{em}}^{\mu }{P}_{\mu }}{4\pi },$$

$$E(U^{\mu }+l^{\mu }){\nabla }_{\mu }\mathcal{I}+\frac{\partial \mathcal{I}}{\partial E}\frac{\mathrm{D}E}{\mathrm{d}\lambda }=\frac{3\mathcal{I}}{E}\frac{\mathrm{D}E}{\mathrm{d}\lambda }-\mathcal{I}{\sigma }_{\mathrm{ab}}{n}_{\mathrm{ab}}E+\frac{{\rho }_{\mathrm{em}}{\epsilon }_{\mathrm{em}}E}{4\pi },$$

$$\frac{\partial \mathcal{I}}{\partial \eta }-\hat{\mathbf{n}}\cdot \nabla \mathcal{I}-\mathcal{H}E\frac{\partial \mathcal{I}}{\partial E}=-3\mathcal{H}\mathcal{I}-\mathcal{I}{\sigma }_{\mathrm{ab}}{n}_{\mathrm{ab}}a+\frac{{\rho }_{\mathrm{em}}{\epsilon }_{\mathrm{em}}a}{4\pi },$$

$$\begin{split}{\cal I}(\eta,\bm{x},E,\hat{\bm{n}})&=\bigg{(}\frac{1+z(\eta)}{1+z_{\ast}}\bigg{)}^{3}e^{-\tau(\eta,\,\bm{x},\,E,\,\hat{\bm{n}})}\,{\cal I}_{\ast}\big{(}\eta_{\ast},\bm{x}+\hat{\bm{n}}(\eta-\eta_{\ast}),E(\eta_{\ast},\eta)\big{)}\\ &=\bigg{(}\frac{1+z(\eta)}{1+z_{\ast}}\bigg{)}^{3}e^{-\tau(\eta,\,\bm{x},\,E,\,\hat{\bm{n}})}\,\frac{\Delta\eta\,a_{\ast}\,{\rho}_{\rm em}\big{(}\eta_{\ast},\bm{x}+\hat{\bm{n}}(\eta-\eta_{\ast})\big{)}\,\varepsilon_{\rm em}\big{(}\eta_{\ast},E(\eta_{\ast},\eta)\big{)}}{4\pi}\,\,,\end{split}$$

$$E({\eta }^{\prime },\eta )=E\frac{1+z({\eta }^{\prime })}{1+z(\eta )},$$

$$\tau(\eta,\bm{x},E,\hat{\bm{n}})=\int_{\eta_{\ast}}^{\eta}\mathrm{d}\eta^{\prime}\,(\sigma_{\rm ab}n_{\rm ab}a)\big{(}\eta^{\prime},\bm{x}+\hat{\bm{n}}(\eta-\eta^{\prime}),E(\eta^{\prime},\eta)\big{)}\,\,.$$

$$\begin{split}&\delta{n}_{g}|_{\rm ion}(\eta,\bm{x})=\int_{\eta_{\ast}}^{\eta}\mathrm{d}\eta^{\prime}\int\mathrm{d}\hat{\bm{n}}\int_{0}^{+\infty}\mathrm{d}E\,{G}^{(1)}_{g}(\eta,\eta^{\prime},E,\hat{\bm{n}})\,\delta{{\cal I}}\big{(}\eta^{\prime},\bm{x}_{\rm fl}(\eta^{\prime}),E,\hat{\bm{n}}\big{)}\\ &\hphantom{\delta{n}_{g}|_{\rm ion}(\eta,\bm{x})}=\int_{\eta_{\ast}}^{\eta}\mathrm{d}\eta^{\prime}\int\mathrm{d}\hat{\bm{n}}\int_{0}^{+\infty}\mathrm{d}E\,{G}^{(1)}_{g}(\eta,\eta^{\prime},E,\hat{\bm{n}})\,\delta{{\cal I}}(\eta^{\prime},\bm{x},E,\hat{\bm{n}})\,\,,\end{split}$$

$$G_g^{(i)}(\eta ,{\eta }^{\prime },E)=\mathcal{G}{\sigma }_{\mathrm{bf}}(E)G_g^{(i)}(\eta ,{\eta }^{\prime })\text{for i=0,1},$$

$$G^{(i)}_{g}(\eta,\eta^{\prime})\to\bigg{(}\frac{1+z_{\ast}}{1+z(\eta^{\prime})}\bigg{)}^{3}\,{G}^{(i)}_{g}(\eta,\eta^{\prime})\quad\text{for $i=0,1$}\,\,.$$

$$\begin{split}&\delta n_{g}|_{\rm ion}(\eta,\bm{x})=\frac{\Delta\eta\,a_{\ast}\,{\cal G}}{4\pi}\int_{\lvert\bm{y}\rvert\,\leq\,\eta\,-\,\eta_{\ast}}\mathrm{d}^{3}y\,\frac{{G}^{(1)}_{g}(\eta,\eta_{\ast}+\lvert\bm{y}\rvert)}{\lvert\bm{y}\rvert^{2}}\,\delta\rho_{\rm em}(\eta_{\ast},\bm{x}+\bm{y})\,\times\\ &\hphantom{\delta n_{g}|_{\rm ion}(\eta,\bm{x})=\frac{\Delta\eta\,a_{\ast}\,{\cal G}}{4\pi}}\int_{0}^{+\infty}\mathrm{d}E\,\sigma_{\rm bf}(E)\,\varepsilon_{\rm em}\big{(}\eta_{\ast},E(\eta_{\ast},\eta_{\ast}+\lvert\bm{y}\rvert)\big{)}\,{e^{-\hat{\tau}(\eta,\,\bm{y},\,E)}}\,\,,\end{split}$$

$$\hat{\tau}(\eta,\bm{y},E)=\hat{\tau}(\eta,\lvert\bm{y}\rvert,E)=\lvert\bm{y}\rvert\int_{0}^{1}\mathrm{d}u\,(n_{\text{HI}}a)(\eta_{\ast}+u\lvert\bm{y}\rvert)\,\sigma_{\rm bf}\big{(}E(\eta_{\ast}+u\lvert\bm{y}\rvert,\eta_{\ast}+\lvert\bm{y}\rvert)\big{)}\,\,.$$

$$\sigma_{\rm bf}(E)=\frac{64\pi\alpha^{3}}{3E_{\infty}^{2}}\bigg{(}\frac{E_{\infty}}{E}\bigg{)}^{\frac{7}{2}}\equiv\sigma_{\infty}\bigg{(}\frac{E_{\infty}}{E}\bigg{)}^{\frac{7}{2}}\,\,,$$

$$\begin{split}\hat{\tau}(\eta,\lvert\bm{y}\rvert,E)&=\lvert\bm{y}\rvert\int_{0}^{1}\mathrm{d}u\,(n_{\text{HI}}a)(\eta_{\ast}+u\lvert\bm{y}\rvert)\,\sigma_{\rm bf}\big{(}E(\eta_{\ast}+u\lvert\bm{y}\rvert,\eta_{\ast}+\lvert\bm{y}\rvert)\big{)}\\ &\approx\lvert\bm{y}\rvert\sigma_{\rm bf}\big{(}E(\eta_{\ast},\eta_{\ast})\big{)}(n_{\text{HI}}a)(\eta_{\ast})=\lvert\bm{y}\rvert\sigma_{\rm bf}(E)(n_{\text{HI}}a)(\eta_{\ast})\equiv\frac{\lvert\bm{y}\rvert}{\lambda_{\rm ion}(\eta_{\ast},E)}\,\,.\end{split}$$

$$\begin{split}&\delta n_{g}|_{\rm ion}(\eta,\bm{x})=\frac{\Delta\eta\,a_{\ast}\,{\cal G}\,\sigma_{\infty}\,\varepsilon_{\infty}}{4\pi}\int_{\lvert\bm{y}\rvert\,\leq\,\eta\,-\,\eta_{\ast}}\mathrm{d}^{3}y\,\frac{{G}^{(1)}_{g}(\eta,\eta_{\ast}+\lvert\bm{y}\rvert)}{\lvert\bm{y}\rvert^{2}}\,\bigg{(}\frac{1+z_{\ast}}{1+z(\eta_{\ast}+\lvert\bm{y}\rvert)}\bigg{)}^{s}\,\times\\ &\hphantom{\delta n_{g}|_{\rm ion}(\eta,\bm{x})=\frac{\Delta\eta\,a_{\ast}\,{\cal G}\,\sigma_{\infty}\,\varepsilon_{\infty}}{4\pi}}\,\delta\rho_{\rm em}(\eta_{\ast},\bm{x}+\bm{y})\int_{E_{\infty}}^{+\infty}\mathrm{d}E\,\bigg{(}\frac{E}{E_{\infty}}\bigg{)}^{s\,-\,\frac{7}{2}}\,{e^{-\frac{\lvert\bm{y}\rvert}{\lambda_{\rm ion}(\eta_{\ast},\,E)}}}\,\,.\end{split}$$

$$\int_{E_{\infty}}^{+\infty}\mathrm{d}E\,\bigg{(}\frac{E}{E_{\infty}}\bigg{)}^{s\,-\,\frac{7}{2}}\,{e^{-\frac{\lvert\bm{y}\rvert}{\lambda_{\rm ion}(\eta_{\ast},\,E)}}}\approx\frac{2E_{\infty}}{5-2s}\,e^{{-\frac{2s\,-\,5}{2(s\,-\,6)}}\lvert\bm{y}\rvert\sigma_{\infty}(n_{\text{HI}}a)(\eta_{\ast})}\,\,,$$

$$\begin{array}{rl}\delta n_g{|}_{\mathrm{ion}}(\eta ,\mathbf{x})& =\underset{\equiv {\mathcal{C}}_{\mathrm{ion}}}{\underbrace{\frac{\Delta \eta a_{\ast }\mathcal{G}{\sigma }_{\infty }{\epsilon }_{\infty }{E}_{\infty }}{2\pi (5-2s)}}}\int {\mathrm{d}}^{3}y\frac{G_{\mathrm{g}}^{(1)}(\eta ,{\eta }_{\ast }+|\mathbf{y}|)}{|\mathbf{y}{|}^2}\delta {\rho }_{\mathrm{em}}({\eta }_{\ast },\mathbf{x}+\mathbf{y})e^{-\frac{|\mathbf{y}|}{{\lambda }_{\mathrm{eff}}({\eta }_{\ast })}}.\end{array}$$

$$G_g^{(1)}(\eta ,{\eta }^{\prime })=\sum _{n=0}^{+\infty }{g}_n(\eta ){\mathcal{H}}^n(\eta )({\eta }^{\prime }-{\eta }_{\ast }{)}^n,$$

$${\delta }_g{|}_{\mathrm{ion}}(\eta ,\mathbf{k})=f_{\mathrm{ion}}(\eta ,k^2){\delta }_{\mathrm{em}}({\eta }_{\ast },\mathbf{k}),$$

$$\int {\mathrm{d}}^{3}y{e}^{-i\mathbf{k}\cdot \mathbf{y}}|\mathbf{y}{|}^{n-2}{e}^{-\frac{|\mathbf{y}|}{{\lambda }_{\mathrm{eff}}}}\equiv {\lambda }_{\mathrm{eff}}^{n+1}(4\pi n!){\mathcal{F}}_{\mathrm{ion}}^{(n)}(k^2{\lambda }_{\mathrm{eff}}^2),$$

$$f_{\mathrm{ion}}(\eta ,k^2)=\sum _{n=0}^{+\infty }{g}_n(\eta )(\mathcal{H}{\lambda }_{\mathrm{eff}}{)}^n{\mathcal{F}}_{\mathrm{ion}}^{(n)}(k^2{\lambda }_{\mathrm{eff}}^2),$$

$${\cal F}^{(n)}_{\rm ion}(x^{2})=\frac{\sin\!\big{(}n\arctan(x)\big{)}}{nx(1+x^{2})^{n/2}}\,\,,$$

$${\delta }_{\mathrm{em}}(\eta ,\mathbf{k})=b_{\mathrm{em}}(\eta )\delta (\eta ,\mathbf{k})+{ϵ}_{\mathrm{em}}(\eta ,\mathbf{k}),$$

$${⟨{ϵ}_{\mathrm{em}}(\eta ,\mathbf{k}){ϵ}_{\mathrm{em}}({\eta }^{\prime },{\mathbf{k}}^{\prime })⟩}^{\prime }=P_{{ϵ}_{\mathrm{em}}}^{\{0\}}(\eta ,{\eta }^{\prime }).$$

Full text

1 Introduction

A crucial step in connecting cosmological scenarios with large-scale-structure (LSS) observations is the bias expansion for LSS tracers (see [1] for a review). Like all effective field theories, the bias expansion is firmly rooted in the idea of locality. The simplest example is that of dark matter halos. Since their dynamics is governed by gravitational interactions only, we know that at any given point in spacetime the halo overdensity $\delta_{h}$ will be a functional of all possible gravitational observables that one can build from the Newtonian potential. Let us focus, for example, on the total matter overdensity $\delta$. At leading order in perturbation theory, we can write the most general functional dependence of the halo overdensity on $\delta$ as the integral

[TABLE]

where the kernel $F_{h}$ can only depend on $\lvert\bm{y}\rvert$ because of statistical isotropy, and cannot depend on $\bm{x}$ because of statistical homogeneity. Even if we do not know the exact shape of the kernel, we know that it is supported only for $\lvert\bm{y}\rvert\lesssim R(M_{h})$, where $R(M_{h})$ is the Lagrangian radius of the halo. This is because the matter within a given halo originates from a region of size $R(M_{h})$ in Lagrangian space. We can then expand $\delta(\eta^{\prime},\bm{x}+\bm{y})$ in a Taylor series inside the integral, and obtain

[TABLE]

Finally, using the fact that, at linear order in perturbations, $\delta$ evolves in a scale-independent way, we can easily carry out the time integration to arrive at the expansion

[TABLE]

where the bias coefficients $b_{1},b_{\nabla^{2}\delta},\dots$, are related to the moments of the kernel $F_{h}$ [2]. Since the nonlocality scale of the kernel is $R(M_{h})$, we see that $b_{\nabla^{2n}\delta}\sim R^{2n}(M_{h})$.

So far, so good: everything seems to be controlled by the single scale $R(M_{h})$. Moreover, this scale is typically of the same order of magnitude as the scale at which the matter density field becomes nonlinear, so that it does not strongly restrict the validity of the perturbative bias expansion. When we move from halos to galaxies and consider their overdensity $\delta_{g}$, we can still write it as in Eq. (1.3). However, now we have to ask what is the nonlocality scale of the higher-derivative terms for galaxies. If the properties of galaxies in a given sample are completely determined by those of their host halos, then this scale is still $R(M_{h})$. However, the real universe is not so simple. For example, baryons are also present. What is their impact on Eq. (1.3)? Pressure forces contribute to the right-hand side of the Euler equation for baryonic matter through a pressure gradient $\bm{\nabla}\delta p_{b}=c^{2}_{\rm s}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{b}\bm{\nabla}\delta_{b}\approx c^{2}_{\rm s}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{b}\bm{\nabla}\delta$ (approximating $\delta_{b}\approx\delta$ on large scales), and then give rise to a baryon-dark matter relative velocity $\bm{v}_{bc}$. Since this relative velocity is a local observable, it can enter in the bias expansion. At leading order in perturbations we can then add the term $\bm{\nabla}\cdot\bm{v}_{bc}\sim\nabla^{2}\delta$ to Eq. (1.3). Thus, we see that these baryonic effects are also captured by higher-derivative terms. However, the length scale suppressing them is not $R(M_{h})$, but the Jeans length $\lambda_{\rm J}=c_{\rm s}(G_{\rm N}\rho)^{-1/2}$, which depends on the average density of gas in the halo and its speed of sound $c^{2}_{\rm s}\sim T/m$ (for gas particles of mass $m$). Fortunately, this length scale is again quite small, and less than $R(M_{h})$ except for very low-mass halos.

Another source of nonlocality is radiative transfer (RT). Indeed, ionizing radiation can affect star-forming galaxies directly: for example, Refs. [3, 4] showed that it can reduce the cooling rate of the gas accreting onto the parent halo (effectively evaporating it from the halo), slowing down the star-formation rate and leading to a suppression of the stellar-mass-to-halo-mass ratio in low-mass galaxies. In this case the nonlocality scale is some “effective” mean free path (m.f.p.) $\lambda_{\rm eff}$ of ionizing radiation, that in the following will always be understood as the comoving one (we use the term “effective” since we will see that what matters for the bias expansion is some average of the mean free path over the photon energy).

Both these effects are expected to become relevant during reionization, when the progenitors of galaxies observed at lower redshifts were actively forming [5]. It is during this epoch (between redshifts $z\approx 15$ and $z\approx 6$) that newly formed radiating objects like metal-poor stars, supernova explosions, accreting black holes, X-ray binaries, mini-quasars, dwarf galaxies, etc. injected photons into the intergalactic medium, and the universe reverted from neutral to ionized. The gas temperature jumped from a few $\rm K$ to several thousands, leading to a Jeans length $\lambda_{\rm J}\approx$0.08${\,h^{-1}\,\mathrm{Mpc}}$ around the end of reionization [5], while the m.f.p. can be of order $$50${\,h^{-1}\,\mathrm{Mpc}}$ at redshift $z\sim 5\leavevmode\nobreak\ {\text{--}}\leavevmode\nobreak\ 6$ [6, 7].

Once we have the expansion of Eq. (1.3), and similar expansions at higher order in perturbation theory, we can predict the statistics of galaxies in terms of those of the underlying matter field at a given order in perturbation theory. In the best case, this prediction holds up to the nonlinear scale, $k_{\rm NL}\approx$0.3${\,h\,\mathrm{Mpc}^{-1}}$ at $z=0$, that controls the perturbation theory for the dark matter distribution. However, it is clear that we must be able to truncate the series in Eq. (1.3) at some order if we want to be predictive. This is easy to see in Fourier space. Let us consider for example the higher-derivative terms coming from RT effects. In Eq. (1.3) we obtain a series of terms that scale as $(k^{2}\lambda_{\rm eff}^{2})^{n}$ times a coefficient that, naively, we expect to be of order $1$: then, the derivative expansion will break down at the very latest at $k\sim 1/\lambda_{\rm eff}$, since at that scale infinitely many terms will become equally relevant. In this paper we will focus on these terms, answering the following questions. Is it possible, within the effective field theory framework, to resum them in a way that allows us to predict the galaxy statistics also at momenta $k\gtrsim 1/\lambda_{\rm eff}$ (but still $k\ll k_{\rm NL}(z)$)? Can we do this with only a finite number of new bias coefficients? Which assumptions are necessary to achieve this?

These questions are important because $\lambda_{\rm eff}$ can be much larger than $\lambda_{\rm J}$: while there is no problem in treating the higher-derivative terms from pressure forces perturbatively (unless their bias coefficients are unnaturally large), doing so with those from RT can lead to a sizeable loss in predictivity. How does this happen? When we marginalize over the free coefficients of the bias expansion, we need to take into account their prior. The coefficients $b_{\nabla^{2n}\delta}$ of the higher-derivative terms are dimensionful: however, once we identify the longest non-locality scale $\lambda$ (be it $R(M_{h})$, or $\lambda_{\rm eff}$, etc.) that we think affects the formation of the tracers under consideration, it is reasonable to assume that after we factor it out the dimensionless coefficients are of order $1$, as we discussed below Eq. (1.3).111Therefore, we would take the priors on the coefficients $b_{\nabla^{2n}\delta}$ to be $[{-{\cal O}}(1),{\cal O}(1)]\times\lambda^{2n}$, where $\lambda^{2n}$ is fixed. Consequently, if we want to keep a finite number of higher-derivative terms (and then a finite number of coefficients to marginalize over), we are forced to stop at a $k_{\rm max}$ which is not larger than $\sim 1/\lambda$. The longer $\lambda$, then, the smaller is the number of modes that we are using, and this would lead to a degradation of the constraints on cosmological parameters (for more details we refer, for example, to the discussions in Sections 3 and 4 (especially Section 4.2) of [8]). More precisely, these higher-derivative terms will strongly modify the shape of the galaxy two-point correlation function around BAO scales, damping the amplitude of the BAO feature and possibly affecting its measured position (see e.g. [9, 10, 11, 12] for a discussion). Finally, they will also have a negative impact on the constraints on parameters like the neutrino mass and the amplitude of equilateral primordial non-Gaussianity: indeed, the scale-dependent effects induced by neutrinos and equilateral non-Gaussianity are controlled by the free-streaming scale $k_{\rm fs}$ and the equality scale $k_{\rm eq}$ respectively, and both these scales could be close to $1/\lambda_{\rm eff}$.

Before proceeding, we emphasize that these RT effects also affect line emission from diffuse gas, like the Lyman-$\alpha$ forest and $\rm 21cm$ intensity mapping (see e.g. [11, 12]). Actually, we expect that the response of galaxies to ionizing radiation is very suppressed if they reside in halos with mass much larger than the Jeans mass at reionization [5], while the line emission depends strongly on the ionization state of the medium and then on the ambient radiation field. This distinction is not important for this paper: since we follow an effective field theory approach, all our conclusions will apply to any physical tracer of the matter field (the differences between various tracers being encoded in their respective bias coefficients).

Our paper is subdivided in four main sections. We formulate in more detail the above questions in Section 2, and we answer them in Sections 3 and 4. We draw our conclusions in Section 5. Technical details on Sections 3 and 4 are collected in Appendix A.

Notation and conventions

We largely follow the notation of [1]. Some differences are that we denote conformal time by $\eta$ ($\tau$ is reserved for the optical depth), and denote the leading stochastic term in the bias expansion for galaxies by $\epsilon_{g}$ (instead of simply $\epsilon$), since we will be dealing with multiple tracers. For simplicity we omit the dependence of the fluid trajectory $\bm{x}_{\rm fl}$ on the Lagrangian coordinate. The metric signature is $(-,+,+,+)$. We work in natural units $c=\hbar=\varepsilon_{0}=1$, where $\varepsilon_{0}$ is the permittivity of free space. Consequently, we freely exchange energy with (angular) frequency and time with length in our discussion of radiative transfer.

2 Setting up the analysis

As one can easily imagine, we can capture the effect of RT on galaxy formation by allowing their number density to be a functional of the incoming flux of ionizing radiation [5]. In principle, galaxy formation is sensitive to this flux in a region of finite size, corresponding to the extent of the gas cloud around the parent halo, which is of order of the Jeans length of the gas. This nonlocality can then be dealt with via a derivative expansion like that of Eq. (1.1), and at leading order it is enough to consider the flux evaluated along the trajectory $\bm{x}_{\rm fl}$ of a Lagrangian patch enclosing the tracer: higher-order corrections involve derivatives of the incoming flux along this fluid trajectory. The additional scales of nonlocality in the problem are the following:

Mean free path of radiation

As we discussed in the introduction, one is $\lambda_{\rm eff}$. Radiation can travel long distances before reaching the galaxy (while in contrast matter and biased tracers typically move $\lesssim\pi/k_{\rm NL}(z=0)\sim$10${\,h^{-1}\,\mathrm{Mpc}}$ over the entire history of the universe). Therefore, since the emissivity of the sources of photoionizing radiation is also biased with respect to matter, we can conclude that the galaxies are now sensitive to the distribution of matter within their whole past light cone. However, radiation is also being absorbed by the intergalactic medium, and then sources farther than $\lambda_{\rm eff}$ from a given galaxy are not able to influence it. This is why we can assume that higher-derivative operators from RT effects are suppressed by $1/\lambda_{\rm eff}$, i.e. $b_{\nabla^{2n}\delta}\sim f_{(2n)}\lambda^{2n}_{\rm eff}$. With this assumption on the scaling of higher-derivative bias coefficients, Ref. [5] obtained a constraint on the bias coefficient of the operator $\nabla^{2}\delta$ for the galaxies of the BOSS DR12 galaxy sample [13, 14]: $\lvert f_{(2)}\rvert\lesssim$0.002$\,(\lambda_{\rm eff}/$50${\,h^{-1}\,\mathrm{Mpc}})^{-2}$ (the smallness of $f_{(2)}$ is consistent with the fact that the galaxies of the BOSS sample reside in massive halos, $M_{h}\approx${10}^{13}${\,h^{-1}M_{\odot}}$).

Response history to incoming flux

It is possible that galaxy formation does not respond instantaneously to the incoming radiation, but keeps some “memory” of it. To see the consequences of this, imagine for a moment that all ionizing radiation is emitted in an interval $\Delta\eta$ around some redshift $z_{\ast}$, with $\Delta\eta\ll 1/{\cal H}(z_{\ast})$. If the response of galaxies to the incoming flux of ionizing radiation was instantaneous, the galaxy number density at any event along the fluid worldline would only depend on the distribution of sources at the intersection of the past light cone of that event with the hypersurface $z=z_{\ast}$ (see Fig. 2 below). In general however, it also knows about sources at $z_{\ast}$ that are inside the past light cone, closer to the spatial position of the galaxies under consideration by an amount controlled by the response time. Clearly this effect goes in the opposite direction with respect to that of a finite mean free path of radiation, and must also be accounted for.

Before discussing how to account for both these scales in our effective field theory approach, let us also briefly think about the alternative to this. Obviously, one could attempt a direct modeling of what would be called “UV physics” in the language of particle physics. This means both a model of the response of galaxies to the ionizing photons, but also of the spatial distribution of the radiation sources within the light cone of the observed galaxies. However, it is clear that we cannot directly test our models for the latter against observations, since the past light cone of the observed galaxies lies within (and not on) our past light cone: while in principle we could imagine to reconstruct the response of galaxies with observations at different redshifts (since by statistical homogeneity the response cannot depend on the spatial position of the galaxies), we cannot receive the light from all the sources inside it. This is easy to see in an FLRW universe (see for example Figs. 1, 2 in the next section) or any conformally flat spacetime, but holds non-perturbatively.222In general, consider an event $P$ and an event $Q$ in the causal past of $P$. Then, any event $R$ in the causal past of $Q$ belongs to the causal past of $P$, since we can always find a timelike or lightlike curve connecting $P$ with $R$. However, it is clear that there are many events in the causal past of $Q$ that are not on the light cone of $P$, like for example any event that is reached by a past-directed timelike curve from $Q$.

Instead of pursuing this route, we attempt to remain as general as possible and only assume that the sources of ionizing radiation can be described by a bias expansion whose nonlocality scale is much shorter than the m.f.p. of ionizing radiation. Further, since in this paper we are only interested in treating the higher-derivative terms controlled by $\lambda_{\rm eff}$, we will drop all higher-derivative terms that are suppressed by these additional short scales and stop at first order in perturbations. The extension to higher orders in perturbations is well understood when spatial nonlocality is expanded following Eq. (1.2) (see Sections 2 and 4 of [1]). Moreover, we expect RT effects to typically be of relatively small amplitude, so that a linear treatment of their contribution is appropriate.

3 Radiative-transfer effects and the bias expansion

We start by presenting the equations of radiative transfer, and discuss in more detail the properties of the emission and absorption coefficients (Section 3.1).

Then, in order to set the stage, in Section 3.2 we consider a “flash” of radiation emitted around $\eta=\eta_{\ast}$, i.e. emitted within an interval $\Delta\eta$ much shorter than a Hubble time, and focus only on inhomogeneities in the radiation field on this hypersurface (i.e. we neglect the effect of inhomogeneities along the photon geodesics leading to the galaxies). We show that, in this idealized scenario, we can capture the RT effects without expanding in derivatives by adding one new function of $k^{2}$ in the bias expansion: this function encodes the response of the galaxies to the incoming flux of ionizing radiation along the past fluid worldline. We also show that, even if we do not know this response exactly, it is possible to predict the shape of this function at every order in an expansion in powers of the radiation m.f.p. divided by the Hubble radius. That is, as long as the m.f.p. is significantly smaller than the horizon, we can successfully resum the power series in $k^{2}\lambda_{\rm eff}^{2}$ into a computable function. Our ignorance of the small-scale physics is then parameterized by the overall coefficients of this expansion, which are generic functions of time.

This statement is not trivial. It is well known that, in the standard bias expansion, all the Green’s functions that describe the nonlocality in time of galaxy formation can, at each order in perturbation theory, be rewritten in terms of a finite number of time-dependent (but not scale-dependent) bias coefficients. At linear order, this is shown explicitly in Section 1 (compare Eq. (1.2) with Eq. (1.3)). In our case, instead, any time dependence along the past fluid trajectory turns into a dependence on the spatial distribution of sources at $\eta=\eta_{\ast}$ (since the integrated incoming flux receives contributions from everything along the past light cone of the galaxies), and this leads to a modification of all the higher-derivative terms in the bias expansion.

However, as long as the m.f.p. is much shorter than the Hubble radius, photons are actually arriving from a small comoving volume of size $\sim\lambda_{\rm eff}^{3}$ around the past fluid trajectory. This is what ultimately allows a resummation of all these higher-derivative terms into specific functions of $k^{2}$, each proportional to an increasing power of ${\cal H}\lambda_{\rm eff}$. Crucially, this resummation allows us to describe galaxy clustering even for $k\gg 1/\lambda_{\rm eff}$. We will later see whether this also holds beyond the instantaneous flash and homogeneous medium approximations.

We compute the impact of the resummed RT effects on the galaxy power spectrum and possible degeneracies with the neutrino mass in Section 3.3. In Sections 3.4 and 3.5 we discuss briefly how to treat tracers that are immune to the very nonlocal effects discussed above (we label these tracers by “$q$” in the following). This provides an illustration for how multiple LSS tracers might help in understanding the RT contributions.

3.1 Radiative-transfer equation

Let us first write down the equation of radiative transfer. We define the phase-space density of emitted radiation as $\mathcal{N}(x^{\mu},P^{\nu})$, where $P^{\mu}$ is the photon four-momentum. Calling $U^{\mu}$ the four-velocity of the observer that follows the fluid worldline $x^{\mu}_{\rm fl}=(\eta,\bm{x}_{\rm fl}(\eta))$, we can decompose $P^{\mu}$ as $E(U^{\mu}+l^{\mu})$, where $E$ and $l^{\mu}$ are, respectively, the photon energy and its direction as measured by the observer defined by $U^{\mu}$. Consequently, we can define the specific intensity of emitted radiation for the observer $U^{\mu}$ as ${\cal I}=E^{3}\mathcal{N}$. Its evolution is dictated by the Boltzmann equation (see e.g. [15])

[TABLE]

where $\lambda$ is an affine parameter along the photon geodesics, normalized such that $P^{\mu}={\mathrm{d}x^{\mu}}/{\mathrm{d}\lambda}$, and the term $3{\cal I}\times({\rm D}\log E/{\rm d}\lambda)$ encodes, e.g., the dilution due to the expansion of the universe. For simplicity of notation we have suppressed the argument of $\cal I$: in general, like $\cal N$, it is also a function of the spacetime position $x^{\mu}$ and the photon four-momentum $P^{\mu}$.

Let us discuss in more detail the second and the third term on the right-hand side of Eq. (3.1), which are directly related to the so-called emission and absorption coefficients. For simplicity, here we have considered a single family of emitter and absorbers, since the generalization to an arbitrary number is straightforward:

• •

the number current of absorbers (e.g. neutral hydrogen) is given by $j^{\mu}_{\rm ab}=n_{\rm ab}U^{\mu}_{\rm ab}$, where $U^{\mu}_{\rm ab}$ is their four-velocity and $n_{\rm ab}$ their number density;

• •

the absorption cross section is called $\sigma_{\rm ab}$. In the case of absorption by neutral hydrogen with emission of an electron in the continuum (photoionization), this would be the bound-free cross section $\sigma_{\rm bf}$;

• •

$U^{\mu}_{\rm em}$ is the four-velocity of the emitting medium and $\rho_{\rm em}$ its mass density. For emitters of fixed mass $m_{\rm em}$ and number current $j^{\mu}_{\rm em}$ we have $\rho_{\rm em}U^{\mu}_{\rm em}=m_{\rm em}j^{\mu}_{\rm em}$,

• •

the dimensionless function $\varepsilon_{\rm em}(x^{\mu},P^{\nu})$ is the emissivity, defined as the energy emitted by the source medium per unit frequency per unit time per unit mass [15].

With this, the emission and absorption coefficients, defined as the energy emitted per unit time per unit solid angle per unit volume and the cross-sectional area presented by absorbers per unit volume, are $(\rho_{\rm em}\varepsilon_{\rm em})/(4\pi)$ and $\sigma_{\rm ab}n_{\rm ab}$, respectively.

What about the absorption cross section and the emissivity? There is no preferred direction towards which the constituents of the source medium can align (at first order in perturbations): therefore the emissivity does not depend on the photon direction $l^{\mu}$. Moreover, here we are also considering the total absorption cross section, so that $\sigma_{\rm ab}$ depends only on the photon energy.333Consider for example the photoionization of hydrogen in the $1s$ state. If the emitted electron is nonrelativistic, the angular dependence of the differential cross section is given by $\lvert\hat{\bm{k}}_{e}\cdot\hat{\bm{\epsilon}}\rvert^{2}$, where $\bm{k}_{e}$ is the momentum of the outgoing electron and $\hat{\bm{\epsilon}}$ is the polarization of the incoming photon. However, in this case we are not interested in the direction of the outgoing electron, so we integrate over it (spin is conserved for nonrelativistic electrons, so tracing over it is trivial). Further, we average over $\hat{\bm{\epsilon}}$ since we consider unpolarized radiation. We emphasize that these assumptions are common to both analytic studies of reionization [16] and radiative-transfer simulations (see e.g. the discussion in [17]). Finally, $\varepsilon_{\rm em}$ does not depend on the fluid trajectory $\bm{x}_{\rm fl}(\eta)$ but only on time, since the emission of radiation is localized at the position of the sources. For some examples of specific source models we refer to [18, 19, 16].

Before proceeding, we also emphasize that in Eq. (3.1) we have neglected scattering. We discuss this issue in more detail in Section 4. While Eq. (3.1) holds in general spacetimes, for the purposes of this paper it is sufficient to consider an FLRW universe, and neglect metric perturbations. More precisely, lensing will be relevant only at second order in perturbations (since a nontrivial lensing effect only arises if there are anisotropies in the radiation field), and gravitational redshift effects (Sachs-Wolfe and Integrated Sachs-Wolfe) can be neglected on sufficiently sub-horizon scales.

Let us now take advantage of the fact that, if we neglect higher-derivative terms suppressed by the halo Lagrangian radius, all of emitters, absorbers, and receivers are comoving with the matter fluid, i.e. there is no velocity bias [20, 21, 1]. Then, we have that $j^{\mu}_{\rm ab}=n_{\rm ab}U^{\mu}$ and $U^{\mu}_{\rm em}=U^{\mu}$. This tells us that $j^{\mu}_{\rm ab}P_{\mu}={-n_{\rm ab}}E$, and ${-P_{\mu}U_{\rm em}^{\mu}}=E$. Therefore Eq. (3.1) becomes

[TABLE]

where we have expanded ${\rm D}/{\rm d}\lambda$ using the fact that ${\cal I}$ does not depend explicitly on $l^{\mu}$, given our assumptions on $\sigma_{\rm ab}$ and $\varepsilon_{\rm em}$. Then, for an FLRW metric in comoving coordinates we have $l^{\mu}=(0,-\hat{\bm{n}}/a)$, where $\hat{\bm{n}}$ is opposite to the photon direction: $\hat{\bm{n}}$ remains constant because there is no lensing, and ${\rm D}E/{\rm d}\lambda=-HE^{2}$. Moreover, since we stop at first order in perturbations in the bias expansion, it is not necessary to consider displacements [20, 21, 1]: that is, the fluid worldline is just given by $x^{\mu}_{\rm fl}=(\eta,\bm{x}_{\rm fl}(\eta))=(\eta,\bm{x})$ and then $U^{\mu}=\delta^{\mu}_{0}/a$. Eq. (3.2) then becomes

[TABLE]

where now and in the following we write ${\cal I}$ in terms of the arguments ${\cal I}(\eta,\bm{x},E,\hat{\bm{n}})$.444Given our assumptions of isotropic emissivity and no scattering the dependence on $\hat{\bm{n}}$ comes only from the photon free-streaming after emission: we refer to Eq. (A.1) for details.

It is straightforward to solve Eq. (3.3) by an integral along the line of sight. This is done in Appendix A.1. Before discussing the most general solution (which will be done in Section 4), let us focus on the case where the emissivity is non-vanishing only for a short interval of time around $\eta_{\ast}$, and there are no inhomogeneities in the number density of absorbers ($n_{\rm ab}(\eta,\bm{x})=\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\rm ab}(\eta)$). Clearly, the time evolution of the emissivity and the fluctuations in $n_{\rm ab}$ must also be considered. However, these assumptions allow us to introduce in a simple way the response of galaxies to the ionizing radiation. Moreover, as discussed above, we will show that it is in this scenario that we are able to predict how the correction to the galaxy bias from radiative transfer depends on $k$.

3.2 Instantaneous flash and homogeneous optical depth

The spacetime diagram that summarizes our setup is shown in Figs. 1, 2. Radiation is emitted in a single flash around $\eta=\eta_{\ast}$, and can affect the galaxies over the past fluid worldline (blue line in Figs. 1, 2).

At leading order in $\Delta\eta\ll 1/{\cal H}(\eta_{\ast})$, and in absence of scattering and of additional sources after $\eta_{\ast}$, the radiation intensity received by the galaxies at some event $(\eta,\bm{x})$ on $x^{\mu}_{\rm fl}$ takes the form (see Appendix A.2 for a derivation)

[TABLE]

where the emitted intensity ${\cal I}_{\ast}$ is equal to $\Delta\eta(\rho_{\rm em}\varepsilon_{\rm em}a)|_{\eta=\eta_{\ast}}/(4\pi)$ (and is independent of $\hat{\bm{n}}$ given that $\varepsilon_{\rm em}$ is isotropic), and the redshift dependence of the photon energy is, for general $\eta^{\prime}\leq\eta$, given by

[TABLE]

so that $E$ is the energy measured by an observer at $(\eta,\bm{x})$. Given the number density of absorbers and the absorption cross section, the optical depth $\tau$ in Eq. (3.4) is equal to

[TABLE]

From now on we take $\sigma_{\rm ab}$ to be the bound-free cross section $\sigma_{\rm bf}$ for photoionization of the $1s$ state, and $n_{\rm ab}=n_{\text{HI}}$ to be the corresponding number density of neutral hydrogen in the $1s$ state. We will discuss what happens if we consider the full cross section and the more general scenario of multiple absorbers (like higher levels of hydrogen or helium, for example) in Section 3.6. Moreover, since in this section we assume the absorbers to be homogeneous, $n_{\text{HI}}=\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}$, $\tau$ does not depend on $\hat{\bm{n}}$. Also, for simplicity of notation we will drop the bar on $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}$.

Let us then assume that there are some local properties of galaxies that depend on the received flux along their whole past history: for example, these could be the heating and cooling rates of the gas accreting onto the parent halo, which are in turn related to the star-formation rate. Therefore, we expect this dependence to be inherited by the galaxy number density $n_{g}$ [5]. More precisely, we can parameterize its response to the intensity of ionizing radiation by means of a Green’s function $G_{g}$ along the fluid worldline (the subscript “$g$” is to differentiate it from the Green’s function for the $q$ tracers, that is discussed in Section 3.4). Such function encodes the “UV physics” of galaxy formation that the bias expansion is oblivious about. At zeroth order, the RT effects on the mean galaxy density can then be written as

[TABLE]

where $G_{g}^{(0)}$ is the zeroth-order Green’s function and we have used the fact that $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ does not depend on $\hat{\bm{n}}$. Similarly, at first order in perturbations, we have

[TABLE]

where again we have used the fact that $\bm{x}_{\rm fl}$ is equal to $\bm{x}$ at the order we are working at.

Before proceeding, let us discuss Eqs. (3.7), (3.8) in more detail. What we have done is to isolate the contribution of ionizing radiation to the evolution of the galaxy number density (hence the subscript “$\rm ion$” on $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}$ and $\delta{n}_{g}$). At the order in perturbations we are working at, there is no loss of generality if we just sum this contribution to that arising from purely gravitational evolution, whose nonlocality scale is the halo Lagrangian radius $R(M_{h})$. That is, we can write the total dimensionless perturbation $\delta_{g}$ to the galaxy number density as the sum $\delta_{g}|_{\rm grav}+\delta_{g}|_{\rm ion}$, where $\delta_{g}|_{\rm ion}$ is obtained from Eqs. (3.7), (3.8) by simply taking the ratio $\delta n_{g}|_{\rm ion}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}$. We also emphasize that we allowed for two different Green’s functions, one for the response of the average galaxy number density and one for its perturbation. Indeed, there is no physical reason why the response should be the same at all orders in perturbations [1]. Besides, as we discussed at the beginning of this section, the sources of ionizing radiation are also expected to be a biased tracer of the underlying matter distribution. Therefore, $\delta{\cal I}$ will have both a deterministic and a stochastic contribution. If these two components have different emission spectra (which also effectively leads to them having different mean free paths), it is possible for galaxies to have a different response to each of them: we will not consider this case in the following, and just comment briefly in Section 3.3 on how this would modify the results we obtain here.

What are then the properties of these Green’s functions, ${G_{g}^{(0)}}$ and ${G_{g}^{(1)}}$? Let us first focus on their dependence on the photon four-momentum. At this order in perturbations we do not expect the galaxy response to depend on the photon arrival direction $\hat{\bm{n}}$: from the point of view of the large-scale bias expansion, this would amount to the presence of a preferred direction in the rest frame of the observer comoving with the fluid, which is forbidden by rotational invariance.555At higher orders in perturbations we can, for example, use vectors constructed from gradients of the matter overdensity. Then, in order to make progress, we assume that the time and energy dependencies can be factorized. Moreover, we assume that the energy dependence of both Green’s functions is proportional to the absorption cross section itself, i.e. we write

[TABLE]

where ${\cal G}$ is a constant with dimensions of an inverse length squared. This assumption is justified if we are thinking of the effect of photoevaporation of the gas accreting onto the halo (in any case, our conclusions will not change even if we consider a more general dependence on $E$, as we will see in Section 3.6).

We now turn to the time dependence of the Green’s functions, to have an idea of what are the different scales involved in the problem. Clearly, in the case of an instantaneous response of galaxy formation to the incoming flux, they would both be a function of $\eta$ times $\delta(\eta-\eta^{\prime})$. In general, there is a finite response time, i.e. we can imagine that the Green’s functions are peaked around $\eta^{\prime}=\eta$, with some finite width for the peak. The incoming ionizing radiation is (obviously) not interacting with the dark matter in the host halo, but with the gas accreting onto it. The time scale of this interaction is set by atomic physics, and then we expect it to be much faster than a Hubble time. Nevertheless, we also expect the star-formation rate to be tied to the accretion rate of baryonic matter on the host halo. This, in turn, is related to the total mass flow, which we know from gravity-only $N$-body simulations to be controlled by Hubble. For this reason, in this section we are going to consider the case in which the galaxy response is varying only on cosmological time scales. A more detailed discussion is presented in Section 5. With this assumption, there is no loss of generality in redefining the two Green’s functions in such a way as to reabsorb the factor $(1+z(\eta^{\prime}))^{3}/(1+z_{\ast})^{3}$ in Eq. (3.4). More precisely, we write

[TABLE]

We will later use the fact that the responses vary on time scales of order Hubble. We can now carry out the integral over $\mathrm{d}\eta^{\prime}$ and $\mathrm{d}\hat{\bm{n}}$ in Eq. (3.8). The full calculation is carried out in Appendix A.2: eventually, we see that $\delta n_{g}|_{\rm ion}$ is equal to

[TABLE]

where $\hat{\tau}$ is given by

[TABLE]

The photoionization of $1s$ hydrogen requires energies higher than the Lyman limit $E_{\infty}$: therefore the integral over $E$ in Eq. (3.11) starts from $E_{\infty}$.666Notice that $E(\eta^{\prime},\eta)\geq E$ for $\eta^{\prime}<\eta$ (photons lose energy as they travel towards the galaxy): therefore this constraint is automatically satisfied in Eq. (3.12). At energies much higher than $E_{\infty}$, but small enough that the emitted electron is still nonrelativistic, $\sigma_{\rm bf}(E)\sim E^{-7/2}$ (see e.g. [22, 15]). More precisely, we have

[TABLE]

where $\alpha=e^{2}/4\pi$ is the fine-structure constant. Corrections to this scaling can be expanded as a series in $E_{\infty}/E$, and are discussed in more detail in Section 3.6. Let us then make three further assumptions. First, we consider a power-law spectrum of emitted radiation, with a spectral index $s$ (which we take to be not too hard, i.e. $s<1$): this is, for example, the parameterization discussed in [18, 19], while [16] studies the cases $s=0$ and $s=-2$. Therefore, in Eq. (3.11) we can write $\varepsilon_{\rm em}(\eta_{\ast},E)$ as $\varepsilon_{\infty}(\eta_{\ast})(E/E_{\infty})^{s}$. Finally, we further approximate $\hat{\tau}(\eta,\lvert\bm{y}\rvert,E)$ as

[TABLE]

Clearly, the time evolution of the mean number density of neutral hydrogen must be taken into account, since it enters in Eq. (3.14) and gives an additional dependence on $\lvert\bm{y}\rvert$. We know that the ionization fraction basically changes from [math] to $1$ during reionization, and it is possible that this change happened very quickly, leading to a strong dependence of $\hat{\tau}$ on $\lvert\bm{y}\rvert$ (however, it is important to stress that we do not have yet a clear picture of the time evolution of the ionization fraction, since the main constraints come from the measurement of CMB anisotropies, that are only sensitive to the optical depth to redshift $z\approx 1100$). Therefore it is important to take this time dependence into account as well: we will come back to this point in Section 3.6, together with the full dependence of the cross section on redshift in Eq. (3.14). For now, the approximation of Eq. (3.14) is sufficient. With these assumptions, then, Eq. (3.11) becomes

[TABLE]

In the following we are going to neglect the restriction of the integral to the interior of the past light cone of $(\eta,\bm{x})$, since the mean free path is much shorter than $\eta-\eta_{\ast}$. Indeed, if we imagine to observe galaxies at $z=$1.5, which is a typical redshift for the upcoming large-scale galaxy redshift surveys, we have that the ratio between $(\eta-\eta_{\ast})$ and the m.f.p. is of order $100$, even taking conservatively $\eta_{\ast}$ corresponding to the end of reionization ($z\approx 6$), and using that the m.f.p. is of order 50${\,h^{-1}\,\mathrm{Mpc}}$ around these redshifts [6, 7]. In any case, including the constraint $\lvert\bm{y}\rvert\leq\eta-\eta_{\ast}$ is straightforward, but it would only complicate the calculations without adding relevant physics: we briefly discuss it in Section 3.6. Moreover, similarly to what we did in Eq. (3.10) and without loss of generality, we reabsorb the factor $(1+z_{\ast})^{s}/(1+z(\eta_{\ast}+\lvert\bm{y}\rvert))^{s}$ coming from the redshift dependence of the emissivity into the Green’s function. The integral over energy in Eq. (3.15) can then be carried out analytically. Let us however approximate it by

[TABLE]

from which we can define an “effective” m.f.p. ${\lambda}_{\rm eff}(\eta_{\ast})$ such that the integral is proportional to $\exp({-\lvert\bm{y}\rvert}/{\lambda}_{\rm eff}(\eta_{\ast}))$. This approximation works well for $\lvert\bm{y}\rvert\ll{\lambda}_{\rm eff}(\eta_{\ast})$, with corrections starting at order $\lvert\bm{y}\rvert^{2}/{\lambda}^{2}_{\rm eff}(\eta_{\ast})$ (we will discuss in Section 3.6 how to go beyond this approximation). To summarize, we have reduced Eq. (3.15) to the form

[TABLE]

For simplicity of notation we will regroup the overall factors in a single dimensionless quantity ${\cal C}_{\rm ion}$ in the following. Moreover, for the rest of this section we will always intend $\lambda_{\rm eff}$ as evaluated at $\eta_{\ast}$ and $\cal H$ as evaluated at $\eta$, unless stated otherwise. We now use the fact that the responses vary on time scales of order Hubble by expanding $G^{(1)}_{g}$ in a Taylor series around $\eta_{\ast}$:777We could have equivalently expanded the Green’s function around $\eta$. We have chosen $\eta_{\ast}$ as the expansion point since this makes calculations easier.

[TABLE]

where, without loss of generality, we have made the functions $g_{n}(\eta)$ dimensionless by factoring out ${\cal H}^{n}(\eta)$. The response for ${\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta)}$ can be expanded in a similar way, but this is unnecessary (as we will see in a moment).

Introducing the dimensionless perturbation $\delta_{\rm em}=\delta\rho_{\rm em}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\rm em}$ of the number density of emitters, Eq. (3.17) becomes

[TABLE]

Finally, after dividing Eq. (3.19) by the average galaxy number density $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}$ we obtain the expression for the dimensionless fluctuation $\delta_{g}|_{\rm ion}$, i.e.

[TABLE]

What does Eq. (3.20) mean for the bias expansion of $\delta_{g}|_{\rm ion}$? We see that Eq. (3.20) is an example of a nonlocal contribution to the bias expansion: only if its kernel is sufficiently localized is it possible to do a (spatial) derivative expansion and end up with local operators. Here it is the m.f.p. that controls the spatial extent of the kernel (galaxies are sensitive only to ionizing photons coming from a comoving volume of size $\sim\lambda_{\rm eff}^{3}$ along the past fluid trajectory), so the expansion will be in $\lambda_{\rm eff}\bm{\nabla}$. Consequently, all the higher-derivative operators become important at momenta $k\sim 1/\lambda_{\rm eff}$. This is the conclusion obtained in [5]. However, we also see that Eq. (3.20) is actually a resummation of all these higher-derivative terms: if we can treat the sum over $n$ perturbatively, we are able to predict the scale dependence of $\delta_{g}|_{\rm ion}$ also for $k\gtrsim 1/\lambda_{\rm eff}$ in terms of a finite number of functions of time. This can be seen more easily if we work in Fourier space. Eq. (3.20) is a convolution in real space, so $\delta_{g}|_{\rm ion}(\eta,\bm{k})$ is

[TABLE]

where $f_{\rm ion}(\eta,k^{2})$ (which can depend only on $k$ because the kernel in Eq. (3.20) depends only on $\lvert\bm{y}\rvert$, and cannot contain odd terms in $k$ because of locality) is given by

[TABLE]

We can simplify this expression for $f_{\rm ion}$ if we define

[TABLE]

where ${\cal F}^{(n)}_{\rm ion}$ are dimensionless functions of $k^{2}\lambda^{2}_{\rm eff}$, and the overall factor of $4\pi n!$ has been chosen so that they are all equal to $1$ for $k=0$. Further, since $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta)$ depends only on time, for each $n$ we can reabsorb the dimensionless factor $(4\pi n!)\times(\lambda_{\rm eff}{\cal C}_{\rm ion}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}(\eta_{\ast}))/(\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta))$ into the corresponding coefficient $g_{n}$. With these redefinitions, $f_{\rm ion}$ takes the form

[TABLE]

and in the following we will often suppress the time dependencies of the coefficients $g_{n}$ for simplicity of notation. For future reference, we also write down the explicit expression of the functions ${\cal F}^{(n)}_{\rm ion}$, i.e.

[TABLE]

with $f^{(0)}_{\rm ion}(x^{2})=4\pi\arctan(x)/x$. A plot of the first two functions $\smash{{\cal F}^{(n)}_{\rm ion}}$ is shown in Fig. 3.

As we discussed at the end of Section 2, we assume that $\delta_{\rm em}$ can be expressed in terms of the matter overdensity via a bias expansion with a nonlocality scale much shorter than the m.f.p. of ionizing radiation. Therefore everything depends on whether or not we are able to truncate the sums in Eq. (3.24) at a finite order: otherwise we would need to know an infinite number of functions of $\eta$ (i.e. the $g_{n}$) to compute $f_{\rm ion}$. Fortunately, this is possible, since the m.f.p. of ionizing radiation is much shorter than ${\cal H}^{-1}$. Therefore we can safely truncate the sums at some finite order in ${\cal H}{\lambda}_{\rm eff}$, and correspondingly introduce only a finite number of new bias coefficients in the bias expansion for galaxies.888We notice that in the limit of the effective mean free path going to zero there is no response of the galaxy number density to ionizing radiation if the Green’s function for $\delta n_{g}|_{\rm ion}$ is zero at $\eta=\eta_{\ast}$. This makes sense since in this limit only radiation emitted from events along the fluid worldline could have reached the galaxies.

We emphasize again how Eq. (3.24) has nothing to do with an expansion of Eq. (3.22) in powers of $k^{2}{\lambda}^{2}_{\rm eff}$. The point here is that the scale dependence at each order in ${\cal H}{\lambda}_{\rm eff}$ can be computed non-perturbatively in terms of the single parameter ${\lambda}_{\rm eff}$. By solving the full RT equation, and by simple locality arguments on the response of galaxies to the ionizing radiation, we have obtained a resummation of the expansion in spatial derivatives. This gives a precise correction to the shape of the power spectrum of galaxies with an amplitude set by a finite number of new bias coefficients, i.e. the $g_{n}$ appearing in Eq. (3.24).

Before proceeding to the next section, where we are going to compute the power spectrum of galaxies, we also recall that, given that we see a suppression of radiative-transfer effects in the clustering of high-mass galaxies [5], the Green’s function for $\delta n_{\rm g}|_{\rm ion}$ must be suppressed by some factor $p(M_{h},M_{\rm J})$ with respect to the one for the average galaxy density, with $p(M_{h},M_{\rm J})\ll 1$ for $M_{h}\gg M_{\rm J}$. This is clear since, as we can see from Eqs. (3.7), (3.22), the effect of radiative transfer on the fractional galaxy density perturbation is controlled by the logarithmic response of the galaxy density to the radiation, that is the ratio of the two Green’s functions in Eqs. (3.7), (3.8).

3.3 Galaxy statistics

We can now compute the contribution to the galaxy power spectrum from Eq. (3.21). This equation relates $\delta_{g}|_{\rm ion}$ at $\eta$ with $\delta_{\rm em}$ at $\eta_{\ast}$: that is, we only need to know the statistics of the sources at $\eta_{\ast}$ to compute the contribution of RT effects to $P_{gg}$ at $\eta$. At leading order in perturbations and derivatives we can write

[TABLE]

where $b_{\rm em}$ is the linear LIMD bias, and the operator $\epsilon_{\rm em}$ (not to be confused with the emissivity) captures the effect of short-scale physics on the evolution of the sources at this order in perturbations (see e.g. Section 2 of [1] for details). The stochastic term $\epsilon_{\rm em}$ is uncorrelated with $\delta$ and has a $k$-independent correlation function (up to terms of order $k^{2}R^{2}(M_{h})$), i.e.999Notice that, similarly to what happens to dark matter halos (or any other tracer), stochasticity is never completely uncorrelated with large-scale density fluctuations: gravitational evolution couples long- and short-wavelength modes, so that at higher orders in perturbations terms like $\epsilon_{{\rm em},\delta}\delta$ (where $\epsilon_{{\rm em},\delta}$ is a different operator than $\epsilon_{\rm em}$) must be included. See, e.g., [1] for a detailed discussion.

[TABLE]

Without loss of generality, we can write $\delta_{\rm em}(\eta_{\ast},\bm{k})$ in terms of $\delta(\eta,\bm{k})$ by reabsorbing the linear growth factor ${D_{1}(\eta_{\ast})}/{D_{1}(\eta)}$ in the bias parameter $b_{\rm em}(\eta_{\ast})$. Consequently, through $\delta_{g}|_{\rm ion}(\eta,\bm{k})$ the bias expansion for galaxies will gain the two new terms

[TABLE]

We emphasize a minor point about this expression for $\delta_{g}|_{\rm ion}$, that will nevertheless become relevant when we consider multiple emission times in Section 4.1. Here $\delta_{g}|_{\rm ion}$ is a sum of a deterministic and a stochastic term, each multiplied by a function of $k^{2}$. Even if the overall amplitude of these two functions is different, their scale dependence is the same (i.e. it is given by the single function $f_{\rm ion}$). Moreover, we can also see what would happen if we considered two different Green’s function for the stochastic and deterministic contributions to $\delta{\cal I}$ in Eq. (3.8). The calculations leading to Eq. (3.28) go through in the same way: the only difference is that instead of a single function $f_{\rm ion}$ we now have two different functions multiplying $\delta$ and $\epsilon_{\rm em}$. These two functions both have an expansion like that of Eq. (3.24), but each of them has its own set of coefficients $g_{n}$: this captures the difference between the Green’s functions for the stochastic and deterministic contributions.

The full galaxy fractional overdensity is simply given by the sum of $\delta_{g}|_{\rm grav}$ and $\delta_{g}|_{\rm ion}$ at this order in perturbations, as we have seen at the beginning of Section 3.2. That is, it takes the form

[TABLE]

where $\delta_{g}|_{\rm grav}$ contains the linear LIMD and stochastic terms ($b_{1}\delta$ and $\epsilon_{g}$, respectively).

Let us then compute the equal-time galaxy power spectrum. For this purpose, it is useful to factor out $g_{0}$ from our expression for $f_{\rm ion}$, since in any case it is degenerate with the bias coefficient $b_{\rm em}$ and the amplitude of the stochastic term $\epsilon_{\rm em}$: that is, in Eq. (3.29) we redefine

[TABLE]

so that now we have

[TABLE]

i.e. $f_{\rm ion}(\eta,k^{2})$ is independent of the galaxy response to ionizing radiation at leading order in ${\cal H}\lambda_{\rm eff}$. If we further reabsorb $g_{0}$ into $b_{\rm em}$ and $\epsilon_{\rm em}$, the equal-time galaxy power spectrum $P_{gg}(\eta,k)$ takes the form

[TABLE]

where $P_{\rm L}$ is the linear matter power spectrum and, following [1], we have defined the cross power spectrum between the two different stochastic terms $\epsilon_{g}$ and $\epsilon_{\rm em}$ as

[TABLE]

We notice that the contributions involving $\epsilon_{\rm em}$ are basically modifying the form of the higher-derivative corrections to the stochastic term $\epsilon_{g}$. These are present also if we consider gravitational interactions only: while in that case they become relevant at scales of order of the halo Lagrangian radius, here they are controlled by the m.f.p. of ionizing radiation.

From Eq. (3.33) we see that, if we stop at zeroth order in the expansion of Eq. (3.31) (so that the function $f_{\rm ion}$ is univocally determined in terms of the mean free path), at fixed $\eta$ we can predict the galaxy power spectrum up to $4$ constants, i.e. $b_{1}$, ${b}_{\rm em}$, and the amplitudes of the two stochastic terms $\epsilon_{g}$, $\epsilon_{\rm em}$.101010Of course we can consider also the effective mean free path as an additional free parameter, with a prior between $\sim$30${\,h\,\mathrm{Mpc}^{-1}}$ and $\sim$100${\,h\,\mathrm{Mpc}^{-1}}$, to account for uncertainties in the evolution of the mean density of neutral hydrogen. We also see that there are no degeneracies between these four parameters. Clearly we can discriminate between $b_{1},{b}_{\rm em}$ and $\smash{P^{\{0\}}_{\epsilon_{g}}},\smash{P^{\{0\}}_{\epsilon_{\rm em}}}$ thanks to the scale dependence of the matter power spectrum. Then, the scale dependence of $f_{\rm ion}$ allows to respectively distinguish $b_{1}$ from ${b}_{\rm em}$, and $\smash{P^{\{0\}}_{\epsilon_{g}}}$ from $\smash{P^{\{0\}}_{\epsilon_{\rm em}}}$.

If we go beyond the zeroth order in ${\cal H}\lambda_{\rm eff}$, we need more free parameters to compute $f_{\rm ion}$. Thanks to the different scale dependence of the functions $\smash{{\cal F}^{(n)}_{\rm ion}}$, however, it is in principle possible to discriminate between them.

We conclude this section with a brief discussion about the possible degeneracy between these RT effects and massive neutrinos. If we neglect the fact that massive neutrinos in general lead to a scale-dependent bias (see [23, 24, 25, 26, 27, 28]), the effect on the galaxy power spectrum is captured by the modification to the transfer function $T(k)$ for the total matter overdensity. More precisely, the relative correction to the deterministic part of the galaxy power spectrum is given by

[TABLE]

This is the green dot-dashed curve in Fig. 4. The corresponding correction due to RT effects, at leading order in $b_{\rm em}$, is given by (see Eq. (3.32))

[TABLE]

where $f_{\rm ion}$ is given by Eq. (3.31) and we recall that $b_{\rm em}$ must indeed be regarded as a small number since we reabsorbed the leading RT bias coefficient $g_{0}$ into it. Can this mimic the effect due to massive neutrinos? First, we see that the amplitude of this scale-dependent correction at $k=0$ is always degenerate with the linear LIMD bias. In other words, $b_{1}$ is fixed by the amplitude of $P_{gg}$ on large scales (obviously neglecting degeneracies with, e.g., $\sigma_{8}$, since they are not relevant for the sake of this discussion). This corresponds to redefining $b_{1}$ in Eq. (3.32) as $b_{1}(\eta)\to b_{1}(\eta)-b_{\rm em}(\eta_{\ast})f_{\rm ion}(\eta,0)$. Correspondingly Eq. (3.35) sees simply $f_{\rm ion}(\eta,k^{2})$ replaced by $f_{\rm ion}(\eta,k^{2})-f_{\rm ion}(\eta,0)$. This is plotted in Fig. 4 at leading order in ${\cal H}\lambda_{\rm eff}$ (blue curves). From the plot it is clear that these scale-dependent corrections are very similar in shape to those from massive neutrinos, i.e. RT effects could give rise to a bias in the constraints on $\sum m_{\nu}$ if not accounted for.111111In principle we could think of getting around this degeneracy by isolating $f_{\rm ion}$ from the corrections to the stochasticity (see Eq. (3.32)). However, we expect that neutrinos affect also the scale dependence of the stochastic term on scales $k\sim k_{\rm fs}$. Besides, as we are going to see in Section 4, when we consider the generic case of emission of radiation over Hubble time scales and add the inhomogeneities in the optical depth, the function of $k^{2}$ multiplying the stochastic term will not be the same as the one multiplying $\delta$.

3.4 Locally-biased tracers

We now move to the $q$ tracers. These tracers are assumed to be locally sensitive to the ionizing radiation emitted by the sources. We capture their memory of the emitted radiation field by allowing their number density to depend on the integral over $\log E$ and $\hat{\bm{n}}$ of the emission coefficient $(\rho_{\rm em}\varepsilon_{\rm em})/(4\pi)$, which is the total amount of ionizing photons emitted per unit time and unit volume. However, this is not enough. While we expect this dependence to be local in space, we must in principle integrate the emission coefficient also along the past fluid worldline, with a Green’s function $G_{q}$ that describes how $n_{q}$ responds to the emissivity of the sources. As in the case of galaxies, we allow for two different Green’s functions at zeroth and first order in perturbations. In principle these Green’s functions can also depend on the photon energy (while they cannot depend on the photon direction $\hat{\bm{n}}$ because of rotational invariance, at the order in perturbations we are working at). Therefore, we write $\delta{n}_{q}|_{\rm ion}$ as

[TABLE]

where we have used our assumptions of an isotropic emissivity. A similar relation holds for $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{q}|_{\rm ion}$. In Eq. (3.36) we have also used the fact that, at the order we are working at, everything is comoving. In principle we should have written our integral over the emission coefficient as, schematically, $n_{q}\supset{-\int\rho_{\rm em}\varepsilon_{\rm em}U^{\mu}_{\rm em}U_{\mu}}$, where ${-\rho_{\rm em}}U^{\mu}_{\rm em}U_{\mu}$ is the energy density of the emitters as seen from the observer comoving with the fluid. Calling $\bm{v}_{\rm rel}$ the relative velocity between the $q$ tracers and the sources we have ${-U^{\mu}_{\rm em}U_{\mu}}\sim 1+\lvert\bm{v}_{\rm rel}\rvert^{2}/{2}$, so that these effects are very suppressed in the nonrelativistic limit. Moreover, in presence of gravitational interactions only, $\bm{v}_{\rm rel}$ is only sourced starting from first order in derivatives, due to the equivalence principle.

From Eq. (3.36) we see that the integral in $\mathrm{d}E$ can be redefined as a new Green’s function that depends on $(\eta,\eta^{\prime})$ only. In the case of instantaneous emission, then, the fractional overdensity $\delta_{q}|_{\rm ion}(\eta,\bm{x})$ is simply proportional to $\delta_{\rm em}(\eta_{\ast},\bm{x})$ through a $\eta$-dependent function (even if this assumption is relaxed, using the bias expansion for $\delta_{\rm em}$ of Eq. (3.26) we can still simply absorb the integral over $\eta^{\prime}$ in the time dependence of the bias coefficients and stochastic terms, similarly to what happened when we moved from Eq. (1.2) to Eq. (1.3)). Therefore, the final expression for $\delta_{q}$ in Fourier space is simply

[TABLE]

where we kept the dependence on both stochastic terms $\epsilon_{q}$ and $\epsilon_{\rm em}$ (since they are different fields), while $b_{\rm em}\delta\subset\delta_{\rm em}$ was absorbed in the LIMD bias. We see that we cannot reabsorb the factor $b_{\epsilon_{\rm em}}$ in the amplitude of the stochastic term $\epsilon_{\rm em}$ without altering our expression for the galaxy power spectrum of Eq. (3.32).

Before proceeding we notice that, as we discussed below Eq. (3.8), also in this case we could allow for two different Green’s function in Eq. (3.36), one for the response to the deterministic part of the emission coefficient and one for its stochastic part. It is straightforward to see, however, that Eq. (3.37) would not change.

3.5 $gq$ cross-correlation

It is straightforward to compute $P_{gq}$ and $P_{qq}$ using Eqs. (3.29), (3.30), (3.37): we find that $P_{gq}$ is equal to

[TABLE]

while $P_{qq}$ is simply

[TABLE]

We notice that even if we stop at zeroth order in ${\cal H}\lambda_{\rm eff}$ we cannot constrain the parameter $b_{\epsilon_{\rm em}}$: a change in $b_{\epsilon_{\rm em}}$ can always be compensated by a change in the amplitude of $\epsilon_{q}$ in both Eq. (3.38) and Eq. (3.39). This can be easily seen also at the level of the fields from Eq. (3.37).

The importance of observations of the $q$ tracers lies in the (well-known: see [29]) fact that, if we combine $P_{gg}$, $P_{gq}$ and $P_{qq}$, it is possible to beat down cosmic variance in the constraints on the different free bias parameters that $f_{\rm ion}$ depends on. In this case we expect this multi-tracer technique to be useful also because we are able to predict exactly the scale dependence of the different functions of $k$ that make up $f_{\rm ion}$, i.e. Eq. (3.25): we are not expanding it in a power series in $k^{2}\lambda^{2}_{\rm eff}$, like we do with higher-derivative terms coming from gravitational dynamics (which are expanded in powers of $k^{2}R^{2}(M_{h})$).

Still, it is also important to keep in mind that the situation is very different from that of constraints on local primordial non-Gaussianity, in the context of which the multi-tracer technique was originally proposed. Indeed, there the equivalent of $f_{\rm ion}$ is a function that scales as $k^{-2}$ at small $k$, and is therefore completely orthogonal to the higher-derivative corrections from gravitational evolution. A precise analysis on the usefulness of the multi-tracer technique is clearly beyond the scope of this paper, so we will not investigate this topic further here.

3.6 More physics that can be captured in the bias expansion

We can now check all the possible effects that we have neglected in the above discussion and see if and how they can affect the scale dependence of $f_{\rm ion}(\eta,k^{2})$.

Integration over the past light cone

The first and simplest one is the restriction of the spatial integral to the past light cone, $\lvert\bm{y}\rvert\leq\eta-\eta_{\ast}$, in Eq. (3.15). Including this restriction would modify the integral in Eq. (3.23) for every $n$. More precisely, it is straightforward to check that at a given $n$ that Fourier transform is now proportional to a dimensionless function of $k/(\eta-\eta_{\ast})$, $k{\lambda}_{\rm eff}$ and $(\eta-\eta_{\ast})/{\lambda}_{\rm eff}$, with a proportionality factor that still scales as ${\lambda}_{\rm eff}^{n\,+\,1}$. Therefore, the same reasoning applies: these functions can be computed without the need of any perturbative expansion in these three dimensionless variables, and increasing orders in $n$ are still suppressed by $({\cal H}{\lambda}_{\rm eff})^{n}$.

General expression for the optical depth

Let us go back to Eqs. (3.14), (3.16): first, $\hat{\tau}(\eta,\lvert\bm{y}\rvert,E)$ takes the form

[TABLE]

With this definition, we can look in more detail to the integral over energy in the expression for $\delta n_{g}|_{\rm ion}$. As we discussed above Eq. (3.16), this integral can be carried out analytically (it can be written in terms of the incomplete Gamma function): while we could have used this function directly in Section 3.2, we have approximated its behavior as that of an exponential. Now, using Eq. (3.40) instead of Eq. (3.14), we can write Eq. (3.16) as

[TABLE]

which again is correct for small $\hat{\tau}_{\rm eff}$ up to ${\cal O}(\hat{\tau}^{2}_{\rm eff})$ (see discussion below and Appendix A.3). Let us then first see how we can treat the corrections that come from $\hat{\tau}_{\rm eff}(\eta,\,\lvert\bm{y}\rvert)$ not being simply proportional to $\lvert\bm{y}\rvert$. These come from the dependence of the photon energy and the average hydrogen density on redshift, i.e. from the time dependence of the mean free path.

Time dependence of the mean free path

It is relatively easy to deal with the time dependence of $n_{\text{HI}}a$ and $E$ (which is reflected by the fact that $\hat{\tau}_{\rm eff}(\eta,\lvert\bm{y}\rvert)$ is an integral over the variable $u$). We can write the exponential term in Eq. (3.41) as

[TABLE]

From this we see how these corrections give rise to the same type of terms that we already discussed in Section 3.2 (see e.g. Eq. (3.22)). Those involving the derivative of $1+z_{\ast}$ scale as ${\cal H}(\eta_{\ast}){\lambda}_{\rm eff}\sim{\cal H}{\lambda}_{\rm eff}$. Therefore, it is clear that at any fixed order in the expansion of $f_{\rm ion}(\eta,k^{2})$ in ${\cal H}{\lambda}_{\rm eff}$ we only need a finite number of these terms: for example, the first of them modifies the scale dependence of the second term on the right-hand side of Eq. (3.24). However, it is important to emphasize that these new terms do not involve any new free functions of time, so they affect $f_{\rm ion}(\eta,k^{2})$ in a way which is completely under control.

The same reasoning can, in principle, be applied to those involving the derivative of ${\lambda}_{\rm eff}$ at $\eta_{\ast}$. However, now we have

[TABLE]

It is then clear that we can safely reabsorb these terms in the expansion of Eq. (3.24) only if ${\partial\log{\lambda}_{\rm eff}(\eta_{\ast})}/{\partial\log\eta_{\ast}}$ is not much larger than $1$. This is not guaranteed to be true, since the mean free path can change rapidly during reionization. For example, if we assume that the neutral fraction evolves in a step-like fashion with a width $\delta\eta$, for $\eta_{\ast}$ close to the time at which reionization is halfway through we have that ${\partial\log{\lambda}_{\rm eff}(\eta_{\ast})}/{\partial\log\eta_{\ast}}\sim 1/({\delta\eta\cal H}(\eta_{\ast}))$, which could be large if the time for the transition is much shorter than an Hubble time. For example, while Strömgren spheres of ionized gas are forming, the m.f.p. is fixed by the mean size of the bubbles, which can be very short until the bubbles coalesce. After this happens, and the universe reionizes, it is controlled by the mean density of the residual neutral hydrogen. In this scenario the size of the time derivatives of the m.f.p. is controlled by how fast the coalescence happens with respect to a Hubble time.

Nevertheless we emphasize that, as before, the new terms of Eq. (3.43) do not add new free functions of time to the expansion of $f_{\rm ion}(\eta,k^{2})$. Even if the time evolution of the mean density of neutral hydrogen is still uncertain, this is different from our uncertainty on the microphysics of galaxy formation that is encoded by the Green’s functions (i.e. by the bias coefficients $g_{n}$). Moreover, we expect that modeling the time dependence of the average hydrogen density $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\rm HI}(\eta)$ is much less difficult than that of, say, the fluctuations $\delta{n}_{\rm HI}(\eta,\bm{x})$ on the past light cone of the galaxies. In other words, we can capture these corrections non-perturbatively in ${\cal H}{\lambda}_{\rm eff}$ once we assume a time evolution of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}$.

Corrections from the integral over energy

We then have to consider our approximation of the integral over energy in Eq. (3.41). We can account for the corrections to the right-hand side of that equation by using the full expression of the integral in terms of the incomplete Gamma function. This is shown in Appendix A.3. The integral over $E$ shows the same qualitative behavior at large $\lvert\bm{y}\rvert$, i.e. it goes to zero for $\lvert\bm{y}\rvert\gg{\lambda}_{\rm eff}$. Therefore, nothing stops us from doing the same expansion of the Green’s function in powers of $\lvert\bm{y}\rvert$ that led to our expression for $f_{\rm ion}(\eta,k^{2})$ as a power series in ${\cal H}{\lambda}_{\rm eff}$. The only difference will now be in the shape of the functions dimensionless functions of $k^{2}\lambda_{\rm eff}^{2}$ of Eqs. (3.23), (3.25). We show these functions in Fig. 5 of Appendix A.3. It is clear that their scale dependence is similar to those of Eq. (3.25): they both go to $1$ at small $k$ and vanish at large $k$, with the turnaround being at $k\sim 1/\lambda_{\rm eff}$.

Energy dependence of $\sigma_{\rm bf}$, $\varepsilon_{\rm em}$ and the Green’s functions, and multiple absorbers

The observation of the previous subsection is also what allows us to include the corrections to the bound-free cross section that we neglected in Eq. (3.13). They do not change the behavior of the integral over energy. The qualitative dependence of the cross section on energy is still the same, i.e. it goes to zero for $E\gg E_{\infty}$ and scales as $(E_{\infty}/E)^{7/2}$ for small $E$. Therefore, if we neglect the time dependencies (that can be treated perturbatively anyway, as discussed above), the final result is still given by $E_{\infty}$ times a function of $\lvert\bm{y}\rvert\sigma_{\infty}(n_{\text{HI}}a)(\eta_{\ast})$, where this function still goes to zero at very large $\lvert\bm{y}\rvert$. The same conclusions then apply.

We stress once again that in order to compute all the effects that we discussed here we require basically only the knowledge of the absorption cross section (that can be straightforwardly computed from first principles once we know what the absorbers are). Our uncertainty on the small-scale physics is only contained in the free functions $g_{n}$, exactly in the same way as the gravity-only bias expansion.

Before proceeding to the next section, we briefly investigate what happens if we consider different species of absorbers, and different assumptions for the energy dependence of the emissivity or of the galaxy response:

• •

including different absorbers with density $n_{\rm ab}$, each with their own absorption cross section $\sigma_{\rm ab}$, clearly affects the m.f.p., since now the optical depth $\tau$ is given by a sum over $\sigma_{\rm ab}n_{\rm ab}$ of terms like the one in Eq. (3.6). Therefore, as long as we know the energy dependence of the different $\sigma_{\rm ab}$, we can deal with this in the same way as we did when we considered only absorption from neutral hydrogen;

• •

things are a bit more complicated if we modify the energy dependence of the emissivity $\varepsilon_{\rm em}$ or the Green’s functions in Eq. (3.9). Let us assume, for example, that galaxy formation is sensitive to radiation at frequencies lower than $E_{\infty}$. Then, forgetting for a moment about redshift,121212Depending on the difference in redshift between emission and absorption it is possible that some photons reach the galaxy with energies lower than the ionization energy, even if they were emitted with energies higher than $E_{\infty}$. The optical depth for these photons is then still controlled by the bound-free cross section for part of their path to the galaxy. the optical depth for photons of such frequencies is controlled by the bound-bound cross section (i.e., mainly by the Rayleigh cross section). In this case, their effect on the bias expansion can still be treated in a similar way as we did above. More precisely, we can split the integral over $E$ in Eqs. (3.7), (3.8) in such a way that in each interval only one absorption line peaks. Then, for each of these sub-integrals it is possible to repeat the procedure above, ending up with a sum of different functions $f_{\rm ion}$: each of these functions will have its own specific m.f.p., which is controlled by the cross section at the corresponding absorption line;

• •

more problematic is instead the case where galaxy formation is sensitive to photons to which the universe is basically transparent, like for example X-rays, and the emissivity is significant at those frequencies (which is expected for active galactic nuclei and microquasars). In this case the effective m.f.p. is of order of the Hubble radius, and we expect that any derivative expansion breaks down. A perturbative bias expansion would clearly be insufficient to describe this situation, so we will not investigate it further in this paper.

4 Full radiative transfer

In this Section we see what happens if we drop basically all the assumptions on radiative transfer that we made in Section 3. We start by the considering a long epoch of emission of radiation on time scales of order a Hubble time and an inhomogeneous optical depth (Section 4.1). We show that now we need to add three functions of $k^{2}$ in the bias expansion (instead of one) in addition to the linear LIMD bias + stochastic term relation coming from gravitational interactions. The jump from one to two functions is due to the fact that now we need to take into account the dependence of the m.f.p. and of the coefficients $g_{n}$ on $\eta_{\ast}$ (the latter comes just from the fact that we are expanding the Green’s function of galaxies around the emission time). Combined with the fact that the source overdensity is not evolving via a single growth factor because of the presence of the stochastic term $\epsilon_{\rm em}$, which evolves independently from the LIMD term $b_{\rm em}\delta$, we end up with the relation $\delta_{g}|_{\rm ion}=f_{\rm em}(k^{2})\delta+f_{\epsilon_{\rm em}}(k^{2})\epsilon_{\rm em}$. In the same way, adding inhomogeneities in the density of absorbers further increases the number of functions from two to three. This is essentially due to the fact that they play the role of sinks of radiation, i.e. of sources with negative emissivity. The most important point, however, is not the increase in the number of functions, but the fact that now it is impossible to predict the shape of these higher-derivative corrections to the bias expansion at all orders in ${\cal H}\lambda_{\rm eff}$, unlike what we found in Section 3.2. It is enough to consider the case of homogeneous optical depth to see this. Indeed, now we need new response parameters for every significant fraction of a Hubble time during which the emissivity is not vanishing. Finally, in Section 4.2 we briefly study what is the impact of scattering on our analysis.

4.1 Adding multiple flashes and inhomogeneities in the optical depth

It is not difficult to see what happens if we move from the case of a very fast flash of radiation around $\eta=\eta_{\ast}$ to that of radiation emitted during an interval that can be of order of a Hubble time. In the case of a single emission, we have seen that at any given order $n$ in the ratio between the m.f.p. and the Hubble radius the scale dependence of $f_{\rm ion}$ can be predicted at the price of $n$ free functions of time. Now, even if we fix $n$, as we increase the number of emission times we need more and more free coefficients $g_{n}$. In the limit of continuous emission that extends on cosmological time scales, we would then need an infinite number of bias parameters. In a sense (to be made more precise below), we are effectively integrating over $\eta_{\ast}$. The function $f_{\rm ion}$ depends on $\eta_{\ast}$ both through the time dependence of the m.f.p. and, most importantly, through the coefficients $g_{n}$: since we do not know how they depend on time, we cannot carry out this integral.

Let us discuss this in more detail. In the interest of clarity, we focus separately on inhomogeneities in the density of emitters and inhomogeneities in the optical depth.

Inhomogeneous emitting medium

We first consider two short bursts of radiation at $\eta_{\ast,1}$ and $\eta_{\ast,2}$, with $\eta_{\ast,1}<\eta_{\ast,2}$ and $\eta_{\ast,2}-\eta_{\ast,1}\gtrsim{\cal H}^{-1}$. Following the same steps that brought to Eq. (3.19), we now have (see also Appendix A.2 for details)

[TABLE]

where we made explicit the dependence of ${\cal C}_{\rm ion}$ and $g_{n}$ on the emission time (we always intend $\cal H$ as evaluated at $\eta$, as in Section 3.2). From the above equation we see that, in Fourier space, $\delta_{g}|_{\rm ion}$ is given by

[TABLE]

where $f_{\rm ion}(\eta,\eta_{\ast,i},k^{2})$ takes the same form as in Eq. (3.22), and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta,\eta_{\ast,i})$ is the background galaxy number density computed for the single flash of radiation at $\eta_{\ast,i}$, so that $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta)$ is equal to $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta,\eta_{\ast,1})+\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}(\eta,\eta_{\ast,2})$.

Then, from Eq. (4.2) we see that, in full generality, we can write $\delta_{g}|_{\rm ion}$ as the sum

[TABLE]

The two functions $f_{\rm em}$ and $f_{\epsilon_{\rm em}}$ are surely different since the matter overdensity and the stochastic term do not evolve with time in the same way. However, it is more interesting to ask if the ratio between them is scale-independent (as was the case when we considered a single emission time in Section 3.3: see Eq. (3.28)). The answer is negative. More precisely, since the time evolution of both $\delta$ and $\epsilon_{\rm em}$ is independent of $k$ at the order we are working at, this is due only to the dependence of the function $f_{\rm ion}$ on the emission time (through, e.g., the mean free path).

Let us then see what happens in the limit of a continuous emission of radiation over an interval of order ${\cal H}^{-1}$. The question now is whether or not we can predict the scale dependence of the two functions multiplying $\delta$ and $\epsilon_{\rm em}$ in Eq. (4.3). The simplest way to see that the answer is negative is to focus on the leading order in the expansion in the mean free path. These two functions, then, depend on all the bias coefficients $g_{0}(\eta,\eta_{\ast,i})$ for $i=1,\dots,n_{\rm flashes}$ (generalizing the result of Eq. (3.24) to $n_{\rm flashes}$ emission times). As we bring $n_{\rm flashes}$ to infinity, we would then need an infinite number of functions of time, losing all predictivity.

Before proceeding to study the inhomogeneities in the number density of absorbers, we notice that the solution for the intensity of radiation received by the galaxy in the limit of continuous emission over an interval $\eta_{f}-\eta_{i}\sim{\cal H}^{-1}$ can be obtained simply by an integral over $\eta_{\ast}$. Looking at the full solution of the RT equation in Appendix A.1, it is straightforward to see that such limit corresponds to the replacement $\sum_{\eta_{\ast,i}}\Delta\eta\to\int_{0}^{\eta}\mathrm{d}\eta_{\ast}$, where $\Delta\eta$ (defined at the beginning of Section 3.2) is the duration of a single flash. However, we cannot obtain the dimensionless perturbation $\delta_{g}|_{\rm ion}$ in terms of an integral over $\eta_{\ast}$ as easily. The simplest way to see this is that $\Delta\eta$ always cancels in the expressions for $\delta_{g}|_{\rm ion}$ (as shown in Eq. (3.24), for example).

Inhomogeneous absorbing medium

We know that there are inhomogeneities in the number density of neutral hydrogen. The ionization fronts propagating into the neutral medium formed Strömgren spheres of ionized, heated gas, so we cannot treat $n_{\text{HI}}$ as homogeneous on scales of order of the typical bubble size, but we have to take into account $\delta_{\text{HI}}$ for scales $k\gtrsim${10}^{-2}\text{,}\mathrm{Mpc}^{-1}$$ [16]. We can imagine that these inhomogeneities in the absorbing medium can also be treated by a bias expansion on sufficiently large scales. Indeed, [17] has recently developed such a bias expansion for the inhomogeneities in the neutral fraction and, with it, a bias expansion for the $\rm 21cm$ signal from reionization (see also [30] for a recent computation of the position-dependent $\rm 21cm$ power spectrum via separate-universe simulations).

Since we have the full solution of the RT equation (Eqs. (A.1), (A.2), (A.3) in Appendix A.1), it should be easy, in principle, to study the effect of these inhomogeneities in the optical depth. However, since $\tau$ itself is given by an integral along the line of sight, the resulting expressions for $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{g}|_{\rm ion}$ and $\delta n_{g}|_{\rm ion}$ are very complicated. For this reason, we study the effect of an inhomogeneous $\tau$ in the case of radiation emitted in a single flash around $\eta_{\ast}$: this is simpler to treat, and the generalization to continuous emission is straightforward.

Let us go back to Eqs. (3.4), (3.6). We can expand $n_{\text{HI}}$ in a background $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}$ and a dimensionless perturbation $\delta_{\text{HI}}$. We see that, at first order in perturbations, the contribution of the $\tau$ inhomogeneities to the perturbations $\delta{\cal I}$ of the specific intensity of ionizing radiation received by the galaxies is equal to

[TABLE]

where the energy-dependent m.f.p. is given, as in Eq. (3.14), by $\lambda_{\rm ion}=1/(\sigma_{\rm bf}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}a)$, and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ is obtained from Eq. (3.6) by taking $n_{\text{HI}}=\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}$. Then, we see that the structure of Eq. (4.4) is similar to that of the full solution of Appendix A.1, if we take a homogeneous optical depth and an emissivity proportional to $\sigma_{\rm bf}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\text{HI}}$. That is, the absorber medium plays a role of an additional “negative” source term (i.e. a sink), whose emissivity is not localized around a single time $\eta_{\ast}$ (unlike that of the emitters).

Generalizing to an arbitrary number of emission times, this tells us that the inhomogeneities in the density of absorbers contribute to $\delta_{g}|_{\rm ion}$ through two additional functions of $(\eta,k^{2})$, which we can call $f_{\text{HI}}$ and $f_{\epsilon_{{\text{HI}}}}$. These two functions multiply $\delta$ and $\epsilon_{\rm HI}$ respectively. There is no loss of generality in summing $f_{\text{HI}}$ and $f_{\rm em}$ into a single function, given that both multiply the matter overdensity $\delta$. We cannot instead reabsorb $f_{\epsilon_{{\text{HI}}}}$ in $f_{\epsilon_{\rm em}}$ since the two stochastic terms $\epsilon_{\text{HI}}$ and $\epsilon_{\rm em}$ are different fields. However, the most important point to emphasize is that, again, the scale dependence of $f_{\text{HI}}$ and $f_{\epsilon_{{\text{HI}}}}$ cannot be computed perturbatively with a finite number of bias parameters.

4.2 Scattering

While a full discussion of the impact of scattering on the radiative-transfer equation is beyond the scope of this paper, in this short section we briefly draw a qualitative picture of how our previous results would be affected if they play an important role.

The most important effect is the redistribution of the direction of the photons.131313The scattering kernel can be only a function of $\hat{\bm{n}}\cdot\hat{\bm{n}}^{\prime}$ at the order in perturbations we are working at. The reason is the same as the one we used to conclude that the emissivity and the galaxy response cannot depend on $\hat{\bm{n}}$. Let us then consider again the case of a single flash to build some intuition. We can imagine that, after photons are emitted isotropically at $\eta_{\ast}$, some of those that are initially directed away from the galaxy are later scattered towards it at some time $\eta_{\rm sc}>\eta_{\ast}$. This effectively mimics a second burst of radiation at $\eta_{\rm sc}$. Therefore, in the limit of scattering happening continuously over time, we expect to roughly go back to the same situation that we have discussed when we have included inhomogeneities in the optical depth.

Finally, what happens if scattering dominates over absorption? In this case we expect that the intensity is made isotropic on scales larger than the (comoving) diffusion length $D\sim\sqrt{1/(\sigma_{\rm sc}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\rm sc}a\cal H)}$ ($\sigma_{\rm sc}$ and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{n}_{\rm sc}$ being the overall amplitude of the scattering cross section and the average number density of scatterers, respectively). On these scales, the monopole of the intensity will follow a (sourced) diffusion equation, so once we compute $\delta_{g}|_{\rm ion}$ using Eqs. (3.7), (3.8) we see that the higher-derivative terms from RT effects can be expanded in powers of $k^{2}D^{2}$.

5 Discussion and conclusions

In this paper we have investigated whether it is possible to find a resummation of the higher-derivative terms in the bias expansion that come from radiative-transfer (RT) effects.

We have shown that, at linear order in perturbations and when mainly absorption and emission play a relevant role in the RT equation, these effects are captured by three functions of $k^{2}\lambda_{\rm eff}^{2}$, where $\lambda_{\rm eff}$ is the effective mean free path of ionizing radiation. Whether or not it is possible to predict the scale dependence of these three functions, instead of simply relying on their expansion in powers of $k^{2}\lambda_{\rm eff}^{2}$, depends on the following factors:

• •

the time dependence of the emissivity of the sources of ionizing radiation;

• •

the time dependence of the response of galaxies to the flux of ionizing radiation;

• •

the presence of inhomogeneities in the optical depth.

This is summarized in Tabs. 1, 2. For example, let us consider the case of the galaxy response varying on cosmological time scales; that is, the observable properties of a given galaxy sample retain a memory of the ionizing flux received at some earlier time. Then, the dependence of these functions on the details of galaxy formation can be absorbed in generalized RT bias coefficients, without needing a detailed “UV” modeling of galaxy formation, if: a) one can neglect inhomogeneities in the optical depth; b) the emission of ionizing radiation only happens for a short period of time $\Delta\eta\ll{\cal H}^{-1}$. If these assumptions do not hold, we cannot predict the dependence on $k$ of these three functions unless we know precisely how the response of galaxies to the ionizing radiation depends on time. It is important to stress that this does not happen in the gravity-only bias expansion. There, the time dependence of the Green’s function can be reabsorbed, order by order in perturbations and spatial derivatives, into the time-dependent bias coefficients. Here we cannot do this because, through the received flux, the galaxies are sensitive to the inhomogeneities in the distribution of sources and sinks of ionizing radiation evaluated along their past light cone, and not only along the past fluid worldline. Consequently, the time dependence of the Green’s function affects also the scale dependence of the bias.

It is however important to emphasize that it is entirely possible for the response to ionizing radiation to happen on time scales much faster than a Hubble time. For example, in the case of line emission from diffuse gas we can imagine that the response to the ambient radiation field is controlled by the recombination rate, which has nothing to do with Hubble. More precisely, if we wanted to be completely general, we should have considered a Fourier transform of the Green’s functions of Section 3.2 with respect to $\eta-\eta^{\prime}$. This Fourier transform can have support for both $\omega\sim{\cal H}$ and for $\omega\gg{\cal H}$, and in this work we have focused on the case where the Fourier coefficients are nonzero at low frequencies. In this sense, our study is conservative: as we can see from Tabs. 1, 2, the contributions to the bias expansion due to the high-frequency part of the response are always under control.

In this paper we have followed an effective field theory approach, i.e. we have not made any specific assumption on the Green’s functions describing the response to the radiation (apart from very general assumptions on their time dependence, as discussed above). In other words, our results apply equally well to any other tracer of the underlying matter distribution that is sensitive to RT effects.

For galaxies, the suppression of the star-formation rate coming from photo-evaporation of the gas accreting onto the parent halo is not very relevant for halos whose mass is much larger than the Jeans mass. We thus expect that the bias coefficients of the higher-derivative terms from these RT effects are actually a strong function of the parent halo mass, and their amplitude is very small for $M_{h}\gg M_{\rm J}$. Since we do expect the RT to typically lead to a smooth contribution to the galaxy power spectrum (see Fig. 3), a contribution that is very small in amplitude can presumably be absorbed by marginalizing over a sufficiently flexible template for the RT effects. However, this might well lead to a degradation of constraints on the neutrino mass $m_{\nu}$ or equilateral non-Gaussianity $f_{\rm NL}^{\rm equil}$. In particular, Fig. 4 illustrates that the RT effects are expected to have a scale dependence that is quite similar to that induced by nonzero neutrino masses.

However, we do not expect RT effects to be necessarily small for tracers like, e.g., the Lyman-$\alpha$ forest. In this case, we can imagine that treating the higher-derivative terms perturbatively (i.e. through an expansion in powers of $k^{2}$) would to lead to a significant loss of constraining power on cosmological parameters in general. Indeed, the bias expansion would stop being predictive already at the large scales where the RT effects leave their imprint on the power spectrum of these tracers.

Future prospects

An explicit modeling of the response to ionizing radiation is surely a way to go around this problem. However, it is clear that any uncertainty on the theoretical errors of this model would reduce the robustness of galaxy clustering as a cosmological probe. A second way relies on the fact that we do not expect RT effects to induce a sizeable velocity bias. Indeed, let us consider the momentum density $\bm{p}_{\cal I}$ carried by the ionizing radiation field, as seen by an observer comoving with the galaxies. It is straightforward to show that, at leading order in derivatives, it is proportional to $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\cal I}\lambda_{\rm eff}\bm{\nabla}\delta$ (see Appendix A.5), where $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\cal I}$ is the average energy density of the radiation field. The momentum transfer from the radiation field to the galaxies, and consequently the contribution from RT to the velocity bias $\bm{v}_{g}-\bm{v}$, is then controlled by the ratio $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\cal I}/\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{b}=\Omega_{\cal I}/\Omega_{b}$. Since the density of radiation is dominated by the CMB, this ratio is much smaller than $\Omega_{r}/\Omega_{b}\approx 5\Omega_{r}/\Omega_{m}\approx 10^{-3}(1+z)$. From this we conclude that the peculiar velocities of galaxies, and with them the higher-multipole contribution to the galaxy overdensity in redshift space, are very weakly affected by radiative transfer (as pointed out by [12]), allowing us to obtain, in principle, unbiased measurements on cosmological parameters.

Comparison with recent work

While this work was being completed, two papers that investigate similar topics have been published on the arXiv [31, 32]. In [31] the authors compute the inhomogeneities $\delta_{\cal I}$ of the intensity $\cal I$ (integrated over angles and averaged over frequencies with some weighting $w(E)$, whose frequency dependence is assumed to be that of the HI photoionization cross section). The spirit of our calculation is similar to theirs (see for example Appendix A.4), with the difference that they consider directly the physically-motivated case of continued emission over Hubble time scales. The work [32], instead, studies the impact of ionizing radiation on the clustering of tracers. The difference with our computation is that the authors assume that the contribution of radiative transfer, i.e. what we called $\delta_{g}|_{\rm ion}$, is given by $b_{\cal I}\delta_{\cal I}$, with $b_{\cal I}$ a time-dependent bias coefficient (see their Eqs. (1), (2)). That is, they assume that the scale of nonlocality in time of the response of tracers to ionizing radiation is much faster than Hubble, so that their Green’s functions are proportional to $b_{\cal I}(\eta)\delta(\eta-\eta^{\prime})$. Consequently, based on this strong simplifying assumption and after assuming a model for the emissivity and the absorption coefficient, they are able to compute the scale dependence of these corrections to galaxy clustering without requiring an infinite number of bias coefficients. This is in agreement with our conclusions (see e.g. Tabs. 1, 2).

Acknowledgements

It is a pleasure to thank Philipp Busch, Chris Byrohl, Jens Chluba, Ildar Khabibullin, Eiichiro Komatsu, Kaloian Lozanov, Matt McQuinn, and Shun Saito for useful discussions. G. C. and F. S. acknowledge support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” from the European Research Council.

Appendix A Review of radiative-transfer physics

In this appendix we review the solution of the radiative-transfer equation in an FLRW universe and give some details on the calculations of Sections 3 and 5.

A.1 Solution of Boltzmann equation

In Section 3.1 we discussed in detail the emission and absorption coefficients, and derived the equation of radiative transfer, i.e. Eq. (3.3) (with Eq. (3.1) being its fully general-relativistic version). This equation can be straightforwardly solved via an integral along the line of sight. More precisely, the phase-space density can be written as

[TABLE]

where we defined the source term $\cal S$ as

[TABLE]

In Eq. (A.1) we integrate starting from $\eta^{\prime}=0$ since the emissivity vanishes for times earlier than some $\eta_{i}$ in any case. The solution for the optical depth $\tau$, which enters in the integral of Eq. (A.1), is

[TABLE]

In both Eqs. (A.1), (A.3), the evolution of the energy is given by

[TABLE]

A.2 Solution for single flash

Let us now study in more detail the solution in the case of a single flash. The assumption that the emissivity is nonzero only for a short interval $\Delta\eta\ll{\cal H}^{-1}$ around $\eta_{\ast}$, which we obviously take to be less than $\eta$ in Eq. (A.1), amounts to writing

[TABLE]

Plugging this into Eq. (A.1) we arrive directly at Eq. (3.4), where ${\cal I}_{\ast}$ is indeed given by

[TABLE]

as we wrote below Eq. (3.4) (noticing that we called $\varepsilon_{{\rm em},\ast}(E)=\varepsilon_{{\rm em}}(\eta_{\ast},E)$ in Section 3.2). The generalization to multiple quick flashes of duration $(\Delta\eta)_{i}$ at different times $\eta_{\ast,i}$, which has been used in Section 4.1, is straightforward: one can just take the emissivity to be a sum of terms like that of Eq. (A.5). Moreover, it is also easy to see that we can recover the full solution of Eq. (A.1) by integrating the solution for a single flash in $\mathrm{d}\eta_{\ast}$.

We are now in the position to check in detail how to go from Eq. (3.8) to Eq. (3.11). Using Eqs. (3.4), (3.9), (3.10), (3.18), (A.6), we see that $\delta n_{g}|_{\rm ion}$ is given by

[TABLE]

where the integral over energy actually starts from $E_{\infty}$ and the optical depth $\tau$ at all orders in perturbations is given by (see Eq. (3.6))

[TABLE]

To see that the integral in $\mathrm{d}\eta^{\prime}\mathrm{d}\hat{\bm{n}}$ is just an integral over $\mathrm{d}^{3}y/\lvert\bm{y}\rvert^{2}$ we redefine $\eta^{\prime}-\eta_{\ast}\equiv y$ (where $y$ will be equal to $\lvert\bm{y}\rvert$ at the end). Then, we have (dropping the overall proportionality factor $({\Delta\eta\,a_{\ast}\,{\cal G}})/({4\pi})$ to simplify the notation)

[TABLE]

Then, in the integral for $\tau$ we define $\eta^{\prime\prime}=\eta_{\ast}+uy$, so that Eq. (A.9) becomes

[TABLE]

Using $\hat{\bm{n}}y=\bm{y}$, and $\mathrm{d}y\mathrm{d}\hat{\bm{n}}=\mathrm{d}^{3}y/\lvert\bm{y}\rvert^{2}$, we finally recognize Eq. (3.11). Moreover, we see that $\hat{\tau}$ is given by

[TABLE]

so that if we assume a homogeneous density of absorbers (i.e. of $1s$ hydrogen, in this case), we indeed find Eq. (3.12).

A.3 Full integral over energy

In this appendix we compute the full integral of Eq. (3.41) and discuss the corresponding corrections to the scale dependence of the bias. In order to do this, we first make the (physically motivated) assumption that the emissivity vanishes above some energy $\widetilde{E}>E_{\infty}$. Further, we assume that its spectral index $s$ is equal to zero, for simplicity. Dropping all corrections coming from the redshift dependence of the photon energy and the time dependence of the mean free path (which are discussed in detail in Section 3.6), Eq. (3.41) (or, equivalently, Eq. (3.16)) becomes

[TABLE]

where $\lambda_{\rm ion}(\eta_{\ast},E)$ is defined by Eqs. (3.13), (3.14), and in the rest of this appendix we define the effective m.f.p. as $\lambda^{-1}_{\rm eff}(\eta_{\ast})\equiv\sigma_{\infty}(n_{\text{HI}}a)(\eta_{\ast})$ (and drop its dependence on $\eta_{\ast}$ to avoid cluttering the notation).

If we then follow Eq. (3.17) and define the dimensionless parameter ${\cal C}_{\rm ion}$ as

[TABLE]

we can read off from Eq. (A.12) the equivalent of the exponential $\exp({-\lvert\bm{y}\rvert/\lambda_{\rm eff}})$. This function has the same behavior of $\exp({-\lvert\bm{y}\rvert/\lambda_{\rm eff}})$: it goes to $1$ for $\lvert\bm{y}\rvert\ll\lambda_{\rm eff}$, and decays exponentially at large $\lvert\bm{y}\rvert$.141414In order to have this exponential decay it is necessary that the integral over energy does not go up to $E=+\infty$. This is why in Eq. (A.12) we have assumed that the emissivity vanishes for $E>\widetilde{E}$. Therefore, following the same steps that lead to Eqs. (3.23), (3.24), we can define the equivalent of the functions $\smash{{\cal F}^{(n)}_{\rm ion}}$ of Section 3. These will still be dimensionless functions of $k^{2}\lambda^{2}_{\rm eff}$ (moreover, they now depend also on the ratio $\widetilde{E}/E_{\infty}$). Their behavior with $k$ is the same as those of Eq. (3.25): they go to $1$ at small $k$ and vanish at large $k$. A plot of the first two functions for different values of the effective m.f.p. is shown in Fig. 5.

A.4 Fourier transform of angle-averaged intensity

In this appendix we compute the dimensionless perturbation to the angle-averaged intensity $\cal I$, in order to make contact with the formalism of [31, 32] and to justify the results of Tabs. 1, 2. Once we have this quantity, if we assume that the response of galaxies to the incoming radiation happens on a time scale $\Delta\eta\ll{\cal H}^{-1}$, we can compute the contribution $\delta_{g}|_{\rm ion}$ to the galaxy overdensity by integrating it over energy and time with a single Green’s function proportional to the photoionization cross section times $b_{\cal I}(\eta)\delta(\eta-\eta^{\prime})$.

The general solution of the RT equation is given in Eq. (A.1). Focusing only on inhomogeneities in the density of emitters, we have

[TABLE]

where we first changed from $\eta^{\prime}$ to $y=\eta-\eta^{\prime}$, and then to $\bm{y}=\hat{\bm{n}}y$ (with $\mathrm{d}y\mathrm{d}\hat{\bm{n}}=\mathrm{d}^{3}y/\lvert\bm{y}\rvert^{2}$). The optical depth in the above equation is given by

[TABLE]

where we defined $\eta^{\prime\prime}=\eta-u(\eta-\eta^{\prime})=\eta-u\lvert\bm{y}\rvert$ to go from the first to the second equality.

Eq. (A.14) is all we need to compute $\delta_{g}|_{\rm ion}$ if the Green’s function of galaxies is proportional to $\delta(\eta-\eta^{\prime})$. After we integrate it over energy with an energy-dependent weight, we can see that its structure is very similar to Eq. (3.8), with the obvious difference that there we have considered only a single flash of radiation. Writing $\delta\rho_{\rm em}=\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\rm em}\delta_{\rm em}$ and expanding the dimensionless fluctuation $\delta_{\rm em}$ as in Eq. (3.26), we can compute $\delta_{g}|_{\rm ion}$ by taking advantage of the fact that $\delta_{\rm em}$ is the sum of two uncorrelated fields both evolving in a scale-independent way, i.e. by means of a growth factor (indeed, if the stochastic term has a power spectrum which is localized in real space at all times, its time evolution cannot depend on $k$). More precisely, we can write $\delta_{\rm em}$ in Eq. (A.14) as

[TABLE]

where $D_{1}$ and $D_{\epsilon_{\rm em}}$ are the linear growth factors for $\delta$ and $\epsilon_{\rm em}$, respectively.

From Eq. (A.16) we see that, technically, $(D_{1}(\eta^{\prime})/D_{1}(\eta))\,b_{\rm em}(\eta^{\prime})$ and $D_{\epsilon_{\rm em}}(\eta^{\prime})/D_{\epsilon_{\rm em}}(\eta)$ play the role of two different Green’s functions for the deterministic and stochastic contributions. In this sense, the case of continuous emission and instantaneous response mirrors that of an extended response and a single flash of radiation discussed in Section 3.2 (compare $\smash{G^{(1)}_{g}(\eta,\eta^{\prime})\propto\delta(\eta-\eta^{\prime})}$ with Eq. (A.5)). However, it is important to emphasize that $(D_{1}(\eta^{\prime})/D_{1}(\eta))\,b_{\rm em}(\eta^{\prime})$ and $D_{\epsilon_{\rm em}}(\eta^{\prime})/D_{\epsilon_{\rm em}}(\eta)$ are very different from the Green’s functions that we have considered in Section 3.2, since they have nothing to do with the response of galaxies to ionizing radiation but only encode how $\delta_{\rm em}$ evolves in time.

These calculations also show that it is still possible to arrive at a resummation of the higher-derivative terms from RT effects even if we include the inhomogeneities in the optical depth. More precisely, expanding $\exp({-\tau})$ at linear order in spatial fluctuations, similarly to what we did in Eq. (4.4), and using the fact that the deterministic and stochastic parts of $\delta_{\rm ab}$ evolve in a scale-independent way (see e.g. Eq. (A.16) above), we end up with an additional “negative” source whose emissivity is not localized around a single time. Thus, we fall in the same situation above. This tells us that the only obstacle to reliably predicting the scale dependence of the bias due to RT effects is having a galaxy response with support also for frequencies comparable to $\cal H$ (see also Tab. 2).

A.5 Momentum density of the radiation field

In this appendix we sketch how to compute the momentum density of the ionizing radiation field (assuming for simplicity an homogeneous optical depth: adding inhomogeneities in $\tau$ will not change these results, as we see for example from the first columns of Tabs. 1, 2 and the discussion at the end of the previous appendix). The energy density and momentum density for the observer $U^{\mu}$ described in Section 3.1 are151515Technically, we should consider the fact that the fluid trajectory $\bm{x}_{\rm fl}$ is different from $\bm{x}$: however, these corrections come in at higher orders in perturbations so we neglect them (see also the discussion in Section 3.4).

[TABLE]

More precisely, paralleling Eq. (A.14), we have

[TABLE]

It is straightforward to see that the same conclusions of Appendix A.4 apply, i.e. it is possible to obtain a resummation of all the higher-derivative contributions to the momentum density $\bm{p}_{\cal I}=\delta\bm{p}_{\cal I}$. However, for the discussion in Section 5 it is sufficient to stop at leading order in derivatives. Expanding the density of emitters in spatial derivatives inside the integral, we see that

[TABLE]

where $\alpha$ is a time-dependent dimensionless coefficient of order $1$ and $\rho_{\cal I}$ is obtained from Eqs. (A.14), (A.17a) by taking $\rho_{\rm em}=\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\rm em}$.