Recall that we had, using Bolometric Radiative Equilibrium, an equation which described the greenhouse effect: $$\sigma T^4=\sigma T_0^4\left[1+{3\over2}\tau\right]$$ Now we want to talk about the effects of the diffusion of photons. For this, we have the general diffusion equation: $$F=-D\nabla n$$ For photons, $F$ is the energy flux, $D$ is $\lambda_{mfp}\cdot c$, and $n\sim{\sigma\over c}T^4$ is the number density of photons. Then: $$F\sim\underbrace{\lambda_{mfp}\over L}_{1\over\tau}c{\sigma\over c}T^4 \sim{\sigma T^4\over\tau}$$ Recall that $F\equiv\sigma T_e^4$, so: $$T^4\sim T_e^4\tau$$ This says that as we go deeper into the atmosphere, the temperature increases, but slowly (as the fourth root).
A while ago, we had the equation: $$\mu{dI_\nu\over d\tau_\nu}=I_\nu-\~\omega_\nu J_\nu-(1-\~\omega_\nu)B_\nu$$ and we decided to simplify our lives, we had set the scattering term $\~\omega_\nu=0$. Now we'd like to solve the compliment problem: we'll set $B_\nu=0$ and look at the equation for pure scattering: $$\mu{dI_\nu\over d\tau_\nu}=I_\nu-\~\omega_\nu J_\nu$$ The first thing we'll do to simplify our integro-differential equation is to say that scattering preserves frequency (which it does, except for Compton scattering). Then we have: \def\wtTemplate:\~\omega $$\mu\ddtau{I_T}=I_T-{\wt\over4\pi}\int{I_Td\Omega}$$ where we say the total intensity $I_T(\tau,\mu)=I_*(\tau,\mu)+I(\tau,\mu)$. $I_*$ is the incident (never scattered) intensity, and $I$ is the diffuse (scattered at least once) intensity. Our equation for the incident component is: $$\mu\ddtau{I_*}=I_*$$ This is a familiar equation with a familiar solution: $$I_*=I_*(\tau=0)e^{\tau\over\mu}$$ You'd expect that $I_*$ should decrease due to scattering, so why is the exponent positive in this equation? Remember by our convention, that $\mu$ points toward the direction of incoming flux, so a positive $\mu$ corresponds to looking upward in the atmosphere, in which direction the intensity of unscattered light should increase. To make our equation look like the solution in Chamberlain, we define: $$I_*(\tau=0)\equiv\pi\mathfrak{F}\delta(\mu+\mu_0)\delta(\phi-\phi_0)$$ $$\boxed{I_*=\pi\mathfrak{F}\delta(\mu+\mu_0)\delta(\phi-\phi_0)}$$ where $\mu_0$ is the angle the incident intensity makes with the plane of the atmosphere, and $\phi$ measures right ascension.\par Now the whole diffusion equation looks like: \begin{aligned}\mu\ddtau{I}&=I-{\wt\over4\pi}\int{I_*+I)d\Omega}\\ &=I-{\wt\over4\pi}\int{Id\Omega}-{\wt\over4\pi}\int{I_*d\Omega}\\ &=I-{\wt\over4\pi}\int{I\overbrace{d\Omega}^{d\phi d\mu}}-{\wt\over4\pi}\pi \mathfrak{F} e^{-\tau\over\mu_0}\\ &=I-{\wt\over4\pi}2\pi\int_{-1}^1{Id\mu}-{\wt\over4}\mathfrak{F} e^{-\tau\over\mu_0}\\ \end{aligned} The key to solving this equation is to get rid of that integral. To do this, we'll do a math trick by turning the integral into a sum: $$\int_{-1}^1{Id\mu}=\sum_{j=1}^n{a_jI(\mu_j)}$$ In general, for some polynomial of order $(2m-1)$, $f(x)=c_0+c_1x+\dots+c_{ 2m-1}x^{2m-1}$: $$\int_{-1}^1{f(x)dx}=\sum_{j=1}^n{a_jf(x_j)}$$ If we want to get an exact answer, each individual term in the polynomial had better give exact equality: $$\int_{-1}^1{c_\ell x^\ell dx}=a_1c_\ell x_1^\ell+a_2c_\ell x_2^\ell+\dots +a_nc_\ell x_n^\ell$$ which is an equation in $2n$ unknowns. Overall, we have $2m$ equations, so to get perfect equality for a $(2m-1)^{th}$ order polynomial, we must demand that we have $n=m$ sample points. For a given $n$, we then need the coefficients $a_i[i=1,\dots,n]$, $x_i[i=1,\dots,n]$. Luckily, values for these are tabulated in books, and in Numerical Recipes. We'll work out a simple case for $n=2$: \def\intmoo#1{\int_{-1}^1{#1}} \begin{aligned}\int_{-1}^1{Id\mu}&=a_1I(\mu_1)+a_2I(\mu_2)\\ \intmoo{\mu^0d\mu}&=2=a_1\mu_1^0+a_2\mu_2^0\\ \intmoo{\mu^1d\mu}&=0=a_1\mu_1^1+a_2\mu_2^1\\ \intmoo{\mu^2d\mu}&={2\over3}=a_1\mu_1^2+a_2\mu_2^2\\ \intmoo{\mu^3d\mu}&=0=a_1\mu_1^3+a_2\mu_2^3\\ \end{aligned} If we stare at these equations for a while, maybe we'll see that the answer is $a_1=a_2=1$, and: $$\begin{matrix}\mu_1=+{1\over\sqrt{3}},&\mu_2=-{1\over\sqrt{3}}\end{matrix}$$ Getting back to our original equation, we now have: \def\aiaio{\left[a_1I(\mu_1)+a_2I(\mu_2)\right]} \begin{aligned}\mu\ddtau{I} &=I-{\wt\over2}\intmoo{Id\mu}-{\wt\mathfrak{F}\over4}e^{-\tau\over\mu_0}\\ (n=2)\,\Rightarrow\,\mu\ddtau{I(\mu)}&=I-{\wt\over2}\aiaio -{\wt\over4}\mathfrak{F} e^{-\tau\over\mu_0}\\ \end{aligned} where $\mu$ is the only remaining parameter we need to solve for: \begin{aligned}\mu_1\ddtau{I(\mu_1)}&=I(\mu_1)-{\wt\over2}\aiaio-{\wt\over4}\mathfrak{F} e^{-\tau\over\mu_0}\\ \mu_2\ddtau{I(\mu_2)}&=I(\mu_2)-{\wt\over2}\aiaio-{\wt\over4}\mathfrak{F} e^{-\tau\over\mu_0}\\ \end{aligned} So we have two coupled, linear, differential equations for $I(\mu_1),I(\mu_2)$. To cut to the chase, these have homogeneous solutions: \begin{aligned}I_{1,hmg}&= {\mathfrak{L}\over1-\mu_1k}e^{k\tau}+{\mathfrak{L}\over1+\mu_1k}e^{-k\tau}\\ I_{2,hmg}&={\mathfrak{L}\over1+\mu_1k}e^{k\tau}+{\mathfrak{L}\over1-\mu_1k}e^{-k\tau}\\ \end{aligned} where $k=\sqrt{1-\wt\over\mu_1^2}$, and $\mathfrak{L}$ is just a number. Compare this to 4.1.17 in the handout, where the only difference is that they pull out $\wt\mathfrak{F}\over4$ in the coefficient. There is also an inhomogeneous solution, which we can figure out by guessing that it looks like the final terms in the coupled differential equations, with some coefficients we need to solve for: $$\begin{matrix} I_{1,part}={\wt\mathfrak{F}\over4}e^{-\tau\over\mu_0}h_1,& I_{2,part}={\wt\mathfrak{F}\over4}e^{-\tau\over\mu_0}h_2\end{matrix}$$ Solving for $h_1, h_2$: $$\begin{matrix} h_1={1-{\mu_1\over\mu_0}\over1-\wt-{\mu_1^2\over\mu_0^2}},& h_2={1+{\mu_1\over\mu_0}\over1-\wt-{\mu_1^2\over\mu_0^2}}\end{matrix}$$