There are two effects it is nice to consider when we work with droplet radiators instead of ordinary panel radiators: inter-absorption and inter-reflection. In the first, light is absorbed, converted to heat, and re-radiated as thermal radiation. In the second, the light just reflects directly.
Already, this is difficult, but the problem is additionally complicated by the fact that when absorption happens, the energy gets routed through the Stefan–Boltzmann Law, which introduces a fourth power of temperature into an otherwise-tractable geometric sum.
To solve this, we exploit a symmetry in the radiometric quantity of radiance: since every droplet is "average", and since radiance is distance-invariant, the incoming radiance to a given droplet must be the same as the radiance that that same droplet is emitting.
By definition, the radiance emitted (\(L_o\), "o" for "out") must be equal to the sum of emitted light (\(L_e\), "e" for "emitted") and reflected light (\(L_r\), "r" for "reflected"):
\[
L_o = L_e + L_r
\]
Meanwhile, \(L_r\) itself is just the fraction (\(1-\epsilon\)) of incoming radiance (\(L_i\), "i" for "incoming") that reflects:
\[
L_r = (1-\epsilon) L_i
\]
But now, the clever part: while our droplet might radiate into another droplet, that other droplet is also radiating back. Because every droplet is "average", both droplets have the same temperature, radiance, etc. In particular, the incoming radiance from an occluding droplet is the same as the outgoing radiance our droplet is sending back—that is, when the incoming direction is from an occluding droplet, \(L_i=L_o\). When it's not, we use the ambient radiance of space (\(L_i=L_s\), "s" for "space").
Call the fraction of occluded directions "\(f\)". In \(f\) of the directions, our droplet is occluded by another droplet emitting \(L_o\). In \((1-f)\) of the directions, we see \(L_s\). Hence, the incoming radiance to our droplet is on average:
\[
L_i = f \cdot L_o + (1-f) L_s
\]
We can substitute all these together and solve for \(L_o\):
\begin{align}
L_o &= L_e + L_r\\
&= L_e + (1-\epsilon) L_i\\
&= L_e + (1-\epsilon) ( f \cdot L_o + (1-f) L_s )\\
(1 - (1-\epsilon) f) L_o &= L_e + (1-\epsilon) (1-f) L_s\\
L_o &= \left( \frac{ L_e + (1-\epsilon) (1-f) L_s }{ 1 - (1-\epsilon) f } \right)
\end{align}
However, what we'll actually be interested in is the net radiance (\(L_n\), "n" for "net"), the difference between the incoming and outgoing radiance:
\begin{align}
L_n &= L_i - L_o\\
&= f \cdot L_o + (1-f) L_s - L_o\\
&= (1-f) (L_s - L_o)\\
&= (1-f) \left( L_s - \frac{ L_e + (1-\epsilon) (1-f) L_s }{ 1 - (1-\epsilon) f } \right)\\
&= \frac{1-f}{1-(1-\epsilon)f} ( \epsilon L_s - L_e )
\end{align}
Recall the aforementioned Stefan–Boltzmann Law (\(\epsilon\) is the absorptivity/emissivity of the material, \(\sigma_{sb}\) is a constant, \(T\) is temperature in degrees Kelvin, \(A_d\) and \(r\) are the surface area and radius of the droplet):
\begin{align}
\Phi_e &= A_d \cdot \epsilon \cdot \sigma_{sb} \cdot T^4\\
&= 4 \pi r^2 \cdot \epsilon \cdot \sigma_{sb} \cdot T^4\\
\end{align}
We also need to relate a droplet's radiant power to its radiance. (This formula applies for incoming, outgoing, or net energy.) Power is the cosine-weighted hemispherical integration of radiance, itself integrated over the droplet's surface area, a function of its radius \(r\):
\begin{align}
\Phi &= \int_{S^2(r)} \int_\Omega L \cdot \cos(\theta) \cdot d \omega_h \cdot d \omega_s\\
&= \int_{S^2(r)} \pi \cdot L \cdot d \omega_s\\
&= 4 \pi^2 r^2 L\\
L &= \frac{\Phi}{4 \pi^2 r^2}
\end{align}
Combining these two formulae with the final formula of the previous section, we can find out the net power \(\Phi_n\) of the droplet; i.e., the rate at which the droplet is losing energy. This takes into account infinite inter-reflection, absorption, and re-radiation!
\begin{align}
L_n &= \frac{1-f}{1-(1-\epsilon)f} ( \epsilon L_s - L_e )\\
\frac{\Phi_n}{4 \pi^2 r^2} &= \frac{1-f}{1-(1-\epsilon)f} \left( \epsilon L_s - \frac{\Phi_e}{4 \pi^2 r^2} \right)\\
\Phi_n &= \frac{1-f}{1-(1-\epsilon)f} \left( 4 \pi^2 r^2 \cdot \epsilon L_s - \Phi_e \right)\\
&= \frac{1-f}{1-(1-\epsilon)f} \left( 4 \pi^2 r^2 \cdot \epsilon L_s - 4 \pi r^2 \cdot \epsilon \cdot \sigma_{sb} \cdot T^4 \right)\\
&= 4 \pi r^2 \frac{(\epsilon) (1-f)}{1-(1-\epsilon)(f)} \left( \pi L_s - \sigma_{sb} T^4 \right)
\end{align}
Thermal energy ("\(J\)") is related to temperature by the specific heat capacity ("c") of the droplet's material. This is considered to be constant (usually accurate, as long as no phase change occurs). The relation is simply:
\[
J(t) = m \cdot c \cdot T(t)
\]
Note that, for a spherical droplet of radius \(r\) and density \(\rho\):
\[
m = \rho \cdot \frac{4}{3} \pi r^3
\]
Because power is the time derivative of energy, we can now combine this equation with the formula from the previous section, and integrate to get energy (or temperature) per time.
Unfortunately, the integration turns out to be horrendous due to the \(L_s\) term. Although it's possible to do in closed form, the result is bad: all logarithms and arctangents—and not even defined in important places. This then needs to be inverted for \(J(t)\). I saw no way to do this, and even various series representations (before or after the integration) are obnoxious.
Therefore, I'll give up, and I'm going to ignore the offending term entirely. It should be very small (\(r\) is small, so \(r^2\) should be tiny, and \(L_s\) should also be small in most cases). However, your improvements would be appreciated.
\begin{align}
\require{cancel}
\frac{d J(t)}{d t} = \Phi_n &=
4 \pi r^2 \frac{(\epsilon)(1-f)}{1-(1-\epsilon)(f)} \left(
\pi L_s - \sigma_{sb} T(t)^4
\right)\\
\int \left( \pi L_s - \sigma_{sb} \left(\frac{J(t)}{m c}\right)^4 \right)^{-1} d J(t)
&=
\int 4 \pi r^2 \frac{(\epsilon)(1-f)}{1-(1-\epsilon)(f)} d t\\
\int \left( \color{red}\cancel{\pi L_s}\color{black} - \sigma_{sb} \left(\frac{J(t)}{m c}\right)^4 \right)^{-1} d J(t)
&\color{red}\,\,=\color{black}
\int 4 \pi r^2 \frac{(\epsilon)(1-f)}{1-(1-\epsilon)(f)} d t\\
\frac{m^4 c^4}{3 \sigma_{sb} J(t)^3} &=
4 \pi r^2 \frac{(\epsilon)(1-f)}{1-(1-\epsilon)(f)} \cdot t + D
\end{align}
Solving the initial value problem for the integration constant \(D\) is trivial using \(t=0\) (note: clearly, we need to know \(J(0)\) or \(T(0)\), the initial energy or temperature of each droplet). Rearranging the result and solving for \(J(t)\):
\begin{align}
\frac{m^4 c^4}{3 \sigma_{sb} J(t)^3} &=
4 \pi r^2 \frac{(\epsilon)(1-f)}{1-(1-\epsilon)(f)} \cdot t +
\frac{m^4 c^4}{3 \sigma_{sb} J(0)^3}\\
\frac{1}{J(t)^3} &=
\frac{12 \pi r^2}{m^4 c^4} \cdot \frac{(\sigma_{sb})(\epsilon)(1-f)}{1-(1-\epsilon)f} \cdot t + \frac{1}{J(0)^3}\\
J(t) &= 1 \bigg/ \sqrt[3]{
\frac{12 \pi r^2}{m^4 c^4} \cdot \frac{(\sigma_{sb})(\epsilon)(1-f)}{1-(1-\epsilon)f} \cdot t + \frac{1}{J(0)^3}
}
\end{align}
One can use the \(J(t) = m \cdot c \cdot T(t)\) relation to get \(T(t)\).
If the average separation between droplet centers is \(q\), and they fly through a region \(A\), then the total number of droplets (from a formula by Eric Rozier quoted here) is:
\[
\text{Number of droplets} = \frac{A}{4 r^2 + 4 r q + q^2}
\]
Since the total energy radiated by a single droplet over one pass taking time \(\Delta t\) is \(J(0)-J(\Delta t)\), the total energy radiated by all droplets during that same \(\Delta t\) is just the product of the droplet energy decrease and the number of droplets. (If this isn't obvious, try thinking about a single droplet in one streamline. Its neighbor droplets don't fly for the whole \(\Delta t\), but droplets that will be emitted during \(\Delta t\) will exactly fill in the part they didn't emit for.) From this, we can calculate the radiator's total power:
\begin{align}
\Phi_{radiator} = \frac{J_{radiator}}{\Delta t} &= \left(\frac{
\left(\frac{A}{4 r^2 + 4 r q + q^2}\right) (J(0) - J(\Delta t))
}{\Delta t}\right)\\
&= \frac{A}{
(4 r^2 + 4 r q + q^2) \Delta t
} \left(
m c T(0) - 1 \bigg/ \sqrt[3]{
\frac{12 \pi r^2}{m^4 c^4} \cdot \frac{(\sigma_{sb})(\epsilon)(1-f)}{1-(1-\epsilon)f} \cdot t + \frac{1}{(m c T(0))^3}
}~
\right)
\end{align}
Some rearranging and real-analysis shows that \(f\) and \(\epsilon\) have a surprisingly small effect on the power, while \(\Delta t\) has the most effect.
One can easily differentiate \(J(t)\) to get the per-droplet power as a function of time: \(\Phi_n(t)\):
\[
\Phi_n(t) =
\frac{-4 \pi r^2}{m^4 c^4} \cdot \frac{(\sigma_{sb})(\epsilon)(1-f)}{1-(1-\epsilon)f} \cdot
\left(
\frac{12 \pi r^2}{m^4 c^4} \cdot \frac{(\sigma_{sb})(\epsilon)(1-f)}{1-(1-\epsilon)f} \cdot t + \frac{1}{J(0)^3}
\right)^{-4/3}
\]
The efficiency of the radiator relative to the case when there is no occlusion can be calculated at \(t=0\) as:
\[
\text{Efficiency} =
\frac{\Phi_{n, f>0}(0)}{\Phi_{n,f=0}(0)} =
\frac{1-f}{1-(1-\epsilon)f}
\]
Note: the original, less-complete and less-correct version of this analysis was posted here.