Dust Self-Absorption `DO WHILE` Loop Can Be Replaced with Analytic Solution

[elijahmathews]

Context

In ADD_DUST, there is a DO WHILE loop (lines 175-192) that iteratively solves for the self-absorption of dust emission. The code simulates the dust absorbing its own radiation, re-emitting it, re-absorbing that re-emission, and so on, until the energy balance converges to within 1%:

        !include dust self-absorption
        !we need to iterate because the dust
        !will re-emit and be re-absorbed, etc.
        tduste = 0.0
        DO WHILE (((lboln-lbold).GT.1E-2).OR.iself.EQ.0)

           oduste = duste
           duste  = duste * diff_dust
           tduste = tduste + duste

           lbold = TSUM(nu,duste)  !after  self-abs
           lboln = TSUM(nu,oduste) !before self-abs

           duste = MAX(mduste/norm*(lboln-lbold),tiny_number)

           iself=1

        ENDDO

Mathematical Derivation

This iterative process is physically a geometric series which converges to an analytic solution.

Let:

• $L_0(\lambda)$ be the initial dust emission spectrum entering the loop.
• $S(\lambda)$ be the normalized shape of the dust emission profile (mduste/norm).
• $T(\lambda)$ be the transmission fraction of the diffuse dust (diff_dust, i.e., ).
• $f_\textrm{self}$ be the fraction of the dust's own emission shape that is re-absorbed by the dust:

$$f_\textrm{self} = \int S(\lambda)\left[1 - T(\lambda)\right]\ d\lambda$$

The total escaping flux $F_\textrm{tot}(\lambda)$ is the sum of the transmitted flux from the initial emission, plus the infinite series of transmitted "echoes" from re-absorption and re-emission.

1. Pass 0 (Initial): The initial spectrum $L_0(\lambda)$ passes through the dust.

• Escaping Flux: $F_0(\lambda) = L_0(\lambda) T(\lambda)$
• Absorbed Energy: $E_{\textrm{abs},0} = \int L_0(\lambda)(1-T(\lambda))\ d\lambda$

2. Pass 1 (First Re-emission): The absorbed energy $E_{\textrm{abs},0}$ is re-emitted with shape $S(\lambda)$.

• Re-emitted Spectrum: $L_1(\lambda) = E_{\textrm{abs},0}S(\lambda)$
• Escaping Flux: $F_1(\lambda) = L_1(\lambda)T(\lambda) = E_{\textrm{abs},0}S(\lambda)$
• Absorbed Energy: $E_{\textrm{abs},1} = \int L_1(\lambda)\left[1-T(\lambda)\right]\ d\lambda=E_{\textrm{abs},0}\int S(\lambda)\left[1-T(\lambda)\right]\ d\lambda = E_{\textrm{abs},0}f_\textrm{self}$

3. Pass (Subsequent Re-emissions): By induction, the energy absorbed in step $k$ is $E_{\textrm{abs},k}=E_{\textrm{abs},0} f_\textrm{self}^k$. The escaping flux at this step is the transmitted portion of that re-emission:

$$F_{k+1} = E_{\textrm{abs},k}S(\lambda)T(\lambda) = E_{\textrm{abs},0}f_\textrm{self}^k S(\lambda)T(\lambda)$$

4. Total Sum: The total escaping dust emission is:

$$\begin{align} F_\textrm{tot} &= F_0 + \sum_{k=0}^{\infty} F_{k+1} \\\ &= L_0(\lambda)T(\lambda) + E_{\textrm{abs},0}S(\lambda)T(\lambda) \sum_{k=0}^\infty f_\textrm{self}^k \end{align}$$

Since $0 \leq f_\textrm{self} < 1$, the geometric series $\sum f_\textrm{self}^k$ converges to $\left(1-f_\textrm{self}\right)^{-1}$.

$$F_\textrm{tot} = L_0(\lambda)T(\lambda) + \left(\frac{E_{\textrm{abs},0}}{1-f_\textrm{self}}\right)S(\lambda)T(\lambda)$$

Proposed Change

We can replace the DO WHILE loop with the exact analytic calculation. This removes the computational overhead of the loop and eliminates the 1% convergence error (1E-2) hard-coded in the current implementation.

Here is a possible Fortran implementation:

    REAL(SP) :: f_self, f_esc, l_abs_initial
    REAL(SP), DIMENSION(nspec) :: dust_shape_norm

    ! Define normalized dust shape S(lambda)
    dust_shape_norm = mduste / norm

    ! 1. Calculate f_self: The fraction of the dust shape that is self-absorbed
    ! f_self = Integral( S * (1 - T) )
    ! It is slightly more efficient to calculate the escaping fraction f_esc = Integral( S * T )
    f_esc = TSUM(nu, dust_shape_norm * diff_dust)
    f_self = 1.0 - f_esc
    
    ! Safety check to prevent division by zero if dust is infinitely opaque
    f_esc = MAX(f_esc, tiny_number)

    ! 2. Calculate the initial pass (Pass 0)
    ! Transmitted part of initial emission: L0 * T
    tduste = duste * diff_dust 
    
    ! Energy absorbed from initial emission: Integral( L0 * (1 - T) )
    ! We calculate this by subtracting transmitted luminosity from total input luminosity
    lboln = TSUM(nu, duste)       ! Total input L
    lbold = TSUM(nu, tduste)      ! Transmitted L
    l_abs_initial = lboln - lbold

    ! 3. Add the infinite series of re-emissions
    ! The total luminosity of all re-emissions is: l_abs_initial * (1 / (1 - f_self))
    ! But we only see the transmitted part of that: Total_Re_Emit * S * T
    ! Note: (1 - f_self) is exactly f_esc.
    
    tduste = tduste + (l_abs_initial / f_esc) * (dust_shape_norm * diff_dust)

This change reduces an $O(N)$ iterative process to an $O(1)$ analytic calculation.