FaRaDiT, The Fargo to Radmc3d Disk Thickening python script

Script overview

Current version of the script

Sources for bc thesis

Pinte et al 2020
Pinte et al 2019

Currently working on…

Adaptive resolution

Create such θ grid so it has few cells in the optically thick and the super optically thin region, but many on the boundary of the two
So far only root grid was implememted
Oct tree (for the future?)

Perturbations

Creating images of the disk with various defects

A ring of higher/lower density

Non-trivial temperature perturbation also required — so far I'm just using a basic gaussian — semi-done ✓

Imaging

proper T slice plot for radmc results
Images at multiple λ for given ncells, nphot, nphotᵢₘ and τ — done ✓
Spectra of few disks

Make

Temperature

Create an additional heating source (vicose) to add to radmc results so midplane temp agrees with fargo (oh boy is this annoying)

I don't understand how to add temperature correctly

creating a procedure to generate gastemper from fargo's temperature outputs

I don't understand how to do this correctly — some papers? termo proff?

loading the radmc3d gastemper file — the result of mctherm — done ✓
analyzing its profile — ~~semi done ✓, missing T = f(z) as opposed to T = f(θ)~~ now also done ✓
creating an interpolation procedure in case of insufficient number of heating photons — done ✓

Results (newest to oldest)

Fourier transform

New approach to opacity distribution

I'm gonna try to observe where the sillicates should evaporate and change the metalicity there. First step suggested by my supervisor is to plot the temp of cells with their density. Final plot is here. Notice the strange structures. Idk what's up with that. For this simulation I used r×θ×φ = (100×40×50)

Attemps at viscous heating

As said by the radmc3d manual (pages 28, 120 and 164), the program 'can' work with spacially distributed heat sources, which is rather useful for viscous heating. It is said multiple times that this doesn't yet work though. I hoped that if I'm brave enough, it will work nonetheless. Heatsource method works by associating power density with each of the cells in the model in units erg·cm⁻³·s¯¹. I tried adding the heatsource to the equatorial plane (as a first approx of viscous heating). Using any power output below ~ 1e-30 shows no real difference above Monte Carlo noise; any higher and radmc throws a "Temp too high" error; the same one as when the disk accidentally reaches the stars surface. I'll perform more tests, but I'm afraid Dullemond knew what he was talking about when he wrote 'NOT YET WORKING'. :(

Perturbations and better grid

I implemented a slightly better θ grid, but I still put too many cells in the midplane region. Like this the cells are distributed as a root distribution. I also generalized the slice plotting procedures for both density and temperature as well as grid plotting. Here is the result

After that I implemented a function that gives gaussian blob at a given r and ran the simulation. The results clearly show an effect of shadowing:

More imaging

I also imaged the disk at multiple λs to see how the disk behaves. I'm using the fargo disks for these ones. All of these are heated by 4·10⁷ photons and imaged by 3·10⁵ photons. Soon I'll also provide the optical depths through the plane of the disk, which are perhaps the most important numbers for such images.

Note the different scales on the colorbars on the left. I chose different "saturate" parameters for each image. This parameter, designed by Attila Juhasz (as the rest of radmc3dPy code) determines the "Highest pixel values to be plotted in terms of the peak value, higher pixel values will be clipped". Not sure how scientific it is, seems to just change the max value for the colorbar, not really the pixels themselves...
I chose these λs to try and see the different behaviour for the opacities of sillicates as in the dustkappa_silicate.inp file (plotted below). One can clearly see the coldness of the midplane as it only shines in far microwave and radio spectrum.

The spectrum of this disk is rather strange...

Temperature shenanigans reloaded

Now finaly I succeeded at implementing a proper vertical temperature plot, here are the results:

Note the smoothness of the curves. I used a 2D cubic interpolation. The linear interpolation gives such result:

The last possibility is the „nearest“ interpolation, which might be more usefull for scatter plots:

Also the temperature plotting function now looks at one of the middle θ indeces (there's always two for no cell can be at exactly the middle) of the 3D T array and just plots out the values at fixed φ. Here are the temperatures as they come out of fargo (i.e. those used to make the disk thick) and the ramdc output (still heated by 4·10⁷ photons, using 10×20×6 r×θ×φ cells)

I plotted out the temperature as calculated by radmc3d of the disk at fixed r and φ as a function of θ. For comparison, in the graph on the right you can see the profila of temperature as an output form fargo. The midplane temperature (i.e. the global minimum of the first graph) is much lower than the expected ~160 K. On the other hand the "constant" maximum temperature is similar, yet not the same as the fargo value. As r increases it gets wrong too though.

Temperature shenanigans

I implemented a function that looks at each ring of temperature at specified r and θ and looks for cells that are colder than cutoff times average of its φ neighbours. Upon finding these, their temperature is replaced by said average. Here are three examples of unaltered, poorly heated disk, one treated with cutoff = 0.7 and one with cutoff = 0.99

Computational load

Below are tables describing the end results of multiple radmc3d simulations of a protoplanetary disk. The parameters of these simulations are

the number of monte carlo simulation photons or heating photons nphot

the number of imaging photons nphotᵢₘ

depth of the disk along a line from the star through the plane of the disk τ

the number of sectors in r, θ and φ directions ncells

the wavelength of the image λ

With these params the simulation is ran always with a 6 core multithread parameter on and the time of the calculation is recorded as well as the end product. Some of the density files of the disks were too large for my RAM so the simulation halted. Radmc doesn't warn of this, so possibly some of the images might be duplicates.

Here are the images of disks with the lowest optical depth (τ≈ 6·10¹² kg/m²):

	nphot = 10⁵	nphot = 10⁶	nphot = 10⁷	nphot = 10⁸
ncells = 10
ncells = 50
ncells = 100

Contrast these with the optically thick (τ≈ 7·10¹³ kg/m²):

	nphot = 10⁵	nphot = 10⁶	nphot = 10⁷	nphot = 10⁸
ncells = 10
ncells = 50
ncells = 100