Statistical theory of photon gas in plasma

TL;DR

This paper explores the thermodynamics of photon gases in plasma, showing how photons acquire effective mass, leading to Bose-Einstein condensation and modified blackbody radiation, with implications for experiments and astrophysics.

Contribution

It provides a statistical physics framework for photon condensation in plasma, highlighting the effects of plasma on blackbody radiation and photon effective mass, and connects theory with recent experiments.

Findings

Photon effective mass in plasma enables Bose-Einstein condensation.

Modified Planck's law with a gap below plasma frequency.

Photon condensate observed in optical microcavity experiments.

Abstract

The thermodynamical properties of the photon-plasma system had been studied using statistical physics approach. Photons develop an effective mass in the medium thus -- as a result of the finite chemical potential -- a photon Bose-Einstein condensation can be achieved by adjusting one of the relevant parameters (temperature, photon density and plasma density) to criticality. Due to the presence of the plasma, Planck's law of blackbody radiation is also modified with the appearance of a gap below the plasma frequency where a condensation peak of coherent radiation arises for the critical system. This is in accordance with recent optical microcavity experiments which are aiming to develop such photon condensate based coherent light sources. The present study is also expected to have applications in other fields of physics such as astronomy and plasma physics.

Tables

Table 1. Table 1 : Exact expressions of the thermodynamical quantities in the x ≫ 1 much-greater-than 𝑥 1 x\gg 1 regime, and their scalings in the critical system for the thermal particles. η ​ ( z ) = Li 5 / 2 ​ ( z ) / Li 3 / 2 ​ ( z ) 𝜂 𝑧 subscript Li 5 2 𝑧 subscript Li 3 2 𝑧 \eta(z)=\text{Li}_{5/2}\left(z\right)/\text{Li}_{3/2}\left(z\right) and κ ​ ( z ) = Li 3 / 2 ​ ( z ) / Li 1 / 2 ​ ( z ) 𝜅 𝑧 subscript Li 3 2 𝑧 subscript Li 1 2 𝑧 \kappa(z)=\text{Li}_{3/2}\left(z\right)/\text{Li}_{1/2}\left(z\right) .

quantity	gaseous phase	scaling at criticality
$\overline{N}$	$\sqrt{\frac{\pi }{2}}\frac{V}{{\lambda }_T^3}{x}^{3/2}{\text{Li}}_{3/2}\left(z\right)$	$\propto T^{3/2}$; $\propto n_p^{3/4}$
$\Phi$	$-\eta \left(z\right)\overline{N}{k}_{B}T$	$\propto -T^{5/2}$; $\propto -n_p^{3/4}$
$\overline{E}$	$\overline{N}{\epsilon }_p+\frac{3}{2}\eta \left(z\right)\overline{N}{k}_{B}T$	$\propto T^{3/2}+\propto T^{5/2}$; $\propto n_p^{3/4}+\propto n_p^{5/4}$
$S$	$\overline{N}{k}_B\left(\frac{5}{2}\eta \left(z\right)-\ln \left(z\right)\right)$	$\propto T^{3/2}$; $\propto n_p^{3/4}$
$P$	$\frac{1}{3}\overline{N}\frac{{\epsilon }_p}{V}+\frac{1}{2}\eta \left(z\right)\overline{N}\frac{k_{B}T}{V}$	$\propto T^{3/2}+\propto T^{5/2}$; $\propto n_p^{3/4}+\propto n_p^{5/4}$
$C_V$	$\frac{3}{2}\overline{N}{k}_B\left(\frac{5}{2}\eta \left(z\right)-\frac{3}{2}\kappa \left(z\right)\right)$	$\propto T^{3/2}$; $\propto n_p^{3/4}$

Equations

$$H$$

$$H\underset{{\mathbf{p}}_i\to 0\forall i}{\overset{C,D}{\to }}\mathcal{H}=\hslash \Omega \left(b^{†}b+\frac{1}{2}\right),$$

$$\rho (\epsilon )=g_s\frac{V}{(2\pi {)}^3}4\pi k^2\frac{dk}{d\epsilon }=\frac{8\pi V}{c^3{h}^3}\epsilon \sqrt{{\epsilon }^2-{\epsilon }_p^2}.$$

$$\Phi =\frac{V}{{\lambda }_T^3}\frac{x^3}{\beta }\underset{1}{\overset{\infty }{\int }}d{\epsilon }^{\prime }{\epsilon }^{\prime }\sqrt{{\epsilon }^{\prime 2}-1}\ln \left(1-e^{-x({\epsilon }^{\prime }-1)}z\right),$$

$$\overline{N}=-{\left({\partial }_{\mu }\Phi \right)}_{V,T},\overline{E}={\left({\partial }_{\beta }\beta \Phi \right)}_{\beta \mu },S=-{\left({\partial }_T\Phi \right)}_{\mu ,V},P=\overline{E}/3V.$$

$$\overline{N}=\frac{V}{{\lambda }_T^3}\underset{1}{\overset{\infty }{\int }}d{\epsilon }^{\prime }\frac{x^3{\epsilon }^{\prime }\sqrt{{\epsilon }^{\prime 2}-1}}{e^{x({\epsilon }^{\prime }-1)}{z}^{-1}-1}.$$

$$\Phi =-\frac{V}{{\lambda }_T^3}\frac{x^2}{\beta }\sum _{j=1}^{\infty }\frac{(e^{x}z{)}^j}{j^2}{K}_2(jx)=\{\begin{array}{ll}-2\frac{V}{{\lambda }_T^3}\frac{1}{\beta }{\text{Li}}_4(e^{\ln z+x}),& \text{for}x\ll 1\\ -\sqrt{\frac{\pi }{2}}\frac{V}{{\lambda }_T^3}\frac{x^{3/2}}{\beta }{\text{Li}}_{5/2}(z),& \text{for}x\gg 1,\end{array}$$

$$\overline{E}=V\underset{{\nu }_p}{\overset{\infty }{\int }}d\nu u(\nu ),\text{with}u(\nu )=\frac{8\pi h}{c^3}\frac{{\nu }^2\sqrt{{\nu }^2-{\nu }_p^2}}{e^{\beta h(\nu -{\nu }_p)}{z}^{-1}-1},$$

$$C=e^{\frac{1}{2}\theta (a^{†}{a}^{†}-aa)},$$

$$D=e^{\sigma b^{†}-{\sigma }^{†}b},$$

$$\mathcal{H}$$

$$\overline{N}=\frac{V}{{\lambda }_{pl}^3}\underset{1}{\overset{\infty }{\int }}d{\epsilon }^{\prime }\frac{{\epsilon }^{\prime }\sqrt{{\epsilon }^{\prime 2}-1}}{e^{x({\epsilon }^{\prime }-1)}{z}^{-1}-1}=\frac{V}{{\lambda }_{pl}^3}{I}_T(x,z),$$

$$\overline{N}=\frac{V}{{\lambda }_T^3}\underset{1}{\overset{\infty }{\int }}d{\epsilon }^{\prime }\frac{x^3{\epsilon }^{\prime }\sqrt{{\epsilon }^{\prime 2}-1}}{e^{x({\epsilon }^{\prime }-1)}{z}^{-1}-1}=\frac{V}{{\lambda }_T^3}{I}_p(x,z).$$

$$P=2\underset{p_x>0}{\int }\frac{d^{3}p}{h^3}2p{v}_{ph}{n}_{BE}co{s}^2\theta =\frac{8\pi }{3h^3}{\int }_0^{\infty }dp{p}^3{v}_{ph}{n}_{BE}.$$

$$v_{ph}=\frac{\Omega }{k}=\frac{1}{k}\sqrt{{\omega }^2+{\omega }_p^2}=\frac{1}{p}\sqrt{(pc{)}^2+{\epsilon }_p^2}.$$

$$P=\frac{8\pi }{3h^3}{\int }_{{\epsilon }_p}^{\infty }d\epsilon \frac{dp}{d\epsilon }{p}^2\epsilon \frac{1}{e^{\beta \left(\epsilon -\mu \right)}-1}=\frac{8\pi }{3h^3{c}^3}{\int }_{{\epsilon }_p}^{\infty }d\epsilon \frac{{\epsilon }^2\sqrt{{\epsilon }^2-{\epsilon }_p^2}}{e^{\beta \left(\epsilon -\mu \right)}-1}.$$

$$\Phi =\frac{V{\beta }^2}{{\lambda }_T^3}{\int }_{{\epsilon }_p}^{\infty }d\epsilon \epsilon \sqrt{{\epsilon }^2-{\epsilon }_p^2}\ln \left(1-e^{-\beta (\epsilon -\mu )}\right)=-\frac{V{\beta }^3}{3{\lambda }_T^3}{\int }_{{\epsilon }_p}^{\infty }d\epsilon {\left({\epsilon }^2-{\epsilon }_p^2\right)}^{3/2}{n}_{BE}(\epsilon )=-\frac{1}{3}\left(\overline{E}-\overline{N}{\epsilon }_p\right),$$

$$\Phi$$

$$\mathcal{L}[\epsilon \sqrt{{\epsilon }^2-1}](jx)$$

$$K_s(x)={\int }_0^{\infty }dt\cosh st{e}^{-x\cosh t},$$

$$\frac{2s}{x}{K}_s(x)=e^{(s-1)\pi i}{K}_{s-1}(x)-e^{(s+1)\pi i}{K}_{s+1}(x).$$

$$\Phi =-\frac{V}{{\lambda }_T^3}\frac{x^2}{\beta }\sum _{j=1}^{\infty }\frac{(e^{x}z{)}^j}{j^2}{K}_2(jx).$$

$$-\frac{V}{{\lambda }_T^3}\frac{x^2}{\beta }\sum _{j=1}^{\infty }\frac{\phi (j)}{j^s},\text{with s=2, and}\phi (j)=(e^{x}z{)}^j{K}_2(jx).$$

$$\overline{N}$$

$$\overline{E}$$

$$S$$

$$P$$

$$C_V$$

$$\Phi =-\frac{V}{{\lambda }_T^3}\frac{2}{\beta }\sum _{j=1}^{\infty }\frac{(e^{x}z{)}^j}{j^4}=-2\frac{V}{{\lambda }_T^3}\frac{1}{\beta }{\text{Li}}_4(e^{\ln z+x}).$$

$$\Phi =-\frac{V}{{\lambda }_T^3}\sqrt{\frac{\pi }{2}}\frac{x^{3/2}}{\beta }\sum _{j=1}^{\infty }\frac{z^j}{j^{5/2}}=-\frac{V}{{\lambda }_T^3}\sqrt{\frac{\pi }{2}}\frac{x^{3/2}}{\beta }{\text{Li}}_{5/2}(z).$$

$$\frac{d}{dy}{\text{Li}}_s(y)=\frac{1}{y}{\text{Li}}_{s-1}(y),$$

$$\overline{N}$$

Full text

I Introduction

Quantum gases are one of the most studied subjects in physics and have well-established statistical descriptions which have been verified by high precision experimental observations. In particular, the theory of the ideal Bose gas (IBG) predicts the phenomenon of the Bose-Einstein condensation (BEC). It was first experimentally realized by using gaseous rubidium and sodium in $1995$ and resulted in a joint Nobel Prize for Cornell, Wiemann and Ketterle, Wiemann ; Ketterle . However, this discovery occurred over seventy years after Bose and Einstein’s prediction Bose ; Einstein , showing the numerous technological challenges that the researchers had overcame in order to demonstrate the existence of BECs. Einstein’s derivation used a modification of Bose’s work on photon gas to describe the statistics of massive particles. In fact, their work showed that only massive gases (following Bose statistics) are able to go through the condensation procedure. Hence, the initially studied photon gas, being massless, does not possess this property despite the bosonic nature of the light quanta. However, condensation is possible by changing experimental conditions, e.g., by considering the photon gas interacting with a plasma. In this scenario, the photon dispersion relation is modified due to the collective oscillation of the charged particles and an effective ”photon mass” can be introduced Anderson ; Varro1 ; Aryeh ; Mendonca1 ; Matip . The mass generation through interaction was first proposed by Anderson using the examples of plasmon theory and superconductivity Anderson . In its more general formulation, this phenomenon became known as the Brout-Englert-Higgs mechanism after the authors who proved the possibility of mass generation thorough spontaneous gauge symmetry breaking in relativistic quantum field theories Higgs ; Englert but others also contributed to the topic Kibble . This idea made it possible to explain the existence of massive vector bosons in the electroweak sector of the Standard Model by Weinberg and Salam Weinberg ; Salam .

In the present paper, the corresponding massive quasiparticle, bosonic in nature, is known as bulk plasmon-polariton 111Throughout the text the terminology ’(bulk) plasmon-polaritons’ and ’photons’ indicate the same physical object while the plasma is present. Photons without interacting with the plasma are referred to as ’free photons’ or ’free photon gas’. It is also worth to emphasize that the plasmon and the plasmon-polariton are not the same: the former is the quanta of electrostatic oscillations, whereas the latter is the quanta of transverse radiation in the plasma. and can be identified as the two transverse mode of the electromagnetic (EM) field in the plasma like in Anderson’s classical argument. It is emerged as the fundamental degrees of freedom in the diagonalized Hamiltonian consisting of charges coupled to a single EM mode Varro1 ; Aryeh ; Matip . The bosonic ensemble built up from plasmon-polaritons must have different statistical properties compared to the free photon gas. Indeed, as a consequence of the effective mass, if the (quasi)particle number is conserved, a finite chemical potential can be introduced, and it makes it possible for the system to undergo the BEC phase transition. The photon-plasma interaction has already been discussed in several papers which even led to the conclusion of the possible existence of the photon BEC. Kompaneets considered a kinetic equation of the photon distribution (Fokker-Planck equation) in plasma where the interaction is mediated exclusively by Compton scatterings, i.e. no emission and absorption processes have been taken into account Kompa . This eventually leads the system to thermal equilibrium. Zel’dovich and Levich, using the results derived in Kompa , showed that the system could possess a BEC phase Zeldo . This topic was then revisited in Levermore where the same result was concluded: it was shown that by including non-number conserving processes the BEC phase cannot develop. However, as no effective ”photon mass” was considered in the approach used in Kompa ; Zeldo ; Levermore , they could only provide valid results in very dilute plasmas where the dispersion relation of the plasmon-polariton has a negligible mass term, i.e. it is the same as the free photon. In Calleb , the effect of the plasma on the thermodynamics of photons has been considered with the assumption that the chemical potential is equal to the ”photon mass”, however, no BEC was found in this case due to the omission of the ground state energy. In Tsint , a kinetic theory of the photon BEC was considered for radiations with sufficiently large intensity (e.g. strong laser pulse) where the photon-photon interactions are the dominant, hence the photon-electron interactions with absorptions are negligible. Different aspects of possible BEC formation are considered in Mendonca2 with the study of the longitudinal mode (plasmon) and four-wave mixing, by using also kinetic theory, where the photon-number conservation is assumed as a consequence of including only Compton scatterings. In both of the latter two studies the bulk plasmon-polariton dispersion relation was used. Despite these findings are correct, as of yet, no comprehensive theory on the statistical description of such system was provided which is the main focus of the present paper. However, there exist other different mechanisms to form a BEC in a photon gas. In Boyce ; Weitz , the photon effective mass comes from the paraxial approximation of field quantization in the Fabry-Pérot microcavity and the thermalization is a consequence of the interaction with dye molecules that fill the cavity Weitz . Experimental observations of photon BEC in the microcavity have been reported in Weitz2 . Further studies of the photon BEC based on variations of the previous ideas can be found in Kruch1 ; Kruch2 ; Kruch3 ; Boi ; Walker1 ; Walker2 .

II Model for the photon-plasma system

A statistical model of the ideal plasmon-polariton gas has been considered using the grand canonical ensemble framework in order to describe the physics of a photon gas in a homogeneous, isotropic plasma in the current paper. The following Hamiltonian defines a system of a $N_{ch}$ charges interacting with a quantized monochromatic EM field,

[TABLE]

where ${\bf{p}}_{i}$’s and $m$ are the momenta and the mass of the charges which have unit charge $e$ and the summation is over the $N_{ch}$ number of these charged particles. The linearly polarized EM mode is represented by the vector potential $\bf{A}$, which in the dipole approximation has the form of ${\bf{A}}=\alpha\,\mathbfcal{E}(a+a^{\dagger})$, where the terms $a^{\dagger}$ and $a$ are the creation and annihilation operators of the mode with angular frequency $\omega$. The parameter $\alpha=c\sqrt{2\pi\hbar/V\omega}$, with the light speed $c$, the quantization volume $V$ and the Planck constant $\hbar=h/2\pi$. The real unit vector $\mathbfcal{E}$ gives the direction of the polarization. In the case of a single charge in (1) ($N_{ch}=1$), the Hamiltonian was exactly diagonalized by Varró and Bergou Varro1 using a displacement ($D$) and a Bogoliubov ($C$) transformation in order to eliminate the linear and quadratic terms of the ladder operators. The resulting Hamiltonian was used to describe nonlinear scattering processes. The solution was generalized to the elliptically polarized cases Matip , as well as to $N_{ch}>1$ charges which was also done in Aryeh for a similar model. In this scenario, the system can be considered as a plasma where the Coulomb interactions are damped between electrons and ions by Debye screening – a free electron gas. The homogeneity of the plasma is achieved by taking the limit ${\bf{p}}_{i}\to\boldsymbol{0}$ for all $i$. This gives a uniformly distributed $N_{ch}e$ net charge in the box (volume $V$) according to the Heisenberg’s uncertainity principle Matip . Applying the above described operations on (1) results in an effective Hamiltonian of a free harmonic quantum oscillator describing a plasmon-polariton system Matip :

[TABLE]

where $b^{(\dagger)}=C^{-1}a^{(\dagger)}C$ are the annihilation and the creation operators of the plasmon-polariton quasiparticles, with $[b,b^{\dagger}]=1$. The effective frequency associated to the quasiparticles is defined as $\Omega=\sqrt{\omega^{2}+\omega_{p}^{2}}$, where $\omega_{p}=\sqrt{4\pi e^{2}n_{p}/m}$ defines the plasma frequency with the plasma density $n_{p}=N_{ch}/V$. Likewise, the ”plasma energy” can be obtained as $\varepsilon_{p}=\hbar\omega_{p}$. By observing the similarity of the dispersion relation to a relativistic massive particle’s this is often considered to be the ”rest energy” of the photon in plasma Anderson . In (2) the $\hbar\Omega/2$ term represents the vacuum energy, and thus will be omitted throughout this analysis. Some details of the diagonalization is presented in Appendix A. The grand potential is defined as $\Phi=\beta^{-1}\sum_{\xi}\ln\left(1-e^{-\beta(\varepsilon_{\xi}-\mu)}\right)$, and hence the quantum statistics of the system follows the Bose-Einstein distribution, $n_{BE}=1/(e^{\beta(\varepsilon-\mu)}-1)$, with $\beta=1/k_{B}T$, where $k_{B}$ is the Boltzmann constant and $T$ is the temperature; $\xi$ indicates a given one-particle state with energy $\varepsilon_{\xi}$ and $\mu$ is the chemical potential associated to the conserving number of quasi-particles. The finite chemical potential makes the properties of the plasmon-polariton system crucially different from the vacuum scenario where the photons are massless, hence it is nonsensical to talk about a definite number of photons in the system as they can be created even with an infinitesimal amount of energy, i.e., $\mu=0$. Thermal equilibrium between the photon gas and the plasma (as well as conserving photon number) is assumed throughout the paper, which can be achieved by one of the previously mentioned mechanisms Kompa ; Zeldo ; Levermore ; Calleb ; Tsint ; Mendonca2 .

III Thermodynamics and Bose-Einstein condensate

Taking the thermodynamic limit, the summation over the one-particle states is replaced by the integral over the phase space: $\sum_{\xi}\longrightarrow\sum_{\text{spin}}\frac{V}{h^{3}}\int d^{3}p$, where the summation goes for the spin degeneracy and gives a factor of $g_{s}$, provided there is no spin dependence. The introduction of the density of states (DOS) enables the integral to be rewritten as $g_{s}\frac{V}{h^{3}}\int d^{3}p=\int d\varepsilon\,\rho(\varepsilon)$ with the DOS

[TABLE]

In the present case $g_{s}=2$ as the degeneracy consists of the two transverse modes. The DOS of the free photon gas, $\rho=8\pi V\varepsilon^{2}/(ch)^{3}$, is naturally recovered in the $\varepsilon_{p}/\varepsilon\to 0$ limit. The grand potential is

[TABLE]

where the following dimensionless quantities are introduced: $\varepsilon^{\prime}=\varepsilon/\varepsilon_{p}$, $x=\beta\varepsilon_{p}$, $\mu_{p}=\mu/\varepsilon_{p}$ and the fugacity $z=\exp{x(\mu_{p}-1)}$. The lower bound of the integral is set to $\varepsilon_{p}$ as there is no real contribution below the ground state energy. The quantity $\lambda_{T}=\beta ch/2\pi^{1/3}$ is the thermal de Broglie wavelength for massless particles. It is possible to introduce another analogous length scale $\lambda_{pl}=ch/2\pi^{1/3}\varepsilon_{p}$ which also will be used throughout the text. The integral is well-defined only for chemical potential values $\mu_{p}\leq 1$, hence $z\in[0,1]$. Using the expression for $\Phi$, all the relevant thermodynamical quantities can be derived; average particle number, energy, entropy and pressure:

[TABLE]

Hence, the average photon number in plasma is

[TABLE]

However, in the thermodynamic limit it only makes sense to talk about the particle density as $V$ and $\overline{N}\to\infty$ while $\overline{n}=\overline{N}/V$ is kept fixed. The same is true for the plasma density $n_{p}$, encapsulated by $\varepsilon_{p}$. These conditions are applied throughout the paper. The integral in (6) is a monotonically increasing function of $z$ and $z\leq 1$ so the maximal value is when $z=1$, where the system is found to be critical, i.e. $\overline{n}_{c}=\overline{n}(z=1)$. As the DOS gives zero weight to the ground state, the formula in (6) provides the particle number only for the thermal states, and above the critical threshold the further particles accumulate in the ground state with energy $\varepsilon_{p}$ by forming a BEC: $\overline{n}_{tot}=\overline{n}_{\varepsilon_{p}}+\overline{n}_{\varepsilon>\varepsilon_{p}}$ (where $\overline{n}_{\varepsilon>\varepsilon_{p}}$ is used for $\overline{n}$ emphasizing its thermal attribute). The same phenomenon happens for massive IBG in three dimensions when reaching the critical temperature or particle density. Besides the photon (plasmon-polariton) density and the temperature, the plasma density is also able to drive the current system to criticality. Thus, the chemical potential depends on all of these parameters and a critical value exists for each of them for $\mu_{p}=1$ or $z=1$, Fig. 1(a). The critical particle density hence can be expressed as a function of $T$ and $n_{p}$, Fig. 1(b) and Fig. 1(c): by fixing the photon density $\overline{n}$ and decreasing $T$ or $n_{p}$ the system reaches criticality at $T_{c}$ and $n_{p}^{c}$, respectively. Decreasing these parameters further a fraction of $\overline{n}_{\varepsilon_{p}}$ photons condense in the ground state while the rest remain in the thermal states with $\overline{n}_{\varepsilon>\varepsilon_{p}}$. However, there is a crucial difference between the two parameter dependences. When the temperature reaches $T_{c}$ the chemical potential becomes $\mu=\varepsilon_{p}$ and preserving this value as $T\to 0$. On the other hand, when the plasma density reaches its critical value the chemical potential acquires the value $\varepsilon_{p}$ again, but in this case $\mu$ must go to zero as $n_{p}\to 0$, since $\varepsilon_{p}\propto\sqrt{n_{p}}$ and the chemical potential cannot exceed the value of $\varepsilon_{p}$. This can be seen in Fig. 1(c): as $n_{p}$ goes to zero the photon number in the thermal states reaches a finite value (dashed line), $\overline{n}=16\pi\zeta(3)/\lambda_{T}^{3}$ – this is exactly the value for the average particle density in a free photon gas. Thus, by lowering the plasma density, a fraction of the initial number of photons in the plasma indeed condense, however, the particle density in the thermal states cannot go below its vacuum value which is determined by the free photon gas. Details of the numerical algorithms for the plots presented in Fig. 1 can be found in Appendix B.

The fingerprint of the phase transition can be observed in the behavior of the thermodynamical quantities usually manifesting as a change in the analytic properties at the given critical parameter. The thermodynamic quantities listed in (5) are shown in Fig. 2. In Fig. 2(a), Fig. 2(b) and Fig. 2(c) the temperature and plasma density dependences of the total energy, the entropy and the heat capacity at fixed particle density are given. The heat capacity is used as the main indicator of phase transitions due to its parameter sensitivity around criticality: in the current case a cusp can be observed at $T_{c}$, which is characteristic to BEC transitions. Similar behavior can be found in case of the IBG and also in interacting Bose systems such as the liquid ${}^{4}He$, where the heat capacity even gets logarithmically singular at the critical temperature, called the ”lambda-transition” Huang ; lambdatr . These quantities are obtained through numerical integration by using the formulae in (4) and (5). In Appendix C the pressure derived from basic principles results in the formula $P=\overline{E}/3V$, which coincides with the free photon case rather than the massive IBG. The pressure dependence on the specific volumes of the photons $\overline{v}=1/\overline{n}$ and plasma $v_{p}=1/n_{p}$ is in Fig. 2(d). It behaves much like for the IBG with respect to $\overline{v}$, at a constant plasma density: below the critical value, $\overline{v}_{c}$, the pressure gets independent of the specific volume of the photons. On the other hand, by keeping $\overline{n}$ fixed, below the critical plasma density, $n_{p}^{c}=1/v_{p}^{c}$, the system enters the BEC phase, whereas for $v_{p}\to 0$ the pressure diverges showing the incompressibility of the plasma. All of these thermodynamical quantities are approaching the values defined by the free photon gas as $n_{p}\to 0$, indicated by dashed lines in Fig. 2.

A better qualitative insight is obtained by applying approximations to enable exact analytical solutions of (4). It is obvious that the energy scale of the system is determined by the temperature and the plasma energy, hence their ratio, $x=\varepsilon_{p}/k_{B}T$, is a crucial dimensionless parameter. The formula for the grand potential in (4) can be recast in the form of an infinite series of Dirichlet-type (Appendix D.1):

[TABLE]

where $K_{s}(y)$ is the Macdonald function with index $s$. The sum can be computed for the limiting cases, as shown in (7), where the function $\text{Li}_{s}(y)$ denotes the polylogarithm with index $s$. In the regime where the temperature dominates the plasma energy ($x\ll 1$), the grand potential $\propto\text{Li}_{4}(e^{x\mu_{p}})$. This approximation fails to study the BEC phenomenon as the domain of the polylogarithm for real values is well defined only in the interval $(-\infty,1]$ and from the previous analysis, the BEC occurs at $\mu_{p}=1$. It only describes the system correctly for $\mu_{p}\leq 0$.

However, in the case when $x\approx 0$, the polylogarithm becomes $\text{Li}_{s}(1)=\zeta(s)$ (where $\zeta$ is the Riemann zeta function), and hence $\Phi=-4\sigma VT^{4}/3c$ with the Stefan-Boltzmann constant $\sigma=2\pi^{5}k_{B}^{4}/15c^{2}h^{3}$, which is exactly the free energy of the free photon gas. Thus, in this extreme regime, by using the formulae in (5), the statistics of the free photon gas is reproduced (see in Appendix E.1).

By considering the region where the temperature is much less than the plasma energy ($x\gg 1$), the argument of the polylogarithm only contains $z$ which equals to unity at $\mu_{p}=1$. Thus, the formation of the BEC can be correctly described. The thermodynamical quantities are shown in Table 1 for both the gaseous and the critical system. The detailed derivation is in Appendix E.2. The critical quantities show power-law behavior both in temperature and plasma density and their scalings in $T$ are the same as in the case of the IBG with a non-zero ground state energy in the BEC phase. The particle density in the condensate fraction can also be expressed in the same way as for the IBG: $\overline{n}_{\varepsilon_{p}}=\overline{n}_{tot}(1-T/T_{c})^{3/2}$, matching the findings in Tsint and Mendonca2 . These results should not be surprising as the dispersion relation of the plasmon-polariton is formally the same as those of the massive relativistic IBG, for which the same functional form of the thermodynamical quantities can be found as for the current system bose1 ; bose2 ; bose3 . The critical value for the temperature and the plasma density can be read off from the critical expression of $\overline{N}$ in Table 1: $T_{c}\propto\overline{n}^{2/3}/\sqrt{n_{p}}$ and $n_{p}^{c}\propto\overline{n}^{4/3}/T^{2}$.

IV Modified blackbody radiation

In different regimes, the system behaves more as a free photon gas ($x\ll 1$) or as a massive Bose gas ($x\gg 1$) but the system describes the behavior of radiation in plasma. Hence, it is sensible to ask whether the thermal radiation is modified in the presence of the medium. The total energy is

[TABLE]

where the frequency $\nu_{(p)}=\omega_{(p)}/2\pi$. This modified form of the blackbody spectral density with changing parameters is shown in Fig. 3.

All curves start at $\nu_{p}$ with a hard cut-off as there is no light propagation below the plasma frequency. In the gaseous phase, the curves are similar to the usual blackbody radiation (with a gap) but a sharp peak appears at $\nu_{p}$ once the particular parameter is tuned to or beyond criticality. This happens when the radiation coming from the ground state start to dominate the radiation associated to the thermal particles. At the critical value (when $z$ reaches unity), the spectrum develops a singularity at the plasma frequency as the particles start to accumulate in a coherent BEC phase. Even though the spectral density becomes singular at $\nu_{p}$ it remains normalizable, hence the overall energy does not diverge. For the general shape comparison of the non-critical, critical and Planckian spectra see Fig. 3(a). In case of tuning $T$, the thermal tail of the distribution vanishes as more particle enters the BEC state, and the normalized distribution eventually tends to a Dirac-delta concentrated at the plasma frequency. Hence, the total energy of the system becomes $\overline{E}=\overline{N}\varepsilon_{p}$ as all the particles joined the condensate, Fig. 3(b). The plasma density dependence is shown in Fig. 3(c). In this case while $n_{p}$ decreases, the gap $\nu_{p}$ shifts towards zero since its value is determined by the plasma density. The condensation peak appears at the critical plasma density and this singularity eventually ends up at $\nu_{p}=0$, for which the spectrum becomes the blackbody radiation in vacuum. Nevertheless, a fraction of photons are still in the condensate but with energy $\varepsilon_{p}\approx 0$ there is no way to observe them in the spectrum. This somewhat resembles the infrared catastrophe in quantum electrodynamics where infinitely many soft photons are radiated during the brehmsstrahlung process, however, no realistic detector can register them due to their finite resolution BN ; JakoMati ; Mati . By adjusting the photon density $\overline{n}$ above $\overline{n}_{c}$ does not change the shape of the distribution further as the system cannot have a higher particle density than its critical threshold, Fig. 3(d). Indeed, any additional particle ends up immediately in the BEC, or by compressing the system, the particles from the normal phase join the BEC in a way that preserves the $n_{c}$ density of the thermal particles. The condensation peak shows up once the critical particle density has been reached. A similar shape of the photon spectrum is derived numerically in Mendonca2 where the kinetic equation, taking into account the processes of Compton and inverse Compton scattering exclusively, drives the plasma-photon system to thermal equilibrium with the condensation peak at the plasma frequency for the right parameters. These results are extremely resembling the findings in Weitz2 , where the peak arises at the cut-off frequency of the cavity, representing the effective photon mass in this case. The idea put forward in Weitz and Weitz2 that the photon BEC in an optical microcavity can be used as coherent light source are applicable for the present system. The spectrum presented in Fig. 3 should be producible with the appropriate experimental setup, for which the present study could serve as a theoretical background.

V Discussion

The study presented in this paper gives a comprehensive theoretical description of the plasma-photon statistical system. It is shown that the BEC formation is possible for a photon gas in a plasma medium by tuning its relevant parameters to their critical values. The radiation spectra shown in Fig. 3 represent the modified blackbody radiation. The distribution develops a sharp condensation peak at the plasma frequency in the BEC state corresponding to coherent radiation which should be detectable in experiment, and perhaps could be used as a new type of source for coherent radiation, likewise in the case of the cavity photon BEC Weitz2 . The thermodynamical properties of the plasma-photon system should also have some potential application in other fields of physics. For instance, non-Planckian modifications of the radiation spectra could be detected in the cosmic microwave background at low frequencies NonPl . Furthermore, the modification of the thermodynamical quantities, e.g. for the radiation pressure and the emitted total energy, might also have some relevance for modeling massive stellar interiors Calleb ; Stellar1 ; Stellar2 and plasma transport properties tokamak1 ; tokamak2 .

Acknowledgments

The author thanks Sándor Varró and Antal Jakovác for useful discussions and comments. The author wishes to thank Piroska Ördög Józsefné and József Ördög for all the kind hospitality shown to him at the Piroska Életmód Vendégház and Andrew Cheesman for thoroughly reading through the manuscript. The ELI-ALPS project (GINOP-2.3.6-15-2015-00001) is supported by the European Union and cofinanced by the European Regional Development Fund.

Appendix A Hamiltonian diagonalization

In order to diagonalize the Hamiltonian in (1) a Bogoliubov transformation must be used to cancel the quadratic terms of $a^{(\dagger)}$:

[TABLE]

where $\theta$ is a real parameter. The quasiparticle ladder operators are defined by its action as $b^{(\dagger)}=C^{-1}a^{(\dagger)}C$, with $[b,b^{\dagger}]=1$. And a shift must be performed on $b^{(\dagger)}$ in order to take care of its linear terms after the action of $C$. This is done by the displacement operator:

[TABLE]

where $\sigma$ is its parameter. By choosing the parameters $\theta$ and $\sigma$ appropriately, the following form of the Hamiltonian can be achieved:

[TABLE]

Further details of the diagonalization can be found in Varro1 and Matip . The resulting Hamiltonian is the sum of $N_{ch}$ free charges and a quantum harmonic oscillator with a shifted particle number by $|\sigma|^{2}$. $\sigma$ depends linearly on the momenta $\{{\bf{p}}_{i}\}$. Thus, by taking the limit ${\bf{p}}_{i}\to{\bf{0}}$ for all $i$ results in a vanishing $\sigma$, and hence (2) remains.

Appendix B Details of the numerical algorithms for the exact case

Numerical evaluations of (4) along with the relations in (5) provide the thermodynamical quantities depicted in Fig. 2. However, this requires the derivation of the parameter dependence of the fugacity, which can be done numerically.

B.1 Temperature dependence of the fugacity

In order to derive the temperature dependence of the fugacity it is instructive to rewrite the formula in (6) as

[TABLE]

where $\lambda_{pl}=ch/2\pi^{1/3}\varepsilon_{p}$, and the temperature dependence only present in $x=\beta\varepsilon_{p}$ in the integral part $I_{T}$. By fixing $\overline{n}$ and solving the equation $\overline{n}\lambda_{pl}^{3}=I_{T}(x,1)$ numerically, the critical temperature is obtained, $T_{c}=T_{c}\left(\overline{n}\right)$, for the given $\overline{n}$. In the next step the equation $\overline{n}\lambda_{pl}^{3}=I_{T}(\chi x_{c}(\overline{n}),z)$ is solved for $z$ for a given $\overline{n}$ numerically, where $\chi=x/x_{c}(\overline{n})$ and $x_{c}=\beta_{c}\varepsilon_{p}$. By interpolating the data points a $z(\tau(\overline{n}))$ function is defined, with $\tau=T/T_{c}$, from which the chemical potential read as $\mu_{p}(\tau)=\ln(z(\tau))\tau T_{c}+1$. The temperature dependence of the fugacity is shown in Fig. 1(a) for fixed $\overline{n}=5/\lambda_{pl}^{3}$.

B.2 Plasma density dependence of the fugacity

For the plasma density dependence it is advisable to use (6)

[TABLE]

In this formula the $\varepsilon_{p}$ dependence is only present in $I_{p}(x,z)$. The plasma energy is expressed with the plasma density as $\varepsilon_{p}=\hbar\sqrt{4\pi e^{2}n_{p}/m}$. In order to obtain the density dependence of the fugacity, the same procedure must be done as for the temperature. The plasma density dependence of the fugacity is shown in Fig. 1(a) for $\overline{n}=5/\lambda_{T}^{3}$.

B.3 Photon number density of the fugacity and the chemical potential

Using (12) with fixed $x$, the critical particle density is $\overline{n}_{c}=I_{T}(x,1)$ in units of $1/\lambda_{T}^{3}$ for that given $x$. Solving the same equation $\overline{n}=I_{T}(x,z)$ now for $z$ (again numerically), and using an interpolation gives the function $z(\overline{n}/\overline{n}_{c})$, which is plotted for $x=1$ in Fig. 1(a).

Appendix C Pressure

It is not trivial if the equation of state is the same as for the massless Bose gas, i.e. $\overline{E}=pV/3$. However, it turns out that this is the case which will be proven in the following by using the argument in Huang . The integral form of the pressure is given as the momentum transfer per photons times the flux of photons. Considering photons reflecting from a wall that is normal to the $x$ axis the photon flux is $v_{ph}\cos\theta$, where $v_{ph}$ is the phase velocity in of the radiation in plasma and $\theta$ is the angle enclosed by the momentum of the photon and the $x$ axis. It is shown in Press1 ; Press2 that the phase velocity is the relevant quantity to use when computing radiation pressure in a dispersive medium. This is in agreement with experiments presented in PressExp . The momentum passed over to the wall by each reflection is $2p\cos\theta$. Hence the total pressure has the form of

[TABLE]

Here, the spherical coordinates were used in the second expression for which the restriction of $p_{x}$ means $0<\theta<\pi/2$ and the factor two comes from the number of polarizations of the EM field. The phase velocity corresponding to the radiation in plasma

[TABLE]

Substituting into (14) and switching to energy variable, $\varepsilon=\sqrt{(pc)^{2}+\varepsilon_{p}^{2}}$,

[TABLE]

This implies that $P=\overline{E}/3V$ as the well-known relation for the photon gas. However, the relationship between the grand potential and the pressure, $\Phi=-PV$ , is no longer true:

[TABLE]

where integration by parts was used. The grand potential does not include the ”effective mass” of the photon. However, as it is acquired by the plasma oscillation it also should contribute to the exerted pressure on the container wall, hence (16) yields the correct formula.

Appendix D Series representations of thermodynamic quantities and their convergence

By expanding the integrand into power series in (4) a series representation of the grand potential and the other thermodynamical quantities can be achieved by applying (5). Even though the results are exact, the sums do not have a closed forms.

D.1 Grand potential

The grand canonical potential of a the system is defined thorough the integral in (4), and hence

[TABLE]

where $\mathcal{L}$ is for the Laplace transform. This Laplace transform can be computed exactly as follows

[TABLE]

Here the integral representation of the Macdonald function is used

[TABLE]

together with the recurrence relation

[TABLE]

Thus the series for the grand potential reads

[TABLE]

Switching the integration with the summation in (18) is well-justified in the present case: after the series expansion the integrand reads $f=\sum_{j}f_{j}=\sum_{j}\varepsilon\sqrt{\varepsilon^{2}-1}[\exp(x)z]^{j}\exp(-jx\varepsilon)/j$, where $\{f_{j}\}_{j}$ is a sequence of positive and integrable functions with a convergent sum. These properties imply that the integration and summation indeed commute and the series representation of the grand potential in (22) converges to the same value as the integral. However, it is not hard to see the convergence of (22) alone by using Abel’s test: the series $\sum_{j}1/j^{2}$ is convergent and $(e^{x}z)^{j}K_{2}(jx)$ monotonically decays to zero which implies the convergence.

The series defined in (22) is a Dirichlet series Apostol of the form

[TABLE]

Moreover, since the sum only with the term $\varphi(1)$ corresponds to the grand potential in the classical limit, $1/\left(\exp\beta(\varepsilon-\mu)-1\right)\approx\exp\beta(\mu-\varepsilon)$, the series in (22) can also be considered as the Dirichlet transform of the classical grand potential. From the general theory of the Dirichlet series with variable $s=\sigma+it$, it is also known that there exists $\sigma_{c}\in\mathbb{R}$ for which if $\sigma>\sigma_{c}$ the series is convergent, hence $\sigma_{c}$ cuts the plane into halves and the region $\sigma>\sigma_{c}$ is called the half-plane of convergence. Since (22) is shown to be convergent, $s=\sigma=2$ is in the half-plane of convergence. Furthermore, if $\sigma$ is in a compact subset of this half-plane then the series is uniformly convergent (Theorem 11.11, Apostol ). Since for (22) this is definitely satisfied, it is uniformly convergent, too. This property is useful when other thermodynamical quantities are computed through the relations (5) where differentiations are involved. The uniform convergence makes it possible to execute the differentiation term by term.

D.2 Other thermodynamical quantities

As it is discussed in D.1 the sum obtained as the series representation for the grand potential in (22) is uniformly convergent. Hence the term-by-term differentiation is possible without altering its result. This property can be exploited when computing the other thermodynamical quantities by using (5):

[TABLE]

For the heat capacity the formula $C_{V}=T(\partial_{T}S)_{V}$ is used as it is easier to handle than the $\beta$ derivative of the energy at fixed particle number. It also can be shown that these expressions are convergent. However, the fugacity $z$ still need to be determined for the non-critical regime which makes the numerical treatment of (4) much more comfortable than the series representation. $z$ simply can be set to unity in the BEC phase.

D.3 Asymptotics

The series given in (22) cannot be summed up to get a closed form, however, it is possible in the limiting cases.

• •

For the case when $x\ll 1$ (and $s>0$), the Macdonald function $K_{s}(y)=2^{s-1}\Gamma(s)y^{-s}+O(y^{2})$, hence

[TABLE]

• •

For the case when $x\gg 1$ (and $s>0$), the Macdonald function $K_{s}(y)=e^{-y}\left[\sqrt{\frac{\pi}{2}}\sqrt{\frac{1}{y}}+O\left(\frac{1}{y^{3/2}}\right)\right]$, hence

[TABLE]

Here, at the summation, the definition of the polylogarithm is used, $\text{Li}_{s}(y)=\sum\limits_{j=1}^{\infty}y^{j}/j^{s}$.

Appendix E Analytical expressions of the thermodynamical quantities in the asymptotic regimes

In D.3 the asymptotic expansions of the Macdonald function (for $x\ll 1$ and $x\gg 1$) made it possible to sum up the series and obtain a closed form for the grand potential. Using the relations in (5) and the differentiation identity of the polylogarithm,

[TABLE]

can provide further analytical expressions for the thermodynamical quantities in these limits.

E.1 High temperature and low plasma density limit ($x\ll 1$)

The thermodynamical quantities in the $x\ll 1$ limit are found to be the following expressions

[TABLE]

The average particle number is kept fixed when computing the heat capacity. Therefore

[TABLE]

And thus $dz/dT$ can be expressed as

[TABLE]

The total energy can be expressed through the average particle number as

[TABLE]

where $\gamma(e^{x}z)=\text{Li}_{4}(e^{x}z)/\text{Li}_{3}(e^{x}z)$. The specific heat at constant volume is obtained by differentiating the energy with respect to the temperature

[TABLE]

where in the third line the expression of $dz/dT$ is used from (30). The argument of the polylogarithms can be rewritten as $\exp(x\mu_{p})$, and thus it shows that this approximation is not suitable for $\mu_{p}>0$, since for these values the polylogarithm becomes complex. As a consequence, these expressions cannot describe the BEC phenomenon as it would require $\mu_{p}=1$. However, at the extreme limit of very high temperatures, i.e. $x\approx 0$, the polyloglogarithms take the form $\text{Li}_{s}(e^{x\mu_{p}})=\text{Li}_{s}(1)=\zeta(s)$ and $dz/dT=0$ from (32). This results that (E.1) and (32) reproduce the thermodynamics of the free photon gas:

[TABLE]

where $\sigma=2\pi^{5}k_{B}^{4}/15c^{2}h^{3}$ is the Stefan-Boltzmann constant, $F$ is the free energy, which in this limit is equal to the grand potential $\Phi$.

E.2 Low temperature and high plasma density limit ($x\gg 1$)

The thermodynamical quantities in the $x\gg 1$ limit are found to be the following expressions

[TABLE]

When computing the specific heat, it must be remembered that the average particle number is kept fixed. Therefore

[TABLE]

And thus $dz/dT$ can be written as

[TABLE]

The total energy can be expressed through the average particle number as

[TABLE]

where $\eta(z)=\text{Li}_{5/2}(z)/\text{Li}_{3/2}(z)$. The specific heat at constant volume is obtained by differentiating the energy with respect to the temperature

[TABLE]

where in the third line the expression of $dz/dT$ is used from (36). Tuning the system to criticality requires $z=1$, for which the polylogarithm $\text{Li}_{s}(1)=\zeta(s)$ when $s\geq 1$ and infinite for $s<1$. Thus the formulae in (E.2) and (38) become

[TABLE]

All the quantities in (E.2) and (38) can be expressed as function of the average particle number which are given in Table 1.

E.3 The fugacity function

In the following the fugacity function in the $x\gg 1$ regime for Appendix E.2 is derived. However, its parameter dependence cannot be obtained in a closed form but can be derived numerically by solving the equation of the average particle density for $z$ from:

[TABLE]

Reduced parameters are used, that are normalized with respect to their critical value. In the following, the temperature dependence is computed, however, the procedure is similar for the $n_{p}$ and $\overline{n}$ dependence. Solving (40) for the temperature gives

[TABLE]

The second formula above gives the critical temperature. From (40), using $\tau=T/T_{c}$ follows

[TABLE]

Hence, the $z$ dependence of $\tau$ is derived from which numerically the $z=z(\tau)$ relation can be obtained, and of course the chemical potential, too, through $\mu=\varepsilon_{p}+\ln(z)k_{B}T$. By using interpolation, the function of $z$ and the chemical potential $\mu$ is depicted in Fig. 4. The full parameter dependence of the thermodynamical quantities is obtained by inserting the numerically defined function $z$ into the expressions of Appendix E.2. An interpolating function for the fugacity $z$ in the case of Appendix E.1 ($x\ll 1$ regime) can be computed in the similar fashion, however, criticality in that case is impossible to achieve, thus reduced parameters, such as $\tau$, can not be defined with the critical parameter.