Advice on using ProDiMo's continuum Radiative Transfer (RT)¶

The continuum radiative transfer in 2D circumstellar disks, including the determination of the internal dust temperature structure, is a challenging numerical problem. Previous benchmark papers (Pascucci et al. 2004, A&A 417, 793) and (Pinte et al. 2009, A&A 498, 967) have shown that modern Monte Carlo techniques, coupled to diffusion-solvers, invoked to the most optically thick parts of the disks, provide the most powerful numerical approach to date, both concerning accuracy and speed, whereas "classical" deterministic approaches, using pre-fixed sets of rays, have difficulties to resolve the complicated angular dependencies, e.g. at sharp edges, and to converge properly in optically thick situations.

ProDiMo nevertheless uses such an approach, calculating new mean intensities Jnu from solving the radiative transfer equation with isotropic scattering along Ntheta x Nphi ~ 200 rays per grid point, then employing a Lambda-iteration technique to evaluate new dust temperatures, then accelerating convergence with the NG-algorthm (see Woitke et al. 2009, A&A 501, 383), with the original idea to provide a realistic UV- and background IR radiation field at any point in the disk for the gas modelling. The continuum results obtained with this technique, however, are inferior to those from highly optimised MC-programs such as MCFOST and MCMAX. Therefore,

TODO: This should not be an advice anymore, we don't even know if this is still working (at best with an very old MCFOST version)

 --> *we advice to first run MCFOST and then transfer the results Jnu(r,z) and Tdust(r,z) to ProDiMo, whenever possible.* <--

If you do want to use ProDiMo anyway for the continuum radiative transfer, in particular for the analysis and interpretation of SEDs, please note the current limitations of ProDiMo

1. isotropic scattering
2. limited angular resolution
3. problems with temperature-convergence in very optically thick cases.

In particular, there is a "trick" implemented to help ProDiMo's RT-convergence in optically thick situations: If the radial and the vertical optical depth > tau_cutoff, the dust opacity is not increased any further in the downward direction from the disc surface to the midplane. We do this for every frequency bin separately. Without this trick, ProDiMo's RT will simply not converge in optically thick cases, i.e. all temperatures in the optically thick core of the disc will stay as they are initialised. Using this trick, however, introduces some uncertainties, especially in the midplane. These deviations rarely have a considerable impact on the calculated SEDs, but under certain circumstances, they may.

Since revision 12f830a5, there is an alternative approach, which I would recommend to use, namely to call a 2D grey diffusion solver to solve

div Flux = div -4pi/(3 kappaRoss) grad J = Heatrate           (1)

in the optically thick region where the "trick" would have reduced the opacities (which is then not done). "Heatrate" is the non-RE dust heating rate [1/cm3/s] which for standard models is =0. We solve this linear system of 2nd order equations based on the frequency-integrated mean intensity J values found just outside of this region (2D boundary value problem). Numerically this leads to N linear equations for N unknowns, namely all J(ix,iz) in the flagged region.

This page offers some advice on how to best choose RT parameters dependent on your problem. The parameters are

tau_cutoff      : threshold optical depth for the "trick" / diffusion solver   (default=10, recommended 10 )
RTitmax         : maximum number of RT-temperature-iterations                  (default=50, recommended 200)
Ntheta          : number of angular grid-points for theta                      (default=17, recommended 19 )
Nphi            : number of angular grid-points for phi                        (default=9,  recommended 13 )
optimise_theta  : adjust theta angles to direction of heating                  (default=T,  recommended T  )
optimise_phi    : adjust phi angles to direction of heating                    (default=T,  recommended T  )
resolve_star    : use 1+Nphi rays to better resolve the star                   (default=T,  recommended T  )
DiffusionApprox : use the 2D diffusion solver                                  (default=T,  recommended T  )

optimise_theta and optimise_phi with adjust the setting of the theta and phi directions (for every point!), during the first few RT iterations, toward the direction where the heating photons are coming from. resolve_star uses 1+Nphi rays to better resolve the star as seen from every point, otherwise the star is treated as point source.

The precision of the results will furthermore depend on the number of frequency bins and the number and position of the grid points, in particular we must make sure that the first few radial layers of the inner rim and subsequent walls are truly optically thin. How to achieve that depends on the disk mass, dust opacity, disk shape parameters, and on the position of the inner rim. For the standard TTauri model, for studying the RT, I would recommend

100         ! NXX
70          ! NZZ
40          ! Ninner
1.5         ! refinement_fac
30          ! NLAM
.true.      ! refine_Spitzer

If you are rather interested in the gas modelling in the inner rim, however, e.g. CO ro-vib, a better choice might be

90          ! NXX
70          ! NZZ
20          ! Ninner
30          ! NLAM
.true.      ! refine_Spitzer

as this setting (without refinement_fac) creates fewer but even finer grid points at the wall to resolve the line optical depths.

We have thoroughly tested the results of ProDiMo's RT against MCFOST and MCMax, using MCcomp.pro in the idl folder. Here are some results for the standard TTauri model which converged after about 200 RT iterations:

As a bonus, the new diffusion solver also copes very well with viscous heating, for example

1.0E-8      ! Mdot [Msol/yr]  : accretion rate
.true.      ! dust_nonRE      : dust not in radiative equilibrium
.false.     ! VHEAT_GAS       : put heat into the dust

What types of artefacts/deviations may occur?¶

1) Due to the "trick", Tdust can be too high in the midplane. If that happens even in the outer parts of the disk, the SED around 100 micron is too high, in this example by about 20%. In the pictures, black is ProDiMo, red is MCFOST (isotropic scattering) and blue is MCFOST (anisotropic scattering). From a model for AB_Aur with tau_cutoff=10.

2) Probably due to the limited angular resolution, the inner rim may be too cold. Consequently, the SED at near-IR wavelengths can be too low, in this example by about 15%. From a tau=10⁴-benchmark model (Pinte et al.2009) with Ntheta=25 and Nphi=15. Please note that the other, quite obvious, temperature deviations visble in the above l.h.s. plot, namely the too warm temperatures in the close midplane, usually have not much of an effect on the SED, because at the short wavelengths this region is situated too deep (optically thick), and for the long wavelengths, the region is not extended enough.

Which types of disks are likely to produce such deviations?¶

Type 1) Very massive disks with vertical A_V>10 even at 200 AU, for example massive disks with a large inner whole, or massive disks with surface density powerlaw index as small as epsilon=0.5.

Type 2) Massive TTauri disks with small Rin, and sharp inner edge.

How to spot a deviation?¶

Type 1) If you see large vertical AV at several 100 AU, in combination with almost isothermal midplane below, your model is likely to suffer from type-1-deviations.

Type 2) I don't know.

What to do?¶

Type 1) Try larger value for tau_cutoff. In the above case, a massive disk model for AB Aur was unexpectedly easy to converge even with tau_cutoff=200. The deviations are gone.

Type 2) Try larger values for NXX, Ntheta and Nphi, but there is no guarantee. I presently do not understand where the temperature-differences to MCFOST at the inner rim come from.

Final remarks¶

ProDiMo creates an output file called RTconv.out, which can be visualised with convergeRT.pro, showing the course of the RT-iterations. Currently, the RT-iterations are stopped if J-change<3.E-3 AND T-change<3.E-4. If your convergence is as nice as in the above case (the cured AB_Aur model depicted above), there is no need to use the "trick" and you should switch it off by using a large value for tau_cutoff, e.g. 200. If you, however, see insufficient convergence you need to use lower value of tau_cutoff, realizing that your results might be affected. You can increase the maximum number of RT-iterations carried out by changing parameter RTitmax.

Here is the log of the standard TTauri model without restart with diffusion solver

This is how it looks when it's converging, but please note that it needed 177 RT-iterations (the plot also counts the NG-iterations) to do so.

In general, if you use ProDiMo for SED analysis, please take some time to investigate these questions. Only looking at the convergence plots is not enough, though, because a small J-change does not necessarily imply that the model has converged. It might just imply that it converges very slowly. Try different RT-parameters and compare the results! If you use automated fitting, it's a good idea to always restart from the last model, because in this way ProDiMo will "add" more RT-iterations with every generation.

More informations about ProDiMo's performance in the RT disk benchmark test (Pinte et al. 2009) can be found here RTbenchmark.