# Monte Carlo Simulations of Radiative Transfer: Overview

FEATURED IMAGE: This is one of the plots that I will be explaining in Post IV. This plot represents the radiation flux emanating from the primary and secondary components of the system considered.

I realize I haven’t posted anything new recently. This is because I have been working on finishing the research project that I worked on during my final semester as an undergraduate.

However, I now intend to share the project on this blog. I will be starting a series of posts in which I will attempt to explain my project and the results I have obtained and open up a dialogue as to any improvements and/or questions that you may have.

Here is the basic overview of the series and projected content:

Post I: Project Overview:

In this first post, I will introduce the project. I will describe the goals of the project and, in particular, I will describe the nature of the system considered.   I will also give a list of references used and mention any and all acknowledgements.

Posts II & III: Radiative Transfer Theory

These posts will most likely be the bulk of the concepts used in the project. I will be defining the quantities used in the research, developing the radiative transfer equation, formulating the problems of radiative transfer, and using the more elementary methods to solving the radiative transfer equation.

Post IV: Monte Carlo Simulations

This post will largely be concerned with the mathematics and physics of Monte Carlo simulations both in the context of probability theory and Ising models. I will describe both and explain how it applies to Monte Carlo simulations of Radiative Transfer.

Post V: Results

I will discuss the results that the simulations produced and explain the data displayed in the posts. I will then discuss the conclusions that I have made and I will also explain the shortcomings of the methods used to produce the results I obtained. I will then open it up to questions or suggestions.

# Solution to Laplace’s Equation

This post deals with the familiar (to the physics student) Laplace’s equation. I am solving this equation in the context of physics, instead of a pure mathematical perspective. This problem is considered most extensively in the context of electrostatics. This equation is usually considered in the spherical polar coordinate system. A lot of finer details are also considered in a mathematical physics course where the topic of spherical harmonics is discussed. Assuming spherical-polar coordinates, Laplace’s equation is

$\displaystyle \frac{1}{r^{2}}\frac{\partial}{\partial r}\bigg\{r^{2}\frac{\partial \psi}{\partial r}\bigg\}+\frac{1}{r^{2}\sin{\theta}}\frac{\partial}{\partial \theta}\bigg\{\sin{\theta}\frac{\partial \psi}{\partial \theta}\bigg\}+\frac{1}{r^{2}\sin{\theta}}\frac{\partial^{2}\psi}{\partial \phi^{2}}=0, (1)$

Suppose that the function $\psi\rightarrow \psi(r,\theta,\phi)$. Furthermore, by the method of separation of variables we suppose that the solution is a product of eigenfunctions of the form $R(r)Y(\theta,\phi)$. Hence, Laplace’s equation becomes

$\displaystyle \frac{Y}{r^{2}}\frac{d}{dr}\bigg\{r^{2}\frac{d^{2}R}{dr^{2}}\bigg\}+\frac{R}{r^{2}\sin{\theta}}\frac{\partial}{\partial \theta}\bigg\{\sin{\theta}\frac{\partial Y}{\partial \theta}\bigg\}+\frac{R}{r^{2}\sin^{\theta}}\bigg\{\frac{\partial^{2}Y}{\partial \phi^{2}}\bigg\}=0. (2)$

Furthermore, we can separate further the term $Y(\theta, \phi)$ into $\Theta(\theta)\Phi(\phi)$. Rewriting (2) and multiplying by $r^{2}\sin^{2}{\theta}R^{-1}Y^{-1}$, we get

$\displaystyle \frac{1}{R}\frac{d}{dr}\bigg\{r^{2}\frac{dR}{dr}\bigg\}+\frac{\sin{\theta}}{\Theta}\frac{d}{d\theta}\bigg\{\sin{\theta}\frac{d\Theta}{d\theta}\bigg\}+\frac{1}{\Phi}\frac{d^{2}\Phi}{d\phi^{2}}=0.(3)$

Bringing the radial and angular component to the other side of the equation and setting the azimuthal component equal to a separation constant $-m^{2}$, yielding

$\displaystyle -\frac{1}{R}\frac{d}{dr}\bigg\{r^{2}\frac{dR}{dr}\bigg\}-\frac{\sin{\theta}}{\Theta}\frac{d}{d\theta}\bigg\{\sin{\theta}\frac{d\Theta}{d\theta}\bigg\}=\frac{1}{\Phi}\frac{d^{2}\Phi}{d\phi^{2}}=-m^{2}.$

Solving the right-hand side of the equation we get

$\displaystyle \Phi(\phi)=A\exp({+im\phi}). (4)$

Now, we set the azimuthal component equal to $m^{2}$ and carry out the derivative in the angular component to get the Associated Legendre equation:

$\displaystyle \frac{d}{d\mu}\bigg\{(1-\mu^{2})\frac{d\Theta}{d\mu}\bigg\}+\bigg\{l(l+1)-\frac{m^{2}}{1-\mu^{2}}\bigg\}\Theta=0, (5)$

where I have let $\mu=\cos{\theta}$. The solutions to this equation are known as the associated Legendre functions. However, instead of solving this difficult equation by brute force methods (i.e. power series method), we consider the case for which $m^{2}=0$. In this case, Eq.(5) simplifies to Legendre’s differential equation discussed previously. Instead of quoting the explicit form of the Legendre polynomials, we equivalently state the Rodrigues formula for these polynomials

$\displaystyle P_{l}(\mu)=\frac{1}{2^{l}l!}\frac{d^{l}}{d\mu^{l}}(\mu^{2}-1)^{l}. (6)$

However, the solutions to Eq.(5) are the associated Legendre functions, not polynomials. Therefore, we use the following to determine the associated Legendre functions from the Legendre polynomials:

$\displaystyle {\Theta_{l}}^{m}(\mu) = (1-\mu^{2})^{\frac{|m|}{2}}\frac{d^{|m|}}{d\mu^{|m|}}P_{l}(\mu). (7)$

For real solutions $|m|\leq l$ and for complex solutions $-l \leq m \leq l$, thus we can write the solutions as series whose indices run from 0 to l for the real-values and from -l to l for the complex-valued solutions.

Now we turn to the radial part which upon substitution of $-l(l+1)$ in place of the angular component, we get

$\displaystyle \frac{d}{dr}\bigg\{r^{2}\frac{dR}{dr}\bigg\}+l(l+1)R=0,$

and if we evaluate the derivatives in the first term we get an Euler equation of the form

$\displaystyle r^{2}\frac{d^{2}R}{dr^{2}}+2r\frac{dR}{dr}+l(l+1)R=0. (8)$

Let us assume that the solution can be represented as

$\displaystyle R(r)=\sum_{j=0}^{\infty}a_{j}r^{j}.$

Taking the necessary derivatives and substituting into Eq.(8) gives

$\displaystyle r^{2}\sum_{j=0}^{\infty}j(j-1)a_{j}r^{j-2}+\sum_{j=0}^{\infty}2rja_{j}r^{j-1}+l(l+1)\sum_{j=0}^{\infty}a_{j}r^{j}=0. (9)$

Simplifying the powers of r we find that

$\displaystyle \sum_{j=0}^{\infty}[j(j-1)+2j+l(l+1)]a_{j}r^{j}=0.$

Now, since $\displaystyle \sum_{j=0}^{\infty}a_{j}r^{j}\neq 0$, this means that

$\displaystyle j(j-1)+2j+l(l+1)=0. (10)$

We can simplify this further to get

$\displaystyle j^{2}+l^{2}+j+l=0.$

Factoring out $(j+l)$, we arrive at

$\displaystyle (j+l)(j+l+1)=0 \implies j=-l ; j=-l-1.$

Since this must be true for all values of l, the solution therefore becomes

$\displaystyle R(r)= \sum_{l=0}^{\infty}Br^{-l}+Cr^{-l-1}. (11)$

Thus, we can write the solution to Laplace’s equation as

$\displaystyle \psi(r,\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=0}^{\infty}[Br^{-l}+Cr^{-l-1}]\Theta_{l}^{m}(\mu)A\exp({im\phi}), (12)$

for real solutions, and

$\displaystyle \psi(r,\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}[Br^{-l}+Cr^{-l-1}]\Theta_{l}^{m}(\mu)A\exp({im\phi}), (13)$

for complex solutions.

# Solution to Legendre’s Differential Equation

Typically covered in a first course on ordinary differential equations, this problem finds applications in the solution of the Schrödinger equation for a one-electron atom (i.e. Hydrogen). In fact, this equation is a smaller problem that results from using separation of variables to solve Laplace’s equation. One finds that the angular equation is satisfied by the Associated Legendre functions. However, if it is assumed that $m=0$ then the equation reduces to Legendre’s equation.

The equation can be stated as

$\displaystyle (1-x^{2})\frac{d^{2}y}{dx^{2}}-2x\frac{dy}{dx}+l(l+1)y(x)=0. (1)$

The power series method starts with the assumption

$\displaystyle y(x)=\sum_{j=0}^{\infty}a_{j}x^{j}. (2)$

Next, we require the first and second order derivatives

$\displaystyle \frac{dy}{dx}=\sum_{j=1}^{\infty}ja_{j}x^{j-1}, (3)$

and

$\displaystyle \frac{d^{2}y}{dx^{2}}=\sum_{j=2}^{\infty}j(j-1)a_{j}x^{j-2}. (4)$

Substitution yields

$\displaystyle (1-x^{2})\sum_{j=2}^{\infty}j(j-1)a_{j}x^{j-2}-2x\sum_{j=1}^{\infty}ja_{j}x^{j-1}+l(l+1)\sum_{j=0}^{\infty}a_{j}x^{j}=0, (5)$

Distribution of the first terms gives

$\displaystyle \sum_{j=2}^{\infty}j(j-1)a_{j}x^{j-2}-\sum_{j=0}^{\infty}j(j-1)a_{j}x^{j}-2ja_{j}x^{j}+l(l+1)a_{j}x^{j}=0, (6)$

where in the second summation term we have rewritten the index to start at zero since the terms for which $j=0,1$ are zero, and hence have no effect on the overall sum. This allows us to write the sum this way. Next, we introduce a dummy variable. Therefore, let $m=j-2\implies j=m+2\implies m+1=j-1$. Thus, the equation becomes

$\displaystyle \sum_{j=0}^{\infty}\bigg\{(j+2)(j+1)a_{j+2}-j(j-1)a_{j}-2ja_{j}+l(l+1)a_{j}\bigg\}x^{j}=0. (7)$

In order for this to be true for all values of j, we require the coefficients of $x^{j}$ equal zero. Solving for $a_{j+2}$ we get

$\displaystyle a_{j+2}=\frac{j(j+1)-l(l+1)}{(j+2)(j+1)}. (8)$

It becomes evident that the terms $a_{2},a_{3},a_{4},...$ are dependent on the terms $a_{0}$ and $a_{1}$. The first term deals with the even solution and the second deals with the odd solution. If we let $p=j+2$ and solve for $a_{p}$, we arrive at the term $a_{p-2}$ and we can obtain the next term $a_{p-4}$. (I am not going to go through the details. The derivation is far too tedious. If one cannot follow there is an excellent video on YouTube that goes through a complete solution of Legendre’s ODE where they discuss all finer details of the problem. I am solving this now so that I can solve more advanced problems later on.) A pattern begins to emerge which we may express generally as:

$\displaystyle a_{p-2n}=\frac{(-1)^{n}(2p-2n)!}{(n!)(2^{p})(p-n)!(p-2n)!}. (9)$

Now, for even terms $0\leq n \leq \frac{p}{2}$ and for odd terms $0 \leq n \leq \frac{p-1}{2}$. Thus for the even solution we have

$\displaystyle P_{p}(x)=\sum_{n=0}^{\frac{p}{2}}\frac{(-1)^{n}(2p-2n)!}{(n!)(2^{p})(p-n)!(p-2n)!}, (10.1)$

and for the odd solution we have

$\displaystyle P_{p}(x)=\sum_{n=0}^{\frac{p-1}{2}}\frac{(-1)^{n}(2p-2n)!}{(n!)(2^{p})(p-n)!(p-2n)!}. (10.2)$

These two equations make up the even and odd solution to Legendre’s equation. They are an explicit general formula for the Legendre polynomials.  Additionally we see that they can readily be used to derive Rodrigues’ formula

$\displaystyle \frac{1}{2^{p}p!}\frac{d^{p}}{dx^{p}}\bigg\{(x^{2}-1)^{p}\bigg\}, (11)$

and that we can relate Legendre polynomials to the Associated Legendre function via the equation

$\displaystyle P_{l}^{m}(x)= (1-x^{2})^{\frac{|m|}{2}}\frac{d^{|m|}P_{l}(x)}{dx^{|m|}}, (12)$

where I have let $p=l$ so as to preserve a more standard notation.  This is the general rule that we will use to solve the associated Legendre differential equation when solving the Schrödinger equation for a one-electron atom.