text:
# Theoretical study of finite temperature spectroscopy in van der Waals clusters. I. Probing phase changes in CaArn

## I Introduction

Large van der Waals clusters offer a convenient medium to study solvation at the molecular level, and also to undertake chemical reactions with well defined thermodynamical conditions.cicr1 ; cicr2 ; cicr3  Reactivity of chromophore-doped rare-gas clusters has been the subject of intense experimental effort.baumfalk00 ; slavicek ; baumfalk01 ; gaveau00 ; liu01 ; briant01  Photodissociation of HBr on argon clusters has been investigated by Baumfalk et al.,baumfalk00 ; slavicek ; baumfalk01  including a theoretical study of the HBr–cluster interaction.slavicek  Chemiluminescence in the bimolecular reaction Ba$`+`$N2O$``$BaO$`{}_{}{}^{}+`$N2 was the subject of a work by Gaveau et al.,gaveau00  and photodetachment in IHI-Arn clusters was used to probe transition state spectroscopy of the I$`+`$HI reaction by Liu and coworkers.liu01  One may also cite the photoinduced reaction dynamics between Ca and HBr towards CaBr$`{}_{}{}^{}+`$H at the surface of a large argon cluster carried out by Briant et al.briant01  Very recently, the reaction between N2O and Li2 mediated by Arn clusters has been investigated by Breckenridge and coworkers.breck02

A practical understanding of the physical and chemical properties of the deposited species relies on a good characterization of their spectroscopy. In turn, spectroscopy gives great insight into the structure of the van der Waals complex, especially in concern with the location of the chromophore.depujo92 ; depujo94 ; epjd  This idea was used to propose a way of detecting structural transformations inside the cluster.hahn ; even  Not long after its was discovered in simulations,berry  the so-called dynamical coexistence phenomenon between solidlike and liquidlike phases received some experimental evidence by Hahn and Whetten who reported changes in the excitation spectrum of benzene-(argon)n clusters with increasing temperature.hahn  Similar experimental results were reported for dichloroanthracene-Arn clusters by Even et al.even  More recently, Curotto and coworkerscurotto1 ; curotto2 ; curotto3  theoretically investigated the infrared spectroscopy in Arn-HF clusters. They observed significant variations in the red shift associated with the two lowest vibrational states of the HF molecule, which they could correlate with the thermodynamical state of the cluster.curotto1 ; curotto2 ; curotto3  Moseler et al.moseler  found that the photoabsorption spectra of Na8 displays some features characteristic of isomerization and structural fluctuations. Lastly, simulations on charged rare-gas clusters have investigated the influence of temperature on the absorption spectrum.grigorov ; galindez  So far, no direct evidence for isomerization in such rare-gas clusters has been found yet.

While spectroscopy has many attractive features for experimental studies in large clusters, several difficulties pave the way for efficient theoretical approaches. Firstly, because of the numerous degrees of freedom, a fully quantum calculation of the absorption spectrum by exact methods (wavepacket or diagonalization in a basis) is not feasible beyond a few atoms and one is usually lead to simplify the problem. A common approximation is the Condon hypothesis, which assumes vertical transitions. Secondly, the equilibrium properties at finite temperature require some sampling of the ground state potential energy surface (PES). This sampling can be hard to get because of the topography of the multidimensional PES (the so-called energy landscape). When the energy landscape is complicated or has multiple funnels, the relaxation to equilibrium may be particularly slow with a glassy-like behavior.walesacp  Broken ergodicity is then likely to occur in conventional molecular dynamics (MD) or Monte Carlo (MC) simulations, and also in quantum simulations. Fortunately some remedies to this problem have recently been suggested, including multicanonical ensemble samplingberg  or parallel tempering.ptmc  Thirdly, most experiments involving chromophore-doped rare-gas clusters take place at low temperatures, at which the classical approximation may be questioned, especially for the lighter elements neon and helium. Path-integral Monte Carlo or centroid MD should then be used instead of classical techniques, however their use with ergodicity restoring strategies remains numerically demanding for systems containing tens or hundreds degrees of freedom.

We have recently proposed an alternative way of calculating thermodynamical observables at canonical equilibrium, using the superposition approximation of quantum oscillators.jcpquantum  This method, originally developed for thermodynamical observables only, was improved phenomenologically by including anharmonic corrections,jcpquantum ; anhardoye  and extended to generally treat other temperature-dependent observables.prop  Its main interest is to circumvent the troubles of broken ergodicity by considering separately the contribution of each isomer to the partition function. The method also naturally describes quantum delocalization in a straightforward way.

One goal of this paper is to extend the quantum superposition method to the realm of calculating photoabsorption spectra. In addition to the ingredients presently available, one needs to calculate the photoabsorption spectrua of each stable isomer. For this, we use a recent theory by Wadi and Pollakwp  that we have extended to the general case of arbitrary potential surfaces. As will be seen below, this Gaussian theory of absorption succeeds in predicting absorption line shapes, widths and shifts for harmonic polyatomic molecules in good agreement with highly accurate Monte Carlo data.

With the quantum superposition method, we are able to carry out systematic studies for a given set of cluster sizes. We can then address our second goal, and try to relate the features of the absorption spectrum to the knowledge of isomerizations and possible phase changes inside the cluster. More specifically, by looking at particular sizes we hope to gain insight in the relationship between structure, temperature, quantum delocalization, and the shape and position of the spectral lines.

This work follows previous efforts aimed at interpreting experimental absorption spectra on large CaArn clusters.epjd  The spectroscopy of the CaAr diatomic molecule has been investigated both experimentally and with ab initio calculations in Ref. abinit, , and a Diatomic-In-Molecules (DIM) Hamiltonian was constructed to model the $`4s`$ ground- and $`4s4p`$ excited-states of CaArn clusters.epjd  In the present context of finite-temperature spectroscopy, the heavy numerical cost of statistical sampling is balanced with this relatively cheap, albeit accurate Hamiltonian. The paper is organized as follows. In the next section, we briefly recall the quantum superposition method in the harmonic approximation. We also include the derivation of approximate perturbative anharmonic quantum corrections to the partition functions. The Gaussian theory of absorption for harmonic molecules, initially proposed by Wadi and Pollak,wp  will be generalized to include possible Duschinskii rotations in Sec. III. These methods will be tested on simple clusters, in both the classical and quantum regimes, by comparing their results to high statistics Monte Carlo data. Their main application to the variations of the photoabsorption spectrum with temperature in some size-specific CaArn clusters will be presented in Sec. IV. The results will be discussed in view of the isomerization and phase changes phenomena. Some concluding remarks will close the paper in Sec. V.

## II Quantum superposition method

In the superposition approach,hoare79 ; stillinger82 ; wales93  also known as the multiple normal modes model,franke  the total partition function $`Z`$ of the $`N`$atom system in equilibrium at temperature $`T=1/k_B\beta `$ is replaced by a sum over stable minima $`\{\alpha \}`$:

$$Z(\beta )=\underset{\alpha }{}n_\alpha Z_\alpha (\beta ).$$
(1)
Each minimum $`\alpha `$, or inherent structure, has a contribution $`Z_\alpha `$ to $`Z`$, further weighted by the degeneracy factor $`n_\alpha `$, which accounts for permutational isomers and global symmetry. In CaArn clusters ($`N=n+1`$), $`n_\alpha `$ is given by $`2n!/h_\alpha `$, where $`h_\alpha `$ is the order of the point group of structure $`\alpha `$. The problem has now been reduced to estimating the individual functions $`Z_\alpha `$. At the simplest level of approximation, the harmonic vibrational frequencies $`\{\omega _{\alpha i}\}`$ provide a zero-th order expression for $`Z_\alpha `$:

$$Z_\alpha ^{(0)}(\beta )=e^{\beta E_\alpha }\underset{i}{}\frac{2}{\mathrm{sinh}\beta \mathrm{}\omega _{\alpha i}/2},$$
(2)
with $`\mathrm{}`$ the reduced Planck’s constant and $`E_\alpha `$ the potential energy minimum of isomer $`\alpha `$. In the classical regime $`\mathrm{}0`$, the decomposition (1) is formally exact, since each point of the configuration space (except saddle points) belong to one and only one basin $`\alpha `$. In the quantum regime (finite $`\mathrm{}`$), tunnelling is neglected.jcpquantum  In principle, the knowledge of all minima of the energy landscape allows us to calculate any thermodynamical observable, including internal energy and heat capacity, from Eq. (1). However, for all but the smallest atomic systems, exhaustive enumeration of these isomers is not possible because the number of local minima is thought to grow exponentially with the number of degrees of freedom.stillexpgrows  To correct for the uncompleteness of the available set of isomers, a reweighting technique can be used.wales93 ; jcpquantum  The idea is to perform a simulation at a reference temperature $`T_0`$, and to compute the probability $`p_\alpha `$ of observing isomer $`\alpha `$. The partition function $`Z`$ is approximated as a partial sum over weighted isomers:

$$Z(\beta )\underset{\alpha }{}^{}g_\alpha n_\alpha Z_\alpha (\beta ),$$
(3)
with the partial character of the sum indicated by the prime symbol. If we assume that the weights $`\{g_\alpha \}`$ are independent of temperature, the probability $`p_\alpha (T_0)`$ of finding isomer $`\alpha `$ at $`T_0`$ is proportional to $`g_\alpha n_\alpha Z_\alpha (\beta _0)`$. Hence, $`Z`$ can be expressed as

$$Z(\beta )\underset{\alpha }{}^{}p_\alpha (\beta _0)\frac{Z_\alpha (\beta )}{Z_\alpha (\beta _0)}.$$
(4)
A classical simulation can be used to sample the minima and their respective weights.jcpquantum  In the above equation, the denominators $`Z_\alpha (\beta _0)`$ should then be replaced by their classical analogue. Even though the resulting partition function (4) should not depend on the classical or quantum nature of the simulation, some differences can be seen when using semiclassical simulations including quantum corrections.prop  This reweighting procedure has been shown to improve significantly the predictions of the superposition method. In particular, Wales was able to obtain a ’S-bend’ in the kinetic temperature of Ar55,wales93  as seen previously in simulations.labastie

The other important limitation of the method is the harmonic approximation used to calculate $`Z_\alpha (\beta )`$. Several phenomenological anharmonic corrections have been proposed in the litterature for classical systems.anhardoye ; ball  In the case of quantum systems, the optimal normal modes method of Cao and Vothvoth  can be used to include anharmonicity to some extent, but the self-consistency required makes it computationally heavy for large sets of isomers. Building upon a recent effort to include systematic corrections from the energy surface topography via classical perturbation theory,jcpanhar  we can attempt to include anharmonicities in the quantum partition functions $`Z_\alpha `$ in a more rigorous fashion. Since the classical corrections are also numerically demanding, we restrict ourselves to low-order corrective terms. In order to make the reading easier, the details of the perturbation expansion are given in the Appendix. The superposition method not only provides the thermodynamical equilibrium properties, and other quantities can be computed provided they can be calculated for each minimum $`\alpha `$.prop  At thermal equilibrium, the average value $`A`$ of the global observable $`A`$ is the weighted sum over minima:

$$A(T)=\frac{_\alpha ^{}p_\alpha (T_0)Z_\alpha (T)A_\alpha (T)/Z_\alpha (T_0)}{_\alpha ^{}p_\alpha (T_0)Z_\alpha (T)/Z_\alpha (T_0)}.$$
(5)
## III Gaussian theory of absorption

We assume that the cluster is in canonical equilibrium at temperature $`T=1/k_B\beta `$, and that it undergoes an excitation by an infinitely narrow pulse at frequency $`\omega `$ from its ground state potential surface $`V_g(𝐑)`$. The normalized excitation spectrum $`(\omega ,\beta )`$ isyan

$`(\omega ,\beta )=`$
$`{\displaystyle \frac{_i_j\delta (\mathrm{\Delta }E+E_{e_j}E_{g_i}\mathrm{}\omega )e^{\beta E_{g_i}}|e_j|\mu |g_i|^2}{_i_je^{\beta E_{g_i}}|e_j|\mu |g_i|^2}},`$ (6)
with $`|g_i`$ and $`|e_j`$ the eigenstates with eigenenergies $`E_{g_i}`$, $`E_{e_j}`$ of the ground- and excited-state Hamiltonians $`H_g`$ and $`H_e`$, respectively. The transition dipole moment $`\mu _{ij}=|e_j|\mu |g_i|^2`$ will be taken as constant in the usual Condon approximation of vertical excitations. $`\mathrm{\Delta }E`$ is the energy difference between the bottom of the ground- and excited-state surfaces. The absorption intensity can be expressed as the Fourier transform of the imaginary-time correlation function:wp

$$(\omega ,\beta )=\frac{\mathrm{}}{2\pi Z(\beta )}_{\mathrm{}}^{\mathrm{}}\chi (\tau ,\beta )e^{i\tau (\mathrm{\Delta }E\mathrm{}\omega )}𝑑\tau ,$$
(7)
with

$$\chi (\tau ,\beta )=\mathrm{Tr}\left[e^{i\tau H_e}e^{(\beta i\tau )H_g}\right],$$
(8)
and where $`Z(\beta )=\mathrm{Tr}e^{\beta H_g}`$ is the partition function. Following Wadi and Pollak,wp  we now perform a short-time expansion in the correlation function. This expansion should be valid at moderate temperatures and for a large number of degrees of freedom. Indeed, in such cases the combination of the numerous vibrational modes will significantly dephase rapidly, and the correlation function will be nonzero only for short times. Using this approximation, the absorption intensity can be calculated analyticallywp  as a Gaussian function:

$$(\omega ,\beta )=\frac{\mathrm{}}{\overline{\omega }\sqrt{2\pi }}\mathrm{exp}\left[\frac{1}{2}\left(\frac{\omega \omega _0\mathrm{\Delta }\omega }{\overline{\omega }}\right)\right],$$
(9)
with $`\omega _0=E_0/\mathrm{}`$, $`E_0`$ being the vertical transition energy. The line shift $`\mathrm{\Delta }\omega `$ and the width $`\overline{\omega }`$ are given as a function of the two first moments of the shifted Hamiltonians difference $`\mathrm{\Delta }H=H_eH_gE_0`$ bywp

$$\mathrm{}\mathrm{\Delta }\omega =\mathrm{\Delta }H,$$
(10)
and

$$\mathrm{}\overline{\omega }=\mathrm{\Delta }H^2\mathrm{\Delta }H^2.$$
(11)
We have noted here $`A`$ the thermal average of observable $`A`$ over the ground state potential surface $`H_g`$. We now write the Hamiltonians $`H_g`$ and $`H_e`$ in harmonic (normal modes) coordinates $`𝐐`$:

$`H_g`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}p_i^2+{\displaystyle \frac{1}{2}}𝐐^{}𝐖_g𝐐`$ (12)
$`H_e`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}p_i^2+𝐆^{}𝐐+{\displaystyle \frac{1}{2}}𝐐^{}𝐖_e𝐐+E_0.`$ (13)
In these equations, $`𝐖_g`$ is the mass-weighted Hessian matrix of the ground state at the lowest energy minimum $`𝐐_0`$, which is assumed to be located at the origin $`\mathrm{𝟎}`$ of the reference frame. $`𝐖_g`$ can be written as $`𝐖_g=\mathrm{diag}(\omega _{g_i}^2)`$, as it is diagonal by construction. $`𝐖_e=(w_e^{ij})`$ is the second derivatives matrix of the excited state surface at $`𝐐_0`$. Lastly, $`𝐆=H_e(𝐐_0)/𝐐`$ is the excited state gradient at $`𝐐_0`$. The matrix $`𝐖_e`$ is not assumed to be diagonal, and may include the so-called Duschinskii rotations.duschinskii  It is then a simple matter of algebra to obtain the expressions for $`\mathrm{\Delta }\omega `$ and $`\overline{\omega }`$. Using the notation $`\rho _i=\omega _{g_i}\mathrm{tanh}\beta \mathrm{}\omega _{g_i}/2`$, we find

$`\mathrm{\Delta }\omega `$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i}{}}{\displaystyle \frac{\omega _{g_i}^2w_e^{ii}}{\rho _i}}`$ (14)
$`\overline{\omega }^2`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{i}{}}\left({\displaystyle \frac{\omega _{g_i}^2w_e^{ii}}{\rho _i}}\right)^2+{\displaystyle \frac{1}{8}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{(w_e^{ij})^2}{\rho _i\rho _j}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}{\displaystyle \frac{G_i^2}{\mathrm{}\rho _i}}.`$ (15)
Apart from the different notations, the above expressions differ from the results by Wadi and Pollak by the second term in the right-hand side of Eq. (15), which was neglected by these authors.wp  Although Eqn. (14) and (15) are only valid for polyatomic, harmonic systems, anharmonic corrections can be carried out in the very same way as the perturbative expansion mentioned in the previous section and detailed in the Appendix. However, the present expressions are already the outcome of an approximation in the correlation function, hence we did not attempt to include such corrections here. In addition, their estimation would require the computation of higher order derivatives of the excited states energy surfaces, which will be a heavy numerical task, even for for DIM potentials.

We have tested the Gaussian theory on small CaArn clusters. The ground- and excited-state Hamiltonians have been fully described previously.epjd  They are based on the Diatomic-In-Molecules (DIM) formalism,dim  and were parameterized on ab initio calculations.abinit  To give a better account of experimental properties of the CaAr diatomics, we have slightly modified the parameters with respect to our previous works.epjd ; prop  The Ca–Ar ground state pairwise potential has been scaled by 0.85, and the excited state surfaces have been shifted by +305 cm-1 to fit the $`4s4p`$ experimental transition in calcium. As a consequence, some global minima may have changed, as in CaAr37.

In the present work, we focus on the $`4s^24s4p`$ transitions, which have received some experimental attention.epjd  The DIM Hamiltonian provides the energy and numerical derivatives, as well as the electronic transition dipole moment, for any atomic configuration $`𝐑`$. It can thus be used to simulate absolute absorption spectra. The ground state vibrational eigenmodes and eigenfrequencies $`\{\omega _{g_i}\}`$ were obtained from the diagonalization of the Hessian matrix whose elements were computed from their analytical expressions. In both the classical and quantum cases, the absolute spectrum is calculated as the sum over the three possible excited states. For the Gaussian theory, the three gaussian curves are weighted by the square of the transition dipole moment at the equilibrium ground state geometry.

First, we have calculated the classical spectrum from the Gaussian theory in the classical limit $`\mathrm{}0`$, and from classical Metropolis Monte Carlo simulations at low temperature, using $`10^6`$ cycles following $`2\times 10^5`$ equilibration cycles. The spectra of CaAr13 and CaAr54 at $`T=5`$ K are displayed in Fig. 1. For comparison, we have also plotted the predictions of the Gaussian theory in the classical limit $`\mathrm{}0`$. In this limit, the line shift and line width tend to their classical values $`(\mathrm{}\mathrm{\Delta }\omega )_c`$ and $`(\mathrm{}\overline{\omega })_c^2`$ given respectively by

$`(\mathrm{}\mathrm{\Delta }\omega )_c`$ $`=`$ $`{\displaystyle \frac{k_BT}{2}}{\displaystyle \underset{i}{}}{\displaystyle \frac{\omega _{g_i}^2w_e^{ii}}{\omega _{g_i}^2}},`$ (16)
$`(\mathrm{}\overline{\omega })_c^2`$ $`=`$ $`{\displaystyle \frac{(k_BT)^2}{2}}{\displaystyle \underset{i}{}}\left({\displaystyle \frac{\omega _{g_i}^2w_e^{ii}}{\omega _{g_i}^2}}\right)^2`$ (17)
$`+{\displaystyle \frac{(k_BT)^2}{2}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{(w_e^{ij})^2}{\omega _{g_i}\omega _{g_j}}}+k_BT{\displaystyle \underset{i}{}}{\displaystyle \frac{G_i^2}{\omega _{g_i}}}.`$
The reliability of our comparison was checked by ensuring that only a single isomer was visited during the Metropolis MC simulation. As can be seen in Fig. 1, the accuracy of the Gaussian theory is very good at such low temperatures. The global minimum geometry of CaAr13 has two spectral lines located at 23658 (2 degenerate lines) and 24034 cm-1, respectively. Due to its finite temperature, the cluster loses its instantaneous symmetry and the degeneracy is lifted. This explains why three distinct lines are seen in Fig. 1(a). In CaAr54 the Gaussian theory of absorption remains in good agreement with Monte Carlo data. This cluster has again 2 degenerate sharp lines at 23715 cm-1, and one broad line at 24080 cm-1. This time, the degeneracy is no longer lifted by temperature, and the two lines still merge into a single peak. In the classical theory the absorption intensity can also be obtained analytically by integrating Eq. (7) under the Condon assumption without performing the short-time expansion. Therefore, in the classical case, the differences between the calculated and simulated spectra in Fig. 1 can be attributed to anharmonicities only.

We now turn to a more complicated situation, namely CaAr37 in the quantum regime, and at moderate temperature $`T=25`$ K. We will show below that this temperature is located just above the melting point of this cluster. Many isomers are present, and this will provide a test of the superposition approximation. In the simulation, quantum effects have been included in a quasiclassical, corrective way by adding the corrective Feynman-Hibbs effective potentialfh  to the ground state energy. This potential was shown to give good results on argon clusters with respect to the more accurate Fourier-path-integral Monte Carlo simulation of Neirotti and coworkers,nfd  especially when compared with the alternative Wigner-Kirkwood corrections.wk  Since it is quasiclassical in nature, it cannot be used in the low temperature regime $`T5`$ K where only few isomers are present, and path-integral simulations would have to be carried out instead. Because the superposition method does not assume any barrier nor pathway between the isomers, it yields supposedly ergodic data within the set of available minima. Therefore a reliable comparison with simulation results must also correct for possible broken ergodicity, and for this we have used parallel tempering Monte Carlo.ptmc  Also, since the effective potential explicitely depends on temperature, the probability of accepting an exchange between adjacent trajectories should be modified accordingly.telec

Concerning the superposition approximation, and to quantify the importance of perturbative corrections in the partition function, we have separately considered zero-th order harmonic as well as second-order, anharmonic-corrected partition functions with the same database of isomers. The absorption spectra obtained with the Gaussian theory in the superposition approach and with quasiclassical MC simulations are represented in Fig. 2. The lowest energy isomer of this cluster, including the zero-point energy contribution, has three distinct spectral lines for the $`4s^24s4p`$ transition, namely at 23840, 23970, and 24260 cm-1. As can be seen from Fig. 2, none of these lines remain visible as peaks. We observe a very good agreement between the present theory and on-the-fly simulation, indicating that the relative weights of the isomers are correct. However, at this point, we cannot tell whether the visible shifts and broadenings are an effect of temperature only, or of the different spectroscopic signatures of the other isomers present. Some shifts in the peaks remain slightly approximate in the superposition approach, and the second-order corrected weights only marginally improve the results. At 25 K, quantum effects turn out to be rather small on the line width and shape, but they may be more important on the weights of the isomers, as will be seen below. Therefore, from Fig 2 we see that anharmonic corrections to the isomers weights have a relatively weak influence for quasiclassical clusters. More information on the importance of these corrections will be gained by looking at the thermodynamical curves in the next section.

The Gaussian theory presented in this section can be extended to purely quantum systems at $`T=0`$, for which delocalization is only due to the quantum nature of the vibrational modes, but not to an actual atomic displacement induced by a finite energy deposit. In this purpose, the Boltzmann weights in Eqn. (6) and (7) are replaced by the probability density of finding the system in its ground state at configuration $`𝐑`$ as the square of the wavefunction. In the harmonic approximation, this probability density is a multidimensional Gaussian function, and the integrated absorption spectrum is found to be Gaussian with line shape and width given by the $`\beta \mathrm{}`$ limits of Eqn. (14) and (15):

$`(\mathrm{\Delta }\omega )_q`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i}{}}{\displaystyle \frac{\omega _{g_i}^2\omega _e^{ii}}{\omega _{g_i}}}`$ (18)
$`(\overline{\omega }^2)_q`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{i}{}}\left({\displaystyle \frac{\omega _{g_i}^2\omega _e^{ii}}{\omega _{g_i}}}\right)^2`$ (19)
$`+{\displaystyle \frac{1}{8}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{\omega _e^{ij}}{\omega _{g_i}\omega _{g_j}}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}{\displaystyle \frac{G_i^2}{\mathrm{}\omega _{g_i}}}.`$
The behavior of the Gaussian theory in the $`T=0`$ quantum case has been tested on CaAr2 and CaAr13 using simple diffusion Monte Carlo (DMC) in internal coordinates (without global rotation and translation), and without importance sampling. The DMC calculations were carried out with 1000 replicas, and the wave functions were averaged over 105 steps after 105 initial equilibration steps. In addition, for CaAr2, the vibrational ground state wavefunction was also computed numerically on grids, using quantum wavepacket relaxation (WPR) method (imaginary time propagation).kos86:223  The propagation was performed for zero total angular momentum, representing the wavepacket by an harmonic oscillator Discrete Variable Representation (DVR) for each of the two Ca–Ar bond coordinates, and by a Legendre DVR for the bending angle $`\widehat{\mathrm{ArCaAr}}`$. The primitive basis for representing one bond coordinate includes 10 harmonic oscillator eigenfunctions. Associated Legendre functions with $`j<60`$ and $`m=0`$ were used for representing the bending mode. In both the DMC and WPR methods, which are both supposed to give essentially exact results, the absorption spectrum was calculated assuming vertical transitions from the quantum sample of geometries.

In Fig. 3 the absorption spectra calculated with these various methods are compared. For the small, triatomic cluster, the three methods give very similar peaks. In its equilibrium geometry, CaAr2 has three nondegenerate spectroscopic lines at 23621, 23670, and 23918 cm-1, respectively. Quantum delocalization of the vibrational modes leads to a rather small shift of these lines, and to a significant broadening. For CaAr13, the simulated peaks are also in good agreement with the Gaussian prediction, but we notice that some asymmetry is visible on the “exact” spectrum from diffusion Monte Carlo. These effects are clearly beyond a Gaussian description, and are due to anharmonic effects. These effects are relatively weak for the rather heavy calcium-doped argon clusters but may become more important for neon or helium clusters.

The above results suggest that the Gaussian theory of photoabsorption developed by Wadi and Pollakwp  upon the formal results of Yan and Mukamelyan  can give quite accurate results in both the classical and quantum regimes at low and moderate temperatures.

## IV Spectroscopic signature of phase changes

The Gaussian theory described in the previous Section will now be used with the superposition approach to investigate the influence of temperature on the absorption spectra in context with isomerization phenomena.

For the two smallest clusters chosen, CaAr6 and CaAr10, we performed extensive sampling of the ground state PES, and we believe that nearly all stable isomers have been identified. To apply the superposition method to larger CaArn clusters, the reweighting procedure from simulations must be used. Quasiclassical parallel tempering Monte Carlo simulations of $`10^7+2\times 10^6`$ equilibration cycles were used, for a set of temperatures equally spaced by 2.5 K in the range 7.5 K$`T60`$ K for $`n54`$, and 15 K$`T75`$ K for CaAr146. The isomers were found by quenching the 50 K trajectory for all sizes except CaAr37, for which the 25 K trajectory was used, and CaAr146, quenched at 60 K. Local minimizations were performed every 500 cycles. To prevent evaporation, a rigid spherical container was used in the MC simulation.lba ; labastie  The radii of the containers and the number of different isomers obtained from the quenches are given in Table 1 for the set of sizes chosen. We have also indicated the point group and energy of the global minima after including the zero-point harmonic contribution. These lowest-energy structures were generated from the basin-hopping or Monte Carlo + minimization algorithm.bh  They are shown in Fig. 4. In most cases, the global minimum of CaArn is very similar to the one of Arn+1, with one argon atom substituted by calcium. Due to the lower binding energy and larger equilibrium distance of CaAr with respect to Ar2,abinit  calcium is always found on the surface of

the remaining argon cluster. The nearly identical masses of calcium and argon atoms also induces relatively small changes in the vibrational properties, and the similarities between CaArn and Arn+1 remain when including the zero-point corrections. However, in the case of CaAr37, we find the somewhat surprising result that the global minimum geometry is decahedral, and remains so after including the harmonic zero point energy contribution. Slightly higher in energy, several icosahedral structures are found, either in the form of uncomplete multilayer icosahedra (Mackay type) or as polyicosahedra (anti-Mackay type). As in Ar38, quantum delocalization favors anti-Mackay structures over Mackay shapes,jcpquantum  yet a Mackay structure is found as the second isomer in both the classical and quantum regimes. This isomer is represented in Fig. 4(f). Additionally, the distorted truncated octahedron geometry lies above, but very close, to these isomers. As a consequence, the energy landscape of this particular cluster reveals interesting features, which will be investigated more deeply in the following article of this series.papdjw  This result is also consistent with the richer thermodynamical behavior generally observed in mixed van der Waals clusters with respect to homogeneous clusters.mixedvdw1 ; mixedvdw2

Therefore, even though CaArn and Arn+1 clusters share many structural features, they also exhibit significant and unexpected differences, which could affect the thermodynamics and dynamics in large extents.

### IV.1 Isomerization in CaAr6 and CaAr10

The lowest energy minimum of CaAr6 is a distorted pentagonal bipyramid. The heat capacity of this cluster, computed from the superposition approach with or without anharmonic perturbative corrections, is represented in Fig. 5(a). For comparison, we have also plotted on the same curves the direct results of parallel tempering Monte Carlo with the Feynman-Hibbs effective potential.fh  As is well known for this potential,jcpquantum  the canonical heat capacity diverges upwards at low temperatures. Due to its small size, this cluster shows a nearly featureless caloric curve, and a flat heat capacity (except the low temperature increase) in Fig. 5(a). The small hump at $`T20`$–25 K is indicative of isomerizations, and is followed by another bump near 45 K in the MC simulation. This latter broad peak is in fact spurious, as it comes from the dissociation of the weakly bound Ca atom. Previous worksnordiek ; calvoc60  have identified this peak as the signature of the liquid-vapor transition, from its dependence on the container radius. Here it would have been difficult to employ a radius smaller than 6Å, because the motion of argon atoms would have been hindered. The 25 K hump is barely visible, because it is hidden and smoothed out by quantum effects. The agreement between the superposition approximation and the quasiclassical MC simulations is very good except for the extra liquid-vapor bump. This feature is obviously missing from the assumptions made in the superposition method. The quantitative effect of perturbative corrections is small for this cluster.

CaAr10 has the ground state of an uncomplete icosahedron. Some isomers involving different locations of the calcium atom lie not far above the global minimum, resulting in a shoulder in the heat capacity near 10–15 K in Fig. 5(b). The main peak near 30 K is due to a more complete rearrangement of argon atoms, and is reproduced quantitatively by the superposition approach after including perturbative corrections in the partition function.

The absorption spectra of these two clusters are represented in Fig. 5(c) and 5(d) as a function of temperature. In each of these figures, the classical spectra obtained from the limit $`\mathrm{}0`$ has been superimposed to the quantum spectra. For the two sizes, the patterns are rather similar: the spectrum varies from a clear three-peak structure at low temperature towards a mostly single peak shape at high temperatures. The transition occurs roughly at 20 K in CaAr6, and near 30 K in CaAr10. In the case of CaAr6 the high temperature main peak appears made of two smaller peaks, in agreement with simulations.isspic  From a general point of view, quantum effects mainly consist in a broadening of the peaks. As expected, they get smaller and smaller as temperature increases, and becomes barely visible above 25 K.

Therefore, in these two small clusters, absorption spectroscopy seems to be an effective mean of detecting isomerization events, as the range of variation in the spectrum approximately matches the width of the heat capacity changes. As cluster size grows, one can thus expect to see sharper transitions, and to possibly use the calcium chromophore as a probe of the argon cluster melting point.

### IV.2 Phase changes in CaArn, $`\mathrm{𝟏𝟑}\mathbf{}𝒏\mathbf{}\mathrm{𝟓𝟒}`$

We have chosen a set of sizes in the range $`13n54`$, spanning between 1 and 2 icosahedral layers, and which are likely to exhibit significant finite size effects. We first discuss the thermodynamical behavior of these clusters. The heat capacities obtained using the superposition method, with harmonic or second-order corrected anharmonic partition functions, are represented in Fig. 6 for the clusters corresponding to $`n=13`$, 30, 37, and 54. In the case of the larger cluster CaAr54, and also for CaAr146, it was not possible to include all the required anharmonic corrections due to the too large computational requirements. Since they are not the main focus of the present paper, these thermodynamical data will be only briefly discussed in order to guide our interpretation of the following absorption spectra.

In general, we find again a rather good agreement between the quantum superposition results and the Monte Carlo simulations for all four sizes. As expected, anharmonic corrections do have a positive quantitative effect by bringing the superposition data closer to the simulation results. However, no new qualitative feature is visible when they are taken into account. All heat capacities exhibit a main peak indicative of the solidlike-liquidlike phase change.labastie  This peak is located near 40–45 K for CaAr13, CaAr30 and CaAr54, and near 20 K for CaAr37. This lower melting point is similar to what is observed in the 38-atom argon cluster described by the Lennard-Jones potential.ar38doye ; ar38calvo  In addition to this main peak, extra humps or even peaks can be observed at lower temperatures in CaAr30 and CaAr37. These features are indicative of preliminary isomerization resulting from solid-solid transitions in these clusters.doyecalvo  Again, they bear similarities with the known finite temperature behavior of LJ clusters with the same size.ar38doye ; ar38calvo ; ar31 ; frantz  In CaAr13, the quantum character of the specific heat hides again a shoulder in the classical curve, as is the case in CaAr6. This cluster happens to have a much larger number of stable isomers than Ar14. In particular the capped icosahedron can be found with the calcium atom either in adatom location (lowest energy structure) or inside the icosahedral shell, resulting in a set of 4 different isomers slightly higher in energy. Structures not based on the capped icosahedron, or with the calcium atom in the center of the icosahedron, are less stable than these five geometries. Therefore the energy spectrum of local minima can be put into three sets, one being the global minimum only, one being the other capped icosahedral isomers with calcium in the icosahedral shell, and all other isomers not build on the capped icosahedron motif. Because of these three sets, the melting process involves essentially two steps. In Fig. 7(a) we plot the probabilities of finding the CaAr13 cluster in its ground state, or in one of the 4 subsequent isomers, or in the remaining, liquid configurations. These probabilities are computed using the anharmonic corrected quantum partition functions. They display a clear two-step isomerization pattern. In the first step, at about 30 K, the calcium atom is able to leave its capping site in favor of one of the twelve icosahedral sites. Then, above 40 K, more disordered structures become thermally available the melting peak takes place.

In CaAr30 the preliminary solid-solid transition near $`2.5`$ K involves polyicosahedral and multilayer icosahedral isomers, as in Ar31.jcpquantum ; ar31  Even though CaAr37 does not have the same global minimum as Ar38, a similar thermodynamic behavior occurs. We have calculated from the harmonic superposition approximation the probabilities of finding the cluster either in its decahedral ground state isomer, in its (Mackay-type) icosahedral isomers, or in any of the remaining isomers. The variation of these probabilities with increasing temperature are displayed in Fig. 7(b). They show that the transition at about 5 K is due to a global decahedral-icosahedral rearrangement, and bears some similarity with the (C60)14 cluster.prlcomment  The anti-Mackay and truncated octahedral structures are only marginally populated at thermal equilibrium, and the main heat capacity peak is associated with the emergence of a large number of disordered configurations. Hence it can be considered as the signature of melting. For comparison, we have also plotted in Fig. 7(b) the corresponding curves in the classical limit. While the global behavior remains essentially unchanged, the preliminary transition occurs at a higher temperature, near 10 K.

The photoabsorption spectra of these clusters are represented in Fig. 8 for increasing temperatures in the range $`0T60`$ K for all sizes except CaAr37, for which $`0T40`$ K. The classical spectra have also been plotted as dashed lines for $`T>0`$. Except for CaAr37, the classical spectrum looks very much like its quantum counterpart with slightly shifted but noticeably narrower excitation peaks. The four clusters display different absorption spectra, characteristic of finite-size effects. In CaAr13, the spectrum gradually evolves from a two-peak structure (at about 23700 cm-1 and 24100 cm-1) to another two-peak structure (at about 23650 and 23800 cm-1). The peak ending at 23800 cm-1 appears first at 30 K, and originates from isomers where the Ca atom lies in the main icosahedral shell. No real signature of melting at 40 K is seen in Fig. 8(a), however the present spectra can still be interpreted in terms of isomerizations within the cluster.

In the case of CaAr30, the low temperature isomerization below 5 K has a strong effect on the absorption peaks, which are red shifted by about 200 wavenumbers. The spectra of this cluster also evolves very continuously, and the two main peaks at 23550 and 24000 cm-1 undergo a slight red shift and thermal broadening. The weak melting phase change is reflected on these spectra, and the most prominent feature remains the sharp isomerization below 5 K.

Near $`T=10`$ K, the multiple peak structure of the spectrum in CaAr37 also reveals that several isomers are coexisting with distinct spectroscopic signatures. The solid-solid transition seen in the heat capacity of CaAr37 at low temperature is reflected on the absorption spectra in Fig. 8(c). The three-peak $`T=0`$ K spectrum is peculiar to the decahedral ground state. At $`T10`$ K, a significant blue shift occurs in the middle peak, providing a signature of the icosahedral isomers. The quantitative difference in the solid-solid transition temperature induced by quantum effect is also seen in the absorption spectra between 5 and 15 K. Near 15 K, the spectrum can be clearly separated into more than three peaks, and suggests that distinct icosahedral isomers are present. At around 20 K, the various peaks merge into a double-peak structure, where the peak at 23800 cm-1 dominates the one at 24750 cm-1. Here we find a clear thermodynamical influence on the spectroscopic properties. Hence, in this cluster, the absorption spectrum contains notable informations about isomerization and phase changes, including some rather fine details on the characteristic temperatures.

The equilibrium geometry of CaAr54 is a double layer Mackay icosahedron where the Ca atom is located on a vertex site, thus keeping a relatively high symmetry (C5v). As in Ar55, the heat capacity near peak 40 K indicates volume melting,labastie ; wales93  while preliminary surface melting has a very small thermal signature. In Fig. 8(d) we observe that the two absorption peaks undergo a significant change near 30 K, and this corresponds to coexisting isomers with the Ca atom located on different substitutional sites of the external icosahedral shell. In this case the photoexcitation spectrum is sensitive to surface melting. At higher temperatures $`T>40`$ K, many isomers are found in which the metal atom is more fully solvated by argon. The spectra show another qualitative change, as the peak near 23600 cm-1 disappears and is replaced by a broad bump near 24300 cm-1.

### IV.3 CaAr146

It is interesting to compare the situation in CaAr54 with the one in the larger Mackay icosahedron, CaAr146. As in CaAr54, the ground state equilibrium geometry has C5v symmetry, with calcium located in a vertex site of the external icosahedral layer. This cluster has a huge number of isomers, and many high lying minima are missed during the quenching process. In the vast majority of structures in our sample the calcium atom remains in the outer argon layers. This is not surprising, since, and as was previously mentioned, the calcium-argon bond is weaker by about 70% than the argon-argon bond.epjd  Despite the relatively small size of our sample with respect to the number of existing structures, the heat capacities in Fig. 9(a) exhibit a reasonably good agreement between simulation and the superposition method, the high peak near 55 K being the signature of the solidlike-liquidlike phase change. Melting is essentially caused by the motion of argon atoms, and is relatively unsensitive to the location of calcium. This explains the shape and the evolution of the absorption spectrum in Fig. 9(b). The two peaks at 23720 and 24150 cm-1 become broader and vanish between 25 K and 30 K, to be replaced by a major peak at 23600 cm-1. This transition is correlated with fluctuations in the location of the calcium atom over the cluster surface. On the other hand, melting at 55 K does not lead to any visible change in the absorption spectrum.

### IV.4 Discussion

The above results may not be immediately amenable to experimental comparison, because of the problem of size selection. Our simulation results show that complex isomerizations involving strong changes in the local environment of calcium could be, in principle, detected in experiments on size-selected clusters initially thermalized. However a typical sample of clusters in molecular beams is made of a broad size distribution, where the width and the center have values close to each other. The common feature of these clusters is the surface location of the calcium atom, and we may expect the absorption spectrum to be a possible probe of surface melting in these clusters, provided that the surface melting temperature does not significantly depend on cluster size. Even though calcium in bulk argon has a very distinct, three-peak spectroscopic signature,epjd  it may be hard to detect volume melting in which calcium would be fully solvated. In fact, such van der Waals clusters evaporate atoms at temperature close to the melting point. From a statistical point of view, the weakly bonded calcium atom is likely to evaporate first. In addition the situations where calcium is fully solvated are less probable. Therefore the true volume melting transition occuring with full solvation of calcium may hardly be observed.

Provided that one could size-select these clusters, we have shown in the present work that spectroscopy can possibly probe their phase changes. More generally, intricate finite size effects in the thermodynamical behavior could be expected to be seen, as in CaAr13, CaAr30 or CaAr37. In these special clusters, coexisting isomers and multiple phase-likes induce complex variations of the absorption properties with increasing temperature.

Even with experimental setups involving selected sizes, another problem may come into play to hinder one from getting the complex spectroscopic information. All the above results are supposed to be at thermal equilibrium, without assuming how this equilibrium state had been reaching. In some cases, the relaxation time required for achieving equilibration may be very long because of multiple funnels in the energy landscape and kinetic trapping. This has been observed in the Ar38 cluster by Miller and coworkers,ar38miller  where typical timescales of 10 s have been estimated for the interfunnel transitions at moderate temperatures. Relaxation to equilibrium and time-dependent absorption spectra will be the subject of the following article.papdjw

## V Conclusion

In this paper, we have investigated some methodological and practical aspects of finite temperature spectroscopy in chromophore-doped van der Waals clusters at thermal equilibrium. Using the harmonic superposition approach improved with anharmonic perturbation corrections to the partition function, the basic thermodynamical observables were calculated from samples of local minima on the energy landscape. The photoabsorption spectrum was obtained using the recent Gaussian theory by Wadi and Pollak wp  in the harmonic approximation. These theories were tested and validated on selected CaArn clusters, in both the classical and quantum regimes, at low or moderate temperatures. Although anharmonic corrections could also be incorporated to the various quantities (line shifts and widths) in this theory, we found them unnecessary in the present work because of the other underlying approximations.

We have studied the possible spectroscopic signatures of the various isomerizations and phase changes in CaAr6, CaAr10, CaAr13, CaAr30, CaAr37, CaAr54, and CaAr146. The two smaller clusters show relatively simple isomerizations leading to a similar spectroscopic signature with a clear three- to one-peak pattern. Both CaAr30 and CaAr37 display strong (non monotonic) finite-size effects and solid-solid transitions. In CaAr13, melting involves the calcium atom taking place of one argon atom in the icosahedral shell. In contrast with most other sizes, it is associated with the appearance of a new absorption peak. In CaAr54 and CaAr146, preliminary surface melting occurs at nearly 30 K, and is reflected in the absorption spectrum. In all these phenomena, the local environment of calcium changes rather abruptly, resulting usually in clear shifts in some absorption peaks. Volume melting towards a calcium atom fully solvated by argon has a potentially significant spectral signature, however calcium is less bound to argon than argon itself, hence it is likely to dissociate before volume melting occurs. One should mention that the main heat capacity peak indicates melting of the argon cluster host, especially for large sizes. While the host cluster is liquidlike the calcium is still free, to some extent, to glide over the cluster surface. Being in this state, the rearrangements inside the argon cluster have little impact on the global spectroscopy of the chromophore, and the system can be expected to behave as a pseudo diatomics Ca–Arn. This point will be examined further in the third paper.papjmm

The techniques used here are valid in both the quantum and classical regimes, which refers to the limit $`\mathrm{}0`$. They could be applied to other chromophore-doped van der Waals clusters, provided that potential energy surfaces are built first. In particular, neon clusters could be profitably investigated using the superposition method and the Gaussian theory of absorption. For helium, the harmonic approximation may break down, and even anharmonic corrections should be unsufficient. In this case, quantum Monte Carlo should be coupled with on-the-fly histogram accumulation to compute the spectrum. More generally, the present methods are relevant to any statistical driven absorption spectrum in heterogeneous or homogeneous polyatomic molecules. For instance, the thermal effects on the absorption spectrum in small sodium clusters moseler ; schmidt  could be investigated using appropriate models.ddtb ; naab  As concerns the approximations within the Gaussian theory, the most questionable point is the short-time expansion in the autocorrelation function. This approximation could be overcome using the formal results of Yan and Mukamel.yan  Unfortunately, this method involves a heavy numerical effort in the context of multiple isomers, and it still assumes harmonic forces. Finally, the superposition approximation provides a framework to study time-dependent properties over long time scales not available to molecular dynamics, by solving appropriate master equations.med  The application of this technique to photoabsorption spectra will be made in the next paper of this series.papdjw

## Acknowledgments

We thank CALMIP for a generous allocation of computer resources. FC also acknowledges interesting discussions with Dr. J. Vigué.

## Appendix: Perturbation expansion of the partition function

The ground state Hamiltonian is written in cartesian coordinates $`𝐏,𝐑`$ as $`H(𝐏,𝐑)=H_a(𝐏,𝐑)+\delta V(𝐑)`$, where $`H_a`$ is the harmonic term and $`\delta V`$ the extra, anharmonic contribution. The Taylor expansion of $`\delta V`$ near the minimum $`𝐑_0`$ contains products of at least 3 coordinates. The partition function $`Z(\beta )=\mathrm{Tr}\mathrm{exp}(\beta H)`$ is now expanded in a semiclassical fashion, in order to yield at least the same terms as the corresponding classical expansion. This approximation is actually correct from the quantum point of view as long as we keep it to first order. However, the rigorous perturbation treatment, which should be performed with path integrals, will not allow us to obtain the simple corrections that can be compared with the classical expressions of Ref. jcpanhar, . Thus we write in a phenomenological way

$$Z(\beta )\mathrm{Tr}\left[e^{\beta H_a}\underset{n=0}{\overset{2}{}}\frac{(\beta )^n}{n!}(\delta V)^n\right],$$
(20)
and we keep only terms containing coordinates up to the power 12, which contribute to $`T^2`$ corrections in the low-temperature classical partition function.jcpanhar  As in the classical case, all odd products have a zero contribution, but it is no longer possible to sort the various nonzero terms according to their power in $`\beta `$. It appears more natural to sort them according to the even power in coordinates. Following our previous notationsjcpanhar  we write the anharmonic term $`\delta V`$ as a function of the normal mode coordinates $`𝐐`$ as

$`\delta V(𝐐)`$ $`=`$ $`{\displaystyle \underset{k3}{}}{\displaystyle \frac{1}{k!}}{\displaystyle \underset{i_1\mathrm{}i_k}{}}{\displaystyle \frac{^kV}{Q_{i_1}\mathrm{}Q_{i_k}}}|_{𝐐=0}Q_{i_1}\times \mathrm{}\times Q_{i_k}`$ (21)
$`=`$ $`{\displaystyle \underset{k3}{}}{\displaystyle \frac{1}{k!}}{\displaystyle \underset{i_1\mathrm{}i_k}{}}\mathrm{\Gamma }_{i_1\mathrm{}i_k}^{(k)}Q_{i_1}\times \mathrm{}\times Q_{i_k}.`$
The corrective terms can be expressed as averages on the harmonic ground state:

$$Z(\beta )Z_0(\beta )[1+\gamma _4+\gamma _6+\gamma _8+\gamma _{10}+\gamma _{12}],$$
(22)
where $`\gamma _i`$ labels the correction due to power $`i`$ in coordinates. These terms are explicitely given by

$`\gamma _4`$ $`=`$ $`{\displaystyle \frac{\beta }{24}}{\displaystyle \underset{i_1\mathrm{}i_4}{}}\mathrm{\Gamma }_{i_1\mathrm{}i_4}^{(4)}Q_{i_1}\times \mathrm{}\times Q_{i_4};`$ (23)
$`\gamma _6`$ $`=`$ $`{\displaystyle \underset{i_1\mathrm{}i_6}{}}\left[{\displaystyle \frac{\beta }{720}}\mathrm{\Gamma }_{i_1\mathrm{}i_6}^{(6)}+{\displaystyle \frac{\beta ^2}{72}}\mathrm{\Gamma }_{i_1i_2i_3}^{(3)}\mathrm{\Gamma }_{i_4i_5i_6}^{(3)}\right]Q_{i_1}\times \mathrm{}\times Q_{i_6};`$ (24)
$`\gamma _8`$ $`=`$ $`{\displaystyle \underset{i_1\mathrm{}i_8}{}}\left[{\displaystyle \frac{\beta ^2}{1152}}\mathrm{\Gamma }_{i_1\mathrm{}i_4}^{(4)}\mathrm{\Gamma }_{i_5\mathrm{}i_8}^{(4)}+{\displaystyle \frac{\beta ^2}{720}}\mathrm{\Gamma }_{i_1i_2i_3}^{(3)}\mathrm{\Gamma }_{i_4\mathrm{}i_8}^{(5)}\right]Q_{i_1}\times \mathrm{}\times Q_{i_8};`$ (25)
$`\gamma _{10}`$ $`=`$ $`{\displaystyle \underset{i_1\mathrm{}i_{10}}{}}[{\displaystyle \frac{\beta ^2}{28800}}\mathrm{\Gamma }_{i_1\mathrm{}i_5}^{(5)}\mathrm{\Gamma }_{i_6\mathrm{}i_{10}}^{(5)}+{\displaystyle \frac{\beta ^2}{17280}}\mathrm{\Gamma }_{i_1\mathrm{}i_4}^{(4)}\mathrm{\Gamma }_{i_5\mathrm{}i_{10}}^{(6)}`$ (26)
$`{\displaystyle \frac{\beta ^3}{1728}}\mathrm{\Gamma }_{i_1i_2i_3}^{(3)}\mathrm{\Gamma }_{i_4i_5i_6}^{(3)}\mathrm{\Gamma }_{i_7\mathrm{}i_{10}}^{(4)}]Q_{i_1}\times \mathrm{}\times Q_{i_{10}};`$
$`\gamma _{12}`$ $`=`$ $`{\displaystyle \underset{i_1\mathrm{}i_{12}}{}}[{\displaystyle \frac{\beta ^2}{1036800}}\mathrm{\Gamma }_{i_1\mathrm{}i_6}^{(6)}\mathrm{\Gamma }_{i_7\mathrm{}i_{12}}^{(6)}{\displaystyle \frac{\beta ^3}{51840}}\mathrm{\Gamma }_{i_1i_2i_3}^{(3)}\mathrm{\Gamma }_{i_4i_5i_6}^{(3)}\mathrm{\Gamma }_{i_7\mathrm{}i_{12}}^{(6)}`$ (27)
$`{\displaystyle \frac{\beta ^3}{17280}}\mathrm{\Gamma }_{i_1i_2i_3}^{(3)}\mathrm{\Gamma }_{i_4\mathrm{}i_7}^{(4)}\mathrm{\Gamma }_{i_8\mathrm{}i_{12}}^{(5)}{\displaystyle \frac{\beta ^3}{82944}}\mathrm{\Gamma }_{i_1\mathrm{}i_4}^{(4)}\mathrm{\Gamma }_{i_5\mathrm{}i_8}^{(4)}\mathrm{\Gamma }_{i_9\mathrm{}i_{12}}^{(4)}`$
$`+{\displaystyle \frac{\beta ^4}{31104}}\mathrm{\Gamma }_{i_1i_2i_3}^{(3)}\mathrm{\Gamma }_{i_4i_5i_6}^{(3)}\mathrm{\Gamma }_{i_7i_8i_9}^{(3)}\mathrm{\Gamma }_{i_{10}i_{11}i_{12}}^{(3)}]Q_{i_1}\times \mathrm{}\times Q_{i_{10}}.`$
The above expressions can be considerably simplified by factorizing the average products $`Q_i\times \mathrm{}\times Q_j`$. This is made possible by Wick’s theorem, which holds for general bosonic ground states:

$$Q_i^{2n}=\frac{1}{2^n}\frac{(2n)!}{n!}Q_i^2^n=(2n1)!!Q_i^2^n,$$
(28)
with

$$Q_i^2=\frac{\mathrm{}}{\omega _i}\left(\frac{1}{2}+\frac{1}{e^{\beta \mathrm{}\omega _i}1}\right).$$
(29)
It is then more convenient to write the $`\gamma _i`$’s as diagrams, in which each $`k`$vertex represents one of the $`\mathrm{\Gamma }_{i_1\mathrm{}i_k}^{(k)}`$ term. For instance we have

$$\gamma _4=\frac{\beta }{8}\text{ }\text{}$$
(30)
with

$$\text{ }\text{}=\underset{ij}{}\mathrm{\Gamma }_{iijj}^{(4)}Q_i^2Q_j^2.$$
(31)
The two diagrams giving $`\gamma _6`$ are

$$\gamma _6=\frac{\beta }{48}\text{ }\text{}+\frac{\beta ^2}{8}\text{ }\text{}+\frac{\beta ^2}{12}\text{ }\text{}.$$
(32)
The other perturbative terms $`\gamma _8`$, $`\gamma _{10}`$, and $`\gamma _{12}`$ contain a larger number of diagrams, including non connected ones. The fewest diagrams are obtained in $`\mathrm{ln}Z`$, which still has 46 different contributions corresponding to Eq. (22):

$`\mathrm{ln}Z=\mathrm{ln}Z_0`$ $`{\displaystyle \frac{\beta }{8}}\text{ }\text{}{\displaystyle \frac{\beta }{48}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{8}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{12}}\text{ }\text{}`$ (33)
$`+{\displaystyle \frac{\beta ^2}{48}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{16}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{16}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{12}}\text{ }\text{}`$
$`+{\displaystyle \frac{\beta ^2}{128}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{48}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{240}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{32}}\text{ }\text{}`$
$`+{\displaystyle \frac{\beta ^2}{48}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{16}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{8}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{8}}\text{ }\text{}`$
$`{\displaystyle \frac{\beta ^3}{8}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{12}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{256}}\text{ }\text{}+{\displaystyle \frac{\beta ^2}{192}}\text{ }\text{}`$
$`+{\displaystyle \frac{\beta ^2}{1440}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{64}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{24}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{32}}\text{ }\text{}`$
$`{\displaystyle \frac{\beta ^3}{32}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{16}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{72}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{32}}\text{ }\text{}`$
$`{\displaystyle \frac{\beta ^3}{8}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{48}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{32}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{8}}\text{ }\text{}`$
$`{\displaystyle \frac{\beta ^3}{12}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{24}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{24}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{16}}\text{ }\text{}`$
$`{\displaystyle \frac{\beta ^3}{48}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{48}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{24}}\text{ }\text{}{\displaystyle \frac{\beta ^3}{32}}\text{ }\text{}`$
$`{\displaystyle \frac{\beta ^3}{48}}\text{ }\text{}+{\displaystyle \frac{\beta ^4}{16}}\text{ }\text{}+{\displaystyle \frac{\beta ^4}{16}}\text{ }\text{}+{\displaystyle \frac{\beta ^4}{48}}\text{ }\text{}`$
$`+{\displaystyle \frac{\beta ^4}{8}}\text{ }\text{}+{\displaystyle \frac{\beta ^4}{24}}\text{ }\text{}.`$
It is important to notice that this perturbative expansion includes all 18 terms up to second-order in the classical expansion,jcpanhar  and that 28 remaining terms would contribute classically to 3rd and higher orders in $`T^n`$. These terms correspond to the $`\mathrm{\Gamma }^{(5)}\times \mathrm{\Gamma }^{(5)}`$, $`\mathrm{\Gamma }^{(4)}\times \mathrm{\Gamma }^{(6)}`$, $`\mathrm{\Gamma }^{(3)}\times \mathrm{\Gamma }^{(4)}\times \mathrm{\Gamma }^{(5)}`$, and $`\mathrm{\Gamma }^{(4)}\times \mathrm{\Gamma }^{(4)}\times \mathrm{\Gamma }^{(4)}`$ products, but cannot be neglected a priori with respect to the other corrections in the quantum regime.

In the superposition method framework, these perturbative corrections act on the individual partition functions, where they are expected to give a more accurate estimate of the relative weights of the isomers as a function of temperature, as well as more accurate thermodynamical observables.


format: text

text:
# 1 Introduction

## 1 Introduction

It has been known for a long time that beyond-leading-order calculations in massless thermal field theories exhibit infrared divergences . In the case of abelian theories, systematic resummation techniques have been developed to cure infrared divergences. This leads to perturbative expansions which are nonanalytic in the coupling. In contrast to the case of abelian theories, no resummation method has so far managed to cure infrared divergences in non-abelian gauge theories like QCD . In QED, infrared divergences mainly arise in Feynman diagrams with multiple insertions of the static photon self-energy. A resummation of all diagrams with an arbitrary number of such insertions (ring diagrams) can be shown to eliminate all infrared divergences. Consequently, quantities like the free energy can, in principle, be expanded to arbitrary order in powers of the coupling. In QCD however, since gluons couple to themselves, the behaviour of the theory is much more singular: diagrams with static spatial gluon loops are increasingly infrared divergent when their order gets higher. This may signal the presence of a magnetic mass which would act as an infrared cut-off to regulate infrared divergent integrals. Although such a magnetic mass is potentially of order $`g^2T`$ since spacelike gluons remain massless to lowest order, it cannot be computed perturbatively. This is because beyond one-loop order, all the diagrams in a self-consistent expansion of the magnetic mass will be of the same order in the coupling . As a result, unresummed perturbation theory does not seem to work any more and all those diagrams that contribute to the same order in the coupling should somehow be resummed in a similar fashion to simpler theories like $`\lambda \varphi ^4`$ or QED. Due to the complexity of QCD, this is unfortunately impossible in practice.
If one uses dimensional regularisation, with $`d`$ the dimension of space-time, the breakdown of perturbation theory from some order signals itself through infrared $`1/(d4)`$ poles in the perturbative series from a given order. It reflects the fact that the exact quantity one tries to compute is nonanalytic in the coupling when $`d=4`$. In this work we therefore examine dimensional regularisation at finite temperature in more detail. Dimensional regularisation was devised as a method to separate the finite and divergent parts of Feynman diagrams arising beyond leading order in perturbation theory . At finite temperature, there are infrared $`1/(d4)`$ divergences in addition to the zero temperature ultraviolet ones so, if one use dimensional regularisation in the usual way, one has to make sure that these IR divergences somehow cancel each other in the expression for physical quantities. Consequently, one could wonder whether making an expansion in powers of $`g`$ and $`d4`$ commute. This is relevant because if we keep $`d4`$ nonzero, the perturbative expansion in powers of $`g`$ actually involves $`g/(d4)`$ and therefore is illegal since we are due to take $`d`$ to $`4`$ at the end of the calculation. As a result, an expansion in powers of $`g`$ is is only valid if $`1/(d4)`$ poles somehow cancel each other in the expression for physical quantities. A case where the previous limits do not commute is the following: consider

$$A=\frac{C}{d4+g^2}$$
(1.1)
If we first expand $`A`$ in power series of the coupling, its terms are divergent in the limit $`d4`$ :

$$A=\frac{C}{d4}\frac{Cg^2}{(d4)^2}+\mathrm{}$$
(1.2)
However, in $`4`$ dimensions, $`A=C/g^2`$ is perfectly finite and is nonanalytic in the coupling.
This work investigates the consequences of this kind of problem in the case of large-$`N`$ scalar theory, which will give us some insight about the way things might work in QCD. In next two sections, we derive and renormalise an exact formula for the thermal pressure in the case of large-$`N`$ scalar field theory. In section 4, we study the singular structure of the expansion of the thermal mass and show that infrared divergences arise from a given order onwards in the renormalised perturbation series. We argue that those divergences crop up because the perturbative expansion has been performed before taking $`d`$ to the physical dimension, in the same fashion as for the example given in (1.1). We will see that the way to rescue the situation is to resumm the perturbative series and that we then obtain the correct expansion, ie the one we get if we go to $`4`$ space-time dimensions first. Finally, we show this problem originates in the fact that the finite temperature contribution to the self-energy is actually nonanalytic in the coupling when $`d=4`$ and involves terms proportional to $`g\mathrm{log}g`$ and half-integer powers of the coupling.

## 2 Unrenormalised formula

In this work, we use large-$`N`$ scalar field theory as a toy model to investigate the commutativity of the analytic continuation to four space-time dimensions and the perturbative expansion. Consider the theory described by the unrenormalised Lagrangian density

$$(\lambda ,x)=\frac{1}{2}(_\mu \varphi )^2\frac{1}{2}m^2\varphi ^2\frac{\lambda }{4!}\frac{3}{N+2}\varphi ^4$$
(2.1)
where $`\varphi `$ is an $`N`$-component field. We are interested in the limit $`N\mathrm{}`$ with $`\lambda /N`$ fixed. The grand canonical partition function reads

$$Z_\beta (\lambda )=𝒟\varphi \mathrm{exp}\left[i_\tau ^{\tau i\beta }𝑑x^0d^{d1}x(\lambda ,x)\right]$$
(2.2)
where $`\beta =1/T`$ and $`\tau `$ is some arbitrary time. The thermodynamic definition of the pressure is, as usual,

$$P=TV^1\mathrm{log}Z_\beta $$
(2.3)
If we perform the change of variables $`\varphi \lambda ^{1/2}\varphi `$ in the path integral, we can make the coupling dependence of the partition function more apparent  and obtain

$$\frac{Z_\beta (\lambda )}{Z_\beta (0)}=\frac{𝒟\varphi \mathrm{exp}\left[\frac{i}{\lambda }_\tau ^{\tau i\beta }𝑑x^0d^{d1}x(1,x)\right]}{𝒟\varphi \mathrm{exp}\left[\frac{i}{\lambda }_\tau ^{\tau i\beta }𝑑x^0d^{d1}x(0,x)\right]}$$
(2.4)
If we define $`\widehat{P}=P(T,\lambda )P(T,0)`$ and use the translation invariance property that

$$_\tau ^{\tau i\beta }𝑑x^0d^{d1}x(1,x)=i\beta V(1,0)$$
(2.5)
it is straightforward to obtain from (2.4) the unrenormalised equation

$$\frac{\widehat{P}(T,\lambda )}{\lambda }=\frac{1}{\lambda }\left((\lambda ,0)(0,0)_{FREE}\right)$$
(2.6)
Because it will simplify the equations in what follows we will hereafter use the notation $`\mathrm{\Delta }A<A><A>_{FREE}`$ for any operator $`A`$.

### 2.1 Calculation of $`\mathrm{\Delta }`$

The main analytic non-perturbative tools available here to compute $`\mathrm{\Delta }`$ are the Schwinger-Dyson equations which can be regarded as the quantum equations of motion for Green’s functions . These equations are the same as at zero temperature; the only difference lies in the boundary conditions : at finite temperature, the time periodicity (for bosons) or anti-periodicity (for fermions) of classical fields (ie the integration variable in the path integral) provides Green’s functions with temperature dependent boundary conditions. In large-$`N`$ scalar field theory, the self-energy is exactly calculable and it is therefore desirable to derive an expression for $`\mathrm{\Delta }`$ that only involves the propagator since, using the Dyson resummation, the latter can in turn be written in terms of the self-energy. To achieve this, we have to express the thermal four-point function in $`\mathrm{\Delta }`$ in terms of the propagator. This can be done using the Schwinger-Dyson equation for the propagator. To derive this equation, consider the generating functional of the thermal Green’s functions :

$$Z_\beta (J)=𝒟\varphi \mathrm{exp}\left[i_Cd^dx\left((\varphi )+\varphi (x)J(x)\right)\right]$$
(2.1)
We use the Keldysh variant of the real-time thermal field theory , so the contour $`C`$ runs along the real axis from $`\mathrm{}`$ to $`\mathrm{}`$, back to $`\mathrm{}`$ and then down to $`\mathrm{}i\beta `$. According to the usual Feynman prescription, we have to add $`iϵ\varphi ^2`$ to the Lagrangian density in the generating functional (2.1) in order to ensure the convergence of the path integral for large $`\varphi `$ (in what follows, this prescription will be implicit). Then we have the following identity

$$𝒟\varphi \frac{\delta }{\delta \varphi (x)}\mathrm{exp}\left[i_Cd^dx^{}\left((\varphi )+\varphi (x^{})J(x^{})\right)\right]=0$$
(2.2)
which reads

$`{\displaystyle }𝒟\varphi [(^2+m^2)\varphi (x))+{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{1}{N+2}}(\varphi (x))^3]`$
$`\mathrm{exp}\left[i{\displaystyle _C}d^dx^{}\left((\varphi )+\varphi (x^{})J(x^{})\right)\right]`$ $`=`$ $`Z_\beta (J)J(x)`$ (2.3)
By applying $`\delta /\delta \varphi (y)`$ to both sides of this equation and setting $`J`$ to zero, we obtain

$$T_c\varphi (y)(^2+m^2)\varphi (x)+\frac{\lambda }{2}\frac{1}{N+2}T_c(\varphi (x))^3\varphi (y)=i\delta _c^{(d)}(xy)$$
(2.4)
and similarly

$$T_c\varphi (y)(^2+m^2)\varphi (x)_{FREE}=i\delta _c^{(d)}(xy)$$
(2.5)
where $`T_C`$ and $`\delta _c`$ denote the time ordering and delta function along the integration contour in complex time. When $`x=y`$, the last two equations lead to

$$\varphi (^2+m^2)\varphi +\frac{\lambda }{2}\frac{1}{N+2}\varphi ^4=\varphi (^2+m^2)\varphi _{FREE}$$
(2.6)
This implies

$$\mathrm{\Delta }=\frac{1}{4}\mathrm{\Delta }\left[\varphi (^2+m^2)\varphi \right]$$
(2.7)
With the Keldysh contour, we get

$$\mathrm{\Delta }=\frac{1}{4}\frac{d^dq}{(2\pi )^d}(q^2m^2)\left(D_{12}(q)D_{12,FREE}(q)\right)$$
(2.8)
where the subscript 12 refers to the element of the $`2\times 2`$ thermal matrix propagator. This propagator has the structure

$$𝐃(q)=𝐌(q^0)\stackrel{~}{𝐃}(q)𝐌(q^0)$$
(2.9)
with

$`𝐌(q^0)`$ $`=`$ $`\sqrt{n(q^0)}\left[\begin{array}{cc}e^{\frac{1}{2}\beta |q^0|}& e^{\frac{1}{2}\beta q^0}\\ e^{\frac{1}{2}\beta q^0}& e^{\frac{1}{2}\beta |q^0|}\end{array}\right]`$
$`\stackrel{~}{𝐃}(q)`$ $`=`$ $`\left[\begin{array}{cc}D(q)& 0\\ 0& D^{}(q)\end{array}\right]`$
$`D(q)`$ $`=`$ $`{\displaystyle \frac{i}{q^2m^2\mathrm{\Pi }(q,m,T)+iϵ}}`$ (2.10)
where $`n(q^0)`$ is the Bose distribution $`\left(e^{\beta |q^0|}1\right)^1`$. This implies

$$D_{12}(q)=2ϵ(q^0)n(q^0)\mathrm{Im}\left(\frac{1}{q^2m^2\mathrm{\Pi }(q,T,m)+iϵ}\right)$$
(2.11)
In the case of large-$`N`$ scalar theory, the self-energy $`\mathrm{\Pi }(q,m,T)`$ is real and independent of the momentum. In consequence, we have

$$\mathrm{Im}\left(\frac{1}{q^2m^2\mathrm{\Pi }+iϵ}\right)=\pi \delta (q^2m^2\mathrm{\Pi }(m,T))$$
(2.12)
Hence we find

$$\mathrm{\Delta }=\frac{\mathrm{\Pi }(m,T)}{2}M_T(m^2+\mathrm{\Pi }(m,T))$$
(2.13)
where $`M`$, $`N_T`$ and $`M_T`$ are defined as :

$$M(m^2)=\frac{1}{2}\frac{d^dq}{(2\pi )^d}\frac{i}{q^2m^2+iϵ}=\frac{1}{2}\frac{d^dq}{(2\pi )^d}2\pi \delta ^+(q^2m^2)$$
(2.14)
$$N_T(m^2)=\frac{d^dq}{(2\pi )^d}\mathrm{\hspace{0.17em}2}\pi \delta ^+(q^2m^2)n(q^0)$$
(2.15)
$$M_T(m^2)=M(m^2)+N_T(m^2)$$
(2.16)
To make use of (2.13), we need to express all the bare parameters of the Lagrangian in terms of their renormalised counterparts.

## 3 Renormalisation

We choose a renormalisation scheme that renders $`\mathrm{\Delta }`$ as simple as possible. Since we are dealing with exact equations, it is suitable to define the renormalised mass through the Dyson equation, which is valid to all orders. At $`T=0`$ and in the on-shell renormalisation scheme, define the renormalised mass as

$$m_R^2=m^2+\mathrm{\Pi }(q^2=m_R^2,m,T=0)$$
(3.1)
For large-$`N`$ scalar theory, the Dyson equation reads

$$\mathrm{\Pi }(m,T=0)=\lambda M(m_R^2)$$
(3.2)
so we have

$$m_R^2=m^2+\lambda M(m_R^2)$$
(3.3)
since $`N_T`$ vanishes at $`T=0`$. The integration in (2.14) gives

$$M(m_R^2)=\frac{\mathrm{\Gamma }(1d/2)}{2(4\pi )^{d/2}}m_R^{d2}$$
(3.4)
At nonzero temperature, it is convenient to define a thermal mass

$$\delta m_R^2=\lambda (M_T(m_R^2+\delta m_R^2)M(m_R^2))$$
(3.5)
so (3.3) implies

$$m_R^2+\delta m_R^2=m^2+\lambda M_T(m_R^2+\delta m_R^2)$$
(3.6)
Comparing (3.6) with (3.3), we see that the definition (3.5) is a natural one. In what follows, we shall be interested in the case of massless theories ($`m=0`$ and therefore $`m_R=0`$) because there one encounters the infrared problems of thermal perturbation theory. Since $`M(0)=0`$, (3.5) reduces to the following integral equation for the thermal mass

$$\delta m_R^2=\lambda M_T(\delta m_R^2)$$
(3.7)
and we obtain

$$\mathrm{\Delta }=\frac{1}{2}\lambda M_T^2(\delta m_R^2)=\frac{1}{2}\frac{\delta m_R^4}{\lambda }$$
(3.8)
Beside mass renormalisation, we also have to renormalise the coupling constant: because $`\lambda `$ has dimension $`M^{4d}`$, we introduce a renormalisation scale $`\mu `$ and define a dimensionless renormalised coupling by

$$\lambda =\mu ^{4d}\lambda _RZ(\lambda _R)$$
(3.9)
where $`Z(\lambda _R)`$ is a combination of wave-function and vertex renormalisation factors. It is a function of $`\lambda _R`$ only because it is dimensionless. In the case of massive theories, one usually uses the on-shell renormalisation scheme

$$\lambda _R=\lambda _R(\mu )=\mu ^{d4}\lambda +\lambda \lambda _RM^{}(m_R^2)$$
(3.10)
where the prime denotes differentiation of $`M`$ with respect to its argument. However, if we consider massless theories, this definition is not useful since $`M^{}(0)`$ is infrared divergent. Instead, define

$$\lambda _R(\mu )=\mu ^{d4}\lambda +\lambda \lambda _R(\mu )M^{}(\mu ^2)$$
(3.11)
With (3.4) and (3.9), this gives

$$Z(\lambda _R)=\left(1\frac{\lambda _R\mathrm{\Gamma }(3d/2)}{(4d)(4\pi )^{d/2}}\right)^1$$
(3.12)
In fact, the renormalisation scheme defined by (3.11) is the modified minimal subtraction scheme ($`\overline{\text{MS}}`$) since (3.12) has the expansion

$$Z=\frac{1}{1+\frac{\lambda _R}{16\pi ^2}\left(\frac{1}{d4}+2(\gamma \mathrm{log}4\pi )\right)}+𝒪(d4)$$
(3.13)
If we introduce the $`\beta `$-function as usual as

$$\beta (\lambda _R)=\mu \frac{\lambda _R}{\mu }|_\lambda $$
(3.14)
and if we differentiate the left and right-hand sides of (3.9) with respect to $`\mu `$ at fixed unrenormalised coupling, we get

$$0=4d+\beta (\lambda _R)\left(\frac{1}{\lambda _R}+\frac{d}{d\lambda _R}\mathrm{log}Z(\lambda _R)\right)$$
(3.15)
From this equation and from (3.12),

$$\beta (\lambda _R)=\lambda _R(d4)+\frac{\lambda _R^2\mathrm{\Gamma }(3d/2)}{(4\pi )^{d/2}}$$
(3.16)
We can also differentiate (3.9) with respect to $`\lambda _R`$ at fixed $`\mu `$ and then use (3.15) to obtain

$$\frac{d\lambda }{\lambda }=d\lambda _R\left[\frac{d}{d\lambda _R}\mathrm{log}Z(\lambda _R)+\frac{1}{\lambda _R}\right]=\frac{d4}{\beta (\lambda _R)}d\lambda _R$$
(3.17)
The pressure being a physical quantity, it is the same before and after renormalisation: $`P_R(T,\lambda _R)=P(T,\lambda )`$. We can therefore combine (3.8), (3.12) and (3.17) to turn (2.6) into the renormalised formula

$$\frac{}{\lambda _R}\widehat{P}_R(T,\lambda _R)=\frac{4d}{\beta (\lambda _R)}\mathrm{\Delta }=\frac{\mu ^{d4}}{2\lambda _R}\delta m_R^4$$
(3.18)
This equation is of course consistent with its unrenormalised counterpart (2.6) since $`\lambda =\mu ^{4d}\lambda _R`$ to lowest order and $`\widehat{P}_R=\widehat{P}`$. It also implies that the question whether we get the same answer for the thermal pressure if use perturbation theory before or after we take the $`d4`$ limit reduces to how $`\delta m_R^4`$ behaves depending on the chosen procedure. In the next section, we will find that if we make an expansion in power series of the coupling (ie use standard diagrammatic Feynman perturbation theory) in $`d`$ space-time dimensions, not only we get uncancelled $`1/(d4)`$ divergences from order $`\lambda _R^4`$ onwards but even the finite second order term is incorrect when we compare it to the exact alternative expansion we obtain in $`4`$ dimensions.

## 4 The thermal mass

We now study the structure of the mass shift when the coupling constant is chosen as expansion parameter either before or after we go to $`4`$ dimensions. Our aim is to see what kind of expression we get in a perturbative calculation of the thermal mass and how resummation works to remove infrared $`1/(d4)`$ terms and obtain the unrenormalised non perturbative equation (3.7).

### 4.1 Naive perturbative expansion

In scalar field theories, it is known that the dominating infrared contributions to the self-energy come from the so-called ring diagrams and their extension called super-daisy or Hartree-Fock diagrams (see figure 1).

First compute the value of a ring diagram with $`L`$ external loops attached using the imaginary-time formalism :

$$\mathrm{\Pi }^{(L)}(d,T)=\frac{\lambda T}{2}m_T^{2L}(1)^L\left(2\underset{n=1}{\overset{\mathrm{}}{}}\frac{d^{d1}p}{(2\pi )^{d1}}\frac{1}{(p^2+\omega _n^2)^{L+1}}+\frac{d^{d1}p}{(2\pi )^{d1}}\frac{1}{p^{2L+2}}\right)$$
(4.1)
where $`\omega _n`$ denote the Matsubara frequencies and $`m_TM_T(0)^{1/2}`$ the one-loop thermal mass. According to the rules of dimensional regularisation, the last term in this expression, the static mode contribution, vanishes since the integrand does not contain any dimensionful quantity. However, the sum over $`L`$ of the static mode contributions is

$$\underset{L=0}{\overset{\mathrm{}}{}}\frac{\lambda T}{2}m_T^{2L}(1)^L\frac{d^{d1}p}{(2\pi )^{d1}}\frac{1}{p^{2L+2}}=\frac{\lambda T}{2}\frac{d^{d1}p}{(2\pi )^{d1}}\frac{1}{p^2+m_T^2}=\frac{\lambda Tm_T}{8\pi }$$
(4.2)
The perturbative expansion of the right hand side of (4.2) is inappropriate since it is nonanalytic in the coupling (it is proportional to $`\lambda ^{3/2}`$) and naive use of dimensional regularisation fails here. We can therefore conclude that static mode contributions only make sense when they have been resummed.
We now turn to the non-static modes. In the spirit of dimensional regularisation, integrals are initially computed in a dimension where they converge and then analytically continued to four space-time dimensions. In $`d`$ dimensions, the non-static part of the self energy is

$$\mathrm{\Pi }_{n.s.}^{(L)}(d,T)=\lambda Tm_T^{2L}(1)^L\frac{(2\pi T)^{d2L3}}{(4\pi )^{\frac{d1}{2}}}\frac{\mathrm{\Gamma }(L+1\frac{d1}{2})}{L!}\zeta (2L+3d)$$
(4.3)
where $`\zeta `$ is the standard Riemann zeta-function. As usual, we get to lowest order

$$m_T^2=\mathrm{\Pi }_{n.s.}^{(0)}(4,T)=\frac{\lambda T^2}{24}$$
(4.4)
Since $`\zeta (z)`$ is analytic in the whole complex plane except for $`z=1`$ where it has a simple pole, $`\mathrm{\Pi }_{n.s.}^{(L)}(d,T)`$ is finite for all values of $`L`$ when $`d4`$, except for $`L=1`$. In the latter case, we have

$$\mathrm{\Pi }_{n.s.}^{(1)}(d,T)=\frac{\lambda m_T^2}{16\pi ^2(d4)}\frac{\lambda m_T^2}{16\pi ^2}\left(\frac{\gamma }{2}\frac{1}{2}\mathrm{log}4\pi \mathrm{log}T\right)+𝒪(d4)$$
(4.5)
It means we still get a $`\lambda ^2/(d4)`$ pole even if we have resummed all the static contributions.
Since it is known that resumming all ring diagrams provides a finite result , the resummation of all the $`\mathrm{\Pi }^{(L)}(d,T)`$ for $`L>1`$ should cancel out the $`1/(d4)`$ pole coming from $`\mathrm{\Pi }^{(1)}(d,T)`$ :

$$\underset{L=2}{\overset{\mathrm{}}{}}\mathrm{\Pi }_{n.s.}^{(L)}(d,T)=\frac{\lambda T}{(4\pi )^{\frac{d1}{2}}}\underset{n=1}{\overset{\mathrm{}}{}}(2\pi nT)^{d3}\underset{L=2}{\overset{\mathrm{}}{}}\left(\frac{m_T}{2\pi nT}\right)^{2L}\frac{(1)^L\mathrm{\Gamma }(L\frac{d3}{2})}{L!}$$
(4.6)
If we use the Taylor expansion

$$(1+x)^{\frac{d3}{2}}=1+\left(\frac{d3}{2}\right)x+\underset{L=2}{\overset{\mathrm{}}{}}\frac{1}{L!}\frac{\mathrm{\Gamma }\left(\frac{d1}{2}\right)}{\mathrm{\Gamma }\left(\frac{d1}{2}L\right)}x^L$$
(4.7)
and the identity

$$\mathrm{\Gamma }(x)\mathrm{\Gamma }(1x)=\frac{\pi }{\mathrm{sin}(\pi x)}$$
(4.8)
we find, by summing over $`n`$,

$`{\displaystyle \underset{L=2}{\overset{\mathrm{}}{}}}\mathrm{\Pi }_{n.s.}^{(L)}(d,T)={\displaystyle \frac{\lambda T\pi }{(4\pi )^{\frac{d1}{2}}\mathrm{\Gamma }\left(\frac{d1}{2}\right)}}[{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(\omega _n^2+m_T^2)^{\frac{d3}{2}}+(2\pi T)^{d3}\zeta (3d)`$
$`+\left({\displaystyle \frac{d3}{2}}\right)m_T^2(2\pi T)^{d5}\zeta (5d)]+𝒪(d4)`$ (4.9)
By expanding this expression around $`d=4`$, we obtain

$`{\displaystyle \underset{L=2}{\overset{\mathrm{}}{}}}\mathrm{\Pi }_{n.s.}^{(L)}(d,T)=\lambda M_T(m_T^2)+{\displaystyle \frac{\lambda T}{8\pi }}m_Tm_T^2{\displaystyle \frac{\lambda m_T^2}{16\pi ^2(d4)}}`$
$`+{\displaystyle \frac{\lambda m_T^2}{16\pi ^2}}\left({\displaystyle \frac{\gamma }{2}}{\displaystyle \frac{1}{2}}\mathrm{log}4\pi \mathrm{log}T\right)+𝒪(d4)`$ (4.10)
As expected, resummation indeed removes the pole: it is cancelled out by the rest of the perturbative series. Hence, the resummed result is merely the one-loop self-energy diagram in which the free propagator has been replaced by the first-order dressed one:

$$\underset{L=0}{\overset{\mathrm{}}{}}\mathrm{\Pi }^{(L)}(d,T)=\lambda M_T(m_T^2)+𝒪(d4)$$
(4.11)
Because it comes from the Matsubara sum, it is not clear whether the $`\lambda ^2/(d4)`$ divergence was an infrared or an ultraviolet one. In order to find out, it is interesting to study $`\mathrm{\Pi }^{(1)}(d,T)`$ in the real-time formalism: this is done in appendix A.1.
In this section, we have so far only been dealing with quantities involving the unrenormalised coupling so we still need to renormalise (3.7) to make sure we have properly taken care of UV divergences before jumping to conclusions about possible remaining divergences. Expressing bare quantities in terms of renormalised ones, (3.7) becomes the following integral equation :

$$\delta m_R^2=\mu ^{4d}\lambda _R\left[\stackrel{~}{M}(\mu ^2,\delta m_R^2)+N_T(\delta m_R^2)\right]$$
(4.12)
where $`\stackrel{~}{M}(\mu ^2,\delta m_R^2)M(\delta m_R^2)\delta m_R^2M^{}(\mu ^2)`$ is finite when $`d=4`$. This equation is exact and therefore corresponds to the full summation of all the Hartree-Fock diagrams. By studying the ring diagram contributions to the thermal mass in the ITF and RTF, we have seen that an expansion of $`N_T(\delta m_R^2)`$ in power series of $`\delta m_R^2`$ gives rise to infrared divergences of the form $`\delta m_R^2/(d4)`$ although $`N_T`$ itself is perfectly finite when $`d=4`$. This is because $`N_T`$ is nonanalytic in $`\delta m_R^2`$ at $`\delta m_R^2=0`$ (namely, it involves terms which are proportional to $`\delta m_R\mathrm{log}\delta m_R`$) and so the Feynman perturbative expansion in powers of $`\lambda _R`$ is illegal since $`\delta m_R\lambda _R`$ to lowest order. In the next section, we investigate the exact form of $`N_T`$ and its nonanalyticity in the coupling. This will lead us to the conclusion that the $`\delta m_R^2/(d4)`$ term in the perturbative expansion of $`N_T`$ is effectively replaced with $`\delta m_R^2\mathrm{log}(\delta m_R^2)`$ in the corresponding exact expression (see (4.18)). In section 4.3, we shall see that the divergent term in a perturbative expansion of $`N_T(\delta m_R^2)`$ gives rise to infrared divergences in the renormalised perturbative series from order $`\lambda _R^4`$ onwards.

### 4.2 Alternative expansion of $`N_T`$

In order to study its dependence in the the coupling, we will now expand $`N_T`$ in another and more sensible way than the naive perturbative expansion :

$`N_T(m^2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^{d1}p}{(2\pi )^{d1}}\frac{p^{d2}}{\sqrt{p^2+m^2}}\frac{1}{\mathrm{exp}(\beta \sqrt{p^2+m^2})1}}`$ (4.13)
$`=`$ $`{\displaystyle \frac{T^2}{(4\pi )^{\frac{d1}{2}}\mathrm{\Gamma }(\frac{d1}{2})}}{\displaystyle _0^{\mathrm{}}}𝑑u{\displaystyle \frac{u^{d2}}{\sqrt{u^2+y^2}}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\mathrm{exp}(k(u^2+y^2)^{1/2})`$
where $`y\beta m`$. We evaluate the sum using the Mellin summation formula (see appendix A.2)

$$\underset{k=1}{\overset{\mathrm{}}{}}f(k)=\frac{1}{2i\pi }_{Ci\mathrm{}}^{C+i\mathrm{}}𝑑z\zeta (z)\varphi (z)$$
(4.14)
where $`\varphi (z)`$ is the Mellin transform of $`f(k)`$ and $`C`$ is such that all the poles of the integrand are to the left of the vertical integration line. We obtain

$$N_T(m^2)=\frac{T^2}{2i\pi (4\pi )^{\frac{d1}{2}}\mathrm{\Gamma }(\frac{d1}{2})}_{Ci\mathrm{}}^{C+i\mathrm{}}𝑑z_0^{\mathrm{}}𝑑uu^{d2}(u^2+y^2)^{\frac{z+1}{2}}\mathrm{\Gamma }(z)\zeta (z)$$
(4.15)
with $`C>d2`$. We can evaluate the integral over $`u`$ and the remaining integral becomes

$$\frac{1}{2}_{Ci\mathrm{}}^{C+i\mathrm{}}𝑑zy^{d2z}\frac{\mathrm{\Gamma }(z)}{\mathrm{\Gamma }(\frac{1+z}{2})}\mathrm{\Gamma }(d5/2)\mathrm{\Gamma }\left(\frac{2+zd}{2}\right)\zeta (z)$$
(4.16)
The integrand has single poles at $`z=0,1`$ and $`z=d2k`$ where $`k=1,2,3,\mathrm{}`$. The contour of the $`z`$ integration in the integral (4.15) may be closed in the left half plane since the integrand vanishes asymptotically on the arc at infinity. We now have to compute all the residues of the poles of the integrand. The contribution to $`N_T`$ from the pole at $`z=2d`$ corresponds to the first order thermal mass $`m_T^2`$. From the pole at $`z=1`$, we get $`Tm/(8\pi )`$ which corresponds to the sum of all the static modes we had in the ITF (see (4.2)) and is nonanalytic in the coupling ($`m_T\lambda _R^{1/2}`$). The sum of the contributions to $`N_T`$ from the poles at $`z=0`$ and $`z=d4`$ is

$$\frac{m^2}{16\pi ^2}\left(\mathrm{log}\frac{m}{4\pi T}+\gamma \frac{1}{2}\right)$$
(4.17)
and this is nonanalytic in the coupling as well due to the term in $`\mathrm{log}m`$. Finally, all the remaining poles have residues of the same form and the complete expression for $`N_T`$ reads (when $`d=4`$)

$`{\displaystyle \frac{4\pi ^2}{T^2}}N_T(m^2)`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{6}}{\displaystyle \frac{\pi m}{2T}}{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{m}{T}}\right)^2\left[\mathrm{log}{\displaystyle \frac{m}{4\pi T}}+\gamma {\displaystyle \frac{1}{2}}\right]`$ (4.18)
$``$ $`{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{m}{T}}\right)^2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n{\displaystyle \frac{(2n)!}{(n+1)!n!}}\zeta (2n+1)\left({\displaystyle \frac{m}{4\pi T}}\right)^{2n}`$
From this expression, we see that $`N_T(\delta m_R^2)`$ contains a term involving $`\mathrm{log}\lambda _R`$ and another one proportional to $`\lambda _R^{1/2}`$ . The fact that an expansion of $`N_T(\delta m_R^2)`$ in power series of the coupling is illegal is quite clear from (4.13) : if we try to expand the integrand in power series of $`m^2`$, the expansion parameter is actually $`m^2/p^2`$ and this obviously does not make sense since the momentum integration runs down to $`0`$. If we do it anyway, we obtain what we got in perturbation theory and we get an isolated $`1/(d4)`$ term at order $`\lambda _R^2`$. Moreover, we also have

$$\stackrel{~}{M}(\mu ^2,\delta m_R^2)=\frac{\delta m_R^2}{32\pi ^2}\left(\mathrm{log}\frac{\delta m_R^2}{\mu ^2}1\right)+𝒪(d4)$$
(4.19)
So we find

$`\delta m_R^2`$ $`=`$ $`\lambda _RT^2\{{\displaystyle \frac{1}{24}}{\displaystyle \frac{\delta m_R}{8\pi T}}{\displaystyle \frac{1}{16\pi ^2}}\left({\displaystyle \frac{\delta m_R}{T}}\right)^2[\mathrm{log}{\displaystyle \frac{\mu }{4\pi T}}+\gamma ]`$ (4.20)
$``$ $`{\displaystyle \frac{1}{16\pi ^2}}\left({\displaystyle \frac{\delta m_R}{T}}\right)^2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1)^n{\displaystyle \frac{(2n)!}{(n+1)!n!}}\zeta (2n+1)\left({\displaystyle \frac{\delta m_R}{4\pi T}}\right)^{2n}\}`$
One can solve this equation iteratively to any desired accuracy in $`\lambda _R`$ and the first few terms in the expansion are

$`{\displaystyle \frac{\delta m_R^2}{T^2}}={\displaystyle \frac{\lambda _R}{24}}`$ $``$ $`{\displaystyle \frac{\lambda _R^{3/2}}{16\pi \sqrt{6}}}+\left(3\gamma \mathrm{log}{\displaystyle \frac{\mu }{4\pi T}}\right){\displaystyle \frac{\lambda _R^2}{384\pi ^2}}`$ (4.21)
$``$ $`\left(12\gamma 2\mathrm{log}{\displaystyle \frac{\mu }{4\pi T}}\right){\displaystyle \frac{\lambda _R^{5/2}}{1024\pi ^3\sqrt{2/3}}}+𝒪(\lambda _R^3)`$
It is evident that $`\delta m_R^2`$ does not have an expansion in powers of $`\lambda _R`$, but rather in powers of $`\lambda _R^{1/2}`$.

### 4.3 Renormalised perturbative series

We express $`\delta m_R^2`$ in terms of the renormalised coupling and the renormalisation scale. From (3.9) and (3.11)

$$\lambda =\lambda _R\mu ^{4d}\left[\underset{k=0}{\overset{\mathrm{}}{}}()^k\mu ^{k(4d)}\lambda _R^k\left(M^{}(\mu ^2)\right)^k\right]$$
(4.22)
so diagrams corresponding to different fixed orders in $`\lambda `$ will yield contributions both to the same order in $`\lambda _R`$ and to higher orders. Evidently, to lowest order, we have $`\lambda =\lambda _R`$ so $`\delta m_R^2=\lambda _RT^2/24`$. To order $`\lambda _R^2`$, both the one and two-loop diagrams contribute and the $`1/(d4)`$ divergence coming from $`\mathrm{\Pi }^{(1)}`$ exactly cancels the one that results from the renormalisation of the coupling in the one-loop diagram $`\mathrm{\Pi }^{(0)}`$. We are therefore left with a finite expression up to order $`\lambda _R^2`$:

$$\delta m_R^2=\frac{\lambda _RT^2}{24}\frac{\lambda _R^2T^2}{384\pi ^2}\left[\mathrm{log}\frac{\mu }{4\pi T}+\gamma \right]$$
(4.23)
Having obtained an exact expression for $`\delta m_R^2`$ in terms of $`\lambda _R`$, we investigate the problems associated with first expanding it as a perturbation series in $`\lambda _R`$. Comparing with (4.21), we see that there are missing terms. This is because we also need to include the contribution from the non-vanishing sum of all the static modes : this contribution, which also yields all the terms proportional to an odd power of $`\lambda _R^{1/2}`$ in (4.21), is missed out in the perturbative expansion and only appears when we resumm the series.
After tedious calculations, divergences cancel out again at order $`\lambda _R^3`$ if we take into account all the diagrams that contribute to the thermal mass at this order (they are represented in figure 2): diagram D in figure 2, which is not a ring diagram, rescues the situation by cancelling divergences coming from diagrams A and B and the renormalisation of the bare coupling in these diagrams. It means that, if we only consider ring diagrams, we get uncancelled divergences from order $`\lambda _R^3`$ onwards in the renormalised perturbation series.

From order $`\lambda _R^4`$ onwards, it is not possible to get rid of divergences any more even if we sum up all diagrams that contribute to the thermal mass at this order : worse and worse divergences crop up in the renormalised perturbative series when going to higher orders. For instance, at order $`\lambda _R^4`$, we have

$$\delta m_R^2|_{𝒪(\lambda _R^4)}=\frac{\lambda _R^4T^2}{49152\pi ^6(d4)^3}+𝒪((d4)^2)$$
(4.24)
In conclusion, we have seen in this section that the nonanalyticity of the thermal mass in the coupling causes the failure of the corresponding Feynman perturbative expansion. This failure has two consequences : firstly, infrared divergences appear in the renormalised series from order $`\lambda _R^4`$ and, secondly, the term proportional to $`\lambda ^{1/2}`$ and its contributions to all order in $`\lambda _R`$ are missed out.

## 5 Conclusion

The equation (3.18) derived for the pressure is exact, or non perturbative. Its right-hand side involves $`\delta m_R^4`$ and we have seen in the previous section that we get uncancelled $`1/(d4)`$ infrared singularities if we expand it in powers of $`\lambda _R`$ before we take the $`d4`$ limit. On the other hand, if we take the limit $`d4`$ on the exact expression (4.13), then we get an expression which cannot be expanded in power series of the coupling any more since a coefficient in the expansion would be infinite. It tells us that we have to use the more involved expansion outlined in section 4. In fact, before we take the limit $`d4`$, nothing prevents us from expanding in power series of the coupling even if this expansion is illegal in four space-time dimensions. We could therefore conclude by saying that the correct procedure is to work in four space-time dimensions before expanding the physical quantity under consideration in a suitable way. In this paper, we studied both the perturbative and the exact expression of the finite temperature contribution to the thermal mass and we saw that the infrared divergent $`\lambda ^2/(d4)`$ term in the unresummed perturbative series is effectively replaced with $`\lambda \mathrm{log}\lambda `$ in the corresponding exact expansion (see section 4.2). The same kind of mechanism may be related to the severe infrared problems encountered in thermal QCD. In the case of the QCD thermal pressure, it is known that diagrams with magnetostatic gluon loops become divergent from order $`g^6`$ in perturbation theory. From what we saw in this work in the case of large-$`N`$ scalar theory, it seems sensible to assume that this $`g^6/(d4)`$ divergence in the perturbative QCD pressure corresponds to a term involving $`g^6\mathrm{log}g`$ in the exact expression for the pressure. This will be the subject of subsequent work.

## Acknowledgments

I am indebted to P.V.Landshoff for reading the manuscript and for invaluable discussions. The work is supported by the Cambridge European Trust, Christ’s College and by the Department of Applied Mathematics and Theoretical Physics.

## Appendix A Appendix

### A.1 Ring diagrams in the real-time formalism

In the real-time formalism, the self-energy also has a matrix structure and the matrix element relevant to the evaluation of $`\mathrm{\Pi }^{(1)}(d,T)`$ is

$$\mathrm{\Pi }_{11}^{(1)}(d,T)=i\frac{\lambda }{2}M_T^2\frac{d^dp}{(2\pi )^d}\left[D_{11}^2(p)D_{12}^2(p)D_{21}^2(p)\right]$$
(A.1)
Using (2.9), we get

$$D_{11}^2(p)D_{12}^2(p)D_{21}^2(p)=\left(\frac{i}{p^2+iϵ}\right)^2n(p^0)\left[\frac{1}{(p^2+iϵ)^2}\frac{1}{(p^2iϵ)^2}\right]$$
(A.2)
The finite temperature part of the latter expression can be written as

$$n(p^0)\frac{}{m^2}\left(\frac{1}{p^2m^2+iϵ}\frac{1}{p^2m^2iϵ}\right)|_{m^2=0}$$
(A.3)
which is the derivative of the one-loop diagram for the thermal mass. We therefore obtain

$$\frac{d^dp}{(2\pi )^d}\left[D_{11}^2(p)D_{12}^2(p)D_{21}^2(p)\right]=i\pi (d3)\frac{d\mathrm{\Omega }_{d1}}{(2\pi )^d}_0^{\mathrm{}}𝑑p\frac{p^{d5}}{e^{\beta p}1}$$
(A.4)
Finally, if we use the representation

$$\zeta (z)=\frac{1}{\mathrm{\Gamma }(z)}_0^+\mathrm{}𝑑u\frac{u^{z1}}{e^u1}$$
(A.5)
and the relation

$$\zeta (1z)=2(2\pi )^z\mathrm{cos}\left(\frac{z\pi }{2}\right)\mathrm{\Gamma }(z)\zeta (z)$$
(A.6)
we obtain the same result (4.5) as with the imaginary-time formalism. Since the momentum integral in the right-hand side of (A.4) is infrared divergent in four space-time dimensions, it tells us that the corresponding $`1/(d4)`$ divergence really is an infrared one. Similarly, we could also rederive the general expression for a ring diagram with an arbitrary number of external loops and check that it agrees with what we found using the imaginary-time formalism : one just has to take higher derivatives with respect to $`m^2`$ in (A.3).

### A.2 The Mellin summation formula

The Mellin transform of a function $`f(x)`$ is usually defined by

$$\varphi (z)=_0^{\mathrm{}}x^{z1}f(x)𝑑x$$
(A.7)
The transform $`\varphi (z)`$ exists if the integral

$$_0^{\mathrm{}}\left|f(x)\right|x^{k1}𝑑x$$
(A.8)
is bounded for some $`k>0`$, in which case the inverse transform is

$$f(x)=\frac{1}{2i\pi }_{Ci\mathrm{}}^{C+i\mathrm{}}x^z\varphi (z)𝑑z$$
(A.9)
where $`C>k`$. Together with the standard definition of the Riemann Zeta-function, this gives the summation formula

$$\underset{n=1}{\overset{\mathrm{}}{}}f(n)=\frac{1}{2i\pi }_{Ci\mathrm{}}^{C+i\mathrm{}}𝑑z\zeta (z)\varphi (z)$$
(A.10)