Category Archives: Uncategorized

Solution to the three-dimensional Heat Equation

After taking a topics course in applied mathematics (partial differential equations), I found that there were equations that I should solve since I would later see those equations embedded into other larger-scale equations. This equation was Laplace’s equation (future post). Once I solved this equation, I realized that it becomes a differential operator when acted upon a function of at least two variables. Thus, I could solve equations such as the Schrödinger equation using a three-dimensional laplacian in spherical-polar coordinates (another future post) and the three-dimensional heat equation. I will be solving the latter.


The heat equation initial-boundary-value-problem is therefore

\frac{\partial u}{\partial t}=\frac{1}{c\rho}\nabla^{2}u+Q(x,y,z), (1.1)

subjected to the boundary conditions

u(0,0,0,t)=0, (1.2)

and the initial condition

u(x,y,z,0)= \xi(x,y,z). (1.8)

Now, this solution is not specific to a single thermodynamic system, but rather it is a more general solution in a mathematical context. However, I will be appropriating certain concepts from physics for reasons that are well understood (i.e. that time exists on the interval I=0\leq t < +\infty ). It is nonphysical or nonsensical to speak of negative time.

To start, consider any rectangular prism in which heat flows through the volume, from the origin to the point (X,Y,Z).At a time t=0, the overall heat of the volume can be regarded to be a function \xi(x,y,z). After a time period \delta t=t has passed, the heat will have traversed to the point (X,Y,Z) from the origin (think of the heat traveling along the diagonal of a cube). The goal is to find the heat as a function of the three spatial components and a single time component.  Furthermore, the boundary conditions maintain that at the origin and the final point,  the heat vanishes. In other words, these points act as heat sinks. (A sink is a point where energy can leave the system.)

To simplify the notation, we define \textbf{r}=x\hat{i}+y\hat{j}+z\hat{k}. Thus the initial-boundary-value-problem becomes

\frac{\partial u}{\partial t}=\frac{1}{c\rho}\frac{\partial^{2}u}{\partial \textbf{r}^{2}}, (2.1)

where u\rightarrow u(\textbf{r},t). Also, this definition also reduces the three-dimensional laplacian to a second-order partial derivative of u. The boundary and initial conditions are then


u(R,t)=0, (2.3)

u(\textbf{r},0)=\xi(\textbf{r}). (2.2)

Now, we assume that the solution is a product of eigenfunctions of the form

u(\textbf{r},t)=\alpha(\textbf{r})\Gamma(t). (3)

Taking the respective derivatives and dividing by the assumed form of the solution, we get

\frac{\Gamma^{\prime}(t)}{\Gamma(t)}=\frac{1}{c\rho}\frac{\alpha^{\prime\prime}(\textbf{r})}{\alpha(\textbf{r})}+Q(\textbf{r}). (4)

Now, Eq.(4) can be equal to three different values, the first of which is zero, but this solution does not help in any way,  nor is it physically significant since it produces a trivial solution. The second is \lambda^{2}>0, \lambda \neq 0, so we can apply to the time dependence equation which then becomes

\frac{\Gamma^{\prime}(t)}{\Gamma(t)}=\lambda^{2}, (5.1)

whose solution is

\Gamma(t)=\Gamma_{0}\exp({\lambda^{2}t}). (5.2)

The third case where \lambda^{2}<0, \lambda \neq 0 allows us to write the spatial equation upon rearrangement as

\alpha^{\prime\prime}(\textbf{r})+c\rho(\lambda^{2}+Q)\alpha(\textbf{r})=0, (6.1)

whose solution is

\alpha(\textbf{r})=c_{1}\cos({\sqrt[]{c\rho(\lambda^{2}+Q)}\textbf{r}})+c_{3}\sin({\sqrt[]{c\rho(\lambda^{2}+Q)}\textbf{r}}), (6.2)

where c_{3}=-c_{2}. Next we apply the boundary conditions. Let \textbf{r}=0 in \alpha(\textbf{r}):

\alpha(0)=c_{1}+0=0 \implies c_{1}=0 \implies \alpha(\textbf{r})=C\sin(\sqrt[]{c\rho(\lambda^{2}+Q)} \textbf{r}). (7.1)

and let \textbf{r}=R in (7.1) to get

\alpha(R)\implies \sin({\sqrt[]{c\rho(\lambda^{2}+Q)}R})=0\implies \sqrt[]{c\rho(\lambda^{2}+Q)}R^{2}=(n\pi)^{2}.

Solving for \lambda^{2} above gives

\lambda^{2}=\frac{1}{c\rho}\bigg\{\frac{n\pi}{R}\bigg\}^{2}-Q. (8)

Substituting into  Eq.(7.1) and simplifying gives

\alpha(\textbf{r})=C\sin\bigg\{\frac{n\pi}{c\rho}\bigg\}^{2}\bigg\{\frac{\textbf{r}}{R}\bigg\}. (9)

To find as many solutions as possible we construct a superposition of solutions of the form

u(\textbf{r},t)=u_{l}(\textbf{r},t)=\sum_{l=0}^{\infty}\alpha_{l}(\textbf{r})\Gamma(t). (10)

Rewriting \alpha(\textbf{r}) as \alpha_{l}(\textbf{r}) and using the solution for the time dependence, and also let the coefficients form a product equivalent to the indexed coefficient A_{l}, we arrive at the solution for the heat equation:

u(\textbf{r},t)=\sum_{l=0}^{\infty}A_{l}\sin\bigg\{\frac{l\pi}{c\rho}\bigg\}^{2}\bigg\{\frac{\textbf{r}}{R}\bigg\}\exp(\lambda^{2}t). (11)

Now suppose that t=0 in Eq.(11). In this case we get the initial heat distribution given by \xi(\textbf{r}):

\xi(\textbf{r})=\sum_{l=0}^{\infty} A_{l}\sin\bigg\{\frac{l\pi}{c\rho}\bigg\}^{2}\bigg\{\frac{\textbf{r}}{R}\bigg\}. (12)


Electron Scattering of a Step Potential

The following post was initially one of my assignments for an independent study in modern physics in my penultimate year as an undergraduate. While studying this problem the text that I used to verify my answer was:

R. Eisberg and R. Resnick, Quantum Physics of Atoms, Molecules, Solids, Nuclei, and Particles. John Wiley & Sons. 1985. 6. 

One of the hallmarks of quantum theory is the Schrödinger equation. There are two forms: the time-dependent and the time-independent equation. The former can be turned into the latter by way of assuming stationary states in which case there is no time evolution (i.e. the time derivative \partial_{t} \Psi(x,t)=0. \partial_{t}\Psi(x,t) denotes the derivative of the wavefunction \Psi(x,t) with respect to time.)

The wavefunction describes the state of a system, and it is found by solving the Schrödinger equation. In this post, I’ll be considering a step potential in which

V(x) = V_{0},  x < 0,  (1.1)

V(x) = 0, x > 0.  (1.2)

After solving for the wavefunction, I will calculate the reflection and transmission coefficients for the case where the energy of the electron is less than that of the step potential E < V_{0}.

First, we assume that we are dealing with stationary states, by doing so we assume that there is no time-dependence. The wavefunction \Psi(x,t) becomes an eigenfunction (eigen– is German for “characteristic” e.g. characteristic function, characteristic value (eigenvalue), and so on), \psi(x). There are requirements for this eigenfunction in the context of quantum mechanics: eigenfunction \psi(x) and its first order spatial derivative \frac{d\psi(x)}{dx} must be finite, single-valued, and continuous. Using the wavefunction, we write the time-independent Schrödinger equation as

-\frac{\hbar^{2}}{2m}\frac{d^{2}\psi(x)}{dx^{2}}+V(x)\psi(x) = E\psi(x), (2)

where \hbar is the reduced Planck’s constant \hbar \equiv h(2\pi)^{-1}, m is the mass of the particle, V(x) represents the potential, and E is the energy.

Now, in electron scattering there are two cases regarding a step potential: the case for which E < V_{0} and the case for which E>V_{0} The focus of this post is the former case. In such a case, the potential is given mathematically by Eqs. (1.1) and (1.2) and can be depicted by the image below:


Image Credit:

The first part of this problem is to solve for \psi(x) when x<0. Then Eq.(2) becomes

\frac{d^{2}\psi_{I}(x)}{dx^{2}}=-\kappa_{I}^{2}\psi_{I}(x), (3)

where \kappa_{I}^{2}\equiv \frac{-2mE}{\hbar^{2}}. Eq.(3) is a second order linear homogeneous ordinary differential equation with constant coefficients and can be solved using a characteristic equation. We assume that the solution is of the general form


which upon taking the first and second order spatial derivatives and substituting into Eq.(3) yields


We can factor out the exponential and recalling that the graph of e^{x} \neq 0 (except in the limit x\rightarrow -\infty)we can then conclude that


Hence, r = \pm i\kappa_{I}. Therefore, we can write the solution Schrödinger equation in the region x < 0 as

\psi_{I}(x) = A\exp{+i\kappa_{I}x}+B\exp{-i\kappa_{I}x}. (4)

This is the eigenfunction for the first region. Coefficients A and B will be determined later.

Now we can use the same logic for the Schrödinger equation in the region x > 0 where V(x)=V_{0}:

\frac{d^{2}\psi_{II}(x)}{dx^{2}}= -\kappa_{II}^{2}\psi_{II}(x), (5)

where \kappa_{II}\equiv \frac{\sqrt[]{2m(E-V_{0})}}{\hbar}. The general solution for this region is

\psi_{II}(x) = C\exp{+i\kappa_{II}x}+D\exp{-i\kappa_{II}x}. (6)

Now the next step is taken using two different approaches: the first using a mathematical argument and the other from a conceptual interpretation of the problem at hand.  The former is this: Suppose we let x\rightarrow \infty. What results is that the first term on the right hand side of Eq.(6) diverges (i.e. becomes arbitrarily large). The second term on the right hand side ends up going to zero (it converges). Therefore, D remains finite, hence D\neq 0. However, in order to suppress the divergence of the first term, we let C=0. Thus we arrive at the eigenfunction for x>0

\psi_{II}(x) = D\exp{-i\kappa_{II}x}. (7)

The latter argument is this: In the region x < 0, the first term of the solution represents a wave propagating in the positive x-direction, while the second denotes a wave traveling in the negative x-direction. Similarly, in the region where x > 0, the first term corresponds to a wave traveling in the positive x-direction. However, this cannot be because the energy of the wave is not sufficient enough to overcome the potential. Therefore, the only term that is relevant here for this region is the second term, for it is a wave propagating in the negative x-direction.

We now determine the coefficients A and B. Recall that the eigenfunction must satisfy the following continuity requirements

\psi_{I}(x)=\psi_{II}(x), (8.1)

\psi^{\prime}_{I}(x)=\psi^{\prime}_{II}(x), (8.2)

evaluated when x=0. Doing so in Eqs. (4) and (7), and equating them we arrive at the continuity condition for \psi(x)

A+B = D. (9)

Taking the derivative of \psi_{I}(x) and \psi_{II}(x) and evaluating them for when x=0 we arrive at

A-B = \frac{i\kappa_{II}}{i\kappa_{I}}D. (10)

If we add Eqs. (9) and (10) we get the value for the coefficient A in terms of the arbitrary constant D

A =\frac{D}{2}(1+\frac{i\kappa_{II}}{i\kappa_{I}}. (11)

Conversely, if we subtract (9) and (10) we get the value for B in terms of D

B = \frac{D}{2}(1-\frac{i\kappa_{II}}{i\kappa_{I}}. (12)


Now the reflection coefficient defined as

R = \frac{B^{*}B}{A^{*}A}, (13)

which is the probability that an incident electron (wave) will be reflected. On the other hand, the transmission coefficient is the probability that an electron will be transmitted through the potential (e.g. barrier potential).  These two also must satisfy the relation

R+T =1. (14)

This means that the probability that the electron will be reflected or transmitted is 100%. Therefore, to evaluate R (the reason why I don’t calculate T will become apparent shortly), we take the complex conjugate of A and B and using them in Eq.(13) we get

R = \frac{\frac{D}{2}(1+\frac{\kappa_{II}}{\kappa_{I}})\frac{D}{2}(1-\frac{\kappa_{II}}{\kappa_{I}})} {\frac{D}{2}(1-\frac{\kappa_{II}}{\kappa_{I}}) \frac{D}{2}(1+\frac{\kappa_{II}}{\kappa_{I}})} = 1. (15)

What this conceptually means is that the probability that the electron is reflected is 100%. This implies that it is impossible for an electron to be transmitted through the potential for this system.


Basic Equations of Ideal One-Fluid Magnetohydrodynamics (Part II)

Continuing with the derivation of the ideal one-fluid MHD equations, the next equation governs the motion of a parcel of fluid (in this case plasma). This momentum equation stems from the Navier-Stokes’ equation. The derivation of this equation will be reserved for a future post. However, the solution of this equation will not be attempted. (Incidentally, the proof of the existence and uniqueness of solutions to the Navier-Stokes’ equation is one of the Millennium problems described by the Clay Institute of Mathematics.)

The purpose of this post is to derive the momentum equation from the Navier-Stokes’ equation.


The Navier-Stokes’ equation has the form

\frac{\partial \textbf{v}}{\partial t} + (\nabla \cdot \textbf{v})\textbf{v} = \textbf{F} - \frac{1}{\rho}\nabla P+\nu \nabla^{2}\textbf{v}, (1)

where \textbf{F} represents a source of external forces, \textbf{v} is the velocity field, \nabla P is the pressure gradient, \rho is the material density, \nu is kinematic viscosity, and \nabla^{2}\textbf{v} is the laplacian of the velocity field. More specifically, it is a consequence of the viscous stress tensor whose components can cause the parcel of fluid to experience stresses and strains.

Defining the magnetic force per unit volume as

\textbf{F}=\frac{\textbf{J}\times \textbf{B}}{\rho}, (2)

where we recall that \textbf{J} is the current density defined by Ohm’s law in a previous post, and also recall from Basic Equations of Ideal One-Fluid Magnetohydrodynamics (Part I) the equation

\nabla \times \textbf{B}=\mu_{0}\textbf{J}, (3)

if we solve for the current density \textbf{J}, we get

\textbf{J}=\frac{1}{\mu_{0}\rho}[(\nabla \times \textbf{B})\times \textbf{B}], (4)

where \mu_{0} is the permeability of free space. Now we invoke the vector identity

[(\nabla \times \textbf{B})\times \textbf{B}]=(\nabla \cdot \textbf{B})\textbf{B}-\nabla \bigg\{\frac{B^{2}}{2}\bigg\}. (5)

At this point, we assume that we are dealing with laminar flows in which case the term \nu=0. The Navier-Stokes’ equation becomes the Euler equation at this point. (Despite great understanding of classical mechanics, one phenomena for which we cannot account is turbulence and sources of friction, so this assumption is made out of necessity as well as simplicity. For processes in which turbulence cannot be neglected the best we can do in this regard is to parameterize turbulence in numerical models.) Using the vector identity as well as making use of our assumption of laminar flows, ideal one-fluid momentum equation is

\frac{\partial \textbf{v}}{\partial t}+ (\nabla \cdot \textbf{v})\textbf{v}=-\frac{1}{\rho}\nabla \bigg\{P + \frac{B^{2}}{2\mu_{0}}\bigg\}+\frac{(\nabla \cdot \textbf{B})\textbf{B}}{\mu_{0}\rho}, (6)

where the additive term to the pressure \frac{B^{2}}{2\mu_{0}} is the magnetic pressure exerted on magnetic field lines and the additional term \frac{(\nabla \cdot \textbf{B})\textbf{B}}{\mu_{0}\rho} is the magnetic tension acting along the magnetic field lines.



Basic Equations of Ideal One-Fluid Magnetohydrodynamics (Part I)

The field in which the interaction of electrically conducting fluids (i.e. plasmas) and magnetic fields are studied is called magnetohydrodynamics (MHD). As an undergraduate, I investigated the MHD processes occurring in the core of Jupiter. The project had two components: calculation of the magnetic field by coding via MATLAB and investigating the magnetohydrodynamics of the metallic hydrogen in the jovian interior. 

The purpose of this post is to derive the the induction equation of MHD from first principles as well as to describe one of my undergraduate research topics.


We start with Maxwell’s equations of electrodynamics:

\nabla \cdot \textbf{E}=\frac{\rho}{\epsilon_{0}}, (1)

\nabla \cdot \textbf{B}=0, (2)

\nabla \times \textbf{E}=-\frac{\partial \textbf{B}}{\partial t}, (3)

\nabla \times \textbf{B}=\mu_{0}\textbf{J}+\mu_{0}\epsilon_{0}\frac{\partial \textbf{E}}{\partial t}. (4)

In the above equations, \textbf{E} is the electric field, \textbf{B} is the magnetic field, \rho is the charge density, \epsilon_{0} and \mu_{0} are the permittivity and permeability of free space, respectively, and \nabla is the del differential operator.  In the context of ideal MHD, we assume that the displacement current (the second term in Eq.(4)) becomes negligible since the velocities considered are believed to be non-relativistic. Therefore, Eq.(4) becomes

\nabla \times \textbf{B}=\mu_{0}\textbf{J}. (5)

A few ancillary expressions are required at this point, namely: Ohm’s law \textbf{J}=\sigma \textbf{E} and the Lorentz force \textbf{F}=q(\textbf{E}+\textbf{v}\times\textbf{B}). Combining these two gives

\textbf{J}=\sigma(\textbf{E}+\textbf{v}\times\textbf{B}). (6)

Solving for \textbf{J} in Eq.(5) and substitution into Eq.(6) and then solving for the electric field \textbf{E} yields

\textbf{E}=\frac{1}{\mu_{0}\sigma}(\nabla \times \textbf{B})-(\textbf{v}\times\textbf{B}), (7)

Taking the curl of the electric field above, while noting the identity \nabla \times (\nabla \times \textbf{A})=\nabla(\nabla \cdot \textbf{A})-\nabla^{2}\textbf{A} we get the induction equation

\frac{\partial \textbf{B}}{\partial t}= (\textbf{v}\times \textbf{B})+\lambda \nabla^{2}\textbf{B}, (8)

where \lambda is the magnetic diffusivity.

The induction equation is the first of the basic equations of ideal one-fluid MHD.





An Introduction…

A bit of background about myself: Since my sophomore year of high school I have been interested in astronomy, physics, and mathematics. I received my Bachelor’s degree in Earth and Planetary Sciences concentrating in astronomy/astrophysics from Western Connecticut State University. My coursework and independent readings ultimately led to minors in mathematics and physics.

My intention for this blog is to serve as a reference in astrophysics and related topics to myself as well as others. I aim to share my own research interests and consider selected problems that have fascinated me. I also hope to communicate recent news in the fields of physics and astronomy and discuss the implications of discoveries made.

DISCLAIMER: I am by no means an expert, and as such the posts that I create are of my opinion and my own logic. I may be wrong sometimes, and I hope that the people who see this (assuming that anyone sees this) will respect that.

That being said… Enjoy!


(ABOVE: An image of the moon taken with a lunar and planetary imaging camera mounted to a Newtonian 130mm reflecting telescope.)