Solution to the Hermite Differential Equation

One typically finds the Hermite differential equation in the context of an infinite square well potential and the consequential solution of the Schrödinger equation. However, I will consider this equation is its “raw” mathematical form viz.

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

First we will consider the more general case, leaving \lambda undefined. The second case will consider in a future post \lambda = 2n, n\in \mathbb{Z}^{+}, where \mathbb{Z}^{+}=\bigg\{x\in\mathbb{Z}|x > 0\bigg\}.

PART I: 

Let us assume the solution has the form

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

Now we take the necessary derivatives

\displaystyle y^{\prime}(x)=\sum_{j=1}^{\infty}ja_{j}x^{j-1}, (3)

\displaystyle y^{\prime \prime}(x)=\sum_{j=2}^{\infty} j(j-1)a_{j}x^{j-2}, (4)

where upon substitution yields the following

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

Introducing the dummy variable m=j-2 and using this and its variants we arrive at

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

Bringing this under one summation sign…

\displaystyle \sum_{j=0}^{\infty}[(j+2)(j+1)a_{j+2}-2ja_{j}+\lambda a_{j}]x^{j}=0. (7)

Since \displaystyle \sum_{j=0}^{\infty}x^{j}\neq 0, we therefore require that

\displaystyle (j+2)(j+1)a_{j+2}=(2j - \lambda)a_{j}, (8)

or

\displaystyle a_{j+2}=\frac{(2j-\lambda)a_{j}}{(j+2)(j+1)}. (9)

This is our recurrence relation. If we let j=0,1,2,3,... we arrive at two linearly independent solutions (one even and one odd) in terms of the fundamental coefficients a_{0} and a_{1} which may be written as

\displaystyle y_{even}(x)= a_{0}\bigg\{1+\sum_{j=0}^{j/2}\frac{(-1)^{j}(\lambda -2j)!}{(j+2)!}x^{j}\bigg\}, (10)

and

\displaystyle y_{odd}(x)=a_{1}\bigg\{\sum_{j=0}^{(j-1)/2}\frac{(-1)^{j}(\lambda-2j)!}{(j+2)!}x^{j}\bigg\}. (11)

Thus, our final solution is the following

\displaystyle y(x)=y_{even}(x)+y_{odd}(x), (12.1)

 

\displaystyle y(x)=a_{0}\bigg\{1+\sum_{j=0}^{j/2}\frac{(-1)^{j}(\lambda-2j)!}{(j+2)!}x^{j+2}\bigg\}+a_{1}\bigg\{x+\sum_{j=1}^{(j-1)/2}\frac{(-1)^{j}(\lambda-2j)!}{(j+2)!}x^{j+2}\bigg\}. (12.2)

 

 

 

Legendre Polynomials

Some time ago, I wrote a post discussing the solution to Legendre’s ODE. In that post, I discussed what is an alternative definition of Legendre polynomials in which I stated Rodriguez’s formula:

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

where

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

and

\displaystyle P_{p}(x)=\sum_{n=0}^{\beta}\frac{(-1)^{n}(2p-2n)!}{2^{p}{n!}(p-n)!(p-2n)!} (0.3)

in which I have let \displaystyle \alpha=p/2 and \displaystyle \beta=(p-1)/2 corresponding to the even and odd expressions for the Legendre polynomials.

However, in this post I shall be using the approach of the generating function. This will be from a purely mathematical perspective, so I am not applying this to any particular topic of physics.

Consider a triangle with sides \displaystyle X,Y.Z and angles \displaystyle \theta, \phi, \lambda. The law of cosines therefore maintains that

\displaystyle Z^{2}=X^{2}+Y^{2}-2XY\cos{(\lambda)}. (1)

We can factor out \displaystyle X^{2} from the left hand side of Eq.(1), take the square root and invert this yielding

\displaystyle \frac{1}{Z}=\frac{1}{X}\bigg\{1+\bigg\{\frac{Y}{X}\bigg\}^{2}-2\bigg\{\frac{Y}{X}\bigg\}\cos{(\lambda)}\bigg\}^{-1/2}. (2)

Now, we can expand this by means of the binomial expansion. Let \displaystyle \kappa \equiv \bigg\{\frac{Y}{X}\bigg\}^{2}-2\bigg\{\frac{Y}{X}\bigg\}\cos{(\lambda)}, therefore the binomial expansion is

\displaystyle \frac{1}{(1+\kappa)^{1/2}}=1-\frac{1}{2}\kappa+\frac{3}{8}\kappa^{2}-\frac{5}{16}\kappa^{3}+... (3)

Hence if we expand this in terms of the sides and angle(s) of the triangle and group by powers of \displaystyle (y/x) we get

\displaystyle \frac{1}{Z}=\frac{1}{X}\bigg\{1+\bigg\{\frac{Y}{X}\bigg\}\cos{(\lambda)}+\bigg\{\frac{Y}{X}\bigg\}^{2}\frac{1}{2}(3\cos^{2}{(\lambda)}-1)+\bigg\{\frac{Y}{X}\bigg\}^{3}\frac{1}{2}(5\cos^{3}{(\lambda)}-3\cos{(\lambda)}\bigg\}.(4)

Notice the coefficients, these are precisely the expressions for the Legendre polynomials. Therefore, we see that

\displaystyle \frac{1}{Z}=\frac{1}{X}\bigg\{\sum_{l=0}^{\infty}\bigg\{\frac{Y}{X}\bigg\}^{l}P_{l}(\cos{(\lambda)}\bigg\}, (5)

or

\displaystyle \frac{1}{Z}=\frac{1}{\sqrt[]{X^{2}+Y^{2}-2XY\cos{(\lambda)}}}=\sum_{l=0}^{\infty}\frac{Y^{l}}{X^{l+1}}P_{l}(\cos{(\lambda)}. (6)

Thus we see that the generating function \displaystyle 1/Z generates the Legendre polynomials. Two prominent uses of these polynomials includes gravity and its application to the theory of potentials of a spherical mass distributions, and the other is that of electrostatics. For example, suppose we have the potential equation

\displaystyle V(r)=\frac{1}{4\pi\epsilon_{0}}\int_{V}\rho(R)\frac{\hat{\mathcal{R}}}{\mathcal{R_{0}}}d\tau. (7.1)

We may use the result of the generating function to get the following result for the electric potential due an arbitrary charge distribution

\displaystyle V(\mathcal{R})=\frac{1}{4\pi\epsilon_{0}}\sum_{l=0}^{\infty}\frac{\mathcal{R}^{l}}{\mathcal{R_{0}}^{l+1}}\int P_{l}(\cos{(\lambda)}). (7.2)

(For more details, see Chapter 3 of Griffith’s text: Introduction to Electrodynamics.)

 

Basics of Tensor Calculus and General Relativity-Vectors and Introduction to Tensors (Part I: Vectors)

SOURCE FOR CONTENT: Neuenschwander D.E.,2015. Tensor Calculus for Physics. Johns Hopkins University Press. 

At some level, we all are aware of scalars and vectors, but typically we don’t think of aspects of everyday experience as being a scalar or a vector. A scalar is something that has only magnitude, that is it only has a numeric value. A typical example of a scalar would be temperature. A vector, on the other hand, is something that has both a magnitude and direction. This could be something as simple as displacement. If we wish to move a certain distance with respect to our current location we must specify how far to go and in which direction to move. Other examples, (and there are a lot of them), including velocity, force, momentum, etc. Now, tensors are something else entirely. In Neuenschwander’s “Tensor Calculus for Physics”, he recites the rather unsatisfying definition of a tensor as being

” ‘A set of quantities T_{s}^{r} associated with a point P are said to be components of a second-order tensor if, under a change of coordinates, from a set of coordinates x^{s} to x^{\prime s}, they transform according to

\displaystyle T^{\prime r}_{s}=\frac{\partial x^{\prime r}}{\partial x^{i}}\frac{\partial x^{j}}{\partial x^{\prime s}}T_{j}^{i}. (1)

where the derivatives are evaluated at the aforementioned point.’ “

Neuenschwander describes his frustration when encountered this so-called definition of a tensor. Like him, I found I had similar frustrations, and as a result, I had even more questions.

We shall start with a discussion of vectors, for an understanding of these quantities are an integral part of tensor analysis.

We define a vector as a quantity that has 3 distinct components and an angle that indicates orientation or direction. There are two types of vectors; those with coordinates and those without. I will be discussing the latter first, but I will consider Neuenschwander’s description in the context of the definition of a vector space. Then I shall move on to the former. Consider a number of arbitrary vectors \displaystyle \textbf{U}, \textbf{V}, \textbf{W},etc. If we consider further more and more vectors, we can therefore imagine a space constituted by these vectors; a vector space. To make it formal, here is the definition:

Def. Vector Space:

A vector space \mathcal{S} is the nonempty set of elements (vectors) that satisfy the following axioms

\displaystyle 1.\textbf{U}+\textbf{V}=\textbf{V}+\textbf{U};  \forall \textbf{U},\textbf{V}\in \mathcal{S},

\displaystyle 2. (\textbf{U}+\textbf{V})+\textbf{W}= \textbf{U}+(\textbf{V}+\textbf{W}); \forall \textbf{U},\textbf{V},\textbf{W}\in \mathcal{S},

\displaystyle 3. \exists  0 \in \mathcal{S}|\textbf{U}+0=\textbf{U},

\displaystyle 4. \forall \textbf{U}\in \mathcal{S}, \exists -\textbf{U}\in \mathcal{S}|\textbf{U}+(-\textbf{U})=0,

\displaystyle 5. \alpha(\textbf{U}+\textbf{V})=\alpha\textbf{U}+\alpha\textbf{V}, \forall \alpha \in \mathbb{R},  \forall \textbf{U},\textbf{V}\in \mathcal{S}.

\displaystyle 6. (\alpha+\beta) \textbf{U}= \alpha\textbf{U}+\beta\textbf{U}, \forall \alpha,\beta \in \mathbb{R},

\displaystyle 7. (\alpha\beta)\textbf{U}= \alpha(\beta\textbf{U}), \forall \alpha,\beta \in \mathbb{R},

\displaystyle 8. 1\textbf{U}=\textbf{U}, \forall \textbf{U}\in \mathcal{S},

and satisfies the following closure properties:

 1. If \textbf{U}\in \mathcal{S}, and \alpha is a scalar, then \alpha\textbf{U}\in \mathcal{S}.

2. If \textbf{U},\textbf{V}\in \mathcal{S}, then \textbf{U}+\textbf{V}\in \mathcal{S}.

The first closure property ensures closure under scalar multiplication while the second ensures closure under addition.

In rectangular coordinates, our arbitrary vector may be represented by basis vectors \hat{i}, \hat{j},\hat{k} in the following manner

\displaystyle \textbf{U}=u\hat{i}+v\hat{j}+w\hat{k}, (1)

where the basis vectors have the properties

\displaystyle \hat{i}\cdot \hat{i}=\hat{j}\cdot \hat{j}=\hat{k}\cdot \hat{k}=1, (2.1)

\displaystyle \hat{i}\cdot \hat{j}=\hat{j}\cdot \hat{k}=\hat{i}\cdot \hat{k}=0.(2.2)

The latter of which implies that these basis vectors are mutually orthogonal. We can therefore write these in a more succinct way via

\displaystyle \hat{e}_{i}\cdot \hat{e}_{j}=\delta_{ij},  (2.3)

where \displaystyle \delta_{ij} denotes the Kronecker delta:

\displaystyle \delta_{ij}=\begin{cases} 1, i=j\\ 0, i \neq j \end{cases}. (2.4)

We may redefine the scalar product by the following argument given in Neuenschwander

\displaystyle \textbf{U}\cdot  \textbf{V}=\sum_{i,j=1}^{3}(U^{i} \hat{e}_{i}) \cdot (V^{j} \hat{e}_{j})=\sum_{i,j=1}^{3}U^{i}V^{j}(\delta_{ij})=\sum_{l=1}^{3}U^{l}V^{l}. (3)

Similarly we may define the cross product to be

\displaystyle \textbf{U}\times \textbf{V}=\det\begin{pmatrix} \hat{i}  & \hat{j} & \hat{k} \\ U^{x} & U^{y} & U^{z} \\ V^{x} & V^{y} & V^{z}  \end{pmatrix}, (4)

whose i-th component is

\displaystyle (\textbf{U}\times \textbf{V})^{i}=\textbf{U}\times \textbf{V}=\sum_{i,j=1}^{3}\epsilon^{ijk}U^{j}V^{k}, (5)

where \epsilon^{ijk} denotes the Levi-Civita symbol. If these indices form an odd permutation \epsilon^{ijk}=-1, if the indices form an even permutation \epsilon^{ijk}=+1, and if any of these indices are equal \epsilon^{ijk}=0.

As a final point, we may relate vectors to relativity by means of defining the four-vector. If we consider the four coordinates x,y,z,t, they collectively describe what is referred to as an event. Formally, an event in spacetime is described by three spatial coordinates and one time coordinate. We may replace these coordinates by x^{\mu}, where \mu \in \mathbb{Z}^{\pm} in which I am defining \mathbb{Z}^{\pm}=\bigg\{x\in \mathbb{Z}|x\geq 0\bigg\}. In words, this means that the index \mu is an integer that is greater than or equal to 0.

Furthermore, the quantity x^{0} corresponds to time, x^{1,2,3} correspond to the x,y, and z coordinates respectively.  Therefore, x^{\mu}=(ct,x,y,z) is called the four-position.

In the next post, I will complete the discussion on vectors and discuss in more detail the definition of a tensor (following Neuenschwander’s approach). I will also introduce a few examples of tensors that physics students will typically encounter.



 

 

 

 

 

Elementary Methods to Solving Integral Equations: Overview of Series

After taking a differential equations course, I reasoned that if there are differential equations, then there must be an “inverse theory” of integral equations. I was able to find a text that discussed the basics, but interestingly videos discussing this topic were scarce. So I thought I would do a (short) series on some basic methods towards solving these types of equations and perhaps apply them to a physical problem.

This series will discuss two types of integral equations:

Fredholm Equations

Volterra Equations

I will also discuss…

Symmetric Kernels and Orthogonal Systems of Functions.

Note & Status: 

This series is meant to be a shorter, supplementary series. My main focus right now is the Monte Carlo Simulations of Radiative Transfer series and once that is finished, I will move on to the Basics of Tensor Calculus series. I am doing this series on integral equations out of pure mathematical curiosity. I should have a couple of new posts up in the next few days.

Basics of Tensor Calculus and General Relativity: Overview of Series

Of the many topics that I have studied in the past, one of the most confusing topics that I have encountered is the concept of a tensor. Every time that I heard the term tensor, it was mentioned for the sake of mentioning it and my professors would move on. Others simply ignored the existence of tensors, but every time they were brought up I was intrigued. So the purpose of this series is to attempt to discover how tensors work and how they relate to our understanding of the universe, specifically in the context of general relativity.

The text I will be following for this will be Dwight E. Neuenschwander’s “Tensor Calculus for Physics”.  I will be documenting my interpretation of the theory discussed in this text, and I will use the text “General Relativity: An Introduction for Physicists” by M.P. Hobson, G. Efstathiou, and A.N. Lasenby to discuss the concepts developed in the context of general relativity. I will list the corresponding chapter titles associated with each post.

Here is an approximate agenda (note that this will take some time as I am learning this as I post, so posts in this series may be irregular).

Post I: Tensor Calculus: Introduction and Vectors

Post II: Tensor Calculus: Vector Calculus and Coordinate Transformations

Post III: General Relativity: Basics of Manifolds and Coordinates

Post IV: General Relativity: Vector Calculus on Manifolds

Post V: Tensor Calculus: Introduction to Tensors and the Metric Tensor

Post VI: Tensor Calculus: Derivatives of Tensors

Post VII: General Relativity: Tensor Calculus on Manifolds

Post VIII: Tensor Calculus: Curvature

Post IX: General Relativity: Equivalence Principle and Spacetime Curvature

Post X: General Relativity: The Gravitational Field Equations

The series will end with the Einstein field equations (EFEs) since they are the crux of the general theory of relativity. The posts in this series will come in their own time as they can be quite difficult and time-consuming, but I will do my best to understand it. I welcome any feedback you may have, but please be respectful.

 

Monte Carlo Simulations of Radiative Transfer: Basics of Radiative Transfer Theory (Part I)

SOURCE FOR CONTENT: Chandrasekhar, S., 1960. Radiative Transfer. 1. 

 

In this post, I will be discussing the basics of radiative transfer theory necessary to understand the methods used in this project. I will start with some definitions, then I will look at the radiative transfer equation and consider two simple cases of scattering.

The first definition we require is the specific intensity, which is the amount of energy associated with a specific frequency dE_{\nu} passing through an area dA constrained to a solid angle d\Omega in a time dt. We may write this mathematically as

dE_{\nu}=I_{\nu}\cos{\theta}d\nu d\Sigma d\Omega dt. (1)

We must also consider the net flux given by

\displaystyle d\nu d\Sigma dt \int I_{\nu}\cos{\theta}d\Omega, (2)

where if we integrate over all solid angles \Omega we get

\pi F_{\nu}=\displaystyle \int I_{\nu}\cos{\theta}d\Omega. (3)

Let d\Lambda be an element of the surface \Lambda in a volume V through which radiation passes. Further let \Theta and \theta denote the angles which form normals with respect to elements d\Lambda and d\Sigma. These surfaces are joined by these normals and hence we have the surface across which energy flows  includes the elements d\Lambda and d\Sigma, given by the following:

I_{\nu}\cos{\Theta}d\Sigma d\Omega^{\prime}d\nu = I_{\nu}d\nu \frac{\cos{\Theta}\cos{\theta}d\Sigma d\Lambda}{r^{2}} (4),

where d\Omega^{\prime}=d\Lambda \cos{\Theta}/r^{2} is the solid angle subtended by the surface element d\Lambda at a point P and volume element dV=ld\Sigma \cos{\theta} is the volume that is intercepted in volume V. If we take this further, and integrate over all V and \Omega we arrive at

\displaystyle \frac{d\nu}{c}\int dV \int I_{\nu} d\Omega=\frac{V}{c}d\nu \int I_{\nu}d\Omega, (5)

where if the radiation travels some distance L in the volume, then we must multiply Eq.(5) by l/c, where c is the speed of light.

We now define the integrated energy density as being

U_{\nu}=\displaystyle \frac{1}{c}\int I_{\nu}d\Omega, (6.1)

while the average intensity is

J_{\nu}=\displaystyle \frac{1}{4\pi}\int I_{\nu}d\nu, (6.2)

and the relation between these two equations is

U_{\nu}=\frac{4\pi}{c}J_{\nu}. (6.3)

I will now introduce the radiative transfer equation. This equation is a balance between the amount of radiation absorbed and the radiation that is emitted. The equation is,

\frac{dI_{\nu}}{ds}=-\epsilon \rho I_{\nu}+h_{\nu}\rho, (7)

where if we divide by \epsilon \rho we get

-\frac{1}{\epsilon_{\nu}\rho}\frac{dI_{\nu}}{ds}=I_{\nu}+U_{\nu}(\theta, \phi), (8)

where U(\theta,\phi) represents the source function given by

U_{\nu}(\theta,\phi)=\displaystyle \frac{1}{4\pi}\int_{0}^{\pi}\int_{0}^{2\pi}p(\theta,\phi;\theta^{\prime},\phi^{\prime})I_{\nu}\sin{\theta^{\prime}}d\theta^{\prime}d\phi^{\prime}. (9)

The source function is typically the ratio between the absorption and emission coefficients. One of the terms in the source function is the phase function which varies according to the specific scattering geometry. In its most general form, we can represent the phase function as an expansion of Legendre polynomials:

p(\theta, \phi; \theta^{\prime},\phi^{\prime})=\displaystyle \sum_{j=0}^{\infty}\gamma_{j}P_{j}(\mu), (10)

where we have let \mu = \cos{\theta} (in keeping with our notation in previous posts).

In Part II, we will discuss a few simple cases of scattering and their corresponding phase functions, as well as obtaining the formal solution of the radiative transfer equation. (DISCLAIMER: While this solution will be consistent in a mathematical sense, it is not exactly an insightful solution since much of the more interesting and complex cases involve the solution of either integro-differential equations or pure integral equations (a possible new topic).)

 

Simple Harmonic Oscillators (SHOs) (Part I)

We all experience or see this happening in our everyday experience: objects moving back and forth. In physics, these objects are called simple harmonic oscillators. While I was taking my undergraduate physics course, one of my favorite topics was SHOs because of the way the mathematics and physics work in tandem to explain something we see everyday. The purpose of this post is to engage followers to get them to think about this phenomenon in a more critical manner.

Every object has a position at which these objects tend to remain at rest, and if they are subjected to some perturbation, that object will oscillate about this equilibrium point until they resume their state of rest. If we pull or push an object with an applied force F_{A} we find that this force is proportional to Hooke’s law of elasticity, that is, F_{A}=-k\textbf{r}. If we consider other forces we also find that there exists a force balance between the restoring force (our applied force), a resistance force, and a forcing function, which we assume to have the form

F=F_{forcing}+F_{A}-F_{R}= -k\textbf{r}-\beta \dot{\textbf{r}}; (1)

note that we are assuming that the resistance force is proportional to the speed of an object. Suppose further that we are inducing these oscillations in a periodic manner by given by

F_{forcing}=F_{0}\cos{\omega t}. (2)

Now, to be more precise, we really should define the position vector. So, \textbf{r}=x\hat{i}+y\hat{j}+z\hat{k}. Therefore, we actually have a system of three second order linear non-homogeneous ordinary differential equations in three variables:

m\ddot{ x}+\beta \dot{x}+kx=F_{0}\cos{\omega t}, (3.1)

m\ddot{y}+\beta \dot{y}+ky=F_{0}\cos{\omega t}, (3.2)

m\ddot{z}+\beta \dot{z}+kz=F_{0}\cos{\omega t}. (3.3)

(QUICK NOTE: In the above equations, I am using the Newtonian notation for derivatives, only for convenience.)  I will just make some simplifications. I will divide both sides by the mass, and I will define the following parameters: \gamma \equiv \beta/m, \omega_{0} \equiv k/m, and \alpha \equiv F_{0}/m. Furthermore, I am only going to consider the y component of this system. Thus, the equation that we seek to solve is

\ddot{y}+\gamma \dot{y}+\omega_{0}y=\alpha\cos{\omega t}. (4)

Now, in order to solve this non-homogeneous equation, we use the method of undetermined coefficients. By this we mean to say that the general solution to the non-homogeneous equation is of the form

y = Ay_{1}(t)+By_{2}(t)+Y(t), (5)

where Y(t) is the particular solution to the non-homogeneous equation and the other two terms are the fundamental solutions of the homogeneous equation:

\ddot{y}_{h}+\gamma \dot{y}_{h}+\omega_{0} y_{h} = 0. (6)

Let y_{h}(t)=D\exp{(\lambda t)}. Taking the first and second time derivatives, we get \dot{y}_{h}(t)=\lambda D\exp{(\lambda t)} and \ddot{y}_{h}(t)=\lambda^{2}D\exp{(\lambda t)}. Therefore, Eq. (6) becomes, after factoring out the exponential term,

D\exp{(\lambda t)}[\lambda^{2}+\gamma \lambda +\omega_{0}]=0.  (7)

Since D\exp{(\lambda t)}\neq 0, it follows that

\lambda^{2}+\gamma \lambda +\omega_{0}=0. (8)

This is just a disguised form of a quadratic equation whose solution is obtained by the quadratic formula:

\lambda =\frac{-\gamma \pm \sqrt[]{\gamma^{2}-4\omega_{0}}}{2}. (9)

Part II of this post will discuss the three distinct cases for which the discriminant \sqrt[]{\gamma^{2}-4\omega_{0}} is greater than, equal to , or less than 0, and the consequent solutions. I will also obtain the solution to the non-homogeneous equation in that post as well.

 

 

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.