Multi-component Dynamics of Particles from Nonequilibrium

Background and Motivation:

Many ordered physical phenomena such as ferromagnetic materials and solid state crystals are formed in nature from a nonequilibrium state of many particles or components, often with various distinct types or properties.  Though studying single particle or two particle systems are essential to understanding physical phenomena at fundamental scales, the complex systems presented by multicomponent, many particle systems present interesting problems, and may come with significant practical applications. For example, ferromagnetism is a widely used and familiar phenomenon involving a large number of particles in an ordered solid. In the first half of the 20th century, physicist Ernst Ising was able to provide a mathematical model for magnetism by employing a simple dynamic two-component system of discrete particles in a lattice. Magnetism arises as an inherently quantum mechanical phenomenon that can be described as a collection of atoms with non-zero electron spin states and their associated magnetic dipole moments. The Ising model uses statistical mechanics and simplifies the system to a graph of particles with two possible states, spin up or spin down, interacting according to nearest neighbor approximations. Though the model is a gross simplification of the physical system, it is actually able to predict phase transitions across a temperature threshold from, say, ferromagnetic states to nonmagnetic states. Furthermore, due to the discrete nature of the model, it can smoothly be implemented in computational simulations without much loss of generality. An implementation of the Ising model in MatLAB, and mapping the associated phase transition, is described below.

Models involving multiple-component dynamics from nonequilibrium states can also be used to describe the self-assembly of molecular structures. Despite being a nonequilibrium process, self-assembly has largely been studied according to purely statistical mechanical models that consider the discrete dynamics of the particles to play no role in the evolution of the system except to convey the concept that systems tend to evolve along paths that minimize the free energy of the system [2]. Though this assumption holds for systems in near equilibrium states or for simple systems under mild nonequilibrium conditions, descriptions that do not account for the dynamics that emerge out of local interactions between particles fail to accurately describe the evolution of systems far from equilibrium. This is particularly true for self-assembly processes that reach thermodynamic equilibrium on long time scales with respect to experimental or computational limitations. For example, a system of particles forced into harsh nonequilibrium conditions, such as a liquid that is rapidly supercooled, may experience binding errors during crystallization, leading to kinetically trapped structures corresponding to “randomly” emergent local minima on the free-energy landscape. It is worth noting that kinetic trapping can also affect self-assembly of multi-component solid state systems under mild nonequilibrium conditions due to the slow mixing of components within a given system.

The Ising Model:

The Ising model was originally developed by German physicists Wilhelm Lenz and Ernst Ising in the 1920s. As stated above, the Ising model simplifies a magnetic solid by assuming that particles are fixed to an evenly spaced lattice of points and all have spin states aligned in either the +z or –z direction ({$ s = +1$} or {$-1$}). Due to interactions between nearby magnetic dipoles, particles have particular energy states associated with them according to a nearest-neighbors approximation. Here, a two-dimensional Ising model is implemented with each particle having four nearest-neighbors, as shown in Figure 1. The interaction energy between two particles with parallel spin components is given by {$\epsilon_s$} (the subscript {$s$} denoting “same”), and the interaction energy between two particles with antiparallel spin components is given by [$\epsilon_d$ ({$d$} denoting “different”). According to the theory of magnetism, magnetic dipoles with parallel alignments are in lower energy states than those with antiparallel alignments. This condition is satisfied by {$\epsilon_s - \epsilon_d > 0$}.

isingmodelfigure

Figure 1. Schematic spin model for a ferromagnet. [1]

The dynamics of the model is achieved by randomly choosing particles at some rate {$c$} to flip to its opposite spin state with a success rate proportional to a Boltzmann factor with respect to the energy of the chosen particle. This probability is given by given by {$min( e^{-\Delta E/k_B T},1) $}, where {$\Delta E = \left( \sum_{<i,j>} \epsilon_{<i,j>} s_{i,j,flipped} s_{<i,j>} \right) - \left(\sum_{<i,j>} \epsilon_{<i,j>} s_{i,j,current} s_{<i,j>} \right) $} and {$<i,j>$} indexes over all four nearest neighbors. It becomes apparent that for flip events with {$\Delta E < 0$} flip with a probability of unity, corresponding to a lowering of energy, and that there is also non-zero probability for a particle to flip to a higher energy state [1]. However, the system should tend toward microstates with lower energy on a statistical whole, eventually approaching some minimum value given a finite number of particles and space.

Results and Analysis:

 

monte_carlo_ising_model_under_critical_temp_no_magnetic_field

Figure 2(a) and 1(b): Monte Carlo implementation of Ising model (2D). Figure 2(a) shows a time series plot of the spin lattice. Lattice was initiated with particles in a checkerboard pattern (ie. with alternating spins). Temperature was below critical value and no magnetic field was applied. Figure 2(b) shows plots of the average energy and magnetization over time iterations. Energy decreases over time, and magnetization tends to approach either +1 or -1, implying total parallel alignment.

As expected, the energy decreases over time, and small localized domains of aligned spins begin to appear. Total alignment (or magnetization {$|M|=1$}), is quickly approached as the free energy “reward” for parallel alignments increases. That is, at low temperatures, the Boltzmann factor is very small for positive values of {\Delta E}, and thus the chance of flipping to a higher energy state is lower.

monte_carlo_ising_model_under_critical_temp_h_applied ising_h_em_1

Figure 3(a) and (b): Monte Carlo implementation of Ising model (2D) with constant magnetic field in positive direction (ignore H=0 label).  3(a) shows a time series plot of spin lattice. Lattice was initiated with particles in a checkerboard pattern (ie. with alternating spins). Temperature was below critical value. Figure 3(b) plots the average energy and magnetization over time iterations. Energy approaches finite negative value and magnetization approaches finite positive value.

As a positive constant magnetic field is applied to the entire lattice, all dipoles begin to align along the direction of the magnetic field. However, once magnetized, if the magnetic field is set to zero, the magnetization will generally fall back down to zero unless the interaction energy for alignment for the material is large enough (for example, in ferromagnets).

evt mvt

Figure 3(a) and (b): Monte Carlo implementation of Ising model (2D). Spin dynamics were calculated for each temperature value over 25,000 Monte Carlo iterations, and the last 1,000 values of energy and magnetization were averaged and plotted against temperature. All calculations were done for a 50×50 lattice, with all other relevant parameters normalized to one.

Figures 3(a) and (b) demonstrate the phase transition from a ferromagnetic state to zero magnetization. One can note some instability in both the magnetization and energy plots around 2.25 (units scaled to {$J/k_B$}), indicating a critical point in the temperature. The system becomes incredibly sensitive to small perturbations as the properties of the material are rapidly changing over temperature. Exact analyical solutions for the critical point how that {$T_c = 2/ln(1+\sqrt{2}) \approx 2.27$} [1]. As T approaches 4, the system crosses the phase transition boundary according to mean field theory.

Future work:

The methods utilized by the Ising model to calculate the dynamics of many-particle systems is not limited to spin systems. The assumptions made by the Ising model, such as that the system dynamics are dominated by nearest neighbor interactions between n-types of particles and obey statistical mechanics, can easily be generalized to other systems as well. For example, consider the slow deposition of molecules with two distinct conformers onto a reactive substrate surface. Supposing that there are significant intermolecular interactions between the molecules and the substrate surface and between the molecules themselves, the dynamics of the system may be simplified to a growing analog of an Ising gas with two-components forced onto an underlying lattice. If this underlying substrate has significant nonhomogenous structure on the scale of the lattice spacings, patterns in the depositing overlayer may emerge. This has actually been observed in certain organic-inorganic interfaces, and the resultant structure may play a key role in the functionality of the material [4].

Conclusion:

A Monte Carlo simulation of the Ising model was developed in MatLAB, along with some extensions to multi-component growth modeling of gaseous particles solidifying on a substrate. Ideally, my implementation of Ising model can be generalized and used to develop other multi-component simulations with more complexity. This program had a number of limitations in its practical usage and theoretical robustness. For one, the order of calculates necessary to compute various results increases rapidly with complexity. My implementation in particular was limited to simulating a 2D lattices containing somewhere between a couple hundred and a thousand points, where as magnetic solids are 3D objects containing on the order of {10^23} atoms. Furthermore, though many unnecessary boundary effects were avoided by implementing toroidal boundary conditions on the lattice edges, which is useful for treating magnetism in the bulk, there is a lack of middle ground where both the bulk and the magnet edges can be accurately simulated in this model. There were also a number of significant approximations and assumptions made during the formalism of the model that affects the accuracy and relevance of the model to real physical systems. However, the model is surprisingly qualitatively powerful despite its simplicity, and has a wide range of applications to other multi-component systems.

References:

[1] Giordano, Nicholas J. Computational Physics. 2nd ed. Upper Saddle River, NJ: Prentice Hall, 2006. Print.

[2] S. Whitelam, L. O. Hedges, and J. D. Schmit, Self-Assembly at a Nonequilibrium Critical Point, Phys. Rev. Let., 112, 155504 (2014).

[3] S. Whitelam, Growth on a rough surface, 2016.

[4] C. Ruggieri, S. Rangan, R. Bartynski, E. Galoppini, Zinc(II) Tetraphenylporphyrin Adsorption on Au(111): An Interplay Between Molecular Self-Assembly and Surface Stress, J. Phys. Chem. C, 2015, 119 (11), pp 6101–6110

Gravity Assist

Overview:

The computational problem that I worked on for my final project was to simulate a satellite undergoing a gravity assist as it passes a planet. The “gravity assist” is a technique used to help spacecrafts reach targets that are “higher up” in the sun’s potential well. This program creates 2-D animations that are not centered around either the satellite or the planet to express the gravitational interactions between the two. What is meant by neither body being centered is that the both the planet and satellite are moving. The planet is used to change the satellite’s trajectory. The program itself generates a velocity graph of the gravity assist of a satellite as it passes by a planet and the velocity graph of a satellite that gets caught and temporarily orbits the planet. The program’s purpose is to show how the initial conditions of a spacecraft as it passes a planet determines its trajectory.

Methods:

The following image shows how the differential equations used in the program were derived. The equations are from a model for the earth’s orbit around a stationary sun. Assuming that the earth has a circular orbit around a stationary sun it allows for the use of constants instead of inputting information for the masses of the two bodies and the gravity between them. What was changed from the earth-sun example is that the larger body is in motion. Also instead of the “r” value being the radius from the sun the “r” value represents the distance between our planet and satellite.

sling-shot-equations

The method specifically used was the Euler-Cromer method (using i+1) because it conserves energy when used with oscillatory problems. So the positions of both the planet and satellite are calculated at each time step and then the velocities are calculated accounting for how the distance between the two change. Now the velocities are based on interactions between the two bodies. In this program it is assumed the size difference between the two is great enough that the satellite doesn’t change the planet’s path.

Figures:

assist-the-og

Figure 1: The gravity assist of a satellite as it passes nearby a planet.

catch-and-relaease

Figure 2: The gravitational capture of a satellite by a planet for two revolutions.

assistog

Figure 3: Graphs of the the change in the Velocity over 5 years

Discussion:

The planet’s initial conditions were an initial position at the origin and an x-velocity of 5 AU/yr. In the gravity assist simulation the satellite began 86 AU away at (-50,-75) with an initial x-velocity of 0.15 AU/yr and y-velocity of 2 AU/yr. The program then simulated the position and velocity for 5 years and after interacting with the planet the the satellite was 4,195 AU (the “r” value)  from the planet and had a final x-velocity of 0.717 AU/yr and a final y-velocity of 1.61 AU/yr. So the total velocity changed from ~2 AU/yr to ~1.75 AU/yr. The left of Figure 3 shows the change in the x-velocity over 5 years, the spike in the graph represent the satellite’s interaction with the planet and then as the two get further away the velocity settles to a constant value.

The second simulation represents a “catch and release” scenario where the satellite is caught by the planet’s gravity but as the two move the satellite is able escape the planet’s gravity. In the simulation it appears that the two touch but the figures are similar to point charges in that that the mass used for the calculations is at the very center. What seems to be the outside of the shapes is just an image representation used in the simulations.

In this scenario the planet still begins at the origin but it has an x-velocity of 1 AU/yr. The satellite now begins ~74 AU away (-55,-50), has an initial x-velocity of -0.45 AU/yr  and y-velocity of 0.45 AU/yr. After 5 years and two complete revolutions the satellite is ~500 AU away (the “ODr” value) and has an x-velocity of 0.036 AU/yr and a final y -velocity of 0.157 AU/yr. The total velocity decreased from ~0.64 AU/yr to ~0.16 AU/yr, the change in the x-velocity over the 5 years is shown on the right of Figure 3.

Conclusion:

The simulations accurately modeled the gravity interactions between a satellite and a planet. If the velocities of the satellite and planet are great enough as the bodies interact (~2 AU/yr for the satellite and 5 AU/yr for the planet) then the two will pass each other with only a change in trajectory of the smaller body. But if the two are moving at low velocities that are close in value (~0.64 AU/yr for the satellite and 1 AU/yr for the planet)  then the smaller body will get caught in the gravity of larger.

The “gravity assist” technique allows for fuel conservations as spacecrafts travel to higher orbits because the trajectory can be changed through interactions with planets. The drawback of the technique is the total velocity decreases after interaction with the planet, in the assist case there was a decrease in velocity by 12.5%. The limitation of the program is that it is unable to simulate changing the velocity at any point in the path which could be done with a spacecraft. Further investigations could look at finding the ideal relationships between travel time, speed, and fuel consumption when trying to reach a specific trajectory.

Resources:

Giordano, Nicholas J., and Hisao Nakanishi. Computational Physics. Upper Saddle River, NJ: Pearson/Prentice Hall, 2006. Print.

“A Gravity Assist Primer.” NASA. NASA, n.d. Web. 11 Dec. 2016.

Propagation Of A Wave Packet In An Infinite Well With Potential Barriers

Auther: Samuel Gilbert

Introduction: 

This project models the propagation of a quantum wave packet in an infinite square well, and introduces between zero and three potential energy barriers. The wave packet used is a “Gaussian” wave packet (Eq. 1), and is a solution to the Schrödinger equation.

ψ(x,t=0) = C*exp[-(x-x0)2/σ2],                                                                                                              (Eq. 1)

C is a constant that represents the magnitude of the wave packet, x is the position variable, x0 is the initial position of the wave packet, and σ is the initial width of the wave packet.

The program is designed to be highly interactive, and therefore requests three inputs from the user – the number of potential barriers, the strength of the potential barriers, and the initial energy of the wave packet. An initial setup with three potential barriers is shown below in Fig. 1:

wave-packet-and-three-barriers

Fig. 1: Image of a sample initial setup. 

The motivation for this program was to create an interactive education tool. To this end, the program provides instant feedback to changing variables such as energy of the wave packet and the strength of the barrier, which is a particularly useful tool in understanding quantum reflection, transmission, tunneling, and furthermore, the probabilistic nature of the quantum world.

Three plots are produced simultaneously – the imaginary part of the wave function, the real part of the wave function, and the probability density. In plotting the probability density, the program gives a visual representation of the fact that the position of a quantum object is not well defined, as it is governed by probability. The user will see the probability density split over a potential barrier, representing the probabilities of the wave packet being reflected or transmitted, respectively.

Methodology: 

The computational method used in this program is the “leap-frog” method, as described in Giordano (1). The leap-frog method was used, because it allows for time dependence, unlike the shooting method studied in class. Our wave packet is then allowed to develop over time from the initial wave packet in Eq. 1., which it does. The leap-frog method alternates between evaluating the real and imaginary parts of the wave function. This, in essence, evaluates the wave function in ∆t/2 increments. In addition, it is notable that the time dependence of the Schrödinger equation is contained in the computational method itself, which is a beautiful way of handling time dependence. The leap-frog method can be thought of as an adapted Euler-Cromer method, as the updated real part of the wave function is then used to calculate the updated imaginary part, and so on throughout the routine. Giordano (1) cites that the leap-frog method has a numerical error of the same order as the 2nd order Runge-Kutta method. A significant portion of the code was dedicated to interactivity, and therefore deals with the presentation of the animated plots. This includes the stop button on the plots, which calls an external function which closes all windows, and also making sure the animation runs for the right amount of time.

Results: 

Fig. 2 below shows three animated gifs for three different energies with one potential barrier.

1-barrier-energy-11-barrier-energy-2 1-barrier-energy-3

Fig. 2: One potential barrier modeled for three different energies. 

These three animations represent what a user would see if the inputs were 1 potential barrier of strength 2*10-6, and a wave packet with initial energy of 1.5*10-33. This is just one of the many options that a user can choose. For example, the user can explore using the same initial energy, but changing the strength of the barrier. In addition, portions of the wave packet can get caught between two potential barriers, and add constructively, as shown below in Fig. 3:

trapped

Fig. 3: Wave packet caught between barriers, and constructive interference. 

The wave functions can get complicated very fast by adding enough potential barriers and letting the program run long enough. These simulations therefore allow for complex systems than could be modeled by hand.

Analysis: 

One of the unexpected results of this program is that the Gaussian wave packet decays over time, as shown below in Fig. 4.

decay

Fig. 4: Spreading of the Gaussian wave packet

Sawada et. al. (2) explain that a Gaussian wave packet is the sum of a few Gaussian curves. These Gaussian curves all have different momenta and velocities. Therefore over time, the Gaussian wave packet develops to show these different velocities, which is the cause of the spreading. Although a pitfall of the model used here, it works well, as this wave packet can be thought of as a superposition of quantum states decaying to the ground state. This is not what is really happening mathematically, but it does still work for the program. It can also be thought of as a formulation of the Heisenberg uncertainty principle, in that we are initially containing the wave packet to a defined space, and therefore there must be uncertainty in its momentum, which is the compositional Gaussian curves (1). Thus, it makes sense that the wave packet should spread out. Nevertheless, the spreading wave packet means that the simulation cannot run for too long without becoming meaningless – the most useful simulation comes in the beginning.

Further examination of the wave packets produced by different energies shows that the wave packets that have less energy have fewer peaks, and wave packets with large energies have many peaks. The nature of the compositional Gaussian curves is shown here, as the high energy wave packets are composed of more Gaussian curves, with each of those curves having its own energy. It is also noted that the high energy wave packets travel faster than the low energy wave packets, as expected.

Conclusion:

This program can be developed much further to include different types of potential barriers. The one used here is essentially a delta function, however it was finite in the program, as MATLAB cannot plot delta functions. However nothing prevents the program from being expanded to include step potentials and the like. There range of possible values for the strength of the potential barrier and the energy of the wave packet could also be expanded.

A natural extension to understanding the spreading of the Gaussian wave packet could be to perform Fourier transformations on the wave packet to parse out the composition of the wave packet. The spreading of the wave packet has been analyzed by Pan (3) using the “Bohmian Model.” The crucial assumption that the Bohmian model makes is that the wave packet has a definite location. This assumption, of course, contradicts Heisenberg’s uncertainty principle. However, the Bohmian model provides new inside into understanding the spreading of the Gaussian wave packet. Thus in the future, this program could potentially not only model a quantum object, but also describe mathematically and visually how the wave packet spreads out using the Bohmian model.

Bibliography: 

(1) Giordano, Nicholas J. Computational Physics. Upper Saddle River, NJ: Prentice Hall, 2006. Print.

(2) Sawada, Shin-Ichi, Robert Heather, Bret Jackson, and Horia Metiu. “A Strategy for Time Dependent Quantum Mechanical Calculations Using a Gaussian Wave Packet Representation of the Wave Function.” The Journal of Chemical Physics 83.6 (1985): 3009. Web.

(3) Pan, Alok Kumar. “Understanding the Spreading of a Gaussian Wave Packet Using the Bohmian Machinery.” Pramana 74.6 (2010): 867-874. Web.

 

4 Body Orbital System

My model of Jupiter, Mars, and the Earth orbiting the Sun can be used as a tool for educating people about planetary motion. The planets have maintained their orbits for billions of years, but this model seeks to prove their stability and test under what conditions their orbits might become unstable. By gradually changing the mass of Jupiter, we can see that it is not until the planet’s mass is 1000 times more than its current weight that it would begin to affect the orbits of Mars and Earth. At this extremely increased mass, The Earth and Mars both exhibit unstable behavior. This model is not only interesting to those who study astronomy, but can be informative and reassuring of the Universe’s stability for those with no knowledge of the subject.

Physics and Astronomy often cause students a lot of anxiety. These subjects are not usually intuitive and can intimidate students, resulting in very few people studying the subject past high school.[1] This unfortunate occurrence must be minimized and graphical teaching tools can be an effective way to engage students. Presenting a visual representation of the often unobservable processes studied by physicians and astronomers is invaluable to a struggling student. Astronomy is inherently a visual subject and it makes sense to model computationally phenomena that may take too long to study physically. With a computational model you can map the progress of a planet for years in a matter of seconds. My model maps the orbits of Jupiter, Mars, and Earth for 20 years and includes a legend and text that makes it easily interpreted. My hope is that this model will help to illuminate the effects of planetary gravitational pull to beginning astronomy students.

fin

Graph 1: This graph shows the orbits of Earth, Mars, and Jupiter after 20 years when Jupiter is at its normal mass.

fin-2

Graph 2: The orbits after the mass of Jupiter has been increased by a factor of 1000.

These graphs are made by building equations from Kepler’s laws that make each planet sensitive to the gravity of all other planets in the system. This means that each planet will have four equations to solve for its velocity and position. The following shows the velocity equation for the x direction of the Earth.

vxE(i+1)=vxE(i)-(4*(pi^2)*xE(i).*dt)/(rE(i).^3)-4*pi^2*(J/S)*(xE(i)- xJ(i))*dt/rEJ(i)^3-4*pi^2*(M/S)*(xE(i)- xM(i))*dt/rEM(i)^3;

Each velocity equation must account for the mass and gravitational pull of the surrounding planets and can be done by manipulating Kepler’s laws. These equations should then be iterated using the Euler-Cromer Method.

The Euler Cromer Method is usually the best method to model oscillatory motion. In my first draft of this code, I modeled the lone Earth around the Sun. Because of an error in my equation, the Earth’s orbit was not elliptical. Upon doing further research, I attempted a 2nd Order Runge-Kutta model which proved too complicated to use when finishing the entire system. [2] I eventually was able to successfully use the Euler-Cromer Method which now perfectly models the orbital motion.

Conclusion:

This project was able to successfully build a teaching tool that would be really easily used in a classroom setting. The animation makes the model interesting to students and the year counter lets you know exactly when Earth and Mars would shoot off into space. While the unstable behavior only occurs under extreme circumstances, it is still an excellent lesson on the stability of the solar system. In the future, I would love to add more planets to the system, but I fear that the length of my equations and the number of matrices I would need to create would slow the computing time considerably. When I added Mars to the system, I had to create another 7 matrices and my equations increased sizably in length. In addition, the perihelion of mercury is more complicated to calculate and might pose a problem. Despite these potential issues, a system in which all 8 planets are represented would extremely enhance the model and its educational value.

Resources:

  1. Grim, Nancy. A Force Concept Correlation Study with Instructional Methods, Anxiety, Perceptions of Difficulty and Student Background Variables. 1999. 11 Dec. 2016 <http://eric.ed.gov/?id=ED438164>.
  2. Giordano, Nicholas J., and Hisao Nakanishi. Computational physics. Pearson Education India, 2006.

Solving Kepler’s Equation

Overview

For my final project I solved Kepler’s Equation numerically using Newton’s Method. Because it is a transcendental equation and thus cannot be used to solve for the eccentric anomaly directly, I created a program that uses iterative methods to solve it and then plots the solutions. Beyond this, I focused in particular on making the program user friendly to make it more usable as a teaching tool. For instance, I included detailed explanations of Kepler’s Equation as well as the ability for the user to change parameters. The reason I chose this project was because I wanted to deepen my own understanding of Kepler’s Equation and creating a program that teaches people about the equation helped me do that.

Program

Kepler’s Equation is as follows: M = E – e*sin(E) where M is mean anomaly, E is eccentric anomaly and e is eccentricity. The program focuses on explaining these three quantities and the roles they have in calculating orbits. The following picture (from www.jgiesen.de) depicts mean and eccentric anomaly:

jgiesen

In the picture, the black ellipse is the actual orbit being calculated (slightly eccentric) and the blue circle is the orbit without any eccentricity. The mean anomaly is the angle from the main body (in the solar system: the sun) to where a planet would be if its orbit were perfectly circular.The eccentric anomaly is the angle from the main body to a point above the actual orbit on a tangent circle.

fig1

The first graph that the program produces is a plot of eccentric anomaly versus mean anomaly at all values of eccentricity (0 to 1) in steps of 0.1. The shaded region represents possible values for E at fixed M or vice versa. As eccentricity is increased from 0 to 1, M and E get further apart. Next, we look at the extremes: e = 0 and e = 1:

fig3 fig2

As you can see, when e = 1, M and E are very different at most points. However, at e = 0, M  and E are the same and the plot is a straight line. This makes sense looking back at the definition of M and E; an e value of 1 means that the orbit is a perfect circle and therefore the mean and eccentric anomaly are equal.

screen-shot-2016-12-11-at-2-00-43-pm

The program also allows for some interactivity so it can be used as a teaching tool; above is how the program explains the graphs of mean anomaly versus eccentric anomaly. After the explanation, it waits for user input and then graphs mean anomaly versus eccentric anomaly with the user-specified eccentricity.

Conclusion

Overall, I’m happy with the outcome of this program. I think that Newton’s Method worked well to solve Kepler’s Equation and it was fairly straightforward to implement in Matlab. I also liked how the educational part of the program came out. However, if I were to redo this program I would possibly use a different programming language because I was unable to find any way to make the interactive part of the program user friendly. For example, I would’ve liked the text explanations to display somewhere besides the command window and would’ve added a button to advance in the program rather than through typed numbers.

References

http://www.jgiesen.de/kepler/img/keplerEllipse.gif

https://www.csun.edu/~hcmth017/master/node16.html

http://mathworld.wolfram.com/KeplersEquation.html

Fourier Analysis of Unemployment Rate in U.S.

Overview:

Throughout the course in computational physics this semester, I learned how to apply Matlab to different concepts and lessons in physics. I wanted to apply what I learned and collaborate them with ideas and materials from my major in Economics. For this project, I analyzed the data of the unemployment rate in U.S. using Fourier analysis. The Fourier transform decomposes a function of time into the frequency components that make up the original function. As a result, one can look at the different cycles and waves that make up the function. I look at the percent changes in the unemployment rate monthly from February 1948 to October 2016. The Fourier transform data has the amplitude graphed against its frequency, so I can look at the different cycles of the unemployment rate. The second part of the project was to apply a filter to the Fourier transformation. By applying the low-pass filter, I filter out the high frequency data. Then, I inverse Fourier transform the low-pass filtered data to compare it to the original data. I expected the data to have less noise, because the high frequency data that affect the volatility of the data in low periods will be filtered.

Fourier Transformation:

The first part of the project was to Fourier transform the unemployment data. I needed reliable data recorded frequently, so I used the unemployment rate recorded monthly from the U.S. Bureau of Labor Statistics. Then I had to apply the Fourier transformation.

fo

Eq. 1

According to Fourier’s theorem, every function can be written in terms of a sum of sines and cosines to arbitrary precision; given function f(t), the Fourier transform F(f(t)) is obtained by taking the integral of f(t)’s sines and cosines function. The exponent factor is just the sum [cos(2*pi*frequency*t)-i*sin(2*pi*frequency*t)]. Discrete Fourier transform can be expressed into the next set of equations from the previous ones.

fouriereq

Eq. 2

The top equation in Eq. 2 is the discrete Fourier transform. The bottom equation is used to get the inverse of the Fourier transform or the original function from the Fourier transform equation. Since the equations in Eq. 2 are just a sum of exponential terms, it appears to be a very straightforward numerical evaluation. In reality, discrete Fourier transforms take a very long time to compute. Each term involves computing the exponential factor, which then have to be multiplied by xn and added to the running total. The total number of operations is of order N^2, because each sum for a given frequency component has N terms and there are N frequency components. However, one can see that the exponential terms in the discrete Fourier transform are multiples of one another. This makes it possible to group up terms in the sums in a way you can store those values and reuse them in evaluating the different components of the discrete Fourier transformation. This is how the fast Fourier transformation came about.

With discrete data set of N ~ 1000, I can perform the fast Fourier transform or the FFT to efficiently compute the discrete Fourier transform. Then, I squared the absolute values of the transformed data. Making this change to the data does not affect the data under Parseval’s Theorem that states, “the squared magnitude of the Fourier Series coefficients indicates power at corresponding frequencies.” The reason behind the squaring manipulation was to change the complex numbers to real components for graphing. Then, I graphed the amplitudes that I get from fft against the frequency. The highest frequency was the data being recorded every month, so I had to set the maximum frequency at .5. Bi-monthly data is the lowest period I can get accurate data from due to Nyquist frequency, which states that minimum rate at which a signal can be sampled without introducing errors is twice the highest frequency present in the signal. Any amplitude that belongs to higher frequencies than .5 would have errors. Many Fourier analysis of signals in physics is shown in varying range of Hertz. However, the periods in economic data have relatively very high periods compared to physics data periods measured in degrees of seconds. As a result, the frequency is relatively low in this graph. Also, the frequency is measured in number of unemployment cycles per month.

fouriertransform

Figure 1.1

Left side of figure 1.1 shows the graph of the Fourier transformation. The right side just has amplitudes of less than 30,000 amplitude filtered out, so that it would be easier to look at the frequencies with higher amplitudes. One of the challenges of the project was getting the most accurate graph that I can. There are so many different ways the economic data can come out, however, the most accurate graph would have to show evidence that align with known theories in economics.

Analysis:

The frequencies with the highest peaks are the most significant components of the original data. The Fourier transformation graph of figure 1.1 has large peaks near the frequencies of .01699 cycles/month and .3653 cycles/month. These frequencies correspond to 58.858 months/cycle and 2.737 months/cycle respectively.

Business cycles, which include unemployment cycles, are supported very strongly to have cycle length of 5 years on average, which is 60 months. The high peaks that happen above and below the frequency of .01699 are the cycle components of fluctuations that happen in the business cycle. The Cycles do have an average period length of 58.858 months, but cycles can also happen in shorter or longer lengths. The other significant peak range is at 2.737 months/cycle. This aligns with the seasonal unemployment rate; the unemployment rate has a cycle for every fiscal quarter, which is about 3 months. Through these evidences, I can say that the unemployment rate data has been accurately Fourier transformed.

http___www-nber-org_cycles_cyclesmain-1

Figure 1.2

The known theories support the numbers I have, but I also wanted to compare my numbers to the actual numbers. Figure 1.2 is a record of U.S. business cycle expansions and contractions from the National Bureau of Economic Research. They also provide the average cycle (duration in months) of business cycles for both trough from previous trough and peak from previous peak. I looked at peak from previous peak, since the difference wasn’t significant. The average duration of business cycles from economic data of 1859-2009 was 56.4 months, which is very close to the value I got. However, the data also had recorded the average business cycles from 1945-2009. This applies more directly to my transformed data, because my data was from 1948-2016. The average for this period was 68.5 months. I wanted to look at why this data’s average was a bit off from my most significant peak.

The average does fit in the range of the cycle fluctuations, but I figured out that the most plausible reason is that because my transformed data was not smoothed out. One can see in Figure 1.1 Fourier Transformation graph that the business cycle frequency peaks at x=.01699 cycles/month (58.858 months/cycle) and drops to the next point where x=.01335 cycles/month (74.906 months/cycle). The average of the two months/cycle points is 66.882 months/cycle, which is much closer to the data presented in Figure 1.2. If my transformed data were to be smoother with more data points, then the average would have been more accurate.

Low-Pass Filter:

fourierfilteredfinal

Figure 1.3

For the second part of the project, I wanted to apply a simple filter on the transformed data and see how the original data would be affected by that. The idea behind it was to filter out high frequencies above a certain point so that only the low frequency components would be left over for the data. If I were to build up the unemployment function again using the filtered components, what would the original data look like?

I expected the unemployment data to have less noise after the filter, because only the longer period cycles would build up the unemployment rates. I used the inverse Fourier transform function ifft to return the Fourier transformed data back to unemployment rate vs. time graph. The graph on the very right in figure 1.3 is the original unemployment rate graph taken from U.S. Bureau of Labor Statistics. I simply used ifft on the transformed data to return it.

The first step to filtering was to filter the power component of the transformed data, shown in the left, 1st graph of figure 1.3. This wasn’t an accurate representation, because the y axis has been affected by Parseval’s theorem. In order to fix this problem, I filtered the data only applied with fft. The most accurate filter graph I got was the middle graph in figure 1.3.

Analysis:

I made components from 300 to 412 (higher frequency components) equal to 0, so that we would still keep the information about those components and not lose resolution. Then, I inverse Fourier transformed that data. Consequently due to the change of the total length with data, the inverse produces transformed data with different unemployment rates from the original graph, including complex numbers. I was able to only plot the real number components. However, the filter helped me see the general trend of the unemployment rate changes at the exchange of some loss of information; the real number components that I do care about in the middle graph of figure 1.3 represents a more accurate graph of unemployment data with less high frequency noise. The evidence is shown in the decrease in fluctuations of the unemployment rate changes when the middle and right graph are compared in figure 1.3. The rates’ general trend converge more along 0. We do see a time (phase) shifting, but this does not impact the Fourier series magnitude and the inverse.

Conclusions:

The numbers I have gotten through the Fourier Transformation looks accurate and they align with the actual data published. I view my attempt at learning the material and applying it to have been successful. Since, I have used the data from 1948-2016, I should focus on the business cycle published between 1945 and 2009 in figure 1.2. The average of the business cycle is longer post 1945. Assuming that the trend will continue, I can look at the lower frequency end of the most significant peak fluctuation in the Fourier transform graph of figure 1.1. The reasonably lowest amplitude value of the fluctuation range has an x value of .009709 cycles/month, which is about 103 months or 8.6 years.

Our economy has been expanding since the end of the last business cycle, 2009, after the 2008 recession. The decrease in unemployment rate since 2009 has been polynomially decreasing, which indicates we are entering the contraction stage of the business cycle. Since we 7.5 years out of the 8.6 year business cycle I predict has passed, using my data, I can predict that the recession and increase in unemployment rate will happen in the next 1 or 2 years.

Additionally, it would have been interesting to develop filtering methods and programs further so I can get more smooth and accurate data on low frequency business cycles. However, given the time, I believe that the simple low-pass filter was sufficient to show how a filter affects the Fourier transform data and its inverse.

Reference:

NBER. “US Business Cycle Expansions and Contractions.” The National Bureau of Economic Research. NBER, 20 Sept. 2010. Web. 09 Dec. 2016.

Series, Nber Working Paper. Matthew D. Shapiro (n.d.): n. pag. NBER. NBER, 1988. Web. 9 Dec. 2016.

US. Bureau of Labor Statistics, Unemployment Level [UNEMPLOY], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/UNEMPLOY, November 29, 2016.

MIT. “The Fourier Series and Fourier Transform.” Circuit Analysis (n.d.): 102-23. MIT. MIT, Spring 2007. Web. 9 Dec. 2016.

James, J. F. A Student’s Guide to Fourier Transforms: With Applications in Physics and Engineering. Cambridge: Cambridge UP, 2011. Print.

Numerical Solution to Schrodinger Equation in a 1D Periodic Potential

Summary:

The following article serves as a discussion of a program that I designed to serve as an educational tool for students of introductory solid state physics.  This program calculates the dispersion relation for a particle in periodic 1D potential and presents related information, such as band gap energies and movies of the time-evolved stationary states directly above and below the first band gap.  I begin by introducing the reader to the underlying physical concepts, such as energy gaps, dispersion curves, Bloch’s theorem, and the nearly-free electron model.  I then proceed with a discussion of the program and how it can be used to demonstrate a wide array of phenomena related to to these theoretical topics.  I conclude with a description of how the program operates, its limitations, and suggestions for how this work can be extended in the future.

Introduction:

In many solids, such as crystals and semiconductors, atoms are arranged in an organized and periodic manner. This periodicity in the organization of the atoms leads to an electric potential that is also periodic in space. Describing the behavior of particles, such as electrons, in such potentials is a primary topic in any introductory solid state physics course, and remains an important area of research in modern electronics.

An interesting phenomenon that occurs in periodic potentials, and a topic which is at the heart semiconductor device operation, is the presence of energy bands. In many solids, there are intervals of allowed energies which particles experiencing that potential can take. That is, if one were to measure the energy of an electron in a periodic potential, they might measure values anywhere between two energies, say, a and b, but nowhere between two other energies, say c and d. The allowed energy intervals are known as energy bands and the energy difference between two bands is known as a band gap [1].

In introductory solid state texts, such as “Solid State Physics,” by Kittel, the origin of the band gap is often presented in the context of the nearly-free electron model [1]. In the nearly-free electron model, electrons are described as free particles with an effective mass that depends on the electron’s wavenumber. The energy gap is then explained conceptually as being due to the existence of standing waves that form when the wavelength of the electron is some harmonic of the periodicity of the potential. Two distinct standing waves form, one in which the electron’s wavefunction is concentrated around regions of high potential energy, and one in which the electon’s wavefunction is concentrated around regions of low potential energy. The energy difference between these two standing wave states is taken as the bang gap energy.

While presenting band gaps in this way avoids grappling with the full complexity of quantum mechanics, band gaps are an essentially quantum mechanical phenomenon. By the time many students begin studying solid-state physics, they have already studied at least one of the archetypal quantum mechanics examples, such as the infinite square well or the hydrogen atom. Furthermore, they have learned that only certain energies have associated normalizable solutions to the Schrodinger equation, and that these are the measurable energies of the system. The only difference in infinite periodic potentials, is that instead of discrete energy levels, there are intervals of allowed energies. The goal of this project is to help students bridge the gap between the conceptual nearly-free electron model and the results of quantum mechanics.

Bloch’s Theorem:

One of the most important theorems involving solutions to the Schrodinger equation in a periodic potential is Bloch’s theorem. This states that the normalizable solutions to the time independent Schrodinger equation in a periodic potential have the form \psi_k=u_k(\vec {r})e^{i \vec{k} \cdot \vec{r}}, where \vec{r} is the position vector, \vec{k} is the wave vector, u_k(\vec {r})=u_k(\vec {r}+\vec{T}), and \vec{T} is a lattice vector (the position vector between two lattice sites) [1]. In one dimension, this reduces to the criterion that \psi_k=u_k(x)e^{ikx}, where u_k(x)=u_k(x+na), a is the periodicity of the potential, and n is an integer. It follows from Bloch’s theorem that if \psi_k is a solution to the Schrodinger equation, then \psi_{k+G} is also a solution, provided that G=\frac{m \pi}{a} for some integer m [1]. It is common to plot the allowed energies against k for \frac{- \pi}{a} < k < \frac{\pi}{a} (known as the first Brillouin zone), as this gives a complete account of the allowed energies of the system [1]. The energy of the allowed wavefunctions as a function of k is known as the dispersion relation or dispersion curve [1].

Project and Results:

I created a program that calculates the dispersion relation for three 1-dimensional periodic potentials, a square, triangle, and sinusoidal wave (figure 1).

Fig 1: Periodic potentials which are available which are available to be solved by simulation. This is understood to represent a single cell of the potential functions. The true potentials infinitely many cells sewn together at the edges.

Fig. 1: Periodic potentials available available in the program. This is understood to represent a single cell of the potential functions. The true potentials repeat infinitely many times in both the positive and negative spacial directions.

This program is intended as an educational tool for students of solid-state physics who are beginning to study the concept of bang gaps. It calculates the dispersion curves by finding the first several allowed energy eigenstates (of the potential of interest) for a discrete selection of k values in the first Brillouin zone. It then plots the energies of these states against the wave number, k, which gives the dispersion curve (figure 2). The first several band gaps are evident as the empty regions between dispersion curves.

tr

tr1

Fig. 2: Dispersion relations calculated by the program. The top image corresponds to the triangle-wave potential, the middle image corresponds to the sinusoidal potential, and the bottom image corresponds to the square-wave potential. Regions between the colored horizontal lines represent band gaps.

The program also produces a movie of the time-evolved eigenstates on either side of the first band gap (figures 3 and 4). In these movies, it is worth noting that the state just above the first energy band is concentrated around high-energy regions of the potential and the state just below the first band is grouped around low-energy regions of the potential. One can also see that the periodicity of the allowed wavefunctions on either side of the gap is indeed a harmonic of the periodicity of the potential (actually a sub-harmonic). These results demonstrate that the qualitative explanation of band gaps provided by the nearly-free electron model is in close agreement with the results of quantum mechanics.

www-gifcreator-me_jyqlj

Fig. 3: Animation of the real and imaginary parts of the wavefunction just below the first band gap (click image to see animation)

Fig. 4: Animation of the real and imaginary parts of the wavefunction just above the first band gap (click image to see animation)

Along with having several types of potentials to choose from, the height and periodicity of the potential are left as variable inputs for students to change and explore. By changing these values, students can begin to build and intuition for how the structure of the dispersion curves depend on the size and shape of the potential. For instance, by decreasing the height of a potential of fixed periodicity, students are able to watch the structure transition from being an insulator (large band gap), to a semiconductor (small band gap), to a metal (little or no band gap) [1].

Fig. 4 Dispersion relation for decreasing magnitude of the potential barriers. The decrease in energy gaps between dispersion curves represent the transition from insulator to conductor.

Fig. 5: Dispersion relations for decreasing magnitude of the potential barriers. The decrease in energy gaps between dispersion curves demonstrate the transition from insulator to conductor.

The program can also be used as a tool for students to build a better understanding of Bloch’s theorem. For instance, by instructing the program to calculate the allowed energies for values of k beyond the first Brillouin zone, we see that the dispersion relation is periodic in k (figure 6). This corresponds to the result from Bloch’s theorem that if \psi_k is a solution to the Schrodinger equation in a periodic potential, then \psi_{k+G} is also a solution [1].

Fig. 5 Dispersion relation, which is periodic in k

Fig. 6: Representative dispersion relation that extends beyond the first Brillouin zone.  Note the periodicity in k, as predicted by Bloch’s theorem.

Verification of Results –  The Kronnig-Penney Model:

In order to demonstrate that this program produces solutions that are in agreement with theory, I considered the square-wave potential. This potential corresponds to the Kronnig-Penney model, which has been studied extensively. There is an analytic solution for the allowed energies of this potential in the form of roots to a transcendental equation [1]. Much like the methods employed by this program, this expression is found using Bloch’s theorem to determine appropriate wavefunction boundary conditions.

This equation for the allowed energies has been used to find the analytic dispersion curve for the Kronnig-Penney model. This is plotted along with the dispersion curve determined by the numerical methods employed in this program (figure 7). Indeed, the two dispersion relations are in excellent agreement, which demonstrates that this program does reliably calculate the allowed energies for a 1D periodic potential.

Figure 1: The top image shows the dispersion curve for the square-wave potential calculated by numerical methods. The bottom image shows the dispersion curve calculated using analytic results obtained from Bloch's theorem.

Fig. 7: The top image shows the dispersion curve for a square-wave potential calculated by numerical methods. The bottom image shows the dispersion curve calculated using analytic results obtained from Bloch’s theorem.

A more in depth discussion of the Kronnig-Penney model can be found in “Solid State Physics,” by Kittel [1].

Program Operation:

In order to determine the allowed energies of a particle in a periodic potential, I employed the numerical methods discussed in “Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions,” by Dr. Constantino Diaz [2]. The theory and operation of the procedure outlined in this paper is structured in the following way. The time independent Schrodinger equation is a second order differential equation in x, so there exists two linearly independent solutions for each value of E, which we will denote as \psi_C(k,E,x) and \psi_S(k,E,x) [2,3]. In other words, any wavefunction with wavenumber k and energy E, is expressible as a linear combination of \psi_C(k,E,x) and \psi_S(k,E,x). That is, \psi_k(E,x)=A\psi_C(k,E,x)+B\psi_S(k,E,x) for some complex numbers A and B [2,3].

\psi_C(k,E,x) and  \psi_S(k,E,x) are determined for a chosen value of k by using the Schrodinger equation to calculate the wavefunction in discrete steps [2]. This is done by expressing the Schrodinger equation in finite difference form:

    \[\frac{\hbar^2}{2m} \frac{d^2}{dx^2} \psi+V\psi=E\psi \approx \frac{\hbar^2}{2m} \frac{\psi_{n+1}+\psi_{n-1}-2\psi_{n}}{\Delta x^2}+V_n \psi_n=E\si_n,\]

where V(x) is the potential function, m is the particle mass, \Delta x is the spacial step size, and n denotes the n^{th} spacial spacial step [3]. Rearranging terms gives:

    \[ \psi_{n+1}=\frac{2m(\Delta x^2)}{\hbar^2}(E-V_n)\psi_n-\psi_{n-1}+2\psi_n,\]

which can be used to calculate the wavefunction for some guess value of E [3].  Since the potential functions are symmetric, initial conditions of \psi_C(\frac{a}{2})=1, \frac{d\psi_C}{dx}(\frac{a}{2})=0 and \psi_S(\frac{a}{2})=0, \frac{d\psi_S}{dx}(\frac{a}{2})=1 can be chosen without loss of generality [2,3].

According to Bloch’s theorem, \psi(a)= u_k(a)e^{ika}=u_k(0)e^{ika}=\psi(0)e^{ika}, and \psi^{'}(a)= u_k^{'}(a)e^{ika}+iku_k(a)e^{ika}=u_k^{'}(0)e^{ika}+iku_k(0)e^{ika}=\psi^{'}(0)e^{ika} [1,2]. These boundary conditions can be expressed in terms of the coefficients A and B as A\Delta C+B\Delta S=0 and A\Delta C^{'}+B\Delta S^{'}=0, where \Delta C=\psi_C(a)e^{\frac{ika}{2}} -\psi_C(0)e^{\frac{-ika}{2}, \Delta S=\psi_S(a)e^{\frac{ika}{2}} -\psi_S(0)e^{\frac{-ika}{2}, \Delta C^{'}=\psi_C^{'}(a)e^{\frac{ika}{2}} -\psi_C^{'}(0)e^{\frac{-ika}{2}, and \Delta S^{'}=\psi_S^{'}(a)e^{\frac{ika}{2}} -\psi_S^{'}(0)e^{\frac{-ika}{2} [2]. For this system of equations to be satisfied, the determinant of the coefficients must be equal to zero. That is, the allowed energies for a given value of k satisfy the equation, G(E)=\Delta C \Delta S^{'} - \Delta S \Delta C^{'} =0. If G(E)=0, then E is taken as an allowed energy [2]. If not, then a new value E+\Delta E is treated as the energy and new values of \Delta C, \Delta S, \Delta C^{'}, and \Delta S^{'} are computed. This process is repeated many times so that several allowed energies are found for the chosen value of k. The whole procedure is then repeated for evenly spaced values of k in the first Brillouin zone.

Program Limitations:

The primary limitation of this program is its sensitivity to choices of step size.  If \Delta E is chosen to be too small, an insufficient number of roots are found for G(E) to create a viable dispersion curve.  If the energy step size is chosen too large, then the program is unable to locate the zeros of G(E), particularly at the ends of the first Brillouin zone.  In both of these instances, the program stops running and produces plot of G(E) for the last calculated value of k. If this graph has less than four roots, then the step size was chosen to be too small.  If the graph has many roots, but the program still failed, then the step size was chosen to be too big. The program seems to be most effective if \Delta E is chosen to be roughly one percent of the maximum of the potential functions.

The program also fails to produce reliable results if the spacial step size is chosen to be too large.  As long as the spacial step size is no more than five percent of the width of the potential barrier and 0.5 percent of the periodicity of the potential, the program seems to operate normally.  The program also produces poor results if the step size in k is chosen to be too large, but this only decreases the resolution of the dispersion curve.  An overly large choice in the step size of k will not result in a failure of the program unless it is chosen to be comparable to the size of the first Brillouin zone. As long as operators are aware of these problems, the program performs well.  Still, I recommend that only small changes in the parameters of the potential be made each time the program is run so that the user is able to adjust the step sizes gradually.

Future Work:

A good continuation of this project would be to alter the program to accommodate for potentials that are not symmetric about the point x=\frac{a}{2}. Suggestions on how to implement this modification can be found in “Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions,” by Dr. Constantino Diaz [2]. An even more interesting progression of this work would be to write a program that can solve the Schrodinger equation in a two or three-dimensional potential lattice. The particular method employed in this program is not generalizable to higher dimensions, but other methods that also apply Bloch’s theorem to provide appropriate boundary conditions can be applied to periodic potentials of higher dimension [2]. In the case of higher dimensional potentials, the dispersion relation is anisotropic with respect to spatial orientation. Qualitatively, this is due to the fact that the spacing between potential barriers is different along distinct spatial directions [1]. A program that could allow students to explore this directional dependence of the dispersion relation and plot the associated wavefunctions could be a helpful tool to solidify the effect of spatial anisotropy on wavefunctions in crystals.

Conclusions:

I created a program that calculates the dispersion relation for a particle in a 1D periodic potential, which can be used as a educational tool for students of solid state physics.  It uses the discrete version of the Schrodinger equation along with boundary conditions provided by Bloch’s theorem in order to find the first several allowed energies for values of k in the first Brillouin zone.  Comparison of the results of this program to the analytic solution for the dispersion relation of Kronnig-Penney demonstrates that this program provides valid results. The program demonstrates several key properties of wavefunctions in a periodic potential, including energy bands, the distribution of wavefunctions on either side of the band gap around high and low regions of the potential function, the dependency of energy gaps on the height of potential function, and the periodicity of the dispersion curve predicted by Bloch’s theorem. A good continuation of this work would be to incorporate non-symmetric potential functions or to construct a procedure for calculating the dispersion relation in higher dimensional potentials.

References:

[1] Kittel, Charles. Solid State Physics. 8th ed. New Caledonia: John Wiley & Sons, 2005. Print.

[2] Uteras Diaz, Constantino A. “Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions.” Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions. ARXIV, 23 Mar. 2005. Web. 10 Dec. 2016. <http://cds.cern.ch/record/828907/files/0503171.pdf>.

[3] Giordano, Nicholas J. Computational Physics. 2nd ed. Upper Saddle River, NJ: Prentice Hall, 2006. Print.

 

Orbital Motion of Our Solar System

Overview:

For my final project, I made two animations that model the orbital motion of planets in our solar system based on data sets that I generated. The first animation shows the orbit of the inner planets and the second shows the full solar system. While the spacing of the planets is to scale, the size of the planets is not. The sun and eight planets are simulated in two-dimensional plots with the sun at the origin (sorry Pluto!). The plots are in the reference frame of the sun, so the sun is at rest. The planets are color-coded based on the dominant color of their surface (i.e. yellow for the sun, red for Mars, etc). The full orbital paths are shown in the plots along with the moving planets for reference. The plots show at least one full revolution for each planet. This program was created for educational purposes. My intentions are that this program could be shown to children learning about the solar system to give them an idea of the scale of the solar system relative to Earth and a visual representation of the difference in orbital periods.

Background, Methods, and Assumptions:

The data sets were generated using the Euler-Cromer method. I made several key simplifying assumptions that led to the model I used. I assumed that the sun is exactly centered in each planet’s orbit. I assumed that the effects of gravity were only between an individual planet and the sun, with planets having no effect on each other. By also assuming the orbits are circular, we equate the centripetal force to the force of gravity to get a zero net force. Using units of AU, we get that the gravitational constant G, multiplied by the mass of the sun, Ms, is equal to 4π2. Solving the equation then gives us the following equations to calculate the x and y components of the position and velocity for a planet of mass M:

screen-shot-2016-12-09-at-10-21-53-pm

As both the centripetal force equation and gravity force equation have the planet’s mass involved, the final equations surprisingly do not end up involving a factor of the mass. The use of the Euler-Cromer method comes in to play for the position equations, where the velocity calculated for time i+1 is used instead of the velocity from time i (the previous time, as would be used in the Euler method).  I initialized the positions of the planets to be at their respective average distance from the sun along the x-axis. The initial velocity at this point of the orbit is completely in the y-direction. It was calculated by dividing the circumference (still assuming a circular orbit) by the period in terms of an Earth year (for Earth, the initial velocity was 2π AU/Earth year). Based on these initial conditions, I looped through the equations over a time window of 300 with a time step of .005 to generate my data. This process was repeated for each planet using essentially the same code but with different initial positions and velocities. Once I had my data calculated, all that was left was to plot. The orbits calculated are not consistent as the planets make successive revolutions about the sun, resulting in an overlapping ellipse pattern when the full data set was plotted. I plotted a section of the data for each planet that included one orbit to show the path. The for loops plotted one position pair at a time successively, resulting in the animated planets orbiting the stationary sun.

Visualizations:

innerplanetsInner Planets Simulation

solarsystem Full Solar System Simulation

 

Results

As I was trying to model an existing system, I had information to compare my results with. The main way I judged the effectiveness of the program is the relative orbital periods of the planets. The real orbital periods compared to Earth are shown in the following table:

Planet Approximate Earth orbits per one orbit
Mercury .24
Venus .62
Earth 1
Mars 1.9
Jupiter 11.9
Saturn 29.4
Uranus 83.7
Neptune 163.7

For one Earth period, Mercury should undergo about four and an eighth periods and Venus should have one and two-thirds periods. Using somewhat qualitative observations by tracking the planets in the animations, I found that my program has very accurate relative orbits as compared to the exact data. My orbits are, however, seem more elliptical than they should be.

 

Discussion:

A number of the simplifying assumptions made are blatantly not correct. The assumption made that the orbits are circular is not true, especially for Mercury, which has a higher eccentricity. I found it quite interesting that, despite creating the model around circular orbits, the calculated orbits ended up being elliptical. Assuming that the sun is centered in the orbits is wrong, as Kepler’s laws state the sun should be at one focus of the orbit. As the orbits are not actually circular, this means that the sun would not be centered (although most planets have a low eccentricity so this would be only a small shift). The time step used must be quite small, or else the orbits will break down (this happens very quickly for Mercury as the time step is increased). I made the decision to have two animations in order for the inner planet orbits to be clearly observable. In the full solar system plot, it is clear that the periods of the smaller planets run much shorter than the outer planets. Apart from that fact, little else can be seen regarding those orbits due to the scale of Neptune’s distance from the sun. The length of Neptune’s period compared to the rest of the planets created other issues as well. The animations themselves take a long time to run, which is somewhat necessary to have a complete cycle for Neptune.  Another limitation is the amount of computation used in the animations. I ended up using 32 vectors with 60000 entries each to hold all the data, along with several other variables to hold enough information to track a full revolution for Neptune. The visualizations lack legends when the code is run. When legends are included in the plots, the time between frames in the animation was reset from what I coded it to be and the animation slowed down substantially. For this reason the code for legends was commented out in the Matlab program (it can be uncommented if legends are desired at the user’s risk of very long runtimes for the animations). For the visualizations, I was able to override the time in between images in the .gif, and so legends are included in the figures shown.

 

Conclusion:

I was able to successfully and accurately model a somewhat simplified version of our solar system. The Euler-Cromer method was appropriate for this task and I would continue to use it in further investigations. In terms of moving forward with this program, my next steps might be making the program interactive. This process would include allowing the user to select the number of planets, the mass, radii, color, and other parameters. It would also be nice to incorporate moons into the simulation. In terms of improving the current program, I would try to figure out a way to include the legends in the program without disrupting the speed of the animations in the code (this would especially be necessary if the program was interactive). Creating an algorithm to optimize the number of data points needed to have one full period for each planet would also be useful to streamline the data.

 

References:

Giordano, Nicholas J., and Hisao Nakanishi. Computational Physics. 2nd ed. Upper Saddle River, NJ: Prentice Hall, 2006. Print.

“Kepler’s Laws.” HyperPhysics. N.p., n.d. Web. 30 Nov. 2016.

Williams, David R. “Planetary Fact Sheet – U.S. Units.” NASA. NASA, n.d. Web. 20 Nov. 2016.

Modeling a Binary Star System

Overview:

The stated goal of this project was to create a set of simulations to be used in an educational environment to better understand the orbital dynamics of a binary star system by plotting, in 2-D, the x and y positions of the celestial objects within the system. This project was split up into four subsections, or subprograms within the program. The first subprogram considered a simple binary star system, the second a system with stars that had different velocities, the third a system with stars that had low velocities, and the fourth a simple system with a planet orbiting one of the stars. We deviated slightly from our project plan because we were unable to explore the effects of mass on the orbits of the stars; this was because of an issue that arose within our model, prohibiting us from being able to accurately assess the effects of mass. Here, we present our results and conclusion of our project. The method section has been omitted, but can be found within our program for each subprogram.

Results:

All of our subprograms self-generated data, based on some initial parameters, and plotted that data on graphs for N iterations. We have attached our graphs here in the form of clickable GIFs. The number of iterations N is responsible for the accuracy of our orbits. In our programs, we use low values for N to obtain semi-accurate orbital data without expending too much computational power or too much time. This is why our orbits may not seem to completely line up; if we increase N to around a thousand or fifteen hundred, then we will see that the orbits will line up quite well. It is important to note that the dots are representative of the velocities of the celestial objects. If there is a cluster of consecutive dots, that indicates that the celestial objects are moving with constant velocities; however, if consecutive dots are far apart, that indicates that the celestial objects are moving with increasing velocity. Additionally, the units used here are not standard units; our position units are in AU, our velocity units are of the form 2πAU/yr, and our mass units are in kg. Below we discuss the results of each subprogram.

I. Simple Binary Star System:

sbss

Figure 1.1: Results of Program 1. Click on image to see GIF simulation.

For our first subprogram, as mentioned previously, we tried to generate a simple binary star system. Our initial parameters were as follows:

For star 1, x = .72, y = .64, V_x = 2.42π, V_y = 0, and M = 5*10^30.

For star 2, x = -.72, y = -.64, V_x = -2.42π, V_y = 0, and M = 5*10^30.

Note that our masses for our stars here are equal (5*10^30 kg). We compared our results, shown in Figure 1.1, to some simulations online, which also had the stars being of equal mass: 1 and 2. Our simulation seemed to be in line with these results. We can also understand why our result makes sense using our knowledge of physics. First, we must consider the following equation:

bc8ad1167d8a4def851ec3df6394bcccThis is the equation for gravitational force between two celestial objects with masses of M_S and M_E. Since our binary star orbit is supposed to be stable, F=F_G=ma. If the objects are close together, then the gravitational force between them is maximized, but because this has to be equal to ma, acceleration must increase since mass cannot be changed. This increase in acceleration is the reason for the increase in velocity. Similarly, if the objects are far apart, then the gravitational force between them is minimized and we expect acceleration to decrease, so the velocities of the objects will approach a constant value.

II. What if the Stars have Different Velocities?

wsdv

Figure 1.2: Results of Program 2. Click on image to see GIF simulation.

Once we were able to generate the simple binary star system, we wanted to see the effects on the orbits due to a change in velocity. For this program, our initial parameters were:

For star 1, x = .72, y = .64, V_x = 2.6π, V_y = 0, and M = 5*10^30.

For star 2, x = -.72, y = -.64, V_x = -2.42π , V_y = 0, and M = 5*10^30.

We increased the V_x of star 1 to 2.6π and kept star 2’s V_x at -2.42π, the same as in subprogram 1. The results can be seen in Figure 1.2. It is, however, important to note a flaw on our method. Namely, we assigned mass to the common barycenter of the stars. This was initially useful to stabilize the orbits of the stars, but as seen here, it makes it so that the orbit star 2 is not affected by that of star 1. We did, however, expect the orbit of star 1 to be more elliptical with a higher velocity because an increase in velocity would mean that the star would be getting closer to its escape velocity, at which point it would move in a straight line.

III. What if the Stars have Low Velocities?

Figure 1.3: Results of Program 3. Click on image to see GIF simulation.

Figure 1.3: Results of Program 3. Click on image to see GIF simulation.

Next, since we were unable to consider changes in mass, we considered the case where the velocities of both stars were very low. For this program, our initial parameters were:

For star 1, x = .72, y = .64, V_x = 1.5π, V_y = 0, and M = 5*10^30.

For star 2, x = -.72, y = -.64, V_x = -1.5π, V_y = 0, and M = 5*10^30.

Note the decrease in velocity from 2.42π and -2.42π to 1.5π and -1.5π for star 1 and star 2, respectively. The result we obtained is shown in Figure 1.3. This result is rather interesting because it seems that the orbits of the stars are moving in a counterclockwise direction through space. However, this orbit is not stable because, after some time, the stars would both escape the system (we know this because we ran this for many iterations before, but did not display it here). We attempted to find a simulation online for comparison, but we were unsuccessful in doing so. However, given that this is built on the program of the simple binary star system, which we know is fairly accurate, we can assume that this may also be a fairly accurate representation.

IV. Simple Binary Star System with a Planet:

Figure 1.4: Results of Program 4. Click on image to see GIF simulation.

Figure 1.4: Results of Program 4. Click on image to see GIF simulation.

This subprogram is the one that we built up to using our previous problems. With a good grasp of the flaws of our system and the orbital dynamics of a simple binary star system, we introduced a planet. Below are our initial parameters:

For star 1, x = .72, y = .64, V_x = 2.42π, V_y = 0, and M = 5*10^30.

For star 1’s planet, x = .3, y = .3, V_x = .09π, V_y = 0, and M = 2*10^24.

For star 2, x = -.72, y = -.64, V_x = -2.42π, V_y = 0, and M = 5*10^30.

The result of our program can be seen in Figure 1.4. It is important to note the initial behavior of the planet. Initially, the planet is assisted by the gravitational force of its host star (since their initial velocities are in the same direction) and moves very quickly outside of the system, but because of the gravitational force of its host star, it is then pulled into an orbit around it. The behavior of the planet is interesting because, as its host star moves, it attempts to stay in orbit by moving along with it. This was to be expected. However, after some time, the planet exits the system. This indicates that this orbit is highly unstable. Although, much time was spent attempting to find a stable orbit with no success. This could explain why we were also unable to find a simulation for comparison – it seems extremely rare or difficult for a planet to be in a stable orbit in a binary star system. Although it must be noted that we had only taken into account the gravitational force on the planet by its host star; previously, we experienced crashes of MATLAB when we attempted to consider the gravitational force of both stars on the planet.

Conclusion:

Referring to our project plans, we were able to achieve most of what we aimed to do. We also believe that this set of simulations can be useful in an educational environment for students. However, because of our model, the only thing we were not able to accurately assess was the effect of mass on a binary star system. Although there are already some interactive simulations online, such as this one, that do so. We experienced a lot of difficulty in building our initial model because we were attempting to replicate an example outlined in Computational Physics. However, we managed to create a good base program to explore further simulations. If we can continue to build on this model, we may consider working out the kinks with taking into account the gravitational force of star 1 and 2 on star 1’s planet. Additionally, it would be interesting to consider multiple planets in different orbits around this binary star system.

References:

“Binary Star.” Wikipedia. Wikimedia Foundation, n.d. Web. 28 Nov. 2016. <https://en.wikipedia.org/wiki/Binary_star#Center_of_mass_animations>.

Fix, John D. “Interactives.” Interactives. N.p., n.d. Web. 28 Nov. 2016. <http://highered.mheducation.com/sites/0072482621/student_view0/interactives.html#>.

Giordano, Nicholas J., and Hisao Nakanishi. Computational Physics. 2nd ed. Upper Saddle River, NJ: Prentice Hall, 2006. Print.

Zingale, Michael. “Introductory Astronomy Animations (and Figures).” N.p., n.d. Web. <https://zingale.github.io/astro_animations/>.