Breakdown of Bose-Einstein Distribution in Photonic Crystals

In the last two decades, considerable advances have been made in the investigation of nano-photonics in photonic crystals. Previous theoretical investigations of photon dynamics were carried out at zero temperature. Here, we investigate micro/nano cavity photonics in photonic crystals at finite temperature. Due to photonic-band-gap-induced localized long-lived photon dynamics, we discover that cavity photons in photonic crystals do not obey Bose-Einstein statistical distribution. Within the photonic band gap and in the vicinity of the band edge, cavity photons combine the long-lived non-Markovain dynamics with thermal fluctuations together to form photon states that memorize the initial cavity state information. As a result, Bose-Einstein distribution is completely broken down in these regimes, even if the thermal energy is larger or much larger than the cavity detuning energy. In this investigation, a crossover phenomenon from equilibrium to nonequilibrium steady states is also revealed.

in Methods. Different spectral density J(v) of the PCs, which requires the knowledge of DOS of the PCs, will provide different photon dissipation and fluctuations. In standard quantum optics, the spectral density J(v) is usually treated as a constant in the weak-coupling limit, where the Weisskopf-Winger approximation or the Markovian master equation is valid so that the photon damping rate is time-independent (i.e. memoryless). As a result, all cavity photon modes have a finite lifetime and the cavity photon states will ultimately evolve into thermal equilibrium with the environment after decoherence takes place 19 , and photons inside the cavity must obey Bose-Einstein distribution [20][21][22] . However, this well-known result is no longer satisfied for micro/nano cavities in PCs, due to the presence of the PBG, as we will show below.
In principle, the DOS for different PCs, denoted by % PC v ð Þ, should be calculated by solving photon eigenfrequencies and eigenfunctions from the Maxwell's equations for different photonic crystal structures 2,23,24 . In the literature, different photonic DOS have been introduced for different PCs [9][10][11][12]14,[25][26][27][28][29] . For examples, for 1D PCs, the corresponding DOS is given by where H(v 2 v e ) is the Heaviside step function and v e is the frequency at the PBE. This DOS is also often used for isotropic 3D PCs, which could predict qualitatively correct behaviors of the non-Weisskopf-Winger decay and the photon-atom bound state in PCs [9][10][11][12]14 , but it may overestimate the spontaneous emission rate of atoms due to the absence of the singularity of DOS in reality. Thus, for 3D PCs, the DOS near the PBE may be modeled by an anisotropic DOS [9][10][11] : % PC v ð Þ! ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi v{v e p H v{v e ð Þ, which could be further modified when the vectorial property of EM field is taken into account 25,26 . For 2D PCs with van Hove saddle point singularity [27][28][29] , the photonic DOS exhibits a logarithmic divergence near the PBE and can be approximated by where v 0 is the center of the logarithmic peak. The spectral density J(v) is microscopically defined as the multiplication of the DOS of the PCs with the photon tunneling amplitude V(v) between the cavity and the PCs, Assume that the cavity mode equally couples to all possible modes near the PBE of PCs, i.e. treating V(v) as a constant, then the corresponding spectral densities J(v) for different dimensional PCs are fully determined by the corresponding DOS, which are summarized in Table 1 and are also plotted in Fig. 1(a). We take the unit h~1 hereafter so that both the frequency v and the spectral density J(v) have the same unit as the energy. Consider the cavity frequency v c lies not too far away from the PBE. The dissipative photon dynamics is described by the cavity field propagating Green function u(t, t 0 ) in photonic crystals through the relation AEa(t)ae 5 u(t, t 0 ) AEa(t 0 )ae, and is determined by the dissipative Table 1 | Characters of different photonic crystal structures. For different DOS of different dimensional PCs, the corresponding different spectral densities J(v) are listed and the reservoir-induced photon self-energy corrections S(v) are calculated, which completely determine the dissipative photon dynamics in PCs. The parameters C, g and x are coupling strengths (with a proper rescaling) between the cavity and PCs for 1D, 2D and 3D PCs, respectively. A smooth high-frequency cutoff V C is introduced for avoiding the divergence of DOS in 3D PCs, and a sharp high-frequency cutoff at V d is set up for maintaining the positivity of DOS in 2D PCs. Here Li 2 (x) is a dilogarithm function and erfc (x) is a complementary error function  Table 1 for different DOS of 1D, 2D and 3D PCs (with different colors) are plotted respectively in the vicinity of photonic band edge v e ; (b) The corresponding localized photon mode frequency v b as a function of the detuning d 5 v c 2 v e ; and (c) The corresponding localized photon mode amplitudes, given in Eq. (2). The localized photon mode shows a crossover for 1D and 2D PCs, and a critical transition for 3D PCs when the cavity frequency v c changes from the PBG region to the PB region. The parameters given in Table 1 take as follows: the coupling strengths for 1D, 2D and 3D PCs are C 2/3 5 0.01v e , g 5 0.001v e and x 5 0.014v e , respectively. Experimentally, the typical photonic band edge frequency v e ranges from a few GHz to a few tens GHz for most of 1D, 2D and 3D PCs 30,31 , then the corresponding coupling strengths used in our calculations are ranged from a few MHz to a few hundreds MHz. The cutoff parameters V d < 3.87v e and V C 5 0.1v e ; and the logarithmic divergence center v 0 5 1.04v e . integro-differential equation (19). The general solution of Eq. (19) is given in Ref. 15, is the reservoir-induced cavity photon self-energy correction, The explicit solution of cavity photon self-energy corrections for different spectral densities is also presented in Table 1. The frequency v b in Eq. (2) is the localized photon mode frequency located inside the PBG (0 , v b , v e ) but not too far away from the PBE, and it is determined by the pole condition ! is a principal-value integral.
Equation (2) provides a general solution of the non-Markovian dissipative photon dynamics for micro/nano cavities in various photonic crystal structures. It shows that the cavity photon dynamics in PCs always contains two parts 15 : a localized photon mode (the first term) plus a non-exponential photon damping (the second term). The localized photon mode is a long-lived non-Markovian effect (dissipationless), induced by the PBG structure in PCs. The corresponding frequency v b lies always within the PBG. The non-exponential photon damping comes from the non-analyticity of the photon self-energy correction, which is determined by the DOS profile of the PCs. This non-exponential damping is a short-time non-Markovian memory effect, and it will become an exponential (Markov) decay in the photonic band (PB) region, as we will show later. The contributions of both the localized photon mode and the non-exponential photon damping strongly rely on the detuning d 5 v c 2 v e .
In Fig. 1b, we plot the localized photon mode frequency as a function of the detuning d for 1D, 2D and 3D PCs, respectively. The values of the localized mode frequencies for three different dimensional PCs are indeed very close to each other. However, a carefully check (see the inset in Fig. 1b) shows that for 3D PCs with an anisotropic DOS, the localized photon mode exists only when dv ffiffiffi p p x, namely, the cavity frequency v c must be tuned into the PBG or near the PBE. Here x is a coupling strength between the cavity and the PCs (see Table 1). For 1D and 2D PCs, the localized photon mode exists for any location of the cavity frequency v c . This is because for 1D and 2D PCs with the DOS given in Table 1, the self-energy corrections are negatively divergent when v R v e 2 0 1 so that the pole condition is always satisfied for the localized photon mode. However, this does not ensure that the localized photon mode in 1D and 2D PCs must have a significant contribution to the photon dynamics for the whole range of the cavity frequency. The importance of the localized photon mode is determined by the localized photon mode amplitude, 1= 1{S' v b ð Þ ½ given in Eq. (2), which is plotted in Fig. 1c as a function of the detuning d for all the three different dimensional PCs. It shows that for 3D PCs with an anisotropic DOS, the localized photon mode amplitude vanishes for dw ffiffiffi p p x, indicating a critical transition for the occurrence of the localized mode. For 1D and 2D PCs, the localized photon mode amplitude decays to almost zero when v c is tuned into the PB region, as a crossover phenomenon. However, the overall effect of these localized photon modes on cavity photon dynamics has the very similar behavior for different PCs given in Table 1, as shown in Fig. 1.
The results presented in Fig. 1b-c provide indeed the full steadystate information of the cavity photon dynamics. This is because the non-exponential photon damping [i.e. the second term in Eq. (2)] will decay to zero so that only the localized photon mode contributes to the steady-state cavity photon field. In other words, the steadystate cavity field amplitude quantifies the contribution of the localized photon mode as a dissipationless effect, given by On the other hand, the non-exponential damping term in Eq. (2) characterizes the time-dependent cavity photon damping rate k(t) in the intermediate dissipation time region through the relation 15,16 The detailed photon dynamics in terms of the cavity field amplitude ju(t, t 0 )j is provided in Fig. 2(a) for 1D, 2D and 3D PCs, respectively, with several different detuning d taking from the PBG region to the PB region. The corresponding damping rate k(t) is given in Fig. 2b, except for the case with d 5 0.1v e which will be discussed later. The results show that the cavity dynamics changes dramatically when v c crosses over from the PBG region to the PB region. Since the range of u(t, t 0 ) is given by 1 $ ju(t, t 0 )j $ 0, we define the crossover region with the condition 0:9 * w u t??,t 0 ð Þ j j * w 0:1 which corresponds to {0:025v e * v d * v 0:025v e , as shown in show that the damping rate k(t) will within the PBG (d , 20.025v e ) and in the vicinity of the PBE ({0:025v e * v d * v 0:025v e ). In other words, the photon damping Dissipation and dissiptionless photon dynamics in terms of the cavity field, AEa(t)ae 5 u(t, t 0 ) AEa(t 0 )ae, and (b) the corresponding decay rate k(t), are plotted for (i) 1D PCs; (ii) 2D PCs and (iii) 3D PCs, with several different detuning d. It shows how the dissipative cavity photons becomes dissipationless when the cavity frequency moves from the PB region into the PBG region for different dimensional PCs. For d 5 0.1v e , the photon dynamics exhibits almost a perfect exponential decay, and by fitting the damping with an exponential function, the resulting damping rate is k . rapidly approach to zero rate will approach to zero after some time, due to the localized photon dynamics. In the PB region (d . 0.025v e where either the localized mode vanishes (for 3D PCs) or becomes negligible (for 1D and 2D PCs), then the cavity photons undergo a full dissipation process, and can be approximately characterized by an exponential damping for d?0:025v e , as shown by the black curves in Fig. 2a for d 5 0.1v e . In this case, the damping rate k(t) will oscillate rapidly in time when u(t, t 0 ) decays to zero, see Eq. (5). This rapidly oscillating damping rate, originated from u(t, t 0 ) approaching zero, has no physical consequence because the photon dissipation is almost completed after this point of time. Thus, alternatively we determine the damping rate k by fitting the dashed-dot black curves in Fig. 2a with an exponential function ,exp(2kt). The resulting cavity photon damping rate for d 5 0.1v e is given by k^3:16 {3 v e , 1.24 3 10 22 v e , and 1.72 3 10 22 v e for 1D, 2D and 3D PCs, respectively. Experimentally, the typical photonic band edge frequency v e ranges from a few GHz to a few tens GHz for most of 1D, 2D and 3D PCs fabricated in the microwave regime 30,31 , this corresponds to the damping rate k ranging from a few tens MHz to hundreds MHz when the cavity frequency v c lies inside the PB region.
Combining the localized photon mode together with the photon dissipation dynamics, as shown in Fig. 1 and Fig. 2, we can see that when the cavity mode is tuned away from the PBE with d . 0.025v e , the localized mode has a negligible effect on cavity photon dynamics for 1D and 2D PCs, and vanishes completely for 3D PCs. The cavity photon field decays almost exponentially for d?0,025v e , which is a Markov decay, as shown in Fig. 2a by the dashed-dot black curves with d 5 0.1v e . In the crossover region {0:025v e * v d * v 0:025v e , i.e. the cavity frequency is tuned into the vicinity of the PBE, the contribution of the localized mode to the photon dynamics increases rapidly, while the contribution of the nonexponential photon damping goes in an opposite way, namely the damping rate k(t) is decreased to zero after some time (see Fig. 2b). When the cavity frequency lies inside the PBG (d , 20.025v e ), the cavity field has almost no damping (the damping rate k(t) quickly decays to zero), and the photon dynamics is dominated by the localized photon mode. Thus light can be confined in the defect of the photonic crystal, providing a high-Q micro/nano cavity. In fact, this nontrivial dissipative cavity photon dynamics reproduces the same result with regard to atomic population trapping and atom-photon bound states in the vicinity of the PBG where an atom is placed in photonic crystals [9][10][11][12] .
The importance of the above results is that the localized photon mode provides dissipationless cavity photon dynamics when the cavity frequency v c lies inside the PBG or in the vicinity of the PBE. This dissipationless localized-mode contribution to the steady-state cavity photon field is universal for different dimensional PCs with different DOS given in Table 1, because it only relies on the presence of PBG. The localized photon mode decreases rather quickly but smoothly in the crossover region ({0:025v e * v d * v 0:025v e ), and becomes negligible for d . 0.025v e for 1D and 2D PCs, where the cavity photon dynamics becomes dissipative. For the 3D PCs, the cavity photon dynamics shows the same dissipationless behavior before reaching the critical point d cf ffiffi p p x^0:025v e , where the localized photon mode dominates the photon dynamics inside the PBG (d , 20.025v e ) and decays rapidly near the PBE ({0:025v e * v d * v 0:025v e ). When d . d c , the localized photon mode vanishes for 3D PCs, and the cavity photon dynamics becomes fully dissipative. With these results, a general picture for dissipation and dissipationless photon dynamics in photonic crystals is provided.
Thermal photon fluctuations. The above exact solution of dissipation and dissipationless cavity photon dynamics in photonic crystals can be used to describe thermal photon fluctuations through the fluctuated photon correlation function v(t, t), which is fully determined by the generalized non-equilibrium fluctuationdissipation theorem 15 In Eq. (6), the two-time correlation functiong t 1 ,t 2 ð Þð is the initial photon distribution in PCs at temperature T. When the cavity approaches the steady state, photon fluctuations are simply determined by the modified steady-state fluctuation-dissipation theorem, v t,t?? where The first term is the localized photon mode contribution that modifies the conventional equilibrium fluctuation-dissipation theorem. This additional contribution is negligible when the cavity frequency is tuned into Fig. 1b-c). Consequently, the solution (7) is reduced to the standard equilibrium fluctuation-dissipation theorem: In the high temperature limit, it reproduces Einstein's fluctuationdissipation relation.
Photon fluctuations are presented in Fig. 3a(i)-(iii) respectively for 1D, 2D and 3D PCs at a given temperature. It shows that photon fluctuations evolve in a similar way for different PCs. When v c is tuned into the PB region (d . 0.025v e ), photons continuously flow into the cavity from the photonic crystal until the cavity reaches its steady state. In this case, the cavity thermally equilibrates with the photonic crystal. When v c lies inside the PBE region ({0:025v e * v d * v 0:025v e ), photons also flow into the cavity at the beginning, but some of photons are then transmitted back into the photonic crystal, due to the environment-induced memory effects and the effect induced by the localized photon mode. When v c is tuned into the PBG region (d , 20.025v e ), photons hardly flow from the photonic crystal into the cavity, because the localized photon mode dominates the cavity photon dynamics in the PBG region where the reservoir-induced dissipation is negligible, even though the thermal energy, k B T~ hv e , is much larger than the detuning energy. Thus, thermal fluctuations are suppressed significantly by the photon localization in the PBG region and also near the PBE.
The physical picture of the above photon fluctuations can be seen clearly by connecting the photon correlations with the measurable average cavity photon number (i.e. field intensity) 16,18,32 : n(t) 5 AEa { (t)a(t)ae 5 ju(t, t 0 )j 2 n(t 0 ) 1 v(t, t), where n(t 0 ) is the initial cavity intensity. It shows that the fluctuated photon correlation function v(t, t) is just the thermal-fluctuation-induced average photon number in the cavity. In Fig. 3b, we present the steady-state photon fluctuations at a given temperature, see the solid-blue curve. The dashed-pink curve is the thermal photon distribution in PCs at the same temperature: n v c ,T ð Þ~1 e hvc=kBT {1 Â Ã . It shows that when d . 0.025v e , the thermal-fluctuation-induced steady-state average photon number in the cavity is identical to the thermal photon distribution: n t??
ð Þ~v t,t?? ð Þ n v c ,T ð Þ. In this regime, the initial cavity photons are totally lost into PCs as u(t, t 0 ) R 0 in the steady-state limit, as a consequence of the Weisskopf-Winger decay.  Fig. 3b. In other words, Bose-Einstein distribution no longer works for cavity photons in the PBE region, due to the localized photon effect. When the cavity frequency v c lies inside the PBG (d , 20.025v e ), the thermal fluctuations approach to zero, and thus the cavity is fully determined by the localized photons, even though the cavity mode equally couples to all possible modes in PCs and the thermal energy is much larger than the detuning, k B T?d. Consequently, Bose-Einstein distribution is broken down. This provides a nontrivial phenomenon of the thermal fluctuation dynamics from equilibrium to nonequilibrium steady states. Here nonequilibrium steady states are defined by these steady states which can memorize the initial state information of the system so that the equilibrium hypothesis cannot be satisfied. Such steady states will be given explicitly in the next subsections.
It is recently shown that non-Markovian dynamics in open systems can induce a critical transition from equilibrium to nonequilibrium steady states, when a localized mode occurs 15,22,33 . Such a critical transition appears for the 3D PCs with an anisotropic DOS, where there is a critical detuning for the occurrence of localized photon mode, d c~ffi ffiffi p p x. However, for the 1D and 2D PCs, our result shows that there is no such a critical transition due to the fact that the localized photon mode always exists. The photon fluctuations then undergoes a crossover from equilibrium to nonequilibrium steady states through the change of the detuning d. No matter it is a cross-over or a critical transition, the results in Fig. 3b shows that thermal photon fluctuations cannot make the cavities embedded in PCs to approach equilibrium with the PCs when the cavity frequency lies inside the PBG or near the PBE. This is a general feature for cavity photons in PCs.
Breakdown of Bose-Einstein distribution through the timeevolution of photonic Fock states. Based on the exact cavity photon dissipation and fluctuation dynamics in PCs, we shall explicitly examine the breakdown of Bose-Einstein distribution through the evolution of cavity photonic states, and also the crossover from equilibrium to nonequilibrium steady states. This must be done by solving the exact master equation (15) given in the Methods. To be more specific, we consider first the cavity to be initially in a Fock state with an arbitrary photon number n 0 , i.e. r(t 0 ) 5 jn 0 ae AEn 0 j, which may be prepared experimentally through the realtime quantum feedback control 34 . By solving the master equation (15), the cavity state at arbitrary time t is given by where ð Þ . This result shows that an initial photon Fock state will evolve into a mixed state of different Fock states, and the probability in each Fock state jnae is P n0 ð Þ n t ð Þ. As we have shown, the breakdown of Bose-Einstein distribution relies on the dissipationless photon dynamics. And different dimensional PCs show almost the same dissipationless photon dynamics, determined by the localized photon mode due to the existence of the PBG. To be specific, we will present in the following the numerical solution of the state evolution for 1D PCs. The 2D and 3D PCs must provide a similar solution, based on the universal dissipation and fluctuation photon dynamics provided by Eqs. (2) and (6), and also the explicit numerical results presented in Fig. 2 and Fig. 3. The timeevolution of the cavity photon distribution P n0 ð Þ n t ð Þ for the initial state jn 0 5 5ae is given in Fig. 4. The steady-state limit, P n0 ð Þ n t?? ð Þ, is shown in Fig. 5, where several different initial states have been considered to demonstrate the initial-state dependence of the steady photon states.
The results show that if we tune v c into the PB region (e.g. d 5 0.1v e ), then u(t, t 0 ) will gradually decay to zero. This indicates that photons in the cavity will gradually be damped into the photonic crystal, and photons in the photonic crystal are transferred into the cavity through thermal fluctuations. In this case, the contribution of the localized photon mode is negligible. The cavity state ultimately reaches thermal equilibrium with the photonic crystal, and Bose-Einstein statistical distribution is produced, where the initial state information is completely washed out, and the cavity photon steady state is independent of the initial states, as shown in Fig. 4a and Fig. 5a. Near the PBE (e.g. d 5 0), the cavity still loses photons into the photonic crystal and gains photons from the photonic crystal through thermal fluctuations. At a low temperature (k B T0 :2 hv e ), the photon loss and photon gain make the cavity become  a mixed state of several Fock states jnae only for n , n 0 , and mainly distributed around n 5 n 0 /2, see Fig. 4b-(i) and also the steady-state limit by Fig. 5b-(i) where the initial state dependence is manifested. The photon distribution deviates significantly from the standard Bose-Einstein distribution. At a relatively high temperature (k B T~ hv e ), the cavity state also becomes a mixed state of several Fock states jnae, distributed mainly among n # n 0 /2, but the distribution is broader, see Fig. 4b-(ii) and also the steady-state limit in Fig. 5b-(ii). In this case, the initial-state dependence of the steady state is still shown up and the derivation of the photon distribution from the standard Bose-Einstein distribution is obvious. As the temperature becomes rather high (k B T~10 hv e ), the time-evolution of the cavity state behaves quite differently because thermal fluctuations play a more important role now, thus the structure of the initial state is quickly destroyed, as shown in Fig. 4b-(iii). However, the cavity does not really approach thermal equilibrium with the photonic ð Þ is small but not negligible so that the distribution given by Eq. (9) is still different from the thermal state distribution. The cavity steady states for different initial states do not differ from each other as much as these in the low temperature cases, but the initial-state dependence is still observable, as shown in Fig. 5b-(iii). When the cavity frequency v c lies inside the PBG, e.g. d 5 20.1v e given in Fig. 4c and Fig. 5c, the steady-state cavity photon distribution strongly depends on the initial state. In particular, if the photonic crystal temperature is not too high (k B T~0:2 hv e and k B T~ hv e ), the cavity remains in its initial Fock state as a photon localization state. It only has a small chance to decay to the Fock state jnae with n , n 0 [see Fig. 4c(i) with k B T~0:2 hv e ], and an even smaller probability to be in the Fock state jnae with n . n 0 [see Fig. 4c(ii) when k B T increases to hv e ]. The steady-state cavity photons are distributed mainly in the regime n , n 0 with the maximum peak at n 5 n 0 , see Fig. 5c(i)-(ii). However, at a higher temperature k B T~10 hv e , the thermal fluctuation becomes strong such that some of the initial state information will be lost and the cavity evolves into a mixed state covering several Fock states around the initial one jn 0 ae [see Fig. 4c(iii)]. The steady-state cavity photons are distributed more broadly, but are still centered around the initial photon number n 0 , as shown in Fig. 5c(iii). This is because the localized photon mode dominates the photon dynamics over thermal photon fluctuations. The overall initial-state dependence of the steady states shown in Fig. 5b-c indicates that the equilibrium hypothesis in statistical mechanics is no longer obeyed and Bose-Einstein statistical distribution is completely broken down, even though the initial thermal energy of the photonic crystals is larger or much larger than the detuning energy.
Breakdown of Bose-Einstein distribution through the timeevolution of coherent states. Because Fock sates are highly nonclassical and may not be easy to prepare in experiments, here we examine the case the cavity is initially in a coherent state ja 0 ae. By solving the master equation (15), the cavity state at an arbitrary later time t is given by where is a displacement operator with a(t) 5 u(t, t 0 )a 0 , and is a thermal-like state with average particle number v(t, t). Equation (11) shows that the initial cavity state will evolve into a displaced thermal-like state 35,36 which is the mixture of displaced number states D a t ð Þ ½ n j i 37 . In the photon number representation, Eq. (11) can be written as where the probability of the cavity containing n photons (i.e. cavity photon distribution) is given by the diagonal terms of Eq. (13). The cavity photon distribution in the steady-state limit is presented in Fig. 6. It shows again that when the cavity frequency v c lies inside the PB region (d 5 0.1v e ), the cavity state (11) will evolve into the thermal state r T n v c ,T ð Þ ½ , which indicates that the cavity  www.nature.com/scientificreports SCIENTIFIC REPORTS | 5 : 9423 | DOI: 10.1038/srep09423 photon number obeys Bose-Einstein distribution, see Fig. 6a. If the cavity mode is tuned at the PBE (d 5 0), the cavity steady state will depend on its initial state, as shown in Fig. 6b. For low temperature (e.g. k B T~0:2 hv e ), the thermal fluctuation is suppressed (v(t, t) R 0), and the cavity will evolve into the coherent state ja(t)ae, see Eq. (11). In this case, the cavity photon distribution in the steady state is Poissonian-like, with average photon number AEnae j t R ' 5 ja(t)j 2 / ja 0 j 2 , see Fig. 6b-(i). As the temperature goes higher (k B T~ hv e ), the cavity photon distribution is deviated away from the Poisson distribution due to thermal fluctuations, see Fig. 6b-(ii). At high temperature (k B T~10 hv e ), since thermal fluctuation plays a more important role, the initial-state dependence is not as significant as that in the low temperature, and the cavity photon distribution is not Poissonian-like, but still deviates from Bose-Einstein distribution, see Fig. 6b-(iii). When the cavity frequency lies inside the PBG (d 5 20.1v e ), the cavity steady state strongly depends on its initial state. The photon distribution in the cavity remains almost the same as in its initial state when the temperature is not very high (k B T~0:2 hv e and k B T~ hv e ), see the first two graphs in Fig. 6c. The initial-state dependence of the photon distribution is still strong when the temperature gets higher (k B T~10 hv e ), see Fig. 6c-(iii). The photon distribution in this case is slightly broadened from the initial states due to thermal fluctuations at high temperature, but still centered around n < ja 0 j 2 , as a manifestation of photon localization. These results obtained for different initial coherent states at different temperature behave indeed in the same way as the results obtained for different initial Fock states under the same temperature.

Discussion
The breakdown of Bose-Einstein distribution in PCs, explored in details through the time-evolution of various cavity photon states for 1D PCs, is also valid for 2D and 3D PCs with different DOS listed in Table 1. This is because the time evolution of the cavity photon states, solved analytically by Eqs. (8) and (11), is fully determined by the photon dissipation and fluctuation dynamics of Eqs. (2) and (6). The steady-state behaviors of Eqs. (2) and (6) for different PCs are almost the same, as shown in Fig. 2 and Fig. 3. There is only a slight difference on the photon dissipation and fluctuation dynamics between 1D and 2D PCs with 3D PCs. For 1D and 2D PCs, varying the detuning leads to a crossover from equilibrium to nonequilibrium steady states for the photon dynamics, while for 3D PCs with an anisotropic DOS, it gives a critical transition, rather than a crossover, due to the existence of a critical detuning for the existence of the localized mode, see explicitly the results presented in Fig. 1 to Fig. 3. Such a critical transition is also found recently in other open systems 22,33 but not for the crossover phenomenon. Thus, photon dissipation and fluctuation dynamics investigated in this work reveal a new nontrivial property, i.e. a crossover phenomenon for photon dynamics in low-dimensional photonic crystals.
In conclusion, we show that when the cavity frequency lies inside the PBG or near the PBE in photonic crystals, Bose-Einstein statistical distribution is broken down for cavity photons. This conclusion is generally valid for various photonic band gap structures in PCs. For the 1D and 2D PCs, the breakdown of Bose-Einstein distribution leads to a crossover from equilibrium to nonequilibrium cavity steady states, while for 3D PCs with an anisotropic DOS, the breakdown of Bose-Einstein distribution corresponds to a critical transition rather than a crossover. No matter whether it is a crossover or a critical transition, the breakdown of Bose-Einstein distribution is a consequence of localization photons due to the presence of PBG structures in PCs. Therefore the conclusion is also valid for other nanomaterials with band gap structures. It could provide a hitherto unexplored challenge on photon statistics. Furthermore, this nontrivial photon dynamics can be examined via quantum non-demolition measurement 38,39 , by sending sequences of circular Rydberg atoms through the photonic crystal micro/nano cavity, which carry information without destroying the cavity photon state. In particular, the cavity photon state can be measured using an experimental setup similar to that given in Ref. 38. Such experiments could be done with microcavities in low-dimensional photonic crystals in the microwave regime.

Methods
To investigate photon dynamics of micro/nano cavities which are coupled each other through waveguides embedded in PCs at finite temperature, we treat both the PCs and waveguides as reservoirs of the cavities. Thus the entire system of the micro/nano cavities (defects) embedded in photonic crystals can be described by the Fano-Anderson model (a model of impurity electrons coupled with continuous states introduced by Anderson 40 in solid-state physics, and discrete states embedded in a continuum proposed by Fano 41 in atomic spectra). The corresponding Fano-Anderson Hamiltonian is given by Refs.16, 40-43: where a i a { Consider initially the photonic crystals are in equilibrium state. Integrating out completely the reservoir degrees of freedom of photonic crystals via the influence functional 44,45 in the coherent state representation 46 , arbitrary photon states for micro/nano cavities in PCs are then governed by the following exact master equation [15][16][17][18]32,43,47,48 Here r (t) is the reduced density matrix for cavity states, H' c t ð Þ~X i,j v' cij t ð Þa { i a j is the renormalized Hamiltonian of cavities with the environment-modified cavity frequencies v' cii t ð Þ:v' ci t ð Þ and the environment-induced couplings between different cavities v' cij t ð Þ, after the environmental degrees of freedom are completely integrated out. The coefficients k ij (t) andk ij t ð Þ characterize photon dissipations and fluctuations in PCs at finite temperature. The renormalized frequency, v' cij , and the time-dependent dissipation and fluctuation coefficients, k ij (t) andk ij t ð Þ, can be exactly and non-perturbatively determined by the dissipation-fluctuation Dyson equation, as shown explicitly in our early works 16,17 , and will also be given later. www.nature.com/scientificreports It might be worth pointing out that the first master equation derived from the original Feynman-Vernon influence functional 44 was obtained by Caldeira and Leggett for the Brownian motion with a high temperature environment 45 which is also called the Caldeira-Leggett master equation in the literature. For the Brownian motion, the system-reservoir coupling Hamiltonian is given by H I~X k l k xq k , where x and q k are the positions of the principal harmonic oscillator (as the Brownian particle) and all other harmonic oscillators in the reservoir, respectively. In terms of the second quantization, this system-reservoir coupling Hamiltonian can be rewritten as which shows that the amplitudes for particle tunneling processes between the system and the reservoir (the first two terms) is equal to the amplitudes of particle pair production and annihilation processes (the last two terms), both are given by V k . However, in quantum optics, it is well-known that the photon tunneling amplitudes cannot be the same as the amplitudes of two-photon pair production and annihilation processes. As a result, the Caldeira-Leggett model is not applicable to photonics. To make the cavity photon dynamics in thermal photonic crystals more specific, we consider a single-mode micro/nano cavity in photonic crystals. The photon dissipation and fluctuations, characterized by the dissipation and fluctuation coefficients k (t) andk t ð Þ (all the sub-indices (i, j) in Eq. (15) can be dropped now), are determined non-perturbatively and exactly by the nonequilibrium Green functions 49,50 through the relations [15][16][17][18] : k t ð Þ~_ v t,t ð Þz2v t,t ð Þk t ð Þ: ð18Þ Here u (t, t 0 ) is the cavity photon field propagating Green's function describing the photon field relaxation, and v (t, t) characterizes the reservoir-induced photon thermal fluctuations. These two Green functions, u (t, t 0 ) and v (t, t), are determined by the following integrodifferential dissipation equation and the nonequilibrium fluctuation-dissipation theorem, respectively 15,16 , v t,t ð Þ~ð where n v,T ð Þ~1 e hv=kB T {1 Â Ã is the initial photon distribution in PCs at temperature T. The spectral density J(v) is microscopically defined as a multiplication of the density of states (DOS) % v ð Þ of PCs with the photon scattering amplitudes V k between the cavity and PCs, In the second equality we have taken the continuous photonic modes of PCs so that V k R V(v), and the index i of V ik in Eq. (14) has also been dropped for single-mode cavity. For arbitrary spectral density J(v), the general solution of Eq. (19) can be obtained exactly 15 , which is given by Eq. (2). By solving Eqs. (16)- (18) and (15) through Eqs. (19)- (20), the complete solution of photon dissipation and fluctuations can be obtained, as we presented in details in the paper.