At Stampede's dedication, Austin TX, April 2013. Each of these towers has 3,000 processors.

My current research interests revolve around the processes taking place in circumstellar disks, with emphasis on their importance for planet formation. As these disks are gaseous, hydrodynamics is a component of paramount importance. Since the central equation of hydrodynamics (the Navier-Stokes equation) is non-linear, we have to resort to numerical simulations in order to progress in our understanding. You can find here some movies of my computer models, along with brief explanations of the research highlights, and links to the respective publications.

Most simulations were done with the Pencil Code, a an open source, collaborative, high order finite difference MHD code that is highly modular and versatile. I have been a co-developer of Pencil since 2005.

The videos are also available in my YouTube channel.



Planet formation

The current paradigm in planet formation theory is developed around a hierarquical growth of solid bodies, from interstellar dust grains to rocky planetary cores. A particularly difficult phase in the process is the growth from meter-size boulders to planetary embryos of the size of our Moon or Mars. Subject to angular momentum loss through nebular drag, peebles and boulders are expected to drift extremely rapid in a protoplanetary disk, so that they would generally fall into the central star well before larger bodies can form.

Ways around this problem have focused on the development of pressure traps in the disk. The way these pressure traps work are shown in the figures below. While a negative pressure gradient slows down the gas, a positive one speeds it up. Particles thus feel a headwind in a negative pressure gradient, and a tailwind in a positive one. The result is that particles drift toward pressure maxima.  


Turbulence is a way to generate several of such pressure maxima. The simulation below shows the saturated state of turbulence brought about by magnetorotational instability (MRI) in a 3D stratified wedge of a protoplanetary disk. The quantity shown in the ratio of thermal to magnetic pressure (known as "plasma beta").  



The movie below shows the behaviour of boulders in MRI turbulence. Structure forms as several transient high pressure regions form as a result of the turbulence, trapping particles in them. The color code refers to the bulk density of solids. The bright (yellow-white) is saturated, meaning that these areas are much denser than they appear.  


Publication: Lyra et al. 2008a, A&A, 479, 883. [arXiv] [ADS]

But turbulence is not present throughout the whole extent of the disk. Because sufficient ionization is a main requirement, a significant portion of the cold, dense disk, will be neutral and therefore unmagnetized. In this region, the disk will be "dead" to the MRI and we can expect laminar motion. At the boundary between active and dead zone, however, the transition in turbulent viscosity will trigger Rossy wave instability (RWI), the form that the Kelvin-Helmholtz instability takes in differentially rotating disks. The saturated state of the RWI are vortices, that are very effective in concentrating particles. The figure below shows the particle concentration that results from vortex trapping in this scenario. It leads, in this model, to the formation of over three hundred bound clumps of particles -- 20 of them more massive than Mars.  


Publications: Lyra et al. 2008b, A&A, 491, L41.[A&A highlight] [arXiv] [ADS]
Publications: Lyra et al. 2009a, A&A, 497, 869. [arXiv] [ADS]  

In the video below, a Jupiter-mass planet is perturbing a planetesimal disk of meter-sized bodies. Particles are trapped in the Lagrangian points of the planet and at the border of the gas gap, where large vortices form. In 50 orbits the particles achieve critical density and collapse. The mass attained by the most massive planet is around 5 Earth masses.  


Publication: Lyra et al 2009b, A&A, 493, 1125.
[A&A cover] [arXiv] [ADS]


Disk vortices

Vortices are hydrodynamical features similar to the Great Red Spot. Long hypothesized to occur in protoplanetary disks, vortices have two properties that make them attractive locations for planet formation. The first is that they are equilibrium solutions to the equations of hydrodynamics, which makes them persistent structures. The second is that under the influence of the drag force, solid bodies inside a vortex experience a radial net force, making them sink to its center. This has the convenient side effect of dramatically enhancing the solids-to-gas ratio locally.

The image below is a high resolution model of 2D vortices in disks, with particles. The upper left panel shows the vorticity in the gas phase. The lower left panels the particles in a model without drag force backreaction; the right panel includes backreaction (which becomes important as St approaches 1). Reproduced from the thesis work of Natalie Raettig, with my co-supervision (Raettig, 2013; also Raettig et al. 2013, 2015).


However, although an inverse cascade makes it simple to form and sustain vortices in two dimensions from small scale noise, in three dimensions the direct cascade dissipates their energy fairly rapidly. Therefore the vorticity lost must be replenished by some other mechanism. Two such mechanisms are known, and described below.  


Rossby wave instability

As accretion in protoplanetary disks is enabled by turbulent viscosity, the border between active and inactive (dead) zones constitutes a location where there is an abrupt change in the accretion flow. The gas accumulation that ensues constitutes a localized pressure bump, that behaves as a resonant cavity. Waves are trapped in this cavity if the density transition is sufficiently sharp. For a smooth jump the waves will tunnel away, but for sharp, high amplitude jumps, bound states exist. In resonance, the wave amplitudes will grow exponentially. The saturated state are anticyclonic vortices.

The movie below show the transition between MRI active (outer disk) and MRI dead (inner disk) zones. The transition is as smooth as 15 scale heights in resistivity (from 10-30 AU), following the smooth decrease in ionization due to attenuation of stellar X-rays. At some point the Lundquist number drops rises above 1 and the MRI is abruptly triggered (irrespectively of how smooth the resistivity jump is). Notice a Rossby vortex forming in the transition (the yellowish non-axisymmetric structure at the transition zone by the end of the simulation).

This simulation is the first to show that a smooth jump in resistivity can still trigger the RWI. Previous alpha disks models (Lyra et al. 2009a) suggested that only jumps sharper than 2H would do the trick.  


Publication: Lyra et al. (2014), A&A, accepted. [arXiv]


The model below was the first model of the RWI with a properly modeled active zone. By that we mean with MRI in a MHD disk instead of viscosity jumps in an alpha disk. It represents the transition from the inner active zone to the dead zone, at 0.1 AU, when temperature drops below 900K and collisions are not energetic enough to ionize the alkali metals.  


Publication: Lyra & Mac Low (2012), ApJ, 756, 62.[arXiv] [ADS]


Convective overstability

A hydrodynamical linear overstability exists in protoplanetary disks, powered by buoyancy in the presence of thermal relaxation. We model the system numerically, reproducing the linear growth rate for all cases studied. The figure shows the linear and nonlinear evolution of the overstability. Upper panels: Convergence study with resolution (left) and initial amplitude of perturbation (right). Resolution of 64 points per scale height is enough to resolve the overstability, and initial amplitudes as small as sound speed (c0) are enough to lead to growth, demonstrating the linear nature of the process. The linear growth rate (black dashed curve) is well reproduced in all cases. Lower panels: With the linear overstability raising the amplitude of the initial fluctuations to nonlinear levels, a large-scale vortex is generated in the course of the next 100 orbits. This was the first simulation that generated a self-sustained 3D vortex from linear amplitude perturbation of a quiescent base state.

Publication: Lyra (2014), ApJ, 789, 77. [arXiv] [ADS]  

Before the linearity of the instability had been recognized (Klahr & Hubbard 2014), however, much about the saturated state had been deduced, yet under the assumption that the instability was nonlinear in nature. We used to call it "baroclinic instability", for reasons that have nothing to do with atmospheric physics. It is a misnomer that was around for long enough to make it to a number of publications. To discriminate between that and the geophysical one, we call the former "subcritical baroclinic instability" (SBI) nowadays, following suggestion by Geoffroy Lesur. So, the saturated state of the COS is the SBI. The movie shows the development of the SBI (BI) in a 2D local shearing box model. The box used was small, and only one vortex remained. In general, vortices do not grow past the sonic scale, as they shock, radiating excess vorticity and thus shrinking. The SBI is shown to provide a hydrodynamical source of turbulence, allowing for some degree of accretion (alpha of the order of 5e-3) in the absense of magnetization.

Publications: Lyra & Klahr (2011), A&A, 527, 138. [arXiv] [ADS]
Publications: Raettig et al. (2013), ApJ, 765, 115. [arXiv] [ADS]


Magneto-elliptic instability

Phenomenologically, turbulence can be described as a series of bifurcations, starting with a primary instability that converts shear into vorticity, creating vortices. This is followed by another bifurcation, a secondary instability, to break these vortices into lesser vortical structures. These in turn shall experience a sequence of inertial instabilities, leading to a cascade. Though the Kelvin-Helmholtz instability and the Rayleigh-Taylor instability are well established as examples of primary instabilities, the highly successful theory of the turbulent cascade put forth by Kolmogorov rested on a heuristic picture of secondary instability, established by early experiments. It was not until the 80s that the elliptic instability was introduced as a mechanism for secondary instability. A fluid in rigid rotation supports a spectrum of stable inertial waves, the simplest case being circularly polarized transverse plane waves oscillating at twice the frequency of the base flow. Strain is introduced when the streamlines pass from circular to elliptical, and some modes find resonance with the strain field, leading to de-stabilization. This is called "elliptic instability".

In this simulation we examined the magnetic properties of vortices produced by the SBI. When magnetic fields are introduced in the problem, the addition of Alfven waves enriches the families of unstable modes. When background rotation is introduced, the same as had occurred to the elliptic instability en- sues. When this background rotation runs opposite to the rotation of the vortex (anti-cyclonic motion), a fluid parcel becomes subject to an intense effective shear. Since the magnetic tension resists shear, leading to instability a powerful unstable in-plane mode appears. This instability is of course the magneto-rotational instability, in generalized form. This provides an interesting unification, explaining the magneto-elliptic and magneto-rotational instabilities as different manifestations of the same magneto-elliptic-rotational instability

We find that vortices survive the hydrodynamical elliptic instability, but they fall prey to the more violent magneto-elliptic instability. The conclusion is that vortices are restricted to the dead zones of accretion disks.

Publications: Lyra & Klahr (2011), A&A, 527, 138. [arXiv] [ADS]
Publications: Mizerski & Lyra (2012), JFM, 698, 358. [arXiv] [ADS]


Migration in evolutionary models

In this series of works, I am interested in the final orbit of planets, as a result of migration and disk evolution. As the gas accretes and photo-evaporates, the density (LHS figure) and temperature (RHS picture) fall steadily.

Planet migration occurs towards equilibrium radii with zero torque. These radii themselves migrate inwards as the disk evolves. We see that as the surface density and temperature fall, the planet orbital migration and disk depletion timescales eventually become comparable, with the precise timing depending on the mass of the planet. When this occurs, the planet decouples from the equilibrium radius. At this time, however, the gas surface density is already too low to drive substantial further migration. A higher mass planet, of 10 Earth masses, can open a gap during the late evolution of the disk, and stops migrating. Low mass planets, with 1 or 0.1 Earth masses, released beyond 1 AU in our models, avoid migrating into the star.

Publication: Lyra et al. (2010), ApJ, 715, L68. [arXiv] [ADS]

In the video below we take this model one step further by implementing multiple interacting planets. We incorporate the torques from an evolving non-isothermal disk into an N-body simulation to examine the behavior of planets or planetary embryos interacting in the convergence zone. We find that mutual interactions do not eject objects from the convergence zone. Small numbers of objects in a laminar disk settle into near resonant orbits that remain stable over the 10 Myr periods that we examine. However, either or both increasing the number of planets or including a correlated, stochastic force to represent turbulence drives orbit crossings and mergers in the convergence zone. These processes can build gas giant cores with masses of order ten Earth masses from sub-Earth mass embryos in 2-3 Myr. This simulation and video were done by Brandon Horn, the graduate student at Columbia University, with my co-supervision.

Publication: Horn et al. (2012), ApJ, 750, 34. [arXiv] [ADS]


Transitional disks

The Atacama Large Millimeter Array (ALMA) has been returning images of transitional disks in which large asymmetries are seen in the distribution of mm-sized dust in the outer disk. The explanation in vogue borrows from the vortex literature by suggesting that these asymmetries are the result of dust trapping in giant vortices, excited via Rossby wave instability (RWI) at planetary gap edges. Due to the drag force, dust trapped in vortices will accumulate in the center, and diffusion is needed to maintain a steady state over the lifetime of the disk. While previous work derived semi-analytical models of the process, in this paper we provide analytical steady-state solutions. Exact solutions exist for certain vortex models. The solution is determined by the vortex rotation profile, the gas scale height, the vortex aspect ratio, and the ratio of dust diffusion to gas-dust friction. In principle, all these quantities can be derived from observations, which would give validation of the model, also giving constrains on the strength of the turbulence inside the vortex core. Based on our solution, we derive quantities such as the gas-dust contrast, the trapped dust mass, and the dust contrast at the same orbital location. We apply our model to the recently imaged Oph IRS 48 system, finding values within the range of the observational uncertainties.

Publication: Lyra & Lin 2013.  


Debris disks

Structure seen in debris disks - sharp rings, eccentricities, asymmetric ansae - are usually taken as evidence for the gravitational influence of unseen planets, as these disks are mostly gas-free systems of dust and leftover planetesimals. In this simulation we show that if the system has gas, hydrodynamical forces at work suffice to explain the structure seen in debris disks. This gas comes from outgassing of planetesimals and dust grains via sublimation, photodesorption, or collisions, generating a system of dust-to-gas ratio close to unity, where hydrodynamics cannot be ignored. The backreaction of the drag force from the gas onto the dust shepherds rings. The eccentricity can be simply explained by a standing wave propagating along the ring. The planet possibility, though thrilling, is not necessarily required to explain these systems.

Publication: Lyra & Kuchner 2013.  


Ice shell convection in Europa

We solve the equations of Stokes flow with ice rheology, in a 2D grid ranging from the top of the subsurface ocean to the surface of the moon. The incompressibility condition is satisfied by defining a streamfunction, and solving the elliptic problem iteratively via successive over-relaxation. Temperature evolution is calculated via a high-order non-conservative finite-difference stencil, as in the Pencil Code.

The parameter A is related to ice grain size and is roughly equivalent to the order of magnitude change in viscosity from top to bottom. The Reynolds number is very small, yet the Rayleigh number is very large (~2e7), so convection is the dominant mechanism of heat transfer. In some models, buoyant warm ice reaches the surface, leading to topography of ~100m, consistent with the terrain's morphology observed by Voyager and Galileo. Video and simulation by Leonardo Sattler, undergraduate student from UFRJ/Berkeley, as a result of his 10 week summer internship with me at the Jet Propulsion Laboratory (JPL). Sattler continues working on the project for his undergraduate thesis.  


Star formation


Ages derived from the low mass stars still contracting onto the main sequence often differ from ages derived from the high mass ones that have already evolved away from it. In this work we investigated the general claim of disagreement between these two independent age determinations by presenting UBVRI photometry of the young galactic open clusters NGC 2232, NGC 2516, NGC 2547 and NGC 4755, spanning the age range ~10-150 Myr. The LHS figure shows NGG 2516, also known as the Southern Pleiades.

We derived reddenings, distances and nuclear ages by fitting ZAMS and isochrones to color-magnitudes and color-color diagrams. To derive contraction ages, we used four different pre-main sequence models, with an empirically calibrated color-temperature relation to match the Pleiades cluster sequence. When using exclusively the V vs. V-I color-magnitude diagram and empirically calibrated isochrones, there is consistency between nuclear and contraction ages for the studied clusters (RHS figure). Although the contraction ages seem systematically underestimated, in none of cases they deviate by more than one standard deviation from the nuclear ages. The work was mentioned in Virginia Trimble's Astrophysics in 2006.

Publication: Lyra et al. (2006), A&A, 453, 101 [arXiv] [ADS]  


Stellar Activity

Chromospheric activity is usually measured by the emission in CaII H & K lines. In this work, we presented a calibration of H-alpha as both chromospheric diagnostic and age indicator. The chromospheric diagnostic was built with a statistically significant sample, including 175 solar neighborhood stars. For the age indicator, we used stars from clusters and kinematic groups, whose ages are better determined. The age vs activity relation in H-alpha is shown in the RHS figure. An exponential decay is seen, leveling at < 2 Gyr.

Publication: Lyra & Porto de Mello (2005), A&A, 431, 329. [arXiv] [ADS]  


Alpha Centauri

The alpha Centauri binary system, owing to its duplicity, proximity and brightness, and its components' likeness to the Sun, is a fundamental calibrating object for the theory of stellar structure and evolution and the determination of stellar atmospheric parameters. This role, however, is hindered by a considerable disagreement in the published analyses of its atmospheric parameters and abundances. In this work we presented a new spectroscopic analysis of both components of the alpha Centauri binary system, and compared results from other published analyses. The LHS figure shows the element abundances in a box-and-whiskers plot of all published analyses up to this work. The spread is of 0.1 dex in metallicity among published works.

The analysis is differential with respect to the Sun, based on high-quality spectra, and employed spectroscopic and photometric methods to obtain as many independent Teff determinations as possible. The atmospheric parameters are also checked for consistency against the results of the dynamical analysis and the positions of the components in a theoretical HR diagram. We discuss possible origins of discrepancies, concluding that the presence of NLTE effects is a probable candidate, but we note that there is as yet no consensus on the existence and cause of an offset between the spectroscopic and photometric Teff scales of cool dwarfs. The spectroscopic surface gravities also agree with those derived from directly measured masses and radii. The abundance pattern can be deemed normal in the context of recent data on metal-rich stars. The position of alpha Cen A in an up-to-date theoretical evolutionary diagram (RHS figure) yields a good match of the evolutionary mass and age with those from the dynamical solution and seismology.

Publication: Porto de Mello et al. (2008) A&A, 488, 653. [arXiv] [ADS]