m |
m |
||
Line 9: | Line 9: | ||
The formula is obtained by integration over the joint probability density function $\phi_\Gamma$ for emitting the photons at times $t_i$, which we show in this paper is given by: | The formula is obtained by integration over the joint probability density function $\phi_\Gamma$ for emitting the photons at times $t_i$, which we show in this paper is given by: | ||
$$\phi_\Gamma(t_1,\dots,t_N)=(-1)^{N}N!\gamma_a^N\left(\frac{\Gamma}{\Gamma_-}\right)^{2N }\prod_{i=1}^{N}\left[{d\over dt_i}({e^{-\gamma_a t_i}\over\gamma_a}+{e^{-\Gamma t_i}\over\Gamma}-4{e^{-\Gamma_+ t_i/2}\over\Gamma_+})\right]{\mathbf{1}}_{[t_{i-1}, t_{i+1}[}(t_{i})$$ | $$\phi_\Gamma(t_1,\dots,t_N)=(-1)^{N}N!\gamma_a^N\left(\frac{\Gamma}{\Gamma_-}\right)^{2N }\prod_{i=1}^{N}\left[{d\over dt_i}({e^{-\gamma_a t_i}\over\gamma_a}+{e^{-\Gamma t_i}\over\Gamma}-4{e^{-\Gamma_+ t_i/2}\over\Gamma_+})\right]{\mathbf{1}}_{[t_{i-1}, t_{i+1}[}(t_{i})$$ | ||
− | where $\Gamma_\pm\equiv\Gamma\pm\gamma_a$. When tracing out all the photons to keep only the $k$th of interest, through some clever math hackery, we arrive to the computation of $\int_0^\infty\left({e^{-\gamma_a t}\over\gamma_a}+{e^{-\Gamma t}\over\Gamma}-4{e^{-\Gamma_+ t/2}\over\Gamma_+}\right)^N\,dt$ which is easy to do with the multinomial theorem, which is where all the mad coefficients come from. The nice thing, of course, is that it covers all cases with a single formula, i.e., it can tell you about the 1st photon from a biphoton or the 7th photon from a 15-photon bundle. Now the direct calculation of these particular cases is a bit tedious. We can ask a formal computer language, like Mathematica, to do it. | + | where $\Gamma_\pm\equiv\Gamma\pm\gamma_a$. When tracing out all the photons to keep only the $k$th of interest, through some clever math hackery, we arrive to the computation of $\int_0^\infty\left({e^{-\gamma_a t}\over\gamma_a}+{e^{-\Gamma t}\over\Gamma}-4{e^{-\Gamma_+ t/2}\over\Gamma_+}\right)^N\,dt$ which is easy to do with the multinomial theorem, which is where all the mad coefficients come from. The nice thing, of course, is that it covers all cases with a single formula, i.e., it can tell you about the 1st photon from a biphoton or the 7th photon from a 15-photon bundle. Now the direct calculation of these particular cases is a bit tedious. We can ask a formal computer language, like Mathematica, to do it. This is the generic expression to be summed over: |
<syntaxhighlight lang="python"> | <syntaxhighlight lang="python"> | ||
Line 21: | Line 21: | ||
2) + (\[CapitalGamma] + \[Gamma]a)/2 (k3 + k9))^2) | 2) + (\[CapitalGamma] + \[Gamma]a)/2 (k3 + k9))^2) | ||
</syntaxhighlight> | </syntaxhighlight> | ||
+ | |||
+ | and defining the following module to provide all the ways to write an integer $s$ as the sum of exactly $k$ terms (possibly zero): | ||
+ | |||
+ | <syntaxhighlight lang="python"> | ||
+ | PartitionForMultinomialSum[n_, k_] := Module[{}, | ||
+ | Flatten[ | ||
+ | Permutations /@ (PadRight[#, k] & /@ | ||
+ | Select[IntegerPartitions[n], Length[#] <= k &]), 1]] | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | for instance, | ||
+ | |||
+ | <syntaxhighlight lang="python"> | ||
+ | In:= PartitionForMultinomialSum[3, 2] | ||
+ | |||
+ | Out= {{3, 0}, {0, 3}, {2, 1}, {1, 2}} | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | we can now compute, e.g., the average time of detection for the <tt>k=2</tt>nd photon from a <tt>N=3</tt>-photon bundle as: | ||
+ | |||
+ | <syntaxhighlight lang="python"> | ||
+ | In:= k = 2; \[CapitalNu] = 3; | ||
+ | In:= listkN = Flatten[Table[Flatten[{lk1, lk2, lk3}], | ||
+ | {lk1, PartitionForMultinomialSum[\[CapitalNu] - k, 3]}, | ||
+ | {lk2, PartitionForMultinomialSum[k - 1, 6]}, | ||
+ | {lk3, PartitionForMultinomialSum[2, 2]} | ||
+ | ], 3 - 1] | ||
+ | Out= {{1, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0}, {1, 0, 0, 1, 0, 0, 0, 0, 0, 0, | ||
+ | 2}, {1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0, 0, 0, | ||
+ | 2, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2}, {1, 0, 0, 0, 1, 0, 0, 0, | ||
+ | 0, 1, 1}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0}, {1, 0, 0, 0, 0, 1, 0, | ||
+ | 0, 0, 0, 2}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 0, | ||
+ | 1, 0, 0, 2, 0}, {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2}, {1, 0, 0, 0, 0, | ||
+ | 0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0}, {1, 0, 0, 0, | ||
+ | 0, 0, 0, 1, 0, 0, 2}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1}, {1, 0, 0, | ||
+ | 0, 0, 0, 0, 0, 1, 2, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2}, {1, 0, | ||
+ | 0, 0, 0, 0, 0, 0, 1, 1, 1}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 2, 0}, {0, | ||
+ | 1, 0, 1, 0, 0, 0, 0, 0, 0, 2}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 1, | ||
+ | 1}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 2, 0}, {0, 1, 0, 0, 1, 0, 0, 0, 0, | ||
+ | 0, 2}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1}, {0, 1, 0, 0, 0, 1, 0, 0, | ||
+ | 0, 2, 0}, {0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 2}, {0, 1, 0, 0, 0, 1, 0, | ||
+ | 0, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 2, 0}, {0, 1, 0, 0, 0, 0, | ||
+ | 1, 0, 0, 0, 2}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1}, {0, 1, 0, 0, 0, | ||
+ | 0, 0, 1, 0, 2, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 2}, {0, 1, 0, 0, | ||
+ | 0, 0, 0, 1, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 2, 0}, {0, 1, 0, | ||
+ | 0, 0, 0, 0, 0, 1, 0, 2}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1}, {0, 0, | ||
+ | 1, 1, 0, 0, 0, 0, 0, 2, 0}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2}, {0, | ||
+ | 0, 1, 1, 0, 0, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 2, | ||
+ | 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 2}, {0, 0, 1, 0, 1, 0, 0, 0, 0, | ||
+ | 1, 1}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 2, 0}, {0, 0, 1, 0, 0, 1, 0, 0, | ||
+ | 0, 0, 2}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 1, | ||
+ | 0, 0, 2, 0}, {0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2}, {0, 0, 1, 0, 0, 0, | ||
+ | 1, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 2, 0}, {0, 0, 1, 0, 0, | ||
+ | 0, 0, 1, 0, 0, 2}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1}, {0, 0, 1, 0, | ||
+ | 0, 0, 0, 0, 1, 2, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2}, {0, 0, 1, | ||
+ | 0, 0, 0, 0, 0, 1, 1, 1}} | ||
+ | |||
+ | In=Simplify[Total[ | ||
+ | tkN[#[[1]], #[[2]], #[[3]], #[[4]], #[[5]], #[[6]], #[[7]], #[[ | ||
+ | 8]], #[[9]], #[[10]], #[[11]]] & /@ listkN]] | ||
+ | Out[200]= (150 \[CapitalGamma]^8 + 2345 \[CapitalGamma]^7 \[Gamma]a + | ||
+ | 14493 \[CapitalGamma]^6 \[Gamma]a^2 + | ||
+ | 41371 \[CapitalGamma]^5 \[Gamma]a^3 + | ||
+ | 58786 \[CapitalGamma]^4 \[Gamma]a^4 + | ||
+ | 41371 \[CapitalGamma]^3 \[Gamma]a^5 + | ||
+ | 14493 \[CapitalGamma]^2 \[Gamma]a^6 + | ||
+ | 2345 \[CapitalGamma] \[Gamma]a^7 + | ||
+ | 150 \[Gamma]a^8)/(6 \[CapitalGamma] \[Gamma]a (\[CapitalGamma] + \ | ||
+ | \[Gamma]a) (2 \[CapitalGamma] + \[Gamma]a) (3 \[CapitalGamma] + \ | ||
+ | \[Gamma]a) (5 \[CapitalGamma] + \[Gamma]a) (\[CapitalGamma] + | ||
+ | 2 \[Gamma]a) (\[CapitalGamma] + 3 \[Gamma]a) (\[CapitalGamma] + | ||
+ | 5 \[Gamma]a)) | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | $$\langle t_2^{(3)} \rangle = | ||
+ | \frac{\Big(150 \gamma_a^8 + 2345 \gamma_a^7 \Gamma +14493 \gamma_a^6 \Gamma^2 + 41371 \gamma_a^5 \Gamma^3 + 58786 \gamma_a^4 \Gamma^4 + 41371 \gamma_a^3 \Gamma^5 + 14493 \gamma_a^2 \Gamma^6 +2345 \gamma_a \Gamma^7 + 150 \Gamma^8\Big)}{6 \Gamma \gamma_a (\Gamma +\gamma_a) (2 \Gamma + \gamma_a) (\Gamma + 2 \gamma_a)(3 \Gamma + \gamma_a)(\Gamma + 3 \gamma_a) (5 \Gamma +\gamma_a) (\Gamma + 5\gamma_a)} | ||
+ | $$ | ||
+ | |||
+ | Maybe not particularly enlightening, but that is the result. This can also be written as a sum of inverse combinations of the two emission rates: | ||
+ | |||
+ | $$\langle t_1^{(1)} \rangle =\frac{2}{\Gamma +\gamma_a}+\frac{1}{\Gamma }+\frac{1}{\gamma_a}$$ | ||
+ | |||
+ | $$\langle t_1^{(2)} \rangle =\frac{1}{8} \left(\frac{8}{\Gamma +\gamma_a}+\frac{9}{3 \Gamma +\gamma_a}+\frac{9}{\Gamma +3 \gamma_a}+\frac{4}{\Gamma }+\frac{4}{\gamma_a}\right)$$ | ||
+ | |||
+ | $$\langle t_2^{(2)} \rangle =\frac{3}{8} \left(\frac{8}{\Gamma +\gamma_a}-\frac{3}{3 \Gamma +\gamma_a}-\frac{3}{\Gamma +3 \gamma_a}+\frac{4}{\Gamma }+\frac{4}{\gamma_a}\right)$$ | ||
+ | |||
+ | $$\langle t_1^{(3)} \rangle =\frac{1}{243} \left(2 \left(\frac{62}{2 \Gamma +\gamma_a}+\frac{125}{5 \Gamma +\gamma_a}+\frac{62}{\Gamma +2\gamma_a}+\frac{125}{\Gamma +5 \gamma_a}+\frac{81}{\Gamma +\gamma_a}\right)+\frac{81}{\Gamma}+\frac{81}{\gamma_a}\right)$$ | ||
+ | |||
+ | $$\langle t_2^{(3)} \rangle =\frac{1}{1944}\Bigg(\frac{3240}{\Gamma +\gamma_a}-\frac{1984}{2 \Gamma +\gamma_a}+\frac{6561}{3 \Gamma +\gamma_a}-\frac{4000}{5 \Gamma +\gamma_a}-\frac{1984}{\Gamma +2 \gamma_a}+\frac{6561}{\Gamma +3 \gamma_a}+{}\\{}-\frac{4000}{\Gamma +5 \gamma_a}+\frac{1620}{\Gamma }+\frac{1620}{\gamma_a}\Bigg)$$ | ||
+ | |||
+ | \begin{multline} | ||
+ | \langle t_3^{(3)} \rangle =\frac{1}{1944}\Bigg(\frac{7128}{\Gamma +\gamma_a}+\frac{992}{2 \Gamma +\gamma_a}-\frac{6561}{3 \Gamma +\gamma_a}+\frac{2000}{5 \Gamma +\gamma_a}+\frac{992}{\Gamma +2 \gamma_a}+{}\\{}-\frac{6561}{\Gamma +3 \gamma_a}+\frac{2000}{\Gamma +5 \gamma_a} | ||
+ | \frac{3564}{\Gamma }+\frac{3564}{\gamma_a}\Bigg) | ||
+ | \end{multline} | ||
This is the formula for the average of the squared time of the $k$th photon from a $N$-photon bundle: | This is the formula for the average of the squared time of the $k$th photon from a $N$-photon bundle: |
Our last paper with Guillermo Díaz Camacho et al. on the dynamics of emission of multiphotons, in particular, of photon bundles (?!) involves a lot of complicated combinatoric formulas, such as the average time of detection for the $k$th photon from a $N$-photon Fock state that is spontaneously emitted with a radiative rate $\gamma_a$ and measured by a detector with bandwidth $\Gamma$: \begin{multline} \langle t_{k}^{(N)} \rangle = 2 N! \left(\frac{\Gamma_+ \Gamma \gamma_a}{\Gamma_-^2}\right)^N \sum_{k_1 + k_2 + k_3\atop = N-k}\sum_{k_4 + \cdots+ k_9\atop = k-1} \sum_{k_{10} + k_{11}\atop = 2}\prod_{j=1}^{11}{1\over k_j!}\times\\ \frac{(-1)^{k_3 + k_6 + k_7 + k_8 + k_{11}}}{\gamma_a^{k_1 +k_4+k_7} \Gamma^{k_2+k_5+k_8} \left( \frac{\Gamma_+}{4} \right)^{k_3 + k_6 + k_9} \big( \gamma_a (k_1 + k_7+{k_{11}\over2}) +\Gamma( k_2 +k_8 +{k_{10}\over2}) + \frac{\Gamma_+}{2} (k_3+k_9)\big)^2 }\,. \end{multline}
Here I provide a Mathematica code to turn this general result into particular and limiting cases, and discuss briefly the whole thing. The formula is obtained by integration over the joint probability density function $\phi_\Gamma$ for emitting the photons at times $t_i$, which we show in this paper is given by: $$\phi_\Gamma(t_1,\dots,t_N)=(-1)^{N}N!\gamma_a^N\left(\frac{\Gamma}{\Gamma_-}\right)^{2N }\prod_{i=1}^{N}\left[{d\over dt_i}({e^{-\gamma_a t_i}\over\gamma_a}+{e^{-\Gamma t_i}\over\Gamma}-4{e^{-\Gamma_+ t_i/2}\over\Gamma_+})\right]{\mathbf{1}}_{[t_{i-1}, t_{i+1}[}(t_{i})$$ where $\Gamma_\pm\equiv\Gamma\pm\gamma_a$. When tracing out all the photons to keep only the $k$th of interest, through some clever math hackery, we arrive to the computation of $\int_0^\infty\left({e^{-\gamma_a t}\over\gamma_a}+{e^{-\Gamma t}\over\Gamma}-4{e^{-\Gamma_+ t/2}\over\Gamma_+}\right)^N\,dt$ which is easy to do with the multinomial theorem, which is where all the mad coefficients come from. The nice thing, of course, is that it covers all cases with a single formula, i.e., it can tell you about the 1st photon from a biphoton or the 7th photon from a 15-photon bundle. Now the direct calculation of these particular cases is a bit tedious. We can ask a formal computer language, like Mathematica, to do it. This is the generic expression to be summed over:
tkN[k1_, k2_, k3_, k4_, k5_, k6_, k7_, k8_, k9_, k10_, k11_] := 2 \[CapitalNu]! (((\[CapitalGamma] + \[Gamma]a) \[CapitalGamma] \ \[Gamma]a)/(\[CapitalGamma] - \[Gamma]a)^2)^\[CapitalNu] 1/( k1! k2! k3! k4! k5! k6! k7! k8! k9! k10! k11!) (-1)^( k3 + k6 + k7 + k8 + k11)/(\[Gamma]a^(k1 + k4 + k7) \[CapitalGamma]^( k2 + k5 + k8) ((\[CapitalGamma] + \[Gamma]a)/4)^( k3 + k6 + k9) (\[Gamma]a (k1 + k7 + k11/2) + \[CapitalGamma] (k2 + k8 + k10/ 2) + (\[CapitalGamma] + \[Gamma]a)/2 (k3 + k9))^2)
and defining the following module to provide all the ways to write an integer $s$ as the sum of exactly $k$ terms (possibly zero):
PartitionForMultinomialSum[n_, k_] := Module[{}, Flatten[ Permutations /@ (PadRight[#, k] & /@ Select[IntegerPartitions[n], Length[#] <= k &]), 1]]
for instance,
In:= PartitionForMultinomialSum[3, 2] Out= {{3, 0}, {0, 3}, {2, 1}, {1, 2}}
we can now compute, e.g., the average time of detection for the k=2nd photon from a N=3-photon bundle as:
In:= k = 2; \[CapitalNu] = 3; In:= listkN = Flatten[Table[Flatten[{lk1, lk2, lk3}], {lk1, PartitionForMultinomialSum[\[CapitalNu] - k, 3]}, {lk2, PartitionForMultinomialSum[k - 1, 6]}, {lk3, PartitionForMultinomialSum[2, 2]} ], 3 - 1] Out= {{1, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0}, {1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2}, {1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2}, {1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0}, {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2}, {1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1}, {1, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2}, {1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 2, 0}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 2}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 2, 0}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 2}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1}, {0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0}, {0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 2}, {0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 2, 0}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 2}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 2, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 2}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 2, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 2}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 0}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 2, 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 2}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 2, 0}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 2}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 1, 0, 0, 2, 0}, {0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2}, {0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 2, 0}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 2}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 0, 0, 1, 2, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2}, {0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1}} In=Simplify[Total[ tkN[#[[1]], #[[2]], #[[3]], #[[4]], #[[5]], #[[6]], #[[7]], #[[ 8]], #[[9]], #[[10]], #[[11]]] & /@ listkN]] Out[200]= (150 \[CapitalGamma]^8 + 2345 \[CapitalGamma]^7 \[Gamma]a + 14493 \[CapitalGamma]^6 \[Gamma]a^2 + 41371 \[CapitalGamma]^5 \[Gamma]a^3 + 58786 \[CapitalGamma]^4 \[Gamma]a^4 + 41371 \[CapitalGamma]^3 \[Gamma]a^5 + 14493 \[CapitalGamma]^2 \[Gamma]a^6 + 2345 \[CapitalGamma] \[Gamma]a^7 + 150 \[Gamma]a^8)/(6 \[CapitalGamma] \[Gamma]a (\[CapitalGamma] + \ \[Gamma]a) (2 \[CapitalGamma] + \[Gamma]a) (3 \[CapitalGamma] + \ \[Gamma]a) (5 \[CapitalGamma] + \[Gamma]a) (\[CapitalGamma] + 2 \[Gamma]a) (\[CapitalGamma] + 3 \[Gamma]a) (\[CapitalGamma] + 5 \[Gamma]a))
$$\langle t_2^{(3)} \rangle = \frac{\Big(150 \gamma_a^8 + 2345 \gamma_a^7 \Gamma +14493 \gamma_a^6 \Gamma^2 + 41371 \gamma_a^5 \Gamma^3 + 58786 \gamma_a^4 \Gamma^4 + 41371 \gamma_a^3 \Gamma^5 + 14493 \gamma_a^2 \Gamma^6 +2345 \gamma_a \Gamma^7 + 150 \Gamma^8\Big)}{6 \Gamma \gamma_a (\Gamma +\gamma_a) (2 \Gamma + \gamma_a) (\Gamma + 2 \gamma_a)(3 \Gamma + \gamma_a)(\Gamma + 3 \gamma_a) (5 \Gamma +\gamma_a) (\Gamma + 5\gamma_a)} $$
Maybe not particularly enlightening, but that is the result. This can also be written as a sum of inverse combinations of the two emission rates:
$$\langle t_1^{(1)} \rangle =\frac{2}{\Gamma +\gamma_a}+\frac{1}{\Gamma }+\frac{1}{\gamma_a}$$
$$\langle t_1^{(2)} \rangle =\frac{1}{8} \left(\frac{8}{\Gamma +\gamma_a}+\frac{9}{3 \Gamma +\gamma_a}+\frac{9}{\Gamma +3 \gamma_a}+\frac{4}{\Gamma }+\frac{4}{\gamma_a}\right)$$
$$\langle t_2^{(2)} \rangle =\frac{3}{8} \left(\frac{8}{\Gamma +\gamma_a}-\frac{3}{3 \Gamma +\gamma_a}-\frac{3}{\Gamma +3 \gamma_a}+\frac{4}{\Gamma }+\frac{4}{\gamma_a}\right)$$
$$\langle t_1^{(3)} \rangle =\frac{1}{243} \left(2 \left(\frac{62}{2 \Gamma +\gamma_a}+\frac{125}{5 \Gamma +\gamma_a}+\frac{62}{\Gamma +2\gamma_a}+\frac{125}{\Gamma +5 \gamma_a}+\frac{81}{\Gamma +\gamma_a}\right)+\frac{81}{\Gamma}+\frac{81}{\gamma_a}\right)$$
$$\langle t_2^{(3)} \rangle =\frac{1}{1944}\Bigg(\frac{3240}{\Gamma +\gamma_a}-\frac{1984}{2 \Gamma +\gamma_a}+\frac{6561}{3 \Gamma +\gamma_a}-\frac{4000}{5 \Gamma +\gamma_a}-\frac{1984}{\Gamma +2 \gamma_a}+\frac{6561}{\Gamma +3 \gamma_a}+{}\\{}-\frac{4000}{\Gamma +5 \gamma_a}+\frac{1620}{\Gamma }+\frac{1620}{\gamma_a}\Bigg)$$
\begin{multline} \langle t_3^{(3)} \rangle =\frac{1}{1944}\Bigg(\frac{7128}{\Gamma +\gamma_a}+\frac{992}{2 \Gamma +\gamma_a}-\frac{6561}{3 \Gamma +\gamma_a}+\frac{2000}{5 \Gamma +\gamma_a}+\frac{992}{\Gamma +2 \gamma_a}+{}\\{}-\frac{6561}{\Gamma +3 \gamma_a}+\frac{2000}{\Gamma +5 \gamma_a} \frac{3564}{\Gamma }+\frac{3564}{\gamma_a}\Bigg) \end{multline}
This is the formula for the average of the squared time of the $k$th photon from a $N$-photon bundle:
\begin{multline} \langle \big(t_k^{(N)}\big)^2\rangle = \frac{2\gamma_a^N}{P(N,N)}\left( \frac{\Gamma}{\Gamma_-} \right)^{2N}k {N\choose k} \sum_{k_1+k_2+k_3=N-k}{N-k\choose k_1, k_2, k_3}\\{}\times \sum_{k_4+\cdots+k_9=k-1}{k-1\choose k_4, k_5, k_6, k_7,k_8,k_9} \sum_{k_{10}+k_{11}=2}{2\choose k_{10}, k_{11}}\\{}\times \frac{(-1)^{k_3 + k_6 + k_7 + k_{8}+k_{11}}}{\gamma_a^{k_1+k_4+k_7} \Gamma^{k_2+k_5+k_8} \left( \frac{\Gamma_+}{4} \right)^{k_3 + k_6 +k_9} ( \gamma_a (k_1+k_7+k_{11}/2) +\Gamma(k_2+k_8 +k_{10}/2) + \frac{\Gamma_+}{2} (k_3+k_9))^3} \end{multline}