Institute for Computational Astrophysics

Previous Images of the Month - 2014

January February March April May June July August September October November

December 2014

click on image to enlarge

For almost a century, the star Vega has played an important role as a photometric "standard". However, in the past couple of decades several studies have suggested Vegas is not as "standard" as was once thought.  Compared to stars of a similar type it has been found to have a higher than usual luminosity and radius, as well as showing an uncommon flat-bottomed shapes on most of its weak lines. These properties can partially be explained by assuming Vega is a rapidly rotating star, and interferometric observations plus high resolution spectroscropy support this hypothesis. Yet despite qualitative agreement, these two distinct types of observations disagree on two key parameters, namely the equatorial rotation velocity and the overall inclination. However, new analysis of interferometric results suggest that their results are particularly sensitive to the way gravity darkening is applied in the models used to match the observations. For rotating stars, gravity darkening makes the equatorial regions less bright than the poles, and small changes in how the effect is applied may be able to reconcile the differences between observational methods. Gravity darkening is usually applied using a single parameter that relates the effective temperature and the local effective gravity at the surface of the star. However, this relation is only approximate for stars that rotate rapidly and the use of more realistic models that relate the effective gravity and effective temperature at the surface of the star appear necessary in order to analyze the observed properties of Vega.

With this in mind, PhD candidate Diego Castaneda working together with former ICA director Dr. Robert Deupree and Embry-Riddle Aeronautical University professor Dr. Jason Aufenberg, used a new scaling method (see Castaneda & Deupree (2014)) to generate surface properties of rotating stars based on more realistic 2D rotating models. The models do not require a gravity darkening law to be fixed beforehand, which may help disentangle the set of parameters that describe what we observe from Vega.  The scaling method requires an initial set of observed parameters to generate a guess on the surface properties of the star in question. In the case of Vega, we set: Teff = 9550K, i=5deg, V sin i = 21.5 km/s, Polar radius = 2.20 Solar Radius. After the surface parameters are known, it is possible to calculate synthetic interferometric observations as well as a synthetic SED to compare with the actual observations. The plot shown here presents comparisons between our models and both the interferometric observations (left) and a sample of two weak lines (right).  There is a good agreement for both sets of data  when compared with previous results but a more robust parameter assumption analysis is currently being done in an attempt to help shed some light on Vega's actual properties.

November 2014

click on image to enlarge

GrayStar is a pedagogical application that models the atmospheric structure of a star and a representative spectral line, computes various standard observables including photometric colours and the line equivalent width, and displays the textual and graphical output, all in response to ad hoc user input.  GrayStar is written entirely in JavaScript and HTML and is guaranteed to run using nothing more than a common Web-browser with no special add-ons or plug-ins running on any common personal computing device with any common OS. It implements the gray solution to the model atmosphere problem, and the "core plus wing" approximation to the Voigt spectral line profile problem, and runs in just a few seconds of wall-clock time, allowing for rapid responsiveness during ad hoc class room demonstrations and for real-time parameter exploration during student laboratory projects.

The output includes the five photometric colour indices of the Johnson-Cousins photometric system and the equivalent width of the spectral line in the text output area (top panel), and in the graphical area, the position of the modelled star in the Hertzsprung-Russell (HR) diagram (upper-left), a rendering of the spatially resolved limb-darkened and -reddened star based on the radiative transfer solution (upper-middle), the temperature structure vs optical depth (upper-right) and vs geometric depth (middle-right), the pressure (including gas and radiation pressure) vs optical depth (middle-left), the monochromatic limb darkening at three representative wavelengths (middle), the flux spectral energy distribution (SED) and the intensity SED at two representative angles of emergence (lower-left), the flux spectral line profile and the intensity line profile at two representative angles of emergence (lower-middle), and a hybrid Grotrian/level-population diagram (lower right) showing the number of particles in key atomic energy levels of the line-absorbing species at optical depth one.  There are additional pedagogical aids throughout both the input (not shown) and output display.

The code is available for both execution and download for local installation at .  I encourage people to customize local installations by editing the HTML and JavaScript code, and posting their experience and advice at .

Prof. C. Ian Short, P.Phys

October 2014

click on image to enlarge

Supermassive black holes in the centres of galaxies range in mass from 10^5 all the way up to 10^9 times the mass of the Sun. However, to reach this size they most grow by consuming material from the nucleus of the galaxy they're in. The fastest mass growth comes from consuming cold gas, in so-called giant molecular clouds, which is the same material that fuels star formation. Given that both star formation and supermassive black holes can inject very large amounts of energy into the galaxy, heating the cold gas in molecular clouds, the precise amount of gas consumed by the black hole and that converted into stars can vary enormously during galaxy formation.

This graph, taken from Thacker et al 2014, demonstrates this inherent complexity by plotting up the evolution in the black hole accretion rate (BHAR) compared to the star formation rate in the galaxy (SFR) as a giant elliptical galaxy is formed. The different lines show different literature models for the black hole accretion and the amount of energy returned to the galaxy. Each arrow represents 20 Myr of evolution, and overall we see a rising BHAR, a peaking SFR and then a decay down to lower values at later times (as the elliptical galaxy forms). Clearly, there are significant differences between models, and understanding galaxy mergers, in particular, should help uncover more of the details necessary to disentangle this remarkably complex physics.

This plot was constructed using data from simulations run by Dr James Wurster, a former ICA member, during his PhD thesis. The ICA computing cluster 'Cerberus' was used to conduct the simulations.

Rob Thacker, Professor

September 2014

click on image to enlarge

Galaxies in the universe are not scattered uniformly. Instead, they cluster into lumpy structures and filaments which trace the underlying distribution of dark matter. This animation shows the distribution of dead galaxies (ones that have stopped forming new stars) at progressively earlier and earlier epochs. The size of the circle representing each galaxy corresponds to its estimated stellar mass, ranging from as low as 10^7 solar masses to as high as 10^12 solar mass. These masses have been obtained by using the observed fluxes in different wavebands to find the closest matching stellar population (known as spectral energy distribution (SED) fitting). This can be computationally expensive, and was thus carried out on ACEnet facilities, as it is an optimization process over many different possible stellar populations. SED fitting also allows us to infer the star formation rates and ages of the galaxies.

The galaxies were observed in the near-infrared (NIR) as part of the ultraVISTA survey (McCracken et al. 2012), and combined with earlier observations in the COSMOS survey to estimate precise redshifts (Muzzin et al. 2013).

Anneya Golob, PhD student

August 2014

click on image to enlarge

During the last decade, a large number of surveys have been constructed in our attempt to comprehend the formation and evolution of galaxies.

When trying to understand clustering of galaxies or performing statistical analysis for their different populations, it is important to understand to which degree the processing of the data can suffer observational or instrumental incompleteness. This is when random object simulation becomes very useful.

In this process, simulated objects are constructed to have the same properties as those expected from the real observations and then they are implanted randomly onto the observed images (in order to reflect the survey geometry, noise properties and the overlapping with foreground sources).  In the end, these simulated objects are subject to the same selection procedure as the one used for the real observations.

What we observe in this plot is the object recovery fraction as a function of simulated Ks magnitude. Solid lines represent completeness functions when object recoveries are only based on their coordinates, whereas dashed lines also include photometry. In the black lines, we see the recovery function when simulated objects are implanted onto a simulated background image (simbg). Different colours represent the number of simulated objects that were implanted in the science image and make it possible to study the effect of overlapping in object recovery.

This technique is especially important in order to build statistical analyses of different populations within a field, such as number counts, luminosity functions and to measure angular correlation functions.

Dr. Taro Sato (former ICA postdoctoral fellow), under the supervision of an ICA faculty member: Dr. Marcin Sawicki, performed these simulations thanks to the high-performance computing environment provided by ACEnet and its staff.

For more information please refer to Sato, Sawicki and Arcila-Osejo (Accepted MNRAS).

Liz Arcila Osejo, PhD student, ICA

July 2014

click on image to enlarge

The graph above displays BH activity curves: the differential of the amount of time spent by a BH above a given accretion rate. Activity curves effectively show how much time is spent by a BH at a given luminosity which is an important factor in the quasar luminosity function, the number of quasars of a given luminosity in a given volume. While the activity function has been inferred from observational data, the resulting curve is considerably different to previously theoretically preferred "light bulb" models which spend time essentially at one luminosity and are essentially either on or off. This graph is the result of simulations ran using a parallel version of Hydra (an adaptive particle-particle, particle-mesh/SPH code; Thacker & Couchman 2006). The blue line shows the expected (observed) BH distribution (Hopkins & Hernquist 2009). The green line represents the simulated activity curve using the WT model (Wurster & Thacker 2013), while the red line represents the simulated activity curve using the PNK model (disc viscous time = 5 Myr & accretion radius = 0.05 times the smallest resolved smoothing length; Power, Nayakshin & King 2011). The two models mainly differ in the they deal with particle accretion and energy feedback: 1) The WT model accretes the nearest particle onto the BH at a modified Bondi accretion rate using a continual-conditional algorithm, while the PNK model accretes particles onto an accretion disc that in turn accretes onto the BH at a rate set by the mass of the disc and a viscous timescale. 2) The WT model returns feedback energy thermally using a top-hat kernel, while the PNK model uses a momentum approach where the momentum is distributed radially and isotropically. The simulated activity curves of these models are more dominant at high luminosities; they do not increase monotonically for lower luminosities as suggested by observations although some models (i.e. PNK) produce an activity function with a constant value at lower luminosities rather than tailing-off like the so-called "light bulb" model. This work was performed by MSc candidate Maan H. Hani, in collaboration with ICA faculty member Dr. Rob Thacker.


click on image to enlarge

The Wide-field Infrared Surveyor for High-redshift (WISH) - shown in the top right - is a proposed space mission that will find extremely rare, luminous first galaxies and quasars in the Epoch of Reionization.  The observatory will have a very large focal plane which forces an unorthodox filter exchanger configuration. Since a standard filter wheel would be too large, heavy, and cumbersome, the solution is, instead, to have a set of flip-type filter exchangers (lower left).  This solution saves weight and size, but introduces additional complexity because of the large number of moving elements (48 for the flip-type mechanism assembly instead of just 3 for an equivalent filter wheel system). This raises questions about the reliability of the system.

Simulations carried out by ICA faculty member Prof. Marcin Sawicki show that the expected failure probability for the flip-type system as a whole is indeed worse than that for the simpler filter-wheel assembly, but the difference is not as dramatic as one might intuitively expect: even though the flip-type system has 48 moving elements compared to the wheel system's three, the probability that it remains fully functional until the end of the mission's 5-year lifetime is only slightly worse than that for the wheel system. The lower-right panel above shows the survival probability for the flip-type system (dashed curve) and that of the wheel system (solid); in the model realization shown the probability that the flip-type mechanism will fail before the nominal end of the mission is 0.8% while that for the wheel mechanism is 0.5%.  Given the small difference between the two probabilities, and the fact that the filter wheel mechanism is impractical for WISH in any case, the flip-type mechanism presents a very good solution to the problem of covering a very large focal plane.

May 2014

click on image to enlarge

The Result

Asteroseismic analysis reveals that Procyon has a convective core that is well beyond that predicted by classical stellar models. The plot shows how much beyond the classically predicted size of the convective core the mixed region is extended versus the probability (Bayesian statistics) that the convective core is extended (the x-axis units are pressure scale heights above the classical radius of the connective core). The data and our computer models imply that there is some unknown process deep inside Procyon that is literally stirring things up and significantly extending the size of the inner mixed region.


If true then other stars may also have similarly extended convective cores. The most dramatic consequence of which is that the ages we calculate for young stars will have to be modified to take this into account because the extended convective core will mix fresh hydrogen into the nuclear burning core. Unfortunately, at this time we do not know precisely what is causing the extension of the convective core. 

How we got the result

The simple graph shown here is the end product of the efforts of two Ph.D. students, Michael Gruberbauer (Saint Mary’s University), studying Bayesian statistics, and Christian Straka (Universität Heidelberg), studying nonlocal convection theory, and their respective supervisors, David Guenther (Saint Mary's University) and Pierre Demarque (Yale University). 10 million detailed stellar models of Procyon based on a variety of different interior convection scenarios were computed. The models were compared to the observations and the probability that a given model and its convection scenario were a likely match was computed. The final result distilled from several months of computations on ICA’s computers is shown here.

April 2014

click on image to enlarge

The image shows surface radius variation for the n=1, ℓ=2, m=0 oscillation mode of 1.875 M8 Zero Age Main Sequence models at various rotation speeds. The surface equatorial velocities are given in km/s. At zero rotation the latitudinal variation is that of an ℓ=2 Legendre polynomial. Rotating models incorporate a range of Legendre polynomials instead of just one, gradually changing the nature of the latitudinal variation. As the rotation rate increases, the node shifts towards the pole and eventually a second node is created near the equator. At the highest rotation velocity a third node has appeared near the pole. It should be noted that when the second node appears the mode does not resemble an ℓ=4 Legendre polynomial. While this case shows an increase in nodes with increasing rotation, there are cases which show a decrease in the number of nodes. These mode calculations were performed by Ph. D. student Mr. Diego Castañeda working with Dr. Robert Deupree.

March 2014

click on image to enlarge

            The interaction of turbulent convection and radial stellar pulsation remains one of the more challenging problems in stellar astrophysics. Dr. Chris Geroux, former Ph. D. student in the ICA and currently a postdoctoral fellow at the University of Exeter, and Dr. Robert Deupree are addressing this issue with 2D and 3D hydrodynamic simulations of RR Lyrae variable stars. Performing the calculations in more than one dimension allows the inclusion of at least the large scale convective flow rather than relying on a formulaic procedure to impose a time dependent behavior of convection. They previously showed that 2D simulations matched RR Lyrae light curves quite well (see November 2012 image of the month). However, the convective energy transport in the turbulent outer layers of these stars is clearly 3D in nature, and the 2D calculations have now been extended to 3D.

            The variation in the pulsation amplitude with location in the instability strip is shown for both the 2D and 3D calculations and observational data from Cacciari et al. (2005) for the RR Lyrae stars in the globular cluster M3. The open symbols refer to the data, the red symbols to the 2D simulations, and the blue symbols to the 3D simulations. Circles refer to fundamental mode pulsators and squares to first overtone pulsators. The 2D and 3D amplitudes agree rather well, although some of the 3D amplitudes are somewhat below their 2D counterparts. Work is in progress to examine this difference. Both sets of simulations agree reasonably well with the data, particularly since no effort was made to fine tune the input physics to the individual cluster. A comparison of the light curves for individual stars and specific 2D and 3D models is under way.

February 2014

click on image to enlarge

The 2D surface in this plot, computed by Ian Short, displays a statistical parameter (chi-square) that expresses how well computational models of stellar spectra match (or "fit") the observed spectrum of a red giant star of spectral class K 3-4 III.  The x-axis shows the surface (or "effective") temperature, and the y-axis shows the logarithmic surface gravity, of the model star for which the synthetic spectrum is being compared to the observed spectrum.  Solid lines are lines of constant surface gravity and dotted lines are lines of constant effective temperature.  The lower the value of the surface, the more closely the model of the corresponding temperature and gravity matches the data.

Blue indicates the lines of best matching effective temperature and surface gravity, and the red cross where the blue lines intersect marks the lowest point on the surface, corresponding to a best match model star of effective temperature equal to 4200 K and surface gravity equal to 10 to the power of 1.5 (or about 32) cm/s/s.

The surface represents the result of computing 120 model spectra covering the entire visible band (325 to 750 nm) with 120 different combinations of stellar effective temperature and surface gravity using Version 15 of the PHOENIX stellar atmosphere and spectrum modelling code.  The spectra were convolved ("blurred") to match the low spectral resolution of the observed spectra and sampled every 1.5 nm, so that about 285 wavelength points were used to calculate each chi-square value. The model spectra were computed with the ACEnet high performance computing facility.

January 2014

click on image to enlarge

Red giant stars rank among the brightest stars in the Galaxy. Located in the upper right portion of the Hertzsprung-Russell (H-R) diagram, red giants are evolved stars in late stages of their lives, useful for tracing old stellar populations and Galactic evolution through their properties.  Photometric colours, in particular, are useful for determining the fundamental parameters of stars: Effective (surface) temperature (Teff), surface gravity (log(g)), and composition.

Left image:  Observer’s H-R diagram of 3182 stars in the Yale Bright Star Catalogue 5th edition with recorded values for visual magnitude MV, B-V colour index, and trigonometric parallax, displaying the large variation of the B-V index even among similar spectral classes. Giant stars of spectral class K or M are marked in red. The large orange marker denotes Arcturus, a typical early K giant star.

In collaboration with ICA faculty member Dr. Ian Short, PhD student Mr. Mitchell Young computed the spectra of red giant stars with the stellar atmospheric modelling and spectrum synthesis code PHOENIX, running on the ACEnet cluster Fundy, and then calculated the photometric colours via integrated pass-band fluxes using his own procedure written in the Python programming language.  These calculations are based on the most realistic treatment of the thermodynamic state of the gas and radiation field in the star's atmosphere (non-local thermodynamic equilibrium, NLTE).

Right image:  B-V colour index values for Arcturus-like model stars having differing degrees of horizontal temperature inhomogeneities (ΔT1.5D), simulated by averaging together two 1D spectra.  The solid black line is the B-V colour of the NLTE 1D Arcturus; the dashed line is B-V colour of the LTE 1D Arcturus.

The model input parameters, for which the above color indices were found, were average Teff = 4250 K, log(g) = 2.0 (log cm s^-2), and 1/3 solar "metal" abundance, representative of the standard star Arcturus (alpha Bootes).  Mr. Young used a library of similar synthetic spectra for a range of red giant stars to study how the NLTE treatment affects our ability to distinguish horizontal temperature variations on the surfaces of such stars.  Generally, NLTE effects mimic the brightening of the blue end of the spectrum that can be caused by "hot spots" on the surfaces of stars.

Hoffleit, D. & Jaschek, C. (1991). The Bright Star Catalogue.