Modelling the dust in the HD100546
protoplanetary disk across all ALMA scales
Michael Stroet
11293284
Bachelor Thesis Physics and Astronomy
Supervisor Prof. Dr. Michiel Hogerheijde
Examiner Prof. Dr. Carsten Dominik
Size 15 EC
Conducted between May 2020 and August 2020
University University of Amsterdam
Vrije Universiteit Amsterdam
Faculty Faculty of Science
Institute Anton Pannekoek Institute for Astronomy
Populairwetenschappelijke samenvatting
Al ver voor de tijd van de oude Grieken en Egyptenaren bestudeerde men de planeten aan de hemel. Maar de enige planeten die zij konden waarnemen bevonden zich in ons eigen zonnes-telsel. Dit maakte het onderzoek doen naar de oorsprong van planeten vrij lastig aangezien onze planeten al gevormd zijn. Pas sinds de laatste jaren van de 20e eeuw, bezitten wij de tech-nologie om planeten te ontdekken die om andere sterren ver buiten ons zonnestelsel draaien. In de moderne tijd weten wij dat ons zonnestelsel acht planeten bevat. Vanaf de zon gezien zijn dat Mercurius, Venus, de Aarde en Mars die allemaal een vaste rotsachtige vorm hebben en de gasreuzen Jupiter, Saturnus, Uranus en Neptunus. De rotsachtige planeten en de gasreuzen zijn vanelkaar gescheiden door de astero¨ıdengordel. Maar astronomen hebben ontdekt dat de exoplaneten van andere sterren deze nette scheiding niet altijd hebben. Sterker nog, ons zon-nestelsel valt op tussen de andere zonzon-nestelsels. Een veel voorkomend verschijnsel is een ”hete Jupiter”, planeten zoals Jupiter die heel dicht bij hun ster staan, dichter dan Mercurius bij onze zon. Om er achter te komen hoe zulke zonnestelsels kunnen ontstaan, moeten wij planeetvorm-ing onderzoeken.
Wanneer er in de ruimte een nieuwe ster wordt geboren door het instorten van een enorme wolk van gas en stof, vormt het overgebleven materiaal een grote schijf om de ster heen. In zo’n schijf kunnen de kleine stofdeeltjes samenklonteren totdat zij de beginselen van planeten vormen, pro-toplaneten. Degelijke schijven worden daarom ookwel protoplanetaire schijven genoemd. Deze protoplaneten kunnen dan verder groeien tot rotsachtige planeten of kunnen gas uit de schijf verzamelen om uit te groeien tot gas reuzen. Oftewel, de verdeling van deze stofdeeltjes is erg belangrijk voor de vorming van planeten.
Daarom richtte ik mij in mijn project op het bestuderen van de verdeling van stofdeeltjes in de protoplanetaire schijf rondom de ster HD100546. Aan de hand van een computermodel maak ik een afbeelding van de schijf die dan wordt vergeleken met echte waarnemingen van de ALMA telescoop hoog in de gebergten van Chili. Een algoritme kijkt naar de vergelijking van het model en de data en maakt kleine aanpassingen aan het model om te kijken of zij dan beter overeenkomen. Uit de uiteindelijke resultaten vind ik dat de schijf van HD100546 een ring van stofdeeltjes bevat tussen ongeveer 20 en 50 astronomische eenheden (AE), waar 1 AE gelijk staat aan de afstand tussen de zon en de aarde. Vooralsnog is het nog niet gelukt om de verdeling van stof in de ring zelf concreet vast te stellen. In de toekomst zal er meer onderzoek gedaan moeten worden om de stofverdeling duidelijker in kaart te brengen.
Abstract
Context. The orderly lineup of planets in our solar system, with four inner terrestrial planets and four outer giants, is not the standard configuration across the universe. Planets form inside the interplanetary disks of newborn stars by growing from small dust particles to larger plan-etesimals. The distribution of the dust particles plays a major role in formation of planets. The ALMA telescope is able to resolve these disk and allows us to observe the planetary formation process directly.
Aims. This work aims to model the dust distribution in the interplanetary disk of Herbig Ae/Be star HD100546 using two archival 1.3 mm continuum images taken by the ALMA telescope. The images have to be observed with different baseline configurations of the radio antennas, which enables the study of both the small and large scale dust structures of the disk.
Methods. Using a Markov Chain Monte Carlo (MCMC) method, a variety of models will be fitted to the data images in order to find the model and fit parameters that best describe the dust distribution. The long baseline image will be modelled by a power law and a Gaussian model, whereas the small baseline image will be modelled by another power law extending from the long baseline power law model.
Results. The small scale images were able to be fitted by a ring of dust extending from either 19.8+3.3−5.5 to 55.0+12.1−8.8 AU for the power law model or from 22.0+2.2−6.6 to 60.5+28.6−12.1 AU for the Gaussian model. The fit results for the distribution of dust inside this ring, however, were in-conclusive, neither was any evidence found for an extended power law of dust beyond the ring.
Conclusions. The results provide values for the extent of a ring of dust surrounding the HD100546 star, but were unable to confidently provide an understanding of the distribution the dust inside the ring. Further research could look into the data for an understanding of the distribution of flux of the u-v visibilities across the different scales. Another way forward is the use of more precise data to try and fit the same or different models to the dust distribution in the disk.
Contents
Populairwetenschappelijke samenvatting 3 Abstract 4 1 Introduction 6 2 HD100546 7 2.1 Observational history . . . 7 2.2 Protoplanetary disks . . . 8 2.2.1 General structure . . . 8 2.2.2 The HD100546 system . . . 9 2.3 Observations . . . 10 2.3.1 ALMA telescope . . . 10 2.3.2 Database images . . . 11 2.3.3 Data processing . . . 123 Fitting the intensity profile 14 3.1 Modelling the dust continuum intensity . . . 14
3.1.1 Intensity model . . . 14
3.1.2 Radial intensity profile . . . 15
3.2 Fitting with a MCMC method . . . 16
3.2.1 Motive for MCMC . . . 16 3.2.2 Bayes’ theorem . . . 17 3.2.3 Likelihood probability . . . 18 3.2.4 Prior probability . . . 18 3.2.5 Initialisation . . . 19 4 Results 21 4.1 Long baseline results . . . 22
4.2 Small baseline results . . . 25
5 Analysis 27 5.1 Small scale ring structure . . . 27
5.2 Large scale extended structure . . . 28
5.3 Future recommendations . . . 28
6 Conclusion 29
Acknowledgements 29
Chapter 1
Introduction
For centuries the only laboratory astronomers had for the study of planets and planetary form-ation was our own solar system. Only since the first detection of a planet in orbit around a star other than the sun (Wolszczan and Frail, 1992) have we been able to study the planetary formation process by observing exoplanets. What we found, however, was that our orderly solar system with four inner terrestrial planets and four outer gas giants is not the standard configuration of solar systems. A popular example of such disparate systems is the existence and abundance of hot jupiters. Hot jupiters are gas giants similar to Jupiter, but in an orbit incredibly close to their star, closer than the orbit of Mercury around the Sun. A hot jupiter was the first exoplanet found orbiting a sun-like star (Mayor and Queloz, 1995). A discovery for which Mayor and Queloz were awarded the 2019 Nobel Prize in Physics. The existence of such planets suggested that they can either form in different places than we thought or are able to migrate to a radically different orbit. To understand the history of our own solar system and planetary formation as a whole, the formation of planets around newborn stars must be studied.
When a molecular cloud of gas and dust collapses, leading to the formation of a star, the leftover material falls into a orbit around the newly born protostar, forming a circumstellar disk. Circumstellar comes from the latin words circum stella meaning ’around a star’. Tiny dust particles inside this disk are able to grow into larger planetesimals, the building blocks of all planets. For that reason, circumstellar disks like these are often called protoplanetary disks. If a large volume of gas is present around the orbit of a protoplanet, it is able to grow even larger into a gas giant. The initial distribution of the dusk particles inside of protoplanetary disks is therefore of profound importance to the formation of planets.
In this project, my goal is to model the dust distribution inside the protoplanetary disk around the Herbig Ae/Be star HD100546. Several models of the distribution of dust will be created and compared to archival images taken by the Atacama Large Millimeter Array (ALMA) telescope in Chile. To be able to study both the small scale and large scale structures of the distribu-tion of dust, the archival images must have been observed with different configuradistribu-tions of the ALMA radio antennas. These models will then be fit to the data using an ensemble based Markov Chain Monte Carlo (MCMC) algorithm as implemented by the Python library emcee (Foreman-Mackey et al., 2013).
This thesis is structured by the following outline:
Chapter 2 portrays the history and characteristics of protoplanetary disks and of HD100546 in particular along with the ALMA observations and required data processing. Chapter 3 describes the dust distribution models and the methods used to fit the models to the data. In chapter 4 the results will be presented for each of the models. The results are then analysed in chapter 5. And finally, chapter 6 will conclude the thesis, summarising all that has been done and learned throughout the project.
Chapter 2
HD100546
In recent years, the protoplanetary disks around newborn stars have been an extensively studied topic (Grady et al. (2001), Dullemond and Dominik (2004), Hogerheijde et al. (2016)). An in-teresting example of a protoplanetary disk is the disk around the Herbig Ae/Be star HD100546. Modelling the disk system of this relatively young star is, as mentioned before, the main goal of this project.
This chapter is structured as follows:
Section 2.1 lays out the observational history of stars like HD100546 up to the discovery of HD100546 itself. Section 2.2 will then describe protoplanetary disks in general and give a more detailed description of the disk system of HD100546. Lastly, section 2.3 will focus on the ALMA telescope and the archival continuum images of the HD100546 disk used to model the dust distribution.
2.1
Observational history
When the Italian astronomer Angelo Secchi was observing the spectra of several stars in 1866, he found the first emission lines in the spectrum of a star, namely γ-Cassiopeiae (Secchi, 1866). Γ-Cassiopeiae is a variable spectral type B star showing strong hydrogen emission lines, espe-cially the Hα and Hβ transitions from the Balmer series. Stars like γ-Cassiopeiae are now known as classical Be stars. In 1931 Otto Struve proposed that the emission lines of these Be stars come from hot material surrounding the star (Struve, 1931). This material, mainly consisting of hydrogen, is ionised by the stellar radiation and re-emits the radiation at the specific Balmer frequencies. In the case of classical Be stars, the surrounding material originates from the star itself. These stars have rotational velocities upwards of hundreds of kilometres per second, producing an equatorial bulge big enough to allow material to escape from the surface of the star.
In later years, stars with emission lines in their spectra were also observed in star-forming clouds or nebulae. The emission lines differ slightly from those found in classical Be stars, suggesting a different origin. These emission lines come not from material of the star itself, but from a disk of gas and dust leftover from the formation of that star. There are several classifications for stars like this, such as T Tauri stars (Joy, 1945) or Herbig Ae/Be stars (Herbig, 1960). The modern criteria for distinguishing Herbig Ae/Be stars from other stars are slightly different than those originally defined by Herbig. The three criteria now are:
1. The star must have a spectral type earlier than F0.
2. The stellar spectrum must show Balmer emission lines.
The spectral type of stars is a measure for the effective temperature of a star. The effective temperature of a star is the temperature a black body needs to have to emit the same amount of light as that star and is a reasonable estimate for the surface temperature of a star. The standard Harvard classification of stellar spectral types is a series of classes from O to M. For historical reasons, the spectral types are said to go from early to later types. The full series goes from the earliest O to B, A, F, G, K and finally M. Herbig Ae/Be stars are, as the name sug-gests, A or B type stars with emission lines. Having spectral types earlier than F0 distinguishes them from T Tauri, which are less massive stars with spectral types later than F0. Criteria number 3 then distinguishes Herbig Ae/Be stars from classical Be stars. The excess infrared radiation originates from tiny dust particles, which are not present in the surrounding material of classical Be stars.
In 1989 J.Y. Hu and others were observing three candidate Herbig Ae/Be stars: HD37411, HD104237 and HD100546. (Hu et al., 1989). They identified all three candidates as important Herbig Ae/Be stars, especially HD104237 and HD100546 due to their brightness, allowing for ultraviolet and high resolution visual observations to be taken.
2.2
Protoplanetary disks
2.2.1 General structure
Protoplanetary disks consist of the gas and dust leftover from the formation of a low-mass star. The newly born star is not yet able to fuse hydrogen in its core, preventing the material in the disk from being blown away by the stellar wind. The disk is made up of tiny mm/cm sized dust particles in the central midplane with even smaller sized particles further out. The dust is encased in an atmosphere of gas, mainly consisting of hydrogen and helium along with some other elements like oxygen and carbon.
Figure 2.1: A schematic illustration of the cross-section of a protoplanetary disk (Andrews, 2020). Gas inside the disk is shown as a greyscale area with the dust particles as coloured circles (size not to scale). The right side highlights the general distribution of gas and dust in the disk, whereas the left side roughly shows the origin of the different emission tracers.
There are three main ways of observing gas and dust in disks. High resolution examples of all three methods are shown in figure 2.2. The first emission tracer is the molecular emission spectra from gas molecules in the outer regions of the gas atmosphere. Gas molecules with a permanent dipole moment, like carbon monoxide (CO) or water vapour (H2O) emit light at
specific frequencies, visible in the spectrum as sharp peaks at those frequencies. However, since the most abundant gas molecule is molecular hydrogen (H2), which has no permanent dipole
moment, most of the gas is undetectable with this method. The second emission tracer is light reflected from the µm-sized dust particles on the surface of the disk. The appearance of the scattered light is constrained by the size, shape, composition and albedo’s of the particles. This method benefits from high resolution imagery, but suffers from contrast with the host star, limiting observation of the inner region of the disk. The last tracer is the continuum emission emitted by the larger dust particles in the midplane. The continuum emission is the thermal radiation emitted by the particles. The continuum intensity is related to the surface density of dust as seen from Earth. This method allows for high resolution images without being limited by the contrasting light from the star, making it the most commonly used method for describing disk structures. The downside of continuum emission, however, is a lack of detailed information on the properties of the particles themselves.
Figure 2.2: Three high resolution observations of the TW Hya protoplanetary disk, each image spans 500 AU on either side and was created by a different emission primer. This collection of images also originates from Andrews (2020) and shows: [a] a 1.6 µm wavelength scattered light image (Van Boekel et al., 2017), [b] a 0.9 mm wavelength continuum emission image (Andrews et al., 2016) and [c] an image from the CO J=3-2 emission lines (Huang et al., 2018).
Just like planets, protoplanetary disks come in all shapes and sizes. The most common struc-tures are inner cavities and rings/gaps. In the continuum image of figure 2.2b, several rings and gaps can be seen. Other structures include spirals and asymmetric arcs, possibly caused by protoplanets dragging material behind them as they orbit the star.
2.2.2 The HD100546 system
The disk system of HD100546 is a famous and well-studied example of a protoplanetary disk. HD100546 is, as previously mentioned, a Herbig Ae/Be star with a mass of 2.2 ± 0.2M (Pineda
et al., 2019) at a distance of 110.02 ± 0.62 pc, or roughly 360 light-years from our solar system (Brown et al., 2018). It is best known for HD100546b, an enormous protoplanet orbiting the star at approximately 53 AU. It is one of the largest exoplanets found thus far with the planet and its accompanying disk together having a radius roughly seven times that of Jupiter (Quanz et al., 2015) and a mass of about 20 Jupiter masses (Acke and Van den Ancker, 2006).
Hd100546b was also the first ever direct observation of an exoplanet that has not yet been fully formed due to still being embedded in a disk of gas and dust of its own. Evidence for the existence of this planet was first given in Acke and Van den Ancker (2006) and was first directly
observed in Quanz et al. (2013). Planetary formation research only had access to computer simulations before this discovery, but now, hypotheses could be directly tested by observing real protoplanets.
The protoplanetary disk of HD100546 contains a cavity of both gas and dust up to 14 AU from the star, with a ring of dust extending from this cavity up to roughly 50 AU (Pineda et al., 2019). The gas in the disk extends several hundred AU and shows signs of asymmetry in the midplane (Miley et al., 2019) and possibly a spiral feature (Boccaletti et al., 2013). Signs point to another possible gas giant forming inside the cavity (Walsh et al., 2014), but its existence has not yet been proven. Another gas giant forming inside the disk would lead to HD100546 being a unique opportunity to study the formation of gas giants.
2.3
Observations
2.3.1 ALMA telescope
Although protoplanetary disks are up to hundreds of AU in diameter, their distance to the solar system causes their angular size to measured in arcseconds. An arcsecond is a very small angle equal to one 3600th of a degree. To better understand an angle of a single arcsecond, imagine looking at a football field from far away. To let the field cover only one arcsecond of your vision, you need to be standing at a distance of 20,000 kilometres away. Human eyes are too small to see such tiny objects, so we need to use a telescope. The diameter D a telescope needs to have in order to be able to resolve objects of a certain angular size θ can be estimated with the Rayleigh criterion: θ = 1.22Dλ. For a telescope observing at a wavelength λ of 1 mm, it needs to have a diameter of approximately 250 m to be able to resolve 1 arcsecond objects. That being impractical, a different method is needed: interferometry.
In astronomy, interferometry is a complicated technique for simulating a larger telescope with several smaller telescopes spread across a large area. A recent and famous example of the use of an interferometer is the Event Horizon Telescope (EHT), a global collection of radio telescopes used to acquire the first direct image of a black hole (Akiyama et al., 2019). One of the radio telescopes used in the EHT is the Atacama Large Millimeter Array (ALMA). The ALMA telescope is a radio interferometer located on the 5 km high Chajnantor plateau in the Atacama Desert of northern Chile. Since its official opening in march 2013, ALMA has been used for researching varying subjects like comets (Cordiner et al., 2014) and planetary formation (Stephens et al., 2014).
ALMA consists of 50 main 12-metre radio antennas with an additional four 12-metre and twelve 7-metre antennas. All antennas are able to be moved from a minimal distance between two an-tennas of 9 metres up to a maximum of 16 kilometres. The distance between anan-tennas is known as the baseline, thus the antennas in figure 2.3 are shown in a small baseline configuration. Using ALMA in a long baseline configuration simulates a much larger telescope, resulting in a much higher resolution image. However, this also means that the actual telescopes cover less of the total area of the simulated telescope than when using a small baseline configuration. Ad-ditionally, the lack of small baselines also means a lack of large scale structure in the resulting image.
How interferometry works exactly is rather complicated and beyond the scope of this bachelor thesis, but the result of the interferometry and the necessary data reduction is an image with a resolution given by the size of the beam. The ALMA beam is an elliptical area in which
Figure 2.3: A 2013 aerial photo of the ALMA radio telescopes close together. Image credit: Clem & Adri Bacri-Normier (wingsforscience.com)/ESO
ALMA is able to observe. The beam effectively distorts the ’real’ shape of a disk in the sky and the resulting image been smeared out. Although all telescopes have such a beam, the effect is usually not significant to the quality of the image. That is not the case here, so the distortion caused by the beam needs to be corrected for. The method used for this will be explained in section 2.3.3.
2.3.2 Database images
To model the dust distribution in the HD100546 disk for all possible ALMA scales, two images of different baselines of the 1.3 mm continuum will be used. The two images consist of a small baseline image for the large scale structure and a long baseline image for the small scale structure. Both images are taken from the ALMA FITS Archive from the Japanese Virtual Observatory project and are shown in figure 2.4.
The small baseline image is taken from Miley et al. (2019), with project ID 2015.1.01600.S, who redid the data reduction and calibration from the observations described in Walsh et al. (2014). The observation included 24 antennas with baselines from 21 to 375 m. The resulting image covers a frequency range from 219 to 235 GHz. The beam of the image has a major axis of 1.11 arcseconds and a minor axis of 0.83 arcseconds. The beam has a position angle (PA) of -21.1 degrees measured east from north. The position angle is a measure for the angle on the sky where an upwards line (north) has angle 0 and increases counter-clockwise (west). The small baseline image, zoomed in to a square of 12 by 12 arcseconds, is shown on the left in figure 2.4 with the size and orientation of the beam indicated by the white ellipse.
The observation of the long baseline image is described in P´erez et al. (2020), with project ID 2016.1.00344.S. The observation used 39-42 antennas with baselines up to 12.2 kilometres and the resulting image covers a frequency range from 217 to 219 GHz. The beam of the image has a major axis of 0.25 arcseconds and a minor axis of 0.13 arcseconds and has a PA of -14.1 degrees measured east from north. The small baseline image, zoomed in to a square of 6 by 6 arcseconds, is shown on the right side of figure 2.4.
Figure 2.4: Left: A small baseline ALMA image of the 1.3 mm continuum of HD100546, project ID 2015.1.01600.S (Miley et al., 2019). The beam size is 1.11” by 0.83” with a PA of −21.1o east from
north. The original image was 60” across on either side and has been zoomed to 12” across either side. Right: A long baseline ALMA image of the 1.3 mm continuum of HD100546, project ID 2016.1.00344.S (P´erez et al., 2020). The beam size is 0.25” by 0.13” with a PA of −14.1o east from north. The original image was 45” across on either side and has been zoomed to 6” across either side.
2.3.3 Data processing
Before comparing a model to the continuum images, the images for the distortion effect of the beam needs to be corrected. To do this, we distort the image a second time with another elliptical area. This area is generated as a 2D Gaussian ellipse that has the same size as the beam but is oriented perpendicular to the beam itself. The perpendicular beam is then applied as a kind of filter to every pixel in the image. The process of applying a filter to an image is called convolution. The perpendicular beam is used as the convolution matrix, also known as the kernel, for the convolution process. Both the small and long baseline images are convolved with their respective perpendicular beam. The result of the convolution process is shown in figure 2.5 for the long baseline image.
The final step is to convert the convolved data image into a one-dimensional graph of the intensity profile of the disk. A one-dimensional graph allows for an quick comparison between the data and a model so the quality of the model can be evaluated quickly. The intensity profile will be generated as the average continuum intensity at different radii from the centre of the disk, where we assume the disk to be circularly symmetric. However, since the disk of HD100546 is not viewed face-on, but with an inclination of roughly 42 degrees (Pineda et al., 2019), the circular disk appears elliptical. The eccentricity of an inclined circle is equal to the sine of the inclination, which in the case of HD100546 is approximately 0.67 and the disk is orientated with a PA of 145.14 ± 0.04 degrees east of north (Pineda et al., 2019). For each radii on the radial intensity profile an ellipse was drawn, centred on the disk, with a major axis equal to that radius and a minor axis calculated with the eccentricity and orientated with the disk PA. The average intensity along the border of the ellipse is calculated and added to the profile. This process is done for both images and the radial profile of the long baseline image is shown in figure 2.6. The profiles are normalised to a peak of 1 since we are only interested in the shape of the curve and not the real intensity values.
Figure 2.5: Left: The long baseline image of P´erez et al. (2020) with the relative declination (DEC) and right ascension (RA) up to 1.5 arcseconds from the centre of the disk. Right: The same image when convolved with a 2D Gaussian kernel with the same full with half maximum as the ALMA beam orientated perpendicularly.
Figure 2.6: The normalised radial intensity profile generated from the convolved long baseline image of P´erez et al. (2020).
Chapter 3
Fitting the intensity profile
This chapter describes the method for fitting different models to the radial intensity profile generated from the ALMA images by making use of an ensemble sampler Markov Chain Monte Carlo (MCMC) method from the emcee Python library.
Section 3.1 will describe the different models and how they are used to create radial intensity profiles. Section 3.2 will then explain the MCMC algorithm and describe how it is used to find the best fit parameters.
3.1
Modelling the dust continuum intensity
3.1.1 Intensity model
In order to model the intensity profiles of the 1.3 mm dust continuum images, a formula for the intensity at a given radius is needed. The model for the continuum intensity is inspired by Hogerheijde et al. (2016), in which they modelled the intensity of the 820 µm continuum of the TW Hya disk. The intensity I as a function of the radial distance from the centre R at frequency ν is modelled by
Iν(R) = Bν(T (R))(1 − e−τ (R)), (3.1)
where Bν is the Planck function, T (R) the dust temperature at radius R and τ (R) the optical
depth of the dust at radius R. The Planck function relates the thermal radiation a black body emits to its temperature for different frequencies of light and is expressed as
Bν(T (R)) =
2hν3
c2 (e
hν/kBT (R))− 1)−1. (3.2)
The temperature of the dust in the disk is modelled by an exponential decay using a power law given by T (R) = T0 R Rin −q , (3.3)
where T0 is the temperature of the disk at the inner disk radius Rin and q the exponential decay
factor. The optical depth τ is a measure for the absorption of light travelling through a medium and how the absorbed light is re-emitted. For dust in a protoplanetary disk, the optical depth is expressed as the product of surface density of dust Σd of the disk, the dust opacity at the
observing wavelength κν and the cosine of the inclination of the disk i,
The dust surface density Σdis a function for the distribution of dust throughout the disk when
seen as a flat 2D surface. Since finding the Σd in the HD100546 disk is the goal of this thesis ,
several different models for Σdwill be used. To model the small scale structure seen on the long
baseline image, Σd will be modelled by either a power law or a Gaussian function between an
inner radius of the ring Rin and an outer radius Rout. The power law model is a simple single
power law similar to equation 3.3 and is given by
Σd(R) = Σ0 R Rin −p , (3.5)
where Σ0 is the dust surface density at the inner radius in kg/m2 and p an exponential decay
factor. The Gaussian model is given by
Σd(R) = Σ0· exp " −1 2 R − Rin c 2# , (3.6)
where Σ0 is the same as in equation 3.5 and c the standard deviation of the Gaussian.
To model the large scale structure beyond the ring, an extended power law will be used in tandem with the results of the power law from equation 3.5. The model is expressed as
Σd(R) =
(
Σ1(R/Rring)−p1, R ≤ Rring
Σ2(R/Rring)−p2, R > Rring
, (3.7)
where the first power law models the intensity from the small scale ring structure and the second power law models the dust beyond the ring. At radius Rring, given by the outer radius of the
ring from the power law model, the model discontinuously switches from surface density Σ1 to
Σ2 and continues on with a different decay factor from p1 to p2.
3.1.2 Radial intensity profile
Using equation 3.1, a model image can be created by calculating the distance from the centre for every pixel and then calculate the intensity at that distance. The size of pixels in this image is an arbitrary choice, effectively causing the model to have an infinite resolution. In order to equate the model resolution to the resolution of the convolved data from section 2.3.3, the model has to be convolved with the same kernel as the data. Unlike the convolution of the data, however, the beam itself must also be included. This is because observing the disk already effectively convolved the resulting image with the beam. By convolving the beam and its per-pendicular form together and then normalising the result, a circular kernel is created with which the model image will be convolved. An example of the result of this process is shown in figure 3.1.
With a model image of the same resolution as the convolved data image, a radial profile is able to be generated for comparison to the data. Since the model is a face-on projection of the disk, where we drew ellipses on the data image, we can now use circles. Similar to the process in section 2.3.3, several concentric circles are drawn on the model image and the average intensity along each circle is calculated and saved as the radial intensity profile. The resulting profile for the long baseline model is shown in figure 3.2 together with the data profile from figure 2.6.
Figure 3.1: Left: A model image of arbitrary resolution created with equation 3.1 using the power law model from equation 3.5. Right: The result of convolving the model image with a circular kernel created from the ALMA beam and its perpendicular form. Both images show the modelled continuum intensity at the relative right ascension (RA) and declination (DEC) from the centre of the disk in arcseconds.
Figure 3.2: The radial intensity profiles generated from the convolved data and the convolved model of figure 3.1. Both profiles are normalised with their peaks at 1.
3.2
Fitting with a MCMC method
3.2.1 Motive for MCMC
The model profile in figure 3.2 obviously does not match the data. So in order to find the parameters that best reproduce the data, we have to fit the model to the data. Say you tried checking every possible combination of parameter values and comparing all the resulting
intens-ity profiles with the data. This guarantees that the best possible solution is found. A simple example will, however, make clear why this course of action is not the wisest choice.
Before a model with some choice of parameters values can be compared to the data, a model image needs to be created, then that images must be convolved with the circular kernel and finally, a radial profile must be generated. This is a quite CPU intensive process for a computer. Say this process takes roughly one second to complete and we want to try and fit a model with three parameters where each parameter has a total of 1000 possible values. Checking every possible combination then takes a total of 10003∗ 1 = 109 or 1 billion seconds, which is nearly
32 years. Trying even more parameters or more possible parameter values increases the time needed considerably. Having only a total of one million parameter values and three parameters already allows for this time to be measured in ’universes’, 13.8 billion years each.
In order to drastically reduce the time needed, a Markov Chain Monte Carlo (MCMC) method will be used. MCMC methods are a collection of different algorithms for finding the most likely parameters that best fit a model to some data through the use of Bayesian statistics. These methods often see use in astronomical research because of the method’s capability of handling many unknown variables. For example, Hadden and Lithwick (2017) used MCMC simulations with N-body integrators to find the planetary characteristics and orbital parameters of multiple exoplanets observed by the Kepler space telescope. For modelling the dust in the HD100546 disk an ensemble sampler based MCMC method will be used as proposed by Goodman and Weare (2010) and implemented in the emcee Python library by Foreman-Mackey et al. (2013).
3.2.2 Bayes’ theorem
The math behind MCMC methods is quite complicated, but a basic understanding of Bayesian statistics is needed in order to understand how and why MCMC methods work.
Bayesian statistics is an interpretation of probability in which the probability of an event is not only determined by some given observations, but also by the prior belief of the probability of that event occurring. Imagine there is a disease that we believe affects one percent of the population. The test for this disease has a 90% chance of correctly identifying the disease. If you were to test positive, what is your probability of actually having the disease? While you might say 90%, that is actually false. In a population of 100 people, there is likely only one person affected by the disease. If we were to test everyone, the diseased person will be correctly identified, but another ten people will receive a false positive. Among the eleven people tested positive, your probability of having the disease is then actually 1 in 11, or 9%.
This interpretation of probability is expressed by Bayes’ theorem, first published in Bayes (1763) by Richard Price, and is commonly written as
P (A|B) = P (B|A) × P (A)
P (B) , (3.8)
where P (A|B) is the probability of an hypothesis A being true given event B. P (A|B) is known as the posterior probability and is calculated as the product of the likelihood of event B being true if A is true (P (B|A)) and the belief, or prior probability, of A being true (P (A)) all divided by the probability of event B (P (B)). Another way of saying this is: ”The posterior probability is proportional to the likelihood times the prior probability”, which can be expressed as
Posterior probability ∝ Likelihood probability × Prior probability. (3.9)
The strength of using Bayes’ theorem and MCMC methods comes from the ability to update the prior beliefs with the calculated posterior probability, allowing the calculation of an even more accurate posterior probability.
The goal of the emcee MCMC algorithm is to acquire a posterior distribution function for each parameter of the model. To do this, emcee uses an ensemble of so-called walkers. Each of these walkers have a value for each of the parameters. Each walker can be viewed as having a position inside a multi-dimensional parameter space in which each coordinate is one of the parameters of the model. Each position in this space is tied to a specific intensity profile able to be generated by the model. The posterior probability at a certain position is used to indicate the quality of the model with the corresponding parameters. For every iteration in the emcee MCMC algorithm the posterior probability is calculated for every walker, which is then used to update the prior probabilities of the walkers and move all of them to a new position in the parameter space. A more detailed description of this process will be given in section 3.2.5 However, before we are able to calculate the posterior probability of equation 3.9, the likelihood and prior probabilities must first be determined.
3.2.3 Likelihood probability
The likelihood is represented as the chance of obtaining a particular intensity value by randomly sampling it from a Gaussian probability function around the intensity value of the data. For a single radius r, the likelihood function Plikely,r is expressed as
Plikely,r = 1 σ√2π · exp " −1 2 Idata,r− Imodel,r σ 2# , (3.10)
where σ2 is the data variance, determined by the average noise level of the data, and I the intensities of the data and model intensity profile at radius r. To get the total likelihood function for all radii from 0 to some radius R along the intensity profile, equation 3.10 must be multiplied for each radius. The likelihood function then becomes
Plikely =
YR
r=0Plikely,r. (3.11)
However, since we are multiplying numerous exponentials, this process can quickly become numerically unstable. To avoid this, we can take the natural logarithm of equation 3.11 to drop the exponentials and turn the multiplications into a sum. After taking the logarithm, the final log-likelihood function is ln (Plikely) = − 1 2 XR r=0 " Idata,r− Imodel,r σ 2 + ln 2πσ2 # . (3.12) 3.2.4 Prior probability
Since we are now working with the natural logarithms of probabilities, the posterior probability of equation 3.9 now has to be written as
ln (Pposterior) ∝ ln (Plikely) + ln (Pprior) . (3.13)
Probability is represented by a value between 0 and 1, so the natural logarithm of probabilities is now represented by values between ln(0) and ln(1). The value of ln(1) is 0 and the limit of ln(x) where x goes to 0 is −∞. Since the prior probabilities of the walkers will be constantly updated, the initial prior probability function can be seen as a range of values in which we believe a parameter is located. If each parameter value is located within its range, we give the prior a probability of 1 and if any of the parameters is outside its range, we give it probability 0. Then the final log-prior probability function becomes
ln (Pprior) =
(
0 all parameters inside range,
−∞ otherwise . (3.14) 3.2.5 Initialisation
With equations 3.13, 3.12 and 3.14 all together, we are now able to run the emcee MCMC algorithm. Before fitting the models to the data, the algorithm must first be initialised.
The first step is to decide which parameters will be fixed and which will be used to fit to the data, the free parameters. Since the intensity profiles are normalised to 1, we only care about the shape of the profile. Thus any parameter that only affects the scaling can be fixed to a single value. Additionally, the dust surface density is the only relevant function in modelling the dust distribution. The free parameters per model are then:
• Rin, Rout and p for the power law model of equation 3.5,
• Rin, Rout and c for the Gaussian model of equation 3.6,
• Rout, SigFrac = Σ2/Σ1 and p2 for the extended power law model of equation 3.7.
The values of the fixed parameters are shown under Fixed parameters in table 4.1. The second step is to define the free parameter ranges for the prior probability and to provide an initial guess of the free parameter values for the initialisation of the walkers. These values are shown in table 3.1.
Ppower law Initial guess Prior range
Rin (arcsec) 0.1 0.01 − 0.29 Rout (arcsec) 0.5 0.3 − 1 p 1 0 − 3 Gaussian Rin (arcsec) 0.1 0.01 − 0.29 Rout (arcsec) 0.5 0.3 − 1 c 1 0 − 5
Extended power law
Rout (arcsec) 5 2 − 6
SigFrac 0.5 0 − 1
p2 1 0 − 3
The third and final step is to initialise an amount of walkers and give the total number of steps the algorithm has to run. All walkers are initialised in a small distance from the initial guesses in table 3.1. From there the walkers explore the parameter space bounded by the prior ranges in the so-called burn-in phase. The amount of burn-in steps should be large enough to allow the walkers to explore the entire parameter space. The burn-in phase is followed by the production phase in which the walkers will try to converge to the areas inside the space with the largest posterior probability. The positions of the walkers throughout the entire production phase are bundled together in order to get a posterior distribution for each free parameter. If the walkers have converged for a parameter the posterior distribution of that parameter shall show a peak around the converged parameter value. The more productions steps the walkers have to converge, the higher the quality of the convergence becomes.
The amount of walkers, burn-in steps and production steps were chosen so that the entire fit process for a single model takes about 8 hours to complete. For all three models a total of 500 walkers will be used. The power law and Gaussian models both will run with 100 burn-in steps and 200 production steps, while the extended power law will use 200 burn-in steps and 300 production steps due to the lower resolution speeding up the process.
Chapter 4
Results
The long baseline continuum image of the disk of HD100546 has been modelled by power law and Gaussian ring models using an ensemble sampler MCMC method with 500 walkers, 100 burn-in steps and 200 production steps. The results of these models are presented in section 4.1. The small baseline continuum image has been modelled by the power law ring with an extended power law of dust using 500 walkers, 200 burn-in steps and 300 production steps and the result is presented in section 4.2.
The results of the emcee MCMC fitting process are shown in special corner plots. Each corner plot consists of several panels arranged in a triangular shape. The panels on the diagonal rep-resent the posterior distribution function of each free parameter as a histogram of the position of the walkers across all production steps. The vertical dashed lines show the locations of the median value and the 1σ deviations. The values of the dashed lines are shown on top in the format: median+1σ−1σ. The vertical blue lines represent the initial guesses for the free parameter values. The inner panels show the relative correlation between the different free parameters, showing the 0.5σ, 1σ, 1.5σ and 2σ contours.
4.1
Long baseline results
Figure 4.1: Resulting corner plot of the emcee MCMC algorithm using 500 walkers, 100 burn-in steps and 200 production steps on the long baseline image modelled by the power law shown in equation 3.5. The inner radius of the ring Rinhas been found to be at 0.18+0.03−0.05arcseconds and the outer radius Rout
Figure 4.2: Resulting corner plot of the emcee MCMC algorithm using 500 walkers, 100 burn-in steps and 200 production steps on the long baseline image modelled by the Gaussian shown in equation 3.6. The inner radius of the ring Rinhas been found to be at 0.20+0.02−0.06arcseconds and the outer radius Rout
Figure 4.3: A comparison of the radial intensity profile generated from the long baseline data and the power law and Gaussian models. The data is represented by a blue line and the model with the median parameters from figures 4.1 and 4.2 by a red line. In orange several other model intensities are shown with parameters chosen between the -1σ to 1σ range of the median values.
A summary of the long baseline fitting results are show in table 4.1 together with the values of the fixed parameters. The inner and outer radii values have been converted from arcseconds to astronomical units (AU) by using a distance of 110 parsec and the small angle approximation.
Parameters Ppower law Gaussian
Fixed parameters ν (GHz) 225 225 κν (m2/kg) 0.21 0.21 i (degrees) 42.46 42.46 T0 (K) 30 30 q 0.25 0.25 Σ0 (kg/m2) 0.25 0.25 Free parameters Rin (AU ) 19.8+3.3−5.5 22.0+2.2−6.6 Rout (AU ) 55.0+12.1−8.8 60.5+28.6−12.1 p 1.54+1.08−1.13 c 2.70+1.54−2.08
Table 4.1: Results of fitting the power law and Gaussian models to the long baseline images. The free parameter values are given by the median value of the posterior distribution functions and their −1 and +1 σ deviation.
4.2
Small baseline results
Figure 4.4: Resulting corner plot of the emcee MCMC algorithm using 500 walkers, 200 burn-in steps and 300 production steps on the long baseline image modelled by a extended power law. The inner radius Rin and outer radius Rout values are given in arcseconds and the c factor is dimensionless.
Figure 4.5: A comparison of the radial intensity profile generated from the small baseline data and extended power law model. The data is represented by a blue line and the model with the median parameters from figure 4.4. While the orange profiles from the 1σ deviations are plotted, they are barely visible. See figure 4.6 for the logarithmic version of this figure.
Figure 4.6: The logarithmic version of figure 4.5. The orange intensity profiles from the 1σ deviations are now visible.
Chapter 5
Analysis
5.1
Small scale ring structure
In the figures of section 4.1 the power law and Gaussian models are reasonably able to converge the inner and outer radii of the small scale ring of dust at roughly 0.2 arcseconds or 20 AU and 0.5 arcseconds or 55 AU. These values are consistent with those found in Pineda et al. (2019). The parameters for the distribution of dust within the ring itself, however, are not able to converge to a single value. When varying the p or c parameters on their own, only a small change in the resulting intensity profiles is observed. This implies that convolving the models removes most of the information about the distribution of dust within the ring. The loss of information during convolution is perhaps explained by the size of the ALMA beam relative to the disk. The major axis of the long baseline beam is longer than the entire width of the ring.
As a small test of this hypothesis, a model image was created using the inner and outer radii from figure 4.1 but with a step function for the intensity. The intensity was set to 1 for radii inside the ring and to 0 on the outside. The image was then also convolved with the original and perpendicular beam and an intensity profile was generated by averaging circles. The resulting intensity profile is shown in figure 5.1.
Figure 5.1: A comparison of the data intensity profile and the profile generated with a step function model.
The intensity profile in figure 5.1 does not accurately follow the data, but it is not far off. The fact that a simple step function is able to be this close suggests that the information loss during convolution is the likely cause for the lack of convergence. As long as you have a model containing a downward slope, it will be able to reasonably fit the data with the correct inner and outer radii.
5.2
Large scale extended structure
When looking at the results of section 4.2, one can clearly see the extended power law model was not able to adequately fit to the data. The model profile in figure 4.5 does not follow the data profile accurately and the Rout parameter in the corner plot of figure 4.4 does not converge
on a single value. The SigF rac and p2 parameters are respectively trying to be as small and large as possible. This is due to the fact that the model intensities are consistently larger than the data intensities. Lowering the SigF rac value and raising p2 both result in a lower intensity. However, as one can see in figures 4.5 and 4.6, the effect of the parameters on the intensity profiles is minuscule and only barely visible when using a logarithmic scale. These findings do not conclusively rule out the existence of an extended region of dust as it likely suffers from the same loss of information discussed in the previous section.
5.3
Future recommendations
There are a few options that can tried to better model the large scale dust distribution. Further research could take a look into the flux distribution in the u-v visibilities of the ALMA data images in order to find out if the large scale data is suppressed in some way. If this is the case, the extended power law model could be tried again to see if the results change.
As for the smaller scale, longer baseline images with beam sizes smaller than the thickness of the ring of dust are recommended. While this will result in higher resolution images with substantially longer MCMC run times, the resulting fits would be less susceptible to the in-formation loss of the convolution process, allowing the distribution of dust within the ring to be determined. Further research can then aim for either finding the correct model or excluding specific models which do not fit the data.
Chapter 6
Conclusion
The dust distribution in the protoplanetary disk of Herbig Ae/Be star HD100546 has been modelled using archival data of a long baseline image and a small baseline image, both of the 1.3 mm dust continuum. Models on the long baseline image were able to fit a ring of dust extending from either 19.8+3.3−5.5 to 55.0+12.1−8.8 AU for the power law model and from 22.0+2.2−6.6 to 60.5+28.6−12.1 AU for the Gaussian model. Both models show no convergence on the parameters dealing with the dust surface density and did not provide a definite answer on the distribution of dust inside the ring itself. The extended power law model I used to find an extended distribution of dust in the small baseline image has not delivered any confident evidence for the existence of such a distribution. Further research can take a in-depth look into the data itself and shed light on the flux distribution for different scales across the u-v visibilities. Additionally, even longer baseline observations can build on my work done during the project in order to provide a clearer understanding of the distribution of dust inside the ring by either using the same models or trying different models.
Acknowledgements
I would like to express my deep gratitude towards my supervisor Michiel Hogerheijde for his patient guidance and overall accessibility during my project. His enthusiastic attitude made working with him an enjoyable experience.
Special thanks to Christian Ginski and Kevin Lange for organising the protoplanetary disk group meetings. It was an honour to be a part of this group during the few months I joined them.
Lastly, I wish to thank my father for proofreading this thesis, providing me with numerous helpful suggestions and corrections.
References
Acke, B. and Van den Ancker, M. (2006). Resolving the disk rotation of hd 97048 and hd 100546 in the [oi] 6300 ˚a line: evidence for a giant planet orbiting hd 100546. Astronomy & Astrophysics, 449(1):267–279.
Akiyama, K., Alberdi, A., Alef, W., Asada, K., Azulay, R., Baczko, A.-K., Ball, D., Balokovi´c, M., Barrett, J., Bintley, D., et al. (2019). First m87 event horizon telescope results. iv. imaging the central supermassive black hole. The Astrophysical Journal Letters, 875(1):L4.
Andrews, S. M. (2020). Observations of protoplanetary disk structures. arXiv preprint arXiv:2001.05007.
Andrews, S. M., Wilner, D. J., Zhu, Z., Birnstiel, T., Carpenter, J. M., P´erez, L. M., Bai, X.-N., ¨
Oberg, K. I., Hughes, A. M., Isella, A., et al. (2016). Ringed substructure and a gap at 1 au in the nearest protoplanetary disk. The Astrophysical Journal Letters, 820(2):L40.
Bayes, T. (1763). Lii. an essay towards solving a problem in the doctrine of chances. by the late rev. mr. bayes, frs communicated by mr. price, in a letter to john canton, amfr s. Philosophical transactions of the Royal Society of London, (53):370–418.
Boccaletti, A., Pantin, E., Lagrange, A.-M., Augereau, J.-C., Meheut, H., and Quanz, S. P. (2013). Multiple spiral patterns in the transitional disk of hd 100546. Astronomy & Astro-physics, 560:A20.
Brown, A., Vallenari, A., Prusti, T., De Bruijne, J., Babusiaux, C., Bailer-Jones, C., Biermann, M., Evans, D. W., Eyer, L., Jansen, F., et al. (2018). Gaia data release 2-summary of the contents and survey properties. Astronomy & astrophysics, 616:A1.
Cordiner, M., Remijan, A., Boissier, J., Milam, S., Mumma, M., Charnley, S., Paganini, L., Villanueva, G., Bockel´ee-Morvan, D., Kuan, Y.-J., et al. (2014). Mapping the release of volatiles in the inner comae of comets c/2012 f6 (lemmon) and c/2012 s1 (ison) using the atacama large millimeter/submillimeter array. The Astrophysical journal letters, 792(1):L2.
Dullemond, C. and Dominik, C. (2004). The effect of dust settling on the appearance of protoplanetary disks. Astronomy & Astrophysics, 421(3):1075–1086.
Foreman-Mackey, D., Hogg, D. W., Lang, D., and Goodman, J. (2013). emcee: the mcmc hammer. Publications of the Astronomical Society of the Pacific, 125(925):306.
Goodman, J. and Weare, J. (2010). Ensemble samplers with affine invariance. Communications in applied mathematics and computational science, 5(1):65–80.
Grady, C., Polomski, E., Henning, T., Stecklum, B., Woodgate, B., Telesco, C., Pina, R., Gull, T., Boggess, A., Bowers, C., et al. (2001). The disk and environment of the herbig be star hd 100546. The Astronomical Journal, 122(6):3396.
Hadden, S. and Lithwick, Y. (2017). Kepler planet masses and eccentricities from ttv analysis. The Astronomical Journal, 154(1):5.
Herbig, G. H. (1960). The spectra of be-and ae-type stars associated with nebulosity. The Astrophysical Journal Supplement Series, 4:337.
Hogerheijde, M. R., Bekkers, D., Pinilla, P., Salinas, V. N., Kama, M., Andrews, S. M., Qi, C., and Wilner, D. J. (2016). Steepening of the 820 µm continuum surface brightness profile signals dust evolution in tw hydrae’s disk. Astronomy & Astrophysics, 586:A99.
Hu, J., De Winter, D., et al. (1989). Photometric and spectroscopic study of three candid-ate herbig ae/be stars-hd 37411, hd 100546, and hd 104237. Astronomy and Astrophysics, 208:213–218.
Huang, J., Andrews, S. M., Cleeves, L. I., ¨Oberg, K. I., Wilner, D. J., Bai, X., Birnstiel, T., Carpenter, J., Hughes, A. M., Isella, A., et al. (2018). Co and dust properties in the tw hya disk from high-resolution alma observations. The Astrophysical Journal, 852(2):122.
Joy, A. H. (1945). T tauri variable stars. The Astrophysical Journal, 102:168.
Mayor, M. and Queloz, D. (1995). A jupiter-mass companion to a solar-type star. Nature, 378(6555):355–359.
Miley, J., Pani´c, O., Haworth, T., Pascucci, I., Wyatt, M., Clarke, C., Richards, A., and Ratzka, T. (2019). Asymmetric mid-plane gas in alma images of hd 100546. Monthly Notices of the Royal Astronomical Society, 485(1):739–752.
P´erez, S., Casassus, S., Hales, A., Marino, S., Cheetham, A., Zurlo, A., Cieza, L., Dong, R., Alarc´on, F., Ben´ıtez-Llambay, P., et al. (2020). Long baseline observations of the hd 100546 protoplanetary disk with alma. The Astrophysical Journal Letters, 889(1):L24.
Pineda, J. E., Szul´agyi, J., Quanz, S. P., Van Dishoeck, E. F., Garufi, A., Meru, F., Mulders, G. D., Testi, L., Meyer, M. R., and Reggiani, M. (2019). High-resolution alma observations of hd 100546: Asymmetric circumstellar ring and circumplanetary disk upper limits. The Astrophysical Journal, 871(1):48.
Quanz, S. P., Amara, A., Meyer, M. R., Girard, J. H., Kenworthy, M. A., and Kasper, M. (2015). Confirmation and characterization of the protoplanet hd 100546 b—direct evidence for gas giant planet formation at 50 au. The Astrophysical Journal, 807(1):64.
Quanz, S. P., Amara, A., Meyer, M. R., Kenworthy, M. A., Kasper, M., and Girard, J. H. (2013). A young protoplanet candidate embedded in the circumstellar disk of hd 100546. The Astrophysical Journal Letters, 766(1):L1.
Secchi, A. (1866). Schreiben des herrn prof. secchi, directors der sternwarte des collegio romano, an den herausgeber. Astronomische Nachrichten, 68:63.
Stephens, I. W., Looney, L. W., Kwon, W., Fern´andez-L´opez, M., Hughes, A. M., Mundy, L. G., Crutcher, R. M., Li, Z.-Y., and Rao, R. (2014). Spatially resolved magnetic field structure in the disk of a t tauri star. Nature, 514(7524):597–599.
Struve, O. (1931). On the origin of bright lines in spectra of stars of class b. The Astrophysical Journal, 73:94.
Van Boekel, R., Henning, T., Menu, J., De Boer, J., Langlois, M., M¨uller, A., Avenhaus, H., Boccaletti, A., Schmid, H. M., Thalmann, C., et al. (2017). Three radial gaps in the disk of tw hydrae imaged with sphere. The Astrophysical Journal, 837(2):132.
Walsh, C., Juh´asz, A., Pinilla, P., Harsono, D., Mathews, G. S., Dent, W. R., Hogerheijde, M. R., Birnstiel, T., Meeus, G., Nomura, H., et al. (2014). Alma hints at the presence of two companions in the disk around hd 100546. The Astrophysical Journal Letters, 791(1):L6.
Wolszczan, A. and Frail, D. A. (1992). A planetary system around the millisecond pulsar psr1257+ 12. Nature, 355(6356):145–147.
