Density Distribution of a Bose-Einstein Condensate of Photons in a Dye-Filled Microcavity

The achievement of Bose-Einstein condensation of photons (phBEC) in a dye-filled microcavity has led to a renewed interest in the density distribution of the ideal Bose gas in a two-dimensional harmonic oscillator. We present measurements of the radial profile of photons inside the microcavity below and above the critical point for phBEC with a good signal-to-noise ratio. We obtain a good agreement with theoretical profiles obtained using exact summation of eigenstates.


I. INTRODUCTION
Since the realization of Bose-Einstein condensation (BEC) in dilute atomic gases [1,2], the density distribution of BECs and their surrounding thermal cloud have been well studied using time-of-flight absorption imaging. The in-situ density distribution of a BEC in a harmonic trap is first observed using phase contrast imaging [3]. In later work, the density distribution of BECs in a harmonic potential is also studied and has been well described theoretically using a local density approximation [4]. This work has recently become even more relevant with the creation of exciton-polariton condensates [5,6], which truly have a two-dimensional nature.
However, for ideal bosons the local density approximation cannot be used as the homogeneous system in that case does not show BEC at non-zero temperature. This problem has recently become relevant with the achievement of phBEC [7], which are expected to behave as ideal bosons. Initial theoretical work on the density distribution of phBECs [8], which relies on the local density approximation, should therefore be used with caution.
In this work we study the density distribution of photons in a dye-filled microcavity, both experimentally and theoretically. Experimentally, we exploit the axial symmetry of the system by computing radial averages of the photon gas yielding a good signal-to-noise ratio of the overall signal, especially in the thermal tail of the distribution. This allows us to accurately study the tail and compare the behavior of the tail below and above the phase transition. Theoretically, we calculate the radial density profiles by carrying out an exact summation of harmonic oscillator states weighed by the appropriate Bose-Einstein distribution function. We find good correspondence between theory and experiment. We also find, both experimentally and theoretically, that below the phase transition the thermal distribution has the expected Gaussian form, but that above the transition the thermal tail becomes strongly non-Gaussian. * These authors contributed equally to this work † Corresponding author: D.vanOosten@uu.nl

II. EXPERIMENTAL SETUP
The optical design of our experimental setup is based on the work of the Klaers et al. [7,9]. The core of our experimental setup is a cavity consisting of two ultrahigh reflecting, spherical mirrors with a radius of curvature of 1 m with a separation on the order of 1 µm, between which a droplet of a Rhodamine 6G dye solution is placed. The dye solution is pumped by laser pulses with a duration of 500 ns and a wavelength of 532 nm, which are created using a CW laser and acousto-optic modulators (AOMs). Using three AOMs in series gives us an extinction ratio of 5 × 10 4 , ensuring a high on/off contrast of the pump pulse.
The density distribution of photons escaping through one of the cavity mirrors is imaged. The photons first pass through a lens with a diameter of 25.4 mm placed at the focal distance f 1 = 75 mm away from the cavity. The consequences of the effective aperture of this collecting lens will be discussed later. After passing through the lens, a second lens with a focal length f 2 = 200 mm is used to created an image plane with a magnification M 1 = f 2 /f 1 = 2.67. Figure 1 shows typical images obtained below and above the BEC phase transition acquired using an RGB color camera in the image plane. During the experiments the RGB camera is removed and the image plane is imaged on a high resolution camera using a second set of lenses, yielding an additional magnification M 2 = 6. The camera and the AOMs are synchronized such that we can image the distribution of photons in the cavity during each individual pump pulse. To ensure reproducibility of our results, we perform our experiment in runs during which thousands of images are taken without user intervention.

A. Automation
The heart of the automation system of the experiment is the Beaglebone Black (BBB), a low-cost development platform built around the ARM Cortex-A8 processor (Sitara), running a Linux operating system. The BBB is particularly suitable for real-time applications, because the Sitara processor includes two so-called programmable real-time units (PRUs), which share memory with the ARM core. We use one such PRU to address a serial ten bit digital-to-analog converter (DAC). The voltage from the DAC is used to control the RF power of one of the AOMs, allowing us to change the pump pulse power during the experimental run. Furthermore, the PRU is programmed to provide trigger pulses for the AOMs and the camera. The PRU is instructed to trigger a shot of the experiment by sending the appropriate commands over shared memory using a Python script running on the ARM core.
The camera is an Andor Zyla 5.5 sCMOS, the data of which are read out by a PC over a USB3 bus. At the start of an experimental run the camera is set up for external triggering and a Python script runs that takes the appropriate number of images. Afterwards, we start a Python script on the BBB that sends the corresponding shot-requests to the PRU. Synchronization is maintained by the fact that the camera is externally triggered by the PRU.
An experimental run typically consist of 50 repetitions of the same sequence containing typically 60 shots each with a different pump power. Additionally, for each shot a background image is taken. During a run we take 8 shots per second, which including the background images leads to a frame rate of 16 frames per second (fps). As a single full frame image of the camera contains approximately 11 MBytes, this corresponds to a data rate of 88 MBytes/s. The total amount of data produced in one run is 66 GBytes. To handle the data flow, we temporarily store the data of 10 experimental runs on a solid-statedrive (SSD), with a capacity of 1 Tbytes. At the end of the day we copy the data from the SSD to other storage media during the following night.
As the experiment runs completely without user intervention during an experimental run, the lab is vacated during runs, which provides additional stability and reproducibility contributing to the quality of our data.

III. THEORY
In the grand canonical ensemble, the total number of identical, ideal bosons in a trap with a set of discrete energy levels E n is given by where f BE denotes the well-known Bose-Einstein distribution function, g n the degeneracy of the n-th energy level, µ the chemical potential, k B the Boltzmann constant, and T the temperature. In our dye-filled microcavity, the energy levels are those of the isotropic twodimensional quantum harmonic oscillator. The principal quantum number is n = n x + n y , with n x,y ∈ N 0 and E nxny = Ω(n x + n y + 1), where Ω denotes the harmonic oscillator frequency. We assume that neither the cavity nor the medium inside the cavity exhibits birefringence, i.e. that the harmonic oscillator frequency does not depend on the polarization. Summation over the polarization thus yields a factor 2.
To determine the density of photons inside the cavity, we substitute the probability distribution of each state of (n x , n y ) into the above summation, i.e.
ρ(x, y) = 2 nx,ny where the functions ψ n (x) are the well-known onedimensional harmonic oscillator wave functions. Because the expression contains a double summation over a potentially large number of states, it is unfeasible to use in a least-square fitting procedure. For that reason, we write the density as where g n = 2(n + 1) denotes the degeneracy and, As the above expression for ρ n (x, y) does not depend on the parameters µ and T , we can precompute the summation over n x and carry out the summation over n during the fitting procedure. An added benefit of this method is that in the isotropic system at hand, ρ n can only depend on the radial coordinate r, i.e. it must be radially symmetric. Thus we can compute ρ n (r) by computing the angular average of ρ n (x, y) and use where we have introduced the detection efficiency η(E n ) that depends on the photon energy, due to for instance the quantum efficiency of the camera used in the experiment.

A. Implementation
Even though the eigenfunctions of the one-dimensional quantum harmonic oscillator are analytically known, their computation presents a numerical problem. For quantum numbers from n ≈ 150 and upwards the computation involves intermediate values that exceed the maximum value of the double format. For this reason, we calculate the harmonic oscillator wave functions using the Numerov method [10]. In order to reach a sufficient accuracy for the large quantum numbers, we use a grid with a resolution of ∆x = 2 × 10 −3 l HO , where l HO denotes the harmonic oscillator length. We are left with an array ψ nm , where n is the quantum number and m labels the position according to x m = −50l HO + m∆x. We subsequently take the absolute value squared of the numbers in this array, apply a low-pass filter along the m-direction and sub-sample the array in the m-direction to account for the finite optical resolution. Using the resulting one-dimensional density distributions, we tabulate the functions ρ n (r m ), where r m = m∆x. To use the tabulated data in the fit procedure, we create interpolating functions that allow us to scale the harmonic oscillator length to the experimental value.
Examples of the resulting functions are shown in Fig. 2. For n = 1 we find the familiar Gaussian form. For larger principle quantum numbers the density distribution rapidly oscillates, as expected. Due to the filtering, the fast oscillations are suppressed. For principle quantum numbers in the order of 100 or 1000, the filtered radial distribution behaves more and more like a step function, with only a fast oscillation around the classical turning point r ctp = √ 2n l HO . Therefore, for principle quantum numbers > 2000, we no longer perform the numerical calculation, but rather use an appropriately normalized step function.

A. Experimental profiles
To quantitatively analyse the results, we exploit the fact that the photon distributions are radially symmetric, as the microcavity is isotropic, allowing us to perform a radial average. We first determine the center of the microcavity by selecting a run with a visible condensate and locate the center pixel of the condensate. Performing the procedure for several randomly selected images, we notice that we always find the same center pixel, showing the stability of the setup. From the center pixel we average each image radially outwards, allowing us to increase the overall signal-to-noise ratio, but especially in the tail of the distribution that is otherwise very noise.
In Fig. 3 we plot examples of experimentally obtained radial profiles, below and above the critical photon number N c . The two profiles overlap except close to the center, where the orange line has an additional peak. The additional peak is caused by the macroscopic occupation of the ground state, i.e. the phBEC. The overlap between the profiles for larger distances, is a direct consequence of Bose-enhancement; the thermal cloud is saturated and thus cannot contain more photons. Additional photons must occupy the ground state.
One can also observe that the profiles deviate from each other close to, but outside the condensate. This is a consequence of the fact that we are dealing with a finite system and that we are thus not formally in the thermodynamic limit.
Lastly, we observe the decrease of the thermal tail of the distribution over more than two orders of magnitude. It is crucial to use a collecting lens with a sufficiently large We can detect the decrease of intensity in the thermal tail over more than two orders of magnitude. The insets show the (false color) camera images, from which these radial profiles are computed.
numerical aperture, as illustrated in Fig. 4. The numerical aperture of the collecting lens is reduced by inserting an aperture in front of the lens. Reducing the aperture drastically reduces the thermal tail, which has a large impact on the analysis. First, the condensate fraction is over-estimated, as most of the thermal photons are not recorded. Second, as the slope of the thermal tail is a measure of the temperature, a smaller aperture produces a steeper slope and thus predicts a lower temperature. As seen from the figure, the difference between an aperture of 19.0 mm and 25.4 mm is minor, which suggests that a 25.4 mm aperture is sufficient.

B. Fitting
The theoretical radial profiles depend on several parameters. First of all, we need to know the harmonic oscillator length, which one can obtain from the cutoff wavelength λ cutoff , using [9] l HO = 1 4π where n 0 denotes the refractive index of the solvent, q the longitudinal mode number of the cavity mode and R the radius of curvature of the mirrors. To relate the theoretical profile to the experimental profile, the detection efficiency needs to be determined which takes into  account the quantum efficiency of the camera, and the reflection coefficients of the optical elements in the experimental setup. To determine the efficiency, we use the fact that at the phase-transition the total number of photons is given by Ref. [9]: where Ω can be expressed in the experimental parameters as Using an iterative method, we find the detection efficiency by dividing the value of N c by the total experimental signal N exp , which we obtain by numerically integrating the experimentally obtained radial profile. However, we first need to obtain the temperature of the photon gas at the critical point. We therefore fit the radial profile using the chemical potential µ and the temperature T , both in units of Ω, as fit parameters. For the fit, we temporarily use the detection efficiency as a free parameter. With the fitted temperature, we can then recalculate the detection efficiency. The fitted and recalculated detection efficiencies agree reasonably well. We repeat the procedure for 50 images at the phase-transition of Bose-Einstein condensation, and average the recalculated detection efficiencies. The average detection efficiency is used and kept fixed while we carry out the fits to the entire data set. Figure 4 shows examples of fits as black dashed lines to the radial profiles. One can note the striking agreement between the fit and the large aperture data over many orders of magnitude in signal, both for the thermal cloud and the Bose-Einstein condensate. The fit parameters we find for the different apertures are given in Table I. We notice that the fitted temperature, when converted to SI units, is not equal to room temperature as we would have expected. This could be due to the aperture of the collecting lens. As mentioned in Sec. II, we image the cloud using a lens with a finite aperture. Limiting the aperture suppresses the thermal tail, but more importantly increases the slope. One can see that this plays a role, as the temperature for the smaller apertures deviates further from room temperature.

V. CONCLUSION
In our system, a Bose-Einstein condensate of photons in a dye-filled microcavity, the signal-to-noise ratio is improved by radially averaging the data. This is especially true for the tails of the thermal distributions; a crucial step as it is the tail that gives a reliable measurement of the temperature of the system. Furthermore, we show that for the same reason, it is crucial to have a large numerical aperture in the imaging system. Finally, we present a theoretical model that we use to fit our data and demonstrate a good correspondence between this model and our data.