Exact Diagonalisation of Photon Bose-Einstein Condensates with Thermo-Optic Interaction

Although photon Bose-Einstein condensates have already been used for studying many interesting effects, the precise role of the photon-photon interaction is not fully clarified up to now. In view of this, it is advantageous that these systems allow measuring both the intensity of the light leaking out of the cavity and its spectrum at the same time. Therefore, the photon-photon interaction strength can be determined once via analysing the condensate broadening and once via examining the interaction-induced modifications of the cavity modes. As the former method depends crucially on the concrete shape of the trapping potential and the spatial resolution of the used camera, interferometric methods promise more precise measurements. To this end, the present paper works out the impact of the photon-photon interaction upon the cavity modes. A quantum mechanical description of the photon-photon interaction, including the thermal cloud, builds the theoretical backbone of the method. An exact diagonalisation approach introduced here exposes how the effective photon-photon interaction modifies both the spectrum and the width of the photon gas. A comparison with a variational approach based on the Gross-Pitaevskii equation quantifies the contribution of the thermal cloud in the respective applications.


Introduction
Observing photon-photon interactions is a demanding task. Classical electrodynamics considers light as a linear phenomenon, such that two crossing light beams only interfere with each other, but do not scatter. However, the advent of quantum electrodynamics changed this view, since in this theory electromagnetic fields can polarise the vacuum.
Halpern was the first to express the idea of light-by-light scattering within the Dirac theory of electrons and positrons [1]. Afterwards, references [2] and [3] formalised the idea already to the modern picture, as figure 1 a) portrays it. Here, two photons interact by polarising the vacuum via the production of a virtual electron-positron pair. The pair then recombines back into two photons with different wave vectors. Subsequently, reference [4] works out the corresponding modification of the Maxwell equations in vacuum, which is nowadays referred to as the Euler-Heisenberg Lagrangian. However, the authors of reference [4] point out, that this only works with photons in the Röntgen-or gamma-ray regime, since the energy needs to be high enough for supporting the electron-positron pair creation. The later works [5,6] support these findings by applying the more modern S-matrix apparatus. Since this process appears in fourthorder perturbation theory, the cross-section for this kind of photon-photon interaction is proportional to the fourth power of Sommerfeld's constant. Therefore, only particle accelerators allow access to this kind of processes. In 2017 the ATLAS experiment at the Large Hadron Collider was able to observe light-by-light scattering in vacuum for the first time [7]. Embedding photons in non-linear materials increases the photon-photon interaction significantly. One prominent example for such a non-linear process is the Kerr effect a) b) Figure 2. Photon-photon interaction in photon BEC setups. a) Mechanism of thermooptic interaction. The photon density |ψ| 2 heats the medium to a temperature ∆T . This shifts the refractive index by ∆n, which effectively introduces a photon-photon interaction. b) Timescales of thermo-optic interaction. The timescales of the photon population and the external pump pulse are the fastest timescales, whereas the scale of the temperature diffusion is the slowest one.
[8], see figure 1 b). Here, the refractive index of the material changes with the photon density, leading to effects like the self-focusing of a light beam or the existence of optical solitons. Such experiments are performed in effectively two-dimensional setups. This is remnant from light propagation in, e.g., glass cylinders, where the coordinate along the optical axis, due to the paraxial approximation, acts as the time coordinate and the remaining two dimensions as true spatial dimensions [9,10]. A second possibility relies on confining the light together with the non-linear medium inside a cavity. As the cavity mirrors impose Dirichlet boundary conditions upon the light field, a standing wave emerges along the optical axis, freezing out the motion along this direction. However, the used photons are usually prepared in a dissipative state, as light leaks out of the apparatus and has to be reinjected by an external light source for achieving a steady state. However, filling a microcavity with a dye solution offers the possibility to create light, the steady-state spectral distribution of which resembles thermal equilibrium [11,12]. In these systems, also a small effective photon-photon interaction emerges, which does not stem from a Kerr effect but from a thermal non-linearity, the mechanism of which figure 2 a) sketches. This originates from a heating of the dye molecules, which do not emit back all absorbed photons into the photon gas, but instead convert them into vibronic excitations. As a consequence, the dye solution heats up, changing its refractive index. As the heat diffuses through the experimental setup, the resulting effective thermo-optic photon-photon interaction is non-local in space and retarded in time, which represents the main difference to the photon-photon interactions sketched in figure 1. Since the temporal retardation is large compared to the remaining experimental timescales, see figure 2 b), the thermo-optic photon-photon interaction effectively corresponds to a potential, which increases linearly in time [13,14]. As the thermo-optic interaction turns out to be weak, only the photon BEC (phBEC) regime allows determining its strength by measuring the condensate width [11]. However, due to the smallness of the effective photon-photon interaction, this method lacks of precision. Instead, spectroscopic measurements promise an increased precision by measuring the small energy shifts, which the thermo-optic photon-photon interaction induces during a single pump pulse. The basis of the theoretical description of this thermo-optic interaction is a detailed model for the diffusion of the temperature produced in a single pump pulse. As the temperature slowly builds up, an adiabatic approximation of the corresponding quantum mechanical modelling grants access to the instantaneous eigenvalues. This allows determining the effective photon-photon interaction strength from the energy difference between two eigenmodes. Whilst reference [14] works out the corresponding basic theory and applies first-order perturbation theory for deriving some initial results, this paper goes beyond and works out in detail both the energy shifts and the condensate broadening by applying exact diagonalisation (ED) [15]. The ED of the underlying second-quantised Hamiltonian represents a powerful method in case of increased effective photon-photon interactions. For instance, reference [16] predicts this situation at the dimensional crossover from 2D to 1D, where a significant trap anisotropy enhances the photon-photon interaction via an increased photon density. Experimentally, microstructured mirrors provide such anisotropic traps [17,18,19]. However, the thermal cloud leads to a stronger increase of the condensate width and to a possible overestimation of the effective photon-photon interaction, provided the theoretical modelling only considers the phBEC ground state. To avoid this, the ED method combines the investigation of the thermodynamic behaviour of the noninteracting phBEC [20] with the influence of the thermo-optic interaction upon the phBEC ground state at the dimensional crossover [16]. The present paper is structured as follows: Section 2 starts with a short summary of the underlying second-quantised Hamiltonian and uses ED for working out the new eigenmodes of the harmonically trapped photon gas. Moreover, the corresponding condensate width is compared to former results and special attention lies on the temperature dependence of both the eigenenergies and the width of the photon gas. Subsequently, section 3 deals with the dimensional crossover and presents the corresponding results. Section 4 summarises the main findings of the paper.

Exact Diagonalisation of Harmonic Potential
This section provides a concise introduction to the model used for describing the thermooptic interaction, which is based on reference [13] and worked out more rigorously in the preceding paper [14]. Afterwards, ED serves as a method for benchmarking both the perturbative results on the energy spectrum derived in [14] and the variational approach performed in [13]. The last part of this section goes beyond these findings and investigates the impact of finite temperatures on both the energy spectrum and the cloud width.

Model
As the phBEC itself varies on timescales, which are much faster than the produced temperature, see figure 2 b), reference [14] breaks down the quantum mechanical description of the thermo-optic interaction to an effective potential, which increases linearly in time. Hence, the second-quantised Hamiltonian, bears an adiabatic time dependency introduced by the temporal retardation of the thermo-optic photon-photon interaction. The bosonic creation and annihilation operatorsâ † n (t),â n (t) belong to the instantaneous eigenmodes of the Hamiltonian (1). The Hamiltonian matrix in (1) has the form with the effective time-dependent thermo-optic interaction strength g(t) increasing linearly in time [14]. At the beginning of the experiment, i.e., at t = 0, no interaction is present, yielding the initial value g(0) = 0. The corresponding eigenvalue problem leads to the eigenenergies E n (0) and the eigenfunctions ψ n (x). Furthermore, the non-diagonal matrix contains the information about the thermo-optic interaction. It uses the mode occupation in the form of a Bose-Einstein distribution at the beginning of the experiment, which stems from the temporal retardation of the thermo-optic interaction: Here β = 1/(k B T ) denotes the inverse temperature and µ(0) represents the initial chemical potential, which is fixed by the conserved total particle number N = l a † l (t)a l (t) .

Harmonic Potential
This paper specialises to a harmonic trapping potential of the form with the effective photon mass m and the trapping frequency Ω x in x-direction. The trap-aspect ratio λ = Ω y /Ω x determines the trapping frequency Ω y in y-direction. Therefore, the eigenenergies take the form Furthermore, the Gauß-Hermite functions determine the eigenfunctions of the second-quantised Hamiltonian (1) and the corresponding Hamiltonian matrix (2), where l i = /(mΩ i ), i = x, y, denote the oscillator lengths in both directions and H n are the Hermite polynomials. This paper numerically determines the new eigenmodes of the Hamiltonian matrix (2) with the potential (5) by applying exact diagonalisation. As this method necessitates using a finite number of Gauß-Hermite eigenmodes (7), the sixty lowest energy subspaces are included. This ensures in the Bose-Einstein condensed regime, that the relative occupation of higher excited states is negligible. With this finite number of modes, the interaction matrix (3) is constructed for a given total particle number N and thermodynamic temperature T . Hence, the presented method not only verifies previously published calculations relying on different methods, but also reveals possible deviations from the latter. To be specific, this section focuses on the isotropic case, i.e., λ = 1.

Spectrum
This paragraph presents the results for the adiabatic time-dependent energy eigenvalues E n (t) of the Hamilton matrix (2) for the potential (5). Figure 3 a) compares the eigenenergies obtained by the ED method with the results of first-order perturbation theory from the preceding paper [14]. It plots the eigenenergies as a function of the dimensionless interaction strengthg(t) = mg(t)/ 2 , the magnitude of which is determined by the duration of the external pump beam. The results of both approaches agree well and, thus, the same observation holds for the corresponding energy differences shown in figure 3 b). However, for larger interaction strengths, a slight deviation of up to 10 −3 between the ED results and the first-order perturbation theory [14] becomes visible in figure 3 c). Hence, the perturbation theory derived in [14] yields a good approximation in this regime of parameters. Therefore, the corresponding formulas, which analytically express the photon-photon interaction strength via the energy differences, are precise enough for being applied to spectroscopic measurements.

Condensate Width
The new eigenstates also show a broadening of the photon gas. The unitary matrix U n,m (t), which rotates into the instantaneous eigenbasis of the Hamiltonian matrix (2), defines creation and annihilation operatorsâ n (t),â † n (t) at a given time t: . a) First few eigenenergies of the second-quantised Hamiltonian matrix (2). The dots stand for the results of the ED, whereas the lines represent the firstorder perturbative results from reference [14]. b) Corresponding energy differences ∆E n,n = E n − E n . c) Relative deviation of energy differences once calculated by the perturbation theory and once by the ED method. In all plots, the indices 2 ± = [(20) ± (02)]/ √ 2 denote the mode hybridisation due to the interaction.
With these at hand, the total photon density, is expressed in terms of the instantaneous eigenstate occupations, and the corresponding eigenstate density, Therefore, the total width σ x of the photon gas is defined by the FWHM Here, the second moment of the density distribution with the adiabatic time-dependent photon density (9), takes the form with the adiabatic time-dependent matrix elements Older calculations and measurements [11,13] consider the deep condensate limit N 0 (t) ≈ N . In this case, (12) reduces to the width of the ground state, i.e., σ x (t) ≈ 2 x 2 0 (t). Appendix A reproduces the corresponding variational calculation of the photon-condensate width, which was first given in reference [13]. Figure 4 compares the results for the condensate width, calculated via ED and via the variational approach from Appendix A. As the results agree well, the ED method developed here is found to be able to reliably reproduce the former theoretical and experimental results. However, a closer look at the inset in figure 4 for larger effective interaction strengths reveals that the results tend to deviate for larger effective photon-photon interaction strengths. Since the width calculated from the ED is larger than the corresponding width from the variational approach, this deviation stems from an increased coupling to the thermal cloud. Thus, the larger the photon-photon interaction strength is, the ground-state-only description is less accurate.

Thermal Cloud
Previous investigations of the effective photon-photon interaction only consider the zerotemperature limit. This section, however, deals with the effects of finite temperature, which the Hamiltonian matrix (2) naturally includes. At the beginning of the experiment, the thermo-optic interaction is not present, and the phBEC behaves as an ideal Bose gas in 2D [20]. As the photon-photon interaction strength increases along the duration of an experiment, this section aims at describing the influence of the thermal cloud upon the phBEC properties. First, the focus lies on the temperature dependence of the energy eigenvalues, which figure 5 a) shows. In the thermal regime, the energies undergo a small shift due to the interaction, whereas in the condensed regime, i.e., for T < T c , a more significant shift builds up. This is due to the small ground-state occupation in the thermal phase.
Since the coupling to the ground state is always the largest, significant energy shifts only occur in the condensed regime, as figure 5 a) clearly indicates. Experiments can validate this finding by measuring the differences between the energy levels, which are plotted in figure 5 b). The time-dependent shifts of the energy eigenstates E l (t) also affect the Bose-Einstein distribution (10) of the photons. Figure 6 compares the photon number distribution (10) at the beginning (dots) and at the end of the experiment (crosses). The horizontal shift of the occupation numbers at the end of the experiment indicates interaction-induced shifts of the eigenenergies. The occupation numbers themselves, however, undergo only minor modifications due to the smallness of the interaction. Finally, the ED also allows calculating the width (12) of the total photon gas, which figure 7 depicts. In the zero-temperature limit, the width of the photon gas approaches the ground-state width. At the beginning of the experiment, where no interaction is present, this width corresponds to the oscillator length, whereas during an experimental run the ground-state width increases and so the width of the total gas increases with the photon-photon interaction. In the thermal regime, however, the interaction does not lead to a noticeable change of the photon-gas width. This finding coincides with the observations above, showing only very small shifts in the energy levels in the thermal regime. The saturation of the width, which occurs for large temperatures, is a remnant from the finite number of considered eigenmodes.

Dimensional Crossover
With the development of new experimental techniques for micro-structuring the cavity mirrors [17,18,19], also large trap anisotropies for phBECs can be manufactured. Therefore, these techniques allow realising a dimensional crossover from 2D to 1D. In this case, the ED method turns out to be essential, since reference [16] predicts a considerable enhancement of the effective photon-photon interaction strength.
The beginning of the section maps the 2D gas in a strong anisotropic trap onto a 1D gas, which determines the effective 1D interaction strength. The latter turns out to depend on the trap-aspect ratio λ. With this at hand, the effective 1D energy spectrum is calculated on the one hand numerically and on the other hand approximated analytically. The latter grants access to an approximation formula determining the effective photon-photon interaction strength from the differences of the eigenenergies in the quasi-1D limit. Finally, the calculation of the condensate width in the quasi-1D limit shows the significance of the ED method, as the width turns out to deviate from the variational approach presented in [13], which is recovered in Appendix A. Note that in this chapter the final photon-photon interaction strengthg(t end ) is assumed to be an order of magnitude smaller than in the previous section, such thatg(t end ) = 10 −5 . This adjustment takes into account the assumed shorter condensate lifetime due to the increased photon intensity, resulting from the strong anisotropy. Nevertheless, this section shows the strong increase of the effective photon-photon interaction strength due to the dimensional crossover.

Effective 1D Interaction Strength
Reference [20] demonstrates that for trap-aspect ratios larger than a critical value λ 1D the squeezed direction freezes out in its ground state and that the systems behaves effectively quasi one-dimensional. In addition, the density in the frozen out direction alters the bare 2D interaction strengthg [16]. The resulting effective quasi onedimensional interaction strengthg 1D (λ) turns out to be a function of the trap-aspect ratio λ. The current paper uses the same reasoning for a single pump pulse, which reference [16] applies to the steady state after several pump pulses. According to definition (5) an increasing trap-aspect ratio squeezes the y-direction and, so, in the quasi-1D case the photons only populate the ground mode in this direction. This suggests the factorisation ψ n (x) = χ nx (x)ϕ 0 (y; λ) (16) for the Gauß-Hermite eigenmodes of the potential (5). Therefore, the interaction matrix (3) reduces to where the weight incorporates the influence of the squeezed direction and, thus, depends on the trapaspect ratio λ. In the real 1D case, however, the structure of the Hamiltonian matrix (2) remains the same; only the multi-indices turn into single indices. Hence, the identificationg of an effective 1D interaction strength allows a mapping of the matrix (3) onto the real 1D case. The explicit calculation of the weight (18) for the harmonic potential (5) leads tog The effective 1D interaction (20) agrees with the corresponding result for the contact Kerr interaction obtained in reference [16], but not with the similar one for the thermooptic interaction, although the current work considers only the latter. The reason is the used short-time approximation in (1). Due to this, the matrix elements (3) effectively resemble the corresponding ones for a contact interaction. Reference [16] investigates the steady state after several pump pulses, such that the temperature diffusion has a significant impact on the behaviour. In the current work, however, this does not happen due to the short-time approximation.

Energy Spectrum
As the effective 1D interaction (20) increases along the dimensional crossover, the deviations observed in the previous sections become significantly large and both the perturbative and the variational calculation are no longer valid in the case λ 1. Figure 8 illustrates this finding. In figure 8 a) the new eigenenergies are plotted with respect to the non-interacting ground-state energy E 0 = Ω x (1 + λ 2 )/2 at the dimensional crossover for a fixed interaction strengthg(t), revealing an increasing effective 1D interaction strengthg 1D (t; λ). Due to the enhanced effective interaction strength, increasing the trap anisotropy shifts the eigenenergies to larger values. Here too, the ground-state energy undergoes the largest change, which is clearly non-linear. Figure 8 b) shows that the energy differences likewise increase with the trap-aspect ratio pointing towards the predicted increased photon-photon interaction at the dimensional crossover, cf., reference [16]. Thus, the ED method also allows determining the effective photon-photon interaction strength along the dimensional crossover in the same way by spectroscopic measurements as outlined above. For instance, one approach relies on interfering different cavity modes and observing the resulting beating signal. Appendix B derives an analytic approximation for the eigenvalues, which are numerically calculated by the ED in the deep condensate limit. This approach considers the first four eigenmodes and diagonalises the corresponding 1D Hamiltonian matrix. This reveals the non-linear character of the eigenvalues, shown in figure 8 a) to stem from avoided crossings with the next-higher eigenmodes with the same symmetry. Due to the same reason, this method only approximates the ground and first excited state well, whereas the ED and the analytical approximation for the higher excited states agree less for larger anisotropies. This can be avoided by taking into account even higher excited states, which is necessary for going to larger trap anisotropies respectively larger interaction strength, but makes the analytical diagonalisation more difficult. However, if the effective photon-photon interaction strength is small enough for a linearisation, the analytical approximation for the ground and first-excited state represents a result for determining it according tõ Figure 9. Dimensional crossover. Both plots are for a total particle number N = 10 4 and a maximal interaction strengthg(t end ) = 10 −5 . a) Comparison of condensate width (12) at the end of the experiment for different trap-aspect ratios at fixed interaction strength. Here, the condensate fraction amounts to N 0 /N ≈ 0.9. Dots: results from ED; lines: variational approach. b) Total width of the photon gas at the beginning of the experiment (blue) and at the end of the experiment (orange), as well as the width of the ground state with interaction (green) for trap-aspect ratio λ = 10. The black line indicates the oscillator length, which is also the width of the ground state without interaction. The drawn lines are a guide to the eye. Figure 9 a) compares the widths from the ED with the results from the variational approach for a fixed interaction strength. Whereas in the squeezed y-direction the results from both methods agree quite well at the whole crossover, the results coincide only for small trap-aspect ratios in the un-squeezed x-direction. At larger trap-aspect ratios, the results deviate up to several percent, which is caused by the coupling to the thermal cloud. Since the effective photon-photon interaction increases along the dimensional crossover, in the sense of equation (20), the thermal cloud becomes more and more important in this scenario. Consequently, the ED method plays a crucial role for measuring the effective photon-photon interaction at the dimensional crossover. Also, for larger trap anisotropies the total width of the photon gas can be calculated from (12), which figure 9 b) pictures. Here, the blue line depicts again the total width of the photon gas at the beginning of the experiment for a certain thermodynamic temperature T . As the experiment runs, the effective photon-photon interaction increases and, thus, also the photon gas width grows during the experiment. Alas, only in the condensed regime this growth is significant and due to the trap anisotropy enhanced in comparison to the situation in the isotropic trap, see figure 7. In the thermal regime, however, still no impact of the effective photon-photon interaction is observable.

Summary
The findings presented in this work are crucial for precisely quantifying the effective photon-photon interaction strength in current and future phBEC experiments. The present paper analyses in detail the theory of reference [14] by using exact diagonalisation for a harmonically trapped photon gas. The aim consists in determining the impact of the thermo-optic interaction upon the cavity modes. The accurate prediction of the shifts of the eigenenergies allows for a precise interferometric measurement of the emerging photon-photon interaction. The method reproduces formerly derived results and extends these systematically by taking the thermal cloud into account. The influence of thermal cloud turns out to be crucial for understanding the dimensional crossover and prevents a possible overestimation of the effective photon-photon interaction strength. Where possible, analytic estimates provide aid for determining the effective photon-photon interaction strength from measurements.

Appendix B. Analytical Approximation in the 1D Case
The aim of this appendix is to work out an analytical approximation for the quasi-1D case of a harmonic potential. Since the trap-aspect ratio determines the interaction strengthg 1D (λ), the latter can reach comparatively large values. Therefore, the first four eigenstates are taken into account and the resulting 1D Hamiltonian matrix is diagonalised. Moreover, the deep condensate limit N 0 ≈ N is assumed. With the harmonic-oscillator eigenfunctions, the Hamiltonian matrix reads where the shift due to the unperturbed ground-state energy has been dropped. The Hamiltonian matrix possesses the eigenvalues