SKIRT / SKIRT9

SKIRT version 9 -- advanced radiative transfer in dusty systems
http://www.skirt.ugent.be
GNU Affero General Public License v3.0
33 stars 28 forks source link

Grain size distribution normalisation for an imported CellMedium with multiple grain populations #217

Open currodri opened 1 month ago

currodri commented 1 month ago

Description I am using SKIRT to generate mock UV-NIR spectra of a range isolated galaxy simulations include 4 dust bins (small carbonaceous grains, large carbonaceous grains, small silicate grains and large silicate grains). We are including each of the dust bins in this way in the SKI file:

`

      </CellMedium>`

Where we assume that the underlying size distribution follows the Zubko model. This results in reasonable dust properties when inspected with the DustGrainPopulationsProbe (similar to the values shown in the SKIRT tutorial here). However, when we want to specify a log-normal distribution for each dust bin like this: <LogNormalGrainSizeDistribution minSize="5.00E-4 micron" maxSize="1.00E-1 micron" centroid="5.00E-03 micron" width="7.5E-01"/> results in the DustGrainPopulationsProbe suggesting a dust per hydrogen atom larger by up to 10 orders of magnitude. This resembles the issue with normalisation raised in the tutorial previously mentioned. The question is the following – how should we obtain this normalisation? If the normalisation factor is one given for the log-normal distribution, this comes from integrating the grain size distribution by the mass of each grain and equating this to the local dust density of that population, which would render this normalisation cell-dependent. If the normalisation is done by using the dust mass per hydrogen atom/mass, is this the mass of a single dust grain with respect to the mass of the hydrogen atom? And if that's the case, do we take the grain mass for the grain size of the centroid of the log-normal distribution (assuming this is the centroid size for the assumed log-normal distribution in the simulation)?

If instead this latter normalisation is the local dust mass per total hydrogen mass, this will again be a cell-dependent normalisation, which should again be a problem...

Thank you!

mbaes commented 1 month ago

Dear Curro,

Andrea Gebek, a PhD student in our team, was dealing with a very similar setup last week. I forwarded your question to him, here is his response.

Kind regards, Maarten

In a case where multiple dust populations (with varying grain types and/or grain size distributions) are modelled in a simulation, the most straightforward way to run SKIRT on the output is by including multiple dust media (as you have presumably done). Each of these dust media then describes the grain type, the grain size distribution, and the normalization for the grain size distribution for one single dust population. Importantly, these normalizations only act relative to other dust populations in the same dust medium. Put another way, when the different dust populations are imported with separate dust media (i.e. with four different CellMedium structures in your ski file which all describe one single dust population), the type of normalization as well as the actual number for the normalization does not affect the SKIRT simulation in any way. This is because the amount of dust for each of the four dust populations are already defined by the input file(s).

In cases where one dust medium contains multiple dust populations, the normalizations for the grain size distributions are important as they determine the relative abundance of each dust population within that dust medium. Note that in this case, these relative abundances are spatially constant, which is not what you want when you have spatially resolved information on the dust properties from your simulation. This is also why the normalizations for the grain size distributions do not depend on the properties of the gas cell (e.g. the local dust density) - grain size distribution normalizations only act relative to other dust populations in the same dust medium, while the overall normalization is controlled by the amount of imported dust in a specific gas cell.

currodri commented 1 month ago

Hi Maarten,

Thanks for sharing Andrea's answer to my question. I'm not fully sure, but I believe I'm following the former option and defining 4 different CellMedium structures, just like in the file attached: frame_npov_Draine_test.txt

This is what I was using originally, but it results in too different results with just adding up all dust bins into a single density and assuming one of the built in size distributions. Using the former configuration, the optical depth convergence probe looks like this:

Spatial grid convergence information
------------------------------------

------ Dust ------

Total mass
  Input:   1.19239e+07 Msun
  Gridded: 1.19215e+07 Msun (-0.02 %)

Optical depth at 0.55 micron along full X-axis
  Input:   35.1626
  Gridded: 35.2008 (0.11 %)

Optical depth at 0.55 micron along full Y-axis
  Input:   30.6565
  Gridded: 30.6627 (0.02 %)

Optical depth at 0.55 micron along full Z-axis
  Input:   0.985419
  Gridded: 0.988615 (0.32 %)

------ Cell statistics ------

Optical depth of cell diagonal at 0.55 micron
  Largest: 0.879238
  Average: 0.13307
  90% <  : 0.32356

------------------------------------

And the probe for the grain properties of the first bin looks like this:

# Medium component 0 -- dust mass per hydrogen atom: 3.917304936e-20 kg
# column 1: grain population index (1)
# column 2: grain material type (-)
# column 3: dust mass as a percentage of total (%)
# column 4: dust mass per hydrogen atom (kg)
# column 5: dust mass per hydrogen mass (1)
# column 6: minimum grain size (micron)
# column 7: maximum grain size (micron)
0 DustEM_aSil 100.0000 3.917304936e-20 2.342014783e+07 5.000000000e-04 1.000000000e-01

When I instead add up all the dust density in all bins and use <MRNDustMix numSilicateSizes="10" numGraphiteSizes="10"/> as the materialMix and the underlying grain size distribution, the optical depth convergence probe gives:

Spatial grid convergence information
------------------------------------

------ Dust ------

Total mass
  Input:   1.19239e+07 Msun
  Gridded: 1.19215e+07 Msun (-0.02 %)

Optical depth at 0.55 micron along full X-axis
  Input:   61.8601
  Gridded: 61.9734 (0.18 %)

Optical depth at 0.55 micron along full Y-axis
  Input:   54.2587
  Gridded: 54.2397 (-0.03 %)

Optical depth at 0.55 micron along full Z-axis
  Input:   1.68741
  Gridded: 1.69393 (0.39 %)

------ Cell statistics ------

Optical depth of cell diagonal at 0.55 micron
  Largest: 1.55878
  Average: 0.254042
  90% <  : 0.545573

------------------------------------

and the grain properties look like this:

# Medium component 0 -- dust mass per hydrogen atom: 1.434721851e-29 kg
# column 1: grain population index (1)
# column 2: grain material type (-)
# column 3: dust mass as a percentage of total (%)
# column 4: dust mass per hydrogen atom (kg)
# column 5: dust mass per hydrogen mass (1)
# column 6: minimum grain size (micron)
# column 7: maximum grain size (micron)
0 Draine_Silicate  58.3750 8.375193537e-30 5.007224967e-03 5.000000000e-03 2.500000000e-01
1 Draine_Graphite  41.6250 5.972024975e-30 3.570457497e-03 5.000000000e-03 2.500000000e-01

Clearly, the dust mass per hydrogen atom shouldn't be 10 orders of magnitude larger than assuming MRN, and the optical depths shouldn't have that large inconsistency...

Maybe I'm understanding your answer incorrectly?

Thanks a lot for your time!

Curro

andreagebek commented 1 month ago

Hi Curro,

thanks a lot for performing this test comparing four different dust media to the built-in MRNDustMix! From the ConvergenceInfoProbe, I noted that the dust input/gridded masses are exactly the same in both SKIRT runs. On the other hand, the dust masses per hydrogen atom are different by almost ten orders of magnitude in the DustGrainPopulationsProbe. The reason for this is a bit technical, essentially this is due to the size distributions in SKIRT being defined as number of dust grains with size a per hydrogen atom. This might sound like a normalized grain size distribution (normalized to the number of hydrogen atoms), but in SKIRT we do not have any information on the number of hydrogen atoms in a given gas cell. Hence, the grain size distributions are effectively not normalized at this stage!

We do not have information on the number of hydrogen atoms, but we do know the mass of dust in each gas cell from the imported data. This dust mass is used to normalize the amount of dust in each gas cell by multiplying the grain size distributions (which are defined per hydrogen atom) by the appropriate number. Hence, when the DustGrainPopulationsProbe shows vastly different numbers for the dust mass per hydrogen atom that will still not affect the SKIRT run because the normalization to the total amount of dust accounts for this. Put another way, if for a given gas cell we have 1 kg of dust (known from the import file), SKIRT would simply "put 1/1.43e-29 hydrogen atoms" in this gas cell when using the MRNDustMix or 1/3.9e-20 when using the DustEM_aSil mix such that the dust mass is the same, independent of whatever dust mix is used. This is why the normalization within the dust grain population does not matter as long as you use only one population per medium (the dust mass per hydrogen atom in the DustGrainPopulationsProbe does change, but this does not have any impact on the SKIRT run).

You also found that the optical depths in the ConvergenceInfoProbe vary by a factor of approximately two. This is expected, since the dust grains (and hence the optical properties) as well as the grain size distributions are different in your custom dust mixes and in the MRNDustMix. The MRNDustMix consists of two types of grains (58.4% of the mass in DraineSilicateGrainComposition, the remaining 42.6% in DraineGraphiteGrainComposition), both modelled with a power-law grain size distribution. In your custom dust mixes you used the "aSil" and "Gra" grains from the DustEmGrainComposition class that have different optical properties, and also the size distributions differ since you use two types of log-normal grain size distributions. You can easily verify that a custom dust mix gives exactly the same SKIRT results as a built-in dust mix as long as you use exactly the same grain types and size distributions.

I hope this helps!

Cheers, Andrea

petercamps commented 1 month ago

Curro,

Just adding one more note. If you'd be more comfortable with more reasonable values in the DustGrainPopulationsProbe, you can easily perform the proper normalization (although in your case it shouldn't affect the end result).

As Andrea said, the grain size distribution defines the number of dust grains of each size in the dust population per hydrogen atom. The normalization simply scales the curve so that the integral of the grain mass over the size distribution equals some given dust mass per hydrogen atom. This is a property of the dust population and is independent of the local dust density, because the hydrogen density is assumed to scale with the dust density (at least in this definition scheme).

That said, instead of calculating the integral yourself to obtain the proper scaling factor, you can let SKIRT do the heavy lifting. The DustGrainPopulationsProbe lists the dust mass per hydrogen atom for each population. Simply take the value that you get with the Zubko size distribution and use that to normalize your own size distribution in the ski file. You can adjust the "factor on size distribution" to match the targeted value, or you can directly specify the dust mass per hydrogen atom (avoiding the extra calculation). This is explained in the section "Normalizing the size distribution" in the tutorial you mentioned.

currodri commented 1 month ago

Hi Andrea, Peter,

Thanks a lot for both your answers, that makes a lot of sense! I guess the DustGrainPopulationsProbe should not be used when one provides spatially varying dust densities from a hydrodynamical simulation?

I've tried doing precisely that, Peter, and indeed you can get to the same value of the dust to hydrogen mass in the DustGrainPopulationProbe (while showing exactly the same dust mass and optical depths at 0.55 micron).

I understand your point about the difference in optical depth due to a different grain size distribution and different grain material properties. I have tried setting up the optical properties to the DraineSilicateGrainComposition and DraineGraphiteGrainComposition, and setting the grain size distribution to a power law resembling the MRN like this:

<PowerLawGrainSizeDistribution minSize="0.005 micron" maxSize="0.25 micron" exponent="3.5"/>

resulting in optical depths within ~10% of the one using the fixed MRN. If I set again the size distribution to follow the log-normal distributions centred in the grain size of each bin (0.005 and 0.1 microns), while maintaining the same grain optical properties as the MRNDustMix, the optical depths go back to the lower values from the original test. I guess this implies that the main difference is from assuming a power-law distribution compared to the log-normal distribution.