Geometrically Confined Thermal Field Theory: Finite Size Corrections and Phase Transitions

Motivated by the recent shocking results from RHIC and LHC that show quark-gluon plasma signatures in small systems, we study a simple model of a massless, noninteracting scalar field confined with Dirichlet boundary conditions. We use this system to investigate the finite size corrections to thermal field theoretically derived quantities compared to the usual Stefan-Boltzmann limit of an ideal gas not confined in any direction. Two equivalent expressions with different numerical convergence properties are found for the free energy in $D$ rectilinear spacetime dimensions with $c\le D-1$ spatial dimensions of finite extent. We find that the First Law of Thermodynamics generalizes such that the pressure depends on direction but that the Third Law is respected. For systems with finite dimension(s) but infinite volumes, such as a field constrained between two parallel plates or a rectangular tube, the relative fluctuations in energy are zero, and hence the canonical and microcanonical ensembles are equivalent. We present precise numerical results for the free energy, total internal energy, pressure, entropy, and heat capacity of our field between parallel plates, in a tube, and in finite volume boxes of various sizes in 4 spacetime dimensions. For temperatures and system sizes relevant for heavy ion phenomenology, we find large deviations from the Stefan-Boltzmann limit for these quantities, especially for the pressure. Further investigation of an isolated system of fields constrained between parallel plates reveals a divergent isoenergetic compressibility at a critical length $L_c\sim1/T$. We have thus discovered a new second order phase transition via a first principles calculation, a transition that is driven by the size of the system.


Motivation and goals
The Relativistic Heavy Ion Collider (RHIC) [1,2] at Brookhaven National Laboratory (BNL) in Upton, New York and the Large Hadron Collider (LHC) [3][4][5] at CERN in Geneva allow scientists to probe the state of matter under unprecedented extreme conditions. Multiple experimental signals confirm the creation of the quark-gluon plasma (QGP) [6][7][8][9][10]; remarkably, these colliders recreate for the first time conditions similar to those existing in the early universe. The future is bright with collider experiments such as the Nuclotron-based Ion Collider fAcility [11] at the Joint Institute for Nuclear Research (JINR) in Dubna and the Facility for Antiproton and Ion Research (FAIR) [12,13] at GSI in Darmstadt that will soon further explore the physics of heavy ion collisions (HIC). This abundance of data has led the HIC field to an era of high precision measurements, which demands a commensurate level of theoretical precision.
Thermodynamics is one of the main avenues of theoretical investigation into the dynamics of the QGP, one that has a connection to the experimental measurements at RHIC and LHC through a number of observables such as identified particle spectra [14]. Understanding the thermodynamics of the QGP created in HIC has been, and still is, the subject of intensive theoretical research for the past three decades. Generally, there are three broad classes of theoretical tools used to explore QGP thermodynamics: weakly coupled thermal field theory [15][16][17][18], Monte Carlo methods [19][20][21][22][23][24][25][26][27][28][29][30][31][32], and AdS/CFT [33]. Despite the long history of investigation into the thermodynamics of the QGP, an important aspect appears to have been overlooked by the HIC community: the effect of finite-sized rather than infinitely-sized systems. The importance of finite-size corrections to the thermodynamics of QGP is of especial significance now as many signals of QGP creation in large ion-ion systems appear to depend on measured particle multiplicity rather than colliding system size. Some examples include collective behavior [34][35][36][37], strangeness enhancement [38][39][40], and quarkonium suppression [41,42], while noticeable absence of appreciable jet quenching demands a thorough interrogation of our understanding of energy loss in small systems [43]. Further, a simple estimate of a quark or gluon mean free path of λ ∼ 1 − 2 fm indicates that even a large colliding Pb+Pb system of radius r ∼ 6 fm is not particularly close to the thermodynamic limit [44].
The present work is motivated by the pressing need for a better understanding of small droplets of QGP, including a quantification of the small system size corrections to the usual approximation of using dynamics derived in systems of infinite size. In terms of QGP phenomenology, our goals are modest. As a first step, we concentrate on the finite size corrections to thermodynamic quantities such as the equation of state (EoS) computed in thermal field theory in the usual Stefan-Boltzmann limit of an ideal gas of infinite size in all directions. As a further simplification, we consider only noninteracting fields. Of the various types of boundary conditions one might use, we focus on Dirichlet boundary conditions (BC) for simplicity and, as we will argue below, because these are the most natural BC when considering a finite-sized QGP.
Despite the simplicity of our model, we will find a suggestive phenomenology. In our derivation of thermodynamic quantities, we use two independent methods that yield analytic formulae with different numerical convergence properties: one which converges exponentially fast when the dimensionless scale of the temperature times the system size is small, T × L 1, and one that converges exponentially fast for T × L 1. A careful application of these two formulae allow us to numerically investigate to arbitrary accuracy the various thermodynamic quantities such as the pressure, energy, entropy, and heat capacity for our model. Of particular note for the heavy ion community, the introduction of Dirichlet instead of periodic BC leads to significant small system corrections to the usual infinite size results of the equation of state even out to fairly large system sizes L ∼ 10/T . However, our work is of potential interest beyond heavy ion physics. In some sense there is a long history of investigating finite size effects in (thermal) field theory. Most famously at T = 0 a noninteracting field between two parallel plates induces an attractive force between the plates, which may be interpreted as a negative pressure, the well-known Casimir effect [45][46][47][48][49]. The original Casimir effect of a free scalar field between two slabs has been extended to systems with nontrivial geometries and to Fermi fields [50]. For T > 0, one encounters a number of interesting fundamental questions, many of which have only begun to be explored recently [51][52][53][54][55]. First, if all spatial dimensions are of finite size, the system may not be in the thermodynamic limit [56]. When not in the thermodynamic limit one is no longer guaranteed that the different statistical ensembles, i.e. microcanonical, canonical, etc., agree [52,[57][58][59]. It is then very important to consider the physical setup for an experimental measurement and use the appropriate ensemble; in particular, isolated systems may behave very differently from systems in contact with a heat bath. Second, phase transitions are notoriously difficult to rigorously identify. This difficulty is only compounded in finite-sized systems [59][60][61][62] with some authors claiming that phase transitions in systems of finite size are impossible. Third, in finite-sized systems, one loses extensivity in thermodynamics [63], where by extensivity we mean that the entropy is positively homogeneous of degree 1, i.e. S(λE, λX 1 , λX 2 , . . .) = λS(E, X 1 , X 2 , . . .). This lack of extensivity then leads to questions about the use of the Tsallis distribution [64,65] and the minimal set of assumptions required for a consistent thermodynamics [66]. Fourth, once a system is of finite size, then thermal fluctuations can be important, which has begun to be explored in the heavy ion hydrodynamic community [67]. Fifth, considering systems with directions of finite extent raises fundamental questions about the applicability of statistical mechanics. It is worth remembering that using the canonical ensemble to describe a thermal system requires only that the system under study be small compared to an approximately infinite heat bath; the system under study cannot be too small for the use of the canonical ensemble. Formally, the thermodynamic limit occurs when N → ∞ and V → ∞ while N/V = constant. While one may think of the thermodynamic limit as one in which all directions are large compared to any other scale in the problem, geometrically confined systems such as a field between parallel plates can satisfy the requirement of the thermodynamic limit while having direction(s) of order or smaller than other scales in the problem. Therefore one can have systems in the thermodynamic limit that nonetheless have large finite size corrections compared to the usual Stefan-Boltzmann limit of a system of infinite size in all directions.
We also find a connection between our work and the study of phase transitions. Phenomenologically, one finds substances for which a phase transition can be driven by a change of pressure at constant temperature, for example the liquid-gas phase transition for water above 0 • C. We are aware of only one example of a phase transition driven by the size of the system [62], in which case the transition is due to the self interactions of the system. What we will show is that for an isolated ideal gas constrained within parallel plates, where by isolated we mean a system with constant energy, the system resists compression until the separation of the plates is on the order of the thermal de Broglie wavelength; at this critical length L c ∼ 1/T the susceptibility diverges and the system collapses. This divergence of the susceptibility indicates a second order phase transition driven by the size of the system, which is conjugate to the pressure on the system. Unlike other works that draw a connection between Bose-Einstein condensation between parallel plates [68] and a first order phase transition in a finite volume box [69], the phase transition found here is second order and also exists for Fermionic fields.

Geometric confinement for HIC
There have been a number of studies of finite size effects using periodic boundary conditions [70][71][72][73][74][75][76][77][78]. Consider, however, that a boundary-less manifold, e.g. a three-dimensional sphere, is entirely decoupled from the rest of the universe: there is no possibility for any signal or information to come through the QGP and reach the detectors. Therefore, spatial boundary conditions-other than periodic ones-should be considered for a more realistic approach to the finite size corrections of the EoS 1 .
Consider now a QGP system inside of which the quarks and gluons are color deconfined and propagate relatively freely while outside of this QGP system the quarks and gluons are color confined into hadrons. Inside the system the quark and gluon fields, then, are weakly coupled while these same quark and gluon fields are strongly coupled outside the system. There is further some very small region of space over which the fields transform from weakly to strongly coupled. For our purposes here, we are most interested in the dynamics inside the QGP system. And if we approximate the small transition region from weak to strong coupling as a decoupling, we must impose a boundary condition that prevents an inside weakly coupled field from propagating outside of the geometric region defining the QGP system. We refer readers to [79][80][81] for related investi-gations that demonstrate the need for such more realistic boundaries. We also refer to [82][83][84][85][86] for somewhat related investigations, as well as to [87] regarding the importance of accounting for finite size effects in the different context of proton-proton collisions. While such a decoupling might seem irrelevant from a perturbative point of view, since this change of degrees of freedom across the two regions must involve nonperturbative physics, our system appears to be analogous to the Quantum Electrodynamics (QED) Casimir effect [49,88], which can be reproduced by imposing Dirichlet Boundary Conditions (DBCs). Such boundary conditions indeed follow the requirement for the fields to vanish at the boundary, and any two-point correlation function connecting the inside part to the boundary would identically vanish, thereby decoupling both sides of the boundary. Thus, DBCs for the QGP can avoid the propagation of QGP particles across the boundary, keeping the relatively free quarks and gluons geometrically confined inside the volume of the QGP system, while at the same time allowing for other fields, e.g. electromagnetic, to freely propagate out from the confined system.
We then introduce the notion of geometric confinement for such a system, by implementing appropriate spatial compactifications (depending on the geometry to be characterized). In each of the compactified directions we impose DBCs 2 . It should be noted that such a boundary is not a material boundary, but rather a thin layer subdividing different regions of space with different degrees of freedom. And since we are only interested in the bulk physics away from the boundary (avoiding then possible technical complications [89,90]), we will not consider the microscopic nature of the boundary. As a result, the physical space is separated into two distinct regions: An inside part of nearly free quarks and gluons characterizing the quark-gluon plasma, and an outside part (of nearly free hadrons composed of strongly coupled quarks and gluons). For the sake of our investigations in this manuscript we ignore the details of color confinement and the microscopic nature of the boundary. Moreover, it should be noted that while DBCs may be implemented in more physically realistic geometries, for the sake of simplicity we choose to work here with simple rectangular geometries, that is planar pairwise parallel spatial boundaries forming a cuboid cavity. In three spatial dimensions, we consider two infinite parallel plates, the infinite rectangular tube, and the finite volume box.
We stress that our concept of geometric confinement is in no way related to actual color confinement, for it barely reproduces only one consequence of it (the fact that quark and gluon fields should not propagate outside of the QGP); we do not address the fundamental mechanism of color confinement. In addition, our concept is different from the MIT bag model [91] since our model is not meant to be relevant for strongly coupled fields, but for weakly coupled ones.

A model for first investigations
Recalling that the EoS for a noninteracting gas of gluons is the same, up to a group theory prefactor, as the EoS given by a massless scalar field, for simplicity we choose to consider a single neutral noninteracting massless scalar field at finite temperature under geometric confinement. Further, we will work in the canonical ensemble: we will compute the partition function using thermal field theoretic methods. In heavy ion collisions two nuclei collide in vacuum. Even though the quark-gluon plasma created in heavy ion collisions is an isolated system, we choose to work in the canonical ensemble for two reasons. First, we may more readily connect our results to lattice QCD calculations [19][20][21][22][23][24][25][26][27][28][29][30][31][32]; and second, for simplicity.
In the following, we will keep the mass of the field nonzero for pure convenience during the calculations. At a later stage, we will, however, take the massless limit which is the present case of interest. We will work in D spacetime dimensions of which an arbitrary number c will be compactified spatial dimensions with DBCs. We will leave the noncompactified dimensions, if any, in the usual infinite Euclidean form. Each such compactified spatial direction will then have a distinct compactification length L i , i = 1, . . . c ≤ D − 1. The L i 's do not necessarily have to be equal, which allows us to investigate systems of asymmetric sizes. Our starting point is therefore the free, Wick rotated Euclidean action where beside the usual trace induced periodic boundary condition along the Euclideanized, Wickrotated τ direction, and the new geometric confinement induced DBCs along the z i spatial directions, the x coordinates will be momentarily compactified around circles of radius R, i.e. with periodic boundary conditions. At the end of the calculation we will take R → ∞. The free La- where we are working in = c = k B = 1 natural units.

Deriving the partition function
In this subsection, we derive the partition function Z(T, {L i }) for a single, neutral, noninteracting, massless scalar field that is geometrically confined within c ≤ D − 1 spatial dimensions. Formally, the partition function in a theory with Hamiltonian operatorĤ and no globally conserved charges is obtained from the trace of the density matrix where β = 1/T is the inverse temperature, V the spatial volume of the system, and the trace represents a summation over all possible physical states. We will employ the Matsubara imaginary time formalism [92] in order to compute the partition function using path integral techniques [93]. For more details on this formalism as well as on thermal field theory, we refer readers to [94][95][96][97][98][99]. The usual procedure for expressing the partition function as a path integral leads to a path integral with a periodic boundary condition on the temporal line which, up to an irrelevant constant, reads We now extend the above to a manifold with c compactified spatial dimensions. The procedure will require DBCs for the compactified spatial dimensions, in addition to the periodic boundary condition required by the trace operation. The derivation of the analogue of eq. (2.2) with compactified spatial dimensions closely follows that of the spatially noncompactified case, which we call the Stefan-Boltzmann limit, the details of which can be seen in the aforementioned textbooks. For pedagogical reasons we will first compactify only one dimension (c = 1), which gives the canonical ensemble description for a noninteracting scalar field at temperature T constrained between two infinite parallel plates separated by a distance L 1 . Extending the result to c compactified dimensions will be straightforward. Following the usual procedure, we start by decomposing our field into the relevant Fourier modes. The dimensions that will ultimately be left noncompactified are, as usual, momentarily compactified onto circles of identical sizes R (which means we impose periodic boundary conditions), while the dimension we wish to permanently compactify is done so along a finite interval with length L 1 . We choose a convenient normalization constant (see appendix (B.1) for more details) and express the Fourier decomposition of our field as where the components of x are the spatial components in the non permanently compactified dimensions, z 1 is the spatial component in the compactified dimension, andφ n, 1 (ω k ) are the dimensionless Fourier modes. The Matsubara, the compactified spatial dimension, and the non permanently compactified spatial dimensions related frequencies are given by ω n = 2πnT , ω 1 = π 1 /L 1 and ω k = 2πn k /R, respectively, and we also refer to [100] for more details on the derivation of the modes given DBCs. The latter frequencies are to be replaced by a D − 2 dimensional continuous momentum k (with corresponding momentum integrals), in the asymptotically large R limit.
Setting c = 1 and employing eq. (1.2) for the Lagrangian, the partition function (2.2) becomes, after an integration by parts, where, among others, the periodic boundary condition from the trace is imposed through setting φ cond. such that φ(β, z 1 , x . Substituting eq. (2.3) into eq. (2.4), and simplifying the argument of the exponential by performing the integrations (see appendix (B.1) for more details) [97][98][99], we obtain As in the Stefan-Boltzmann case, we are faced with an issue of double counting when performing explicitly the path integral. This problem can be accounted for, in the noninteracting case, by separating the Fourier modes into real and imaginary parts and since the field is required to be real, we obtain the following restrictions from which one may choose a set of physically relevant independent φ-variables over which to integrate. Following the standard procedure [97][98][99], one may then perform the infinite set of Gaussian integrals (dropping any T -and L i -independent factors [101]) to obtain the partition function of a free scalar field in D−1 spatial dimensions with c = 1 geometrically confined dimension.
However, since thermodynamic quantities are straightforwardly related to the logarithm of the partition function, we will find the more useful quantity to be Formally, eq. (2.8) is our final result for the logarithm of the partition function of a single neutral noninteracting scalar field in between two parallel plates. Extending (2.8) to arbitrary c ≤ D − 1 compactified spatial dimensions is relatively straightforward, and the corresponding logarithm of the partition function with c geometrically confined dimensions is therefore given by (2.9) where we still have to send R to infinity.
In Appendix (A.1) and (A.2) we present two independent methods of evaluating eq. (2.9). We have confirmed numerically that both methods yield the same results. The two methods are mutually complementary as they naturally yield results with different numerical convergence properties. More precisely, the usual method yields a result that is always better suited for small values of the dimensionless parameters (T × L i ). The alternative method yields a result better suited for high values of these dimensionless parameters and thus to recover the usual Stefan-Boltzmann limits. Moreover, the usual result explicitly isolates the T -independent contributions, and one can pass from one compactified dimension case to the next in an iterative manner; therefore the usual method explicitly yields the known zero temperature Casimir pressure.

Evaluating the free energy
We now present our final results for the free energy in the massless limit using both eq. (A. 19) and results from eq. (A. 21). Recall that the free energy density is given by f ≡ F/V , F being the total free energy, and therefore reads

Case I: Two infinite parallel plates
In the usual approach, and after considering D − 1 = 3 − 2 , simplifying the summations with Bessel functions, expanding about = 0, and performing summations using analytic continuation of the Epstein-zeta functions, eq. (A.19) leads to the further refined expression of the free energy density, f (c=1) (T, L 1 ), of a system geometrically confined in between two infinite parallel plates separated by a distance L 1 where the Li n are the usual polylogarithm functions. We can also follow the alternative approach, as presented in appendix (A.2), to compute the same free energy density. Employing similar techniques and making use of eq. (A.20) in order to analytically continue the Epstein-zeta functions without formally manipulating a divergence, we end up with a similar expression which we further resum using the contour integral representation of the polylogarithm function. We find where in the last summation, we kept the −1 explicit, even though there exists a simple closed form for coth(x) − 1, since the less simple form improves the convergence properties of the expression. As previously mentioned, we see that eq. (2.11) converges exponentially fast for low values of the dimensionless variable T L 1 , while eq. (2.12) converges exponentially fast for high values of this dimensionless variable. The latter even has enhanced convergence properties (due to the additional resummation which we performed), for nearly all T L 1 down to T L 1 ∼ 10 −6 .
A similar resummation could have been performed on eq. (2.11), but we will refrain in doing so since the alternative result covers nearly all T L 1 values.

Case II: Infinite rectangular tube
Using the usual approach together with the same type of procedure as in the parallel plates case, we may then set c = 2 in eq. (A. 19), perform a few more refinements, and obtain the free energy density, f (c=2) (T, L 1 , L 2 ), of a system geometrically confined within an infinite rectangular tube of section where K n are the usual modified Bessel function of the second kind. We may also use our alternative approach, in order to compute the very same free energy density, doing so we obtain Case III: Finite volume box Using, again, the same type of procedure as for the two previously cases, we may set c = 3 in eq. (A. 19), perform a few more refinements, and then obtain within the usual approach, the free energy density, f (c=3) (T, L 1 , L 2 , L 3 ), of a system geometrically confined in a finite volume box (2.15) Introducing a new notation, in order to shorten our next representation, namely the set σ of permutations of L i , 3 as well as Ω = {L 1 , L 2 , L 3 }, we may also use our alternative approach to compute the same free energy density, 3 Containing six elements. For instance, if the third element of σ is σ 3 = (L 3 , L 1 , L 2 ), then σ 3 (2) = L 1 .

Thermodynamic expressions, first and third laws of thermodynamics, statistical fluctuations, and nonextensivity
In this section we 1) examine how the First Law of Thermodynamics is altered by the compactified directions, in particular showing how the pressure is no longer a scalar but depends on direction, 2) provide some formal manipulations to arrive at the formulae for the usual thermodynamic quantities such as the EoS that we ultimately evaluate numerically in section (4), 3) comment on the effects of geometric confinement on the Third Law of Thermodynamics, 4) quantify the thermal fluctuations of our geometrically confined system, and 5) show explicitly that a thermal system between parallel plates is not extensive.

Modification of the first law
Recall that the fundamental object that is computed when using the canonical ensemble is the partition function, Z(T, {L i }). The partition function is indeed fundamental in the sense that it allows one to compute any thermodynamic potential, and therefore any thermodynamic quantity of interest. Recall also that the natural potential to use in the canonical ensemble is the Helmholtz free energy F (T, {L i }), as previously derived in section (2). Besides being related to the partition function, the free energy is also defined as the Legendre transform [102] of the total energy E(S, {L i }) with respect to the total entropy S provides a generalization of the definition for the usual Stefan-Boltzmann thermodynamic temperature [103,104]. Recall also that the Legendre transform in eq. (3.1) is well defined only if E is a convex function of S, at constant lengths L i , i.e. E (S, . . .) > 0. We will see below that E is convex for our geometrically confined systems. Furthermore, the total energy as a function of the total entropy S and the lengths L i , i.e. the entropy as a function of only extensive parameters 4 , is a thermodynamically complete function: S(E, {L i }) contains the full thermodynamic information.
Other relations such as, e.g., the pressure as a function of T and L i , p(T, L i ), are equations of state, which may only contain partial information about the thermodynamics of the system. See, e.g., [103,104] for more details on these general concepts. From both the variable dependence of the total energy and eq. (3.2), we can readily write Let us then comment on the subsequent modification of the first law of thermodynamics, recalling that for infinitesimal changes the most general form [103,104] of the First Law is where dQ and dW are respectively the infinitesimal transfers of heat and work. Considering quasistatic processes, without any loss of generality, and since the only work which can be done is the one coming from a displacement of the boundaries in any of the three directions, we obtain a generalization of the First Law of Thermodynamics to include the possibility of asymmetric pressures: where the P i are the pressures along each (i) of the three directions and are defined as Eq. (3.6) shows explicitly that for asymmetric systems, the pressure may be anisotropic. The above manipulations also allow us to make the temperature and pressure functions of the entropy, S, and the system side lengths, Note that in the limit in which all lengths L i become asymptotically large, the second term of eq. (3.5) reduces to the usual −P dV . Thus, in this limit, we recover that the pressure P is thermodynamically complete, an obvious consequence of the relation P V = −F .

Thermodynamic expressions
Let us now give some details concerning the practical computation of various thermodynamic quantities. We start from the free energy as derived in section (2).
Starting from eq. (3.1), we already have that and thus In addition, given the proof in section (B.2), we know that eq. (3.6) is equivalent to Furthermore, using eqs. (3.8) and (3.9), we obtain the trace of the energy-momentum tensor where ε ≡ E/V is the energy density. Not surprisingly, given that our noninteracting model is indeed conformal, this quantity clearly vanishes. Notice further that it is indeed a generalization of the Stefan-Boltzmann limit result, T µ µ = ε−3P , and encodes the aforementioned system anisotropy. Coming back now to eq. (3.4), and using the explicit form for the total amount of infinitesimal work that can be done from eq. (3.5), we obtain a heat function Q which admits an exact differential, when all the lengths L i are kept fixed 5 . More precisely, we have which leads to the following definition of the heat capacity at constant sizes (hence volume): We keep the usual V index on our heat capacity to emphasize the connection with the usual heat capacity at fixed volume. It is worth noting that given eqs. (3.7) and (3.8), the heat capacity can be re-written as (3.13) 5 Hence when the volume is fixed too, even if the latter constraint does not imply the former.
We will see in section (4) that our finite-size heat capacity is always positive. The positivity of the finite size heat capacity means the connection between the pressure, temperature, and entropy is straightforward: the relation T = T (S, {L i }) is invertible and hence equivalent to S = S(T, {L i }). Thus we may equivalently use Further, this statement together with eq. (3.13) imply that E (S, . . .) = T /C v > 0, thus making the Legendre transform connecting the free energy with the total internal energy well defined.

On the third law
The Third Law of Thermodynamics states that the entropy reduces to a constant value (usually 0) as the temperature of a system goes to 0 [103,104]. Some implementations of the Lifschitz theory of the Casimir effect in QED claim that in their system the Third Law of Thermodynamics is violated [106]. We first recall that in our system the spatial compactification(s) required for establishing a geometric confinement lead to a zero temperature Casimir-type geometrical contribution to the free energy F (T, {L i }). However, given the definition of the entropy (3.7), we expect that this temperature independent contribution vanishes from the expression of the entropy. The numerical value of the entropy can then be assessed by appropriately rescaling S by the dimensionless combination T 3 V . Numerical evaluation of this quantity (see section (4)) clearly shows thats → 0 as the temperature vanishes. We refer to [107] for another study of the Third Law in the context of QED Casimir systems.

Fluctuations of the energy
Recall that in the micro-canonical ensemble the total energy is fixed. However, in the canonical ensemble the system is in contact with a heat bath of infinite heat capacity and all energies of the system are accessible through a probability distribution. For the canonical ensemble, then, there is a mean total energy together with a standard deviation (the square root of the statistical variance . It is precisely the variance of the energy of the system that we take to provide us with insight into the energy fluctuations of our canonical ensemble system. One may compute the mean total energy E ≡ E in the usual way with partition functions , (3.14) where the angular brackets denote the usual ensemble average The standard deviation, σ E = ∆E 2 , is then .
Therefore, after a little manipulation, the standard deviation may be simply written as where we recall that the heat capacity C v is defined in eq. (3.12).
Recall that both the heat capacity and the total energy of a system scale with the volume V of the system. Therefore as the system volume increases, the relative size of the energy fluctuations decreases, σ E /E ∼ 1/V 1/2 → 0. Systems that are not compactified in all directions, such as the case of two infinite parallel plates or a rectangular tube of infinite length, have infinite volume and therefore experience no fluctuations in their total energy. Crucially, then, for these systems with no fluctuations in total system energy, we conclude that the canonical ensemble must exactly reproduce the results of the microcanonical ensemble. We will exploit this equivalence of ensembles later to demonstrate a phase transition at a critical length (instead of temperature) for isolated systems of finite extent in some (but not all directions) in section (5).
In figure (1), we plot the mean total energy and its standard deviation as a function of temperature T (in units of 1/L) for the case of a massless, noninteracting scalar field geometrically confined in a finite-sized box. The upper left panel shows the results for a symmetric box, with side ratios 1:1:1; the upper right panel shows the same for an asymmetric box of side ratios 1:1:3, a finite volume symmetric tube; the lower panel shows the same for an asymmetric box of side ratios 1:3:3, a set of two finite area parallel plates. One can see that as T grows large, the system approaches the Stefan-Boltzmann limit, which we denote by "SB". Even out to relatively large T × L ∼ 20, energy fluctuations are on the order of 10%. In the limit that T ×L → 0 the energy density becomes negative, which is the usual case for Casimir-like systems 6 .

Nonextensivity of finite size systems
We show in figure (2) the deviation from extensivity for the parallel plates case. The plot shows the difference between the entropy for a massless, noninteracting scalar field between two parallel plates a distance 2L apart and twice the entropy for plates a distance L apart scaled by T 3 V . For an extensive system, this difference is 0. The nonextensivity goes to zero as T → 0 trivially: the entropy for the parallel plates case falls faster than T 3 V for small T , which we show in detail below in the left panel of figure (7). A system with all lengths infinite is extensive. Hence we may understand why the system approaches extensivity as T increases in the figure as follows. As T increases, the thermal de Broglie wavelength decreases, and the effective size of the system becomes larger.

Thermodynamic properties of a geometrically confined scalar field
In this section, we present a number of thermodynamic quantities, ranging from the free energy to the specific heat capacity at constant lengths, relevant to various geometrically confined systems (infinite parallel plates, infinite tube, and finite volume box). The results are all computed from the canonical ensemble, which is to say for systems in contact with an infinite heat bath held at constant temperature T , for a massless, noninteracting scalar field. Each of the plots has been rescaled by the Stefan-Boltzmann result. Recall that in the Stefan-Boltzmann limit of zero coupling and infinite size in all directions, various thermodynamic expressions we wish to evaluate in the finite size case are In figure (3), we show the free and total internal energies in the case of two infinite parallel plates as a function of the temperature T in units of 1/L, where L is the distance between plates. Recall that for a plasma of temperature ∼ 400 MeV and width of ∼ 2 fm, relevant for a high multiplicity pp or pA collision at RHIC or LHC, T × L ∼ 4. For a mid-central AA collision resulting in a plasma of temperature T ∼ 400 MeV and width ∼ 10 fm, T × L ∼ 20. One can see in the figure that both results tend towards the thermodynamic limit as T × L → ∞. However, both the energy and the free energy are 5% different from their Stefan-Boltzmann limits even at the relatively large value of T × L ∼ 20. The total energy of a system geometrically confined in between two infinite parallel plates separated by a distance L, and in contact with a thermal bath at temperature T , is thus noticeably affected by its finite size.   In figure (4), we display the free and total energies in the three cases of infinite parallel plates (blue lines), infinite tube (yellow lines), and finite volume box (red lines). All results tend towards the Stefan-Boltzmann ones as T ×L → ∞. We see again the large finite size corrections to the Stefan-Boltzmann limits, even for T × L ∼ 20, with the size of the corrections increasing with increasing number of compactified directions. Thus the total and free energies of a system geometrically confined is noticeably affected by its finite sized directions. We note that the peculiar behavior of the tube case, that its total energy reaches a T = 0 limit which is positive unlike the plates and box cases, is not a surprise. It is indeed due to the dimensionality of the space-time [108].

SB Limit
We now turn towards the different pressures. In figure (5) we plot the perpendicular (left) and parallel (right) pressures as a function of the temperature T in units of length 1/L for a massless, noninteracting scalar gas held between two infinite parallel plates that are a perpendicular distance L apart. The plots have a number of interesting properties. First, the pressure in the perpendicular direction (p 1 , left panel), i.e. between the plates, is different from the pressure in the parallel direction (p 2 or p 3 , right panel), i.e. along the plates. Thus no single scalar pressure can be used to describe the finite sized parallel plates case: there is an intrinsic anisotropy due to the The difference in entropy for a massless, noninteracting scalar field between parallel plates a distance 2L apart and twice the entropy for a system with plates a distance L apart, at a temperature T measured in units of 1/L scaled by the quantity T 3 V .  compactification. Furthermore, the compactification of one direction appears to have a greater effect on the pressure in the noncompactified parallel direction than on the pressure in the compactified perpendicular direction: the parallel pressure is nonmonotonic, unlike in the perpendicular case, and approaches the Stefan-Boltzmann limit much slower than the perpendicular pressure. Note that the pressure for our geometrically confined system dimensionally reduces in the correct way. For our infinite parallel plates case in 4D spacetime where vacuum depends only on L. Note that p is related to the parallel pressures p (1) 2/3 between parallel plates discussed earlier (3.9) by p = p (1) 2/3 /L. One can see that in the limit T × L 1 eq. (4.1) reproduces the pressure of a massive, noninteracting scalar theory in 3 spacetime dimen- The pressures, perpendicular (left, p1) and parallel (right, p2 or p3), for a massless, noninteracting scalar field restricted between two infinite parallel plates as a function of temperature T in units of the perpendicular distance L between the two plates. Both quantities are relative to the usual Stefan-Boltzmann pressure.
sions in the Stefan-Boltzmann limit for m ≡ π/L, where we have dropped the temperature independent, infinite zero point pressure from the last expression. One sees then that the inverse of the compactified length L in 4D acts as an effective mass in the 3D theory. We show in figure (6) the pressure in one of the compactified directions for the two infinite parallel plates, the infinite (symmetric) tube and the finite volume (symmetric) box cases, respectively represented by the blue, yellow, and red lines. The pressures are rescaled by their Stefan-Boltzmann limits and are plotted as a function of the temperature T in units of 1/L, where L is the length of the compactified direction.
Notice that the effect of further compactifications on the pressure is drastic. In the finite box case, for T × L ∼ 4, which is relevant for a pp collision resulting in a ∼ 400 MeV QGP, the pressure sees a ∼ 40% correction. Even for T × L ∼ 20, there are ∼ 10% corrections to the pressure of the finite volume box. Note that while it might appear that the pressures diverge at low T , this apparent divergence is an artifact of plotting the ratio of our finite sized results with the Stefan-Boltzmann limit; recall that the Stefan-Boltzmann case scales as T 4 . One may use the un-rescaled expression eq. (2.11) to investigate the asymptotic T = 0 behavior, in which case one recovers the usual Casimir pressure in the perpendicular direction in each compactification case. (This comment also applies to the other pressures, and the free and total energies.) We plot in figure (7) the entropy (left) and specific heat (right) for the infinite plates case as a function of the temperature T in units of 1/L, where L measures the length of the finite direction of the system, and scaled by the Stefan-Boltzmann limit. The inset in the left plot shows that the entropy of our finite sized system indeed goes to zero as the temperature vanishes, as dictated by the Third Law of Thermodynamics, and as opposed to some implementations of Lifschitz's theory for the QED Casimir effect [106]. The inset in the right plot shows that the specific heat for our finite sized system always remains positive.
We plot in figure (8) the entropy (left) and specific heat (right) for the infinite plates (blue), infinite tube (yellow), and finite box (red) cases as a function of the temperature T in units of 1/L, where L measures the length of the finite direction(s) of the systems, and scaled by the Stefan-Boltzmann limit. Although we do not provide insets in these plots, the entropy for all our cases again goes to 0 as T → 0, in accordance with the Third Law of Thermodynamics and the specific heat remains strictly positive. Notice again that the size of the deviations from the Stefan-Boltzmann limit increase as the number of compactified dimensions increases. In particular the deviation from the Stefan-Boltzmann limit is significant (∼ 10 − 15%) for the finite box case even out to T × L ∼ 20.
We would like to understand a bit better the importance of the type of boundary condition imposed on our systems on the size of the finite size corrections to the Stefan-Boltzmann limit. Intuitively, one might expect that periodic boundary conditions cause the least difference since the system has a finite size but is boundaryless. We show in figure (9) the free energy of a massless, noninteracting scalar field between two infinite parallel plates rescaled to its Stefan-Boltzmann  limit as a function of the temperature T in units of 1/L, where L is the distance between the plates. The plot shows the results for Dirichlet (solid blue), Neumann (dotted brown), and periodic (dashed purple) boundary conditions. One can clearly see from the figure that the system with periodic boundary conditions reaches the Stefan-Boltzmann limit at much smaller T × L than either the Neumann or Dirichlet cases. One should perhaps then hesitate to conclude that small finite sized corrections computed for a system with periodic BCs [70][71][72][73][74][75][76][77][78] will remain small for heavy ion phenomenology.
Let us now come back to our geometric confinement boundary conditions. Since our noninteracting scalar field theory is conformal, the trace of the energy momentum tensor should vanish identically. We performed a nontrivial check of our numerics and confirmed that our results do respect T µ µ ≡ ε − (p 1 + p 2 + p 3 ) = 0.  Figure 9. The free energy of a massless, noninteracting scalar field in between two infinite parallel plates rescaled to its Stefan-Boltzmann limit for Dirichlet (blue), Neumann (dotted brown), and periodic (dashed purple) boundary conditions.

A novel geometric phase transition
In this section we consider how our results depend on the physical setup of our system. In the previous section we considered our system to be in contact with an infinite thermal heat bath. We now focus in particular on the infinite parallel plates case and consider the possibility of the system existing in isolation: instead of a system at constant temperature T we rather consider a system with constant energy E. All results shown in this section are derived from the partition function calculated in the canonical ensemble. We are justified in this approach as argued in section (3.4): the ratio of the energy fluctuations away from the mean energy divided by the mean energy computed in the canonical ensemble for the Stefan-Boltzmann parallel plates case is zero; hence we will have ensemble equivalence between the canonical and microcanonical ensembles. Thus the canonical ensemble provides exact results for all thermodynamic quantities for a set of parallel plates kept in isolation, i.e. with fixed energy E.
As we explore the different physics associated with a noninteracting scalar field between two infinite parallel plates, it is worth keeping in mind two related thermodynamic concepts. First, a system is thermodynamically stable if, when squeezed and all other parameters are held fixed, the pressure increases. I.e. our parallel plates system separated by a length L is thermodynamically stable so long as for fixed E and plate area V 2 , which is to say that the pressure in the system decreases with increasing system size; see, e.g., [69]. Second, a system undergoes a second order phase transition should a susceptibility diverge; see, e.g., [109]. Recall that the susceptibility is the derivative of an extensive parameter with respect to its conjugate intensive parameter. The susceptibility is then nothing more than the multiplicative inverse of the second derivative of the entropy with respect to some parameter. For our case of two parallel plates separated by a distance L, the susceptibility is the length-scaled negative of the compressibility: Comparing the condition for stability, eq. (5.1), to the definition of susceptibility above, one can see that a second order phase transition occurs when the system goes from being thermodynamically stable to unstable. Further, the order parameter associated with the compressibility is the size of the system L.
We can further understand physically what happens at a phase transition by examining the relationship between the entropy and an intensive variable such as the pressure. There are three common, and equivalent, definitions of pressure: The first expression is the quantity that must be equal for two systems in thermodynamic equilibrium separated by a moveable wall. The second expression is the generalized force conjugate to the volume. And the third is the generalized thermodynamic intensive variable conjugate to the extensive volume variable. In the parallel plates case of current interest, the derivatives with respect to volume V become derivatives with respect to separation length L. One can then see that a phase transition occurs precisely when the second derivative of the entropy changes sign; i.e. when the entropy goes from a convex function of L to a concave function of L. While all three expressions of the pressure in eq. (5.3) are equivalent, the different expressions are naturally a function of different variables: p(E, L), p(S, L), and p(β, L), respectively, where we have already switched over to using the plate separation distance L instead of the volume V for our system. For simplicity in this section, we denote the previously defined perpendicular pressure p (1) 1 by p. One may freely switch between different definitions of pressure by using relations between the various independent variables. For example, equipped with an expression that relates the energy to the entropy and volume, one can equivalently use the first definition of pressure in the same way as the second definition with p E(S, L), L = p(S, L).
In figure (10) we compare the pressure of a massless, noninteracting scalar field between parallel plates of fixed area V 2 as a function of the plate separation length L for fixed temperature T (left) and for fixed energy E (right); i.e. the left plot shows the pressure as a function of L for a system in contact with a heat bath whereas the right plot shows the pressure as a function of L for an isolated system. We computed p(L) from the partition function eq. (2.11) by numerically inverting E(T, V 2 , L) to find p T (E, V 2 , L), V 2 , L . It is perhaps not so easy to see in the figure, but for the system in contact with a heat bath, the pressure always decreases for decreasing L; i.e. ∂p/∂L > 0 and the system is always unstable: the system always wants to collapse. The isolated system, however, resists collapse as the system size is decreased-i.e. ∂p/∂L < 0 and the system is thermodynamically stable-so long as the system starts off large enough, which is to say the length L is greater than some critical length L c . As soon as L ≤ L c , the system is unstable and will collapse, shrinking until the plates are no longer separated at all. We show explicitly the susceptibilities as a function of plate separation L for these two systems in figure (11). As claimed, the isothermal compressibility is purely negative while the isoenergetic compressibility clearly exhibits a divergence. We thus conclude that our massless, noninteracting scalar field theory constrained between two parallel plates exhibits a phase transition at a critical length L c . This novel phase transition is interesting for two main reasons. First, we believe this is the first example of the explicit derivation of a phase transition in an ideal gas induced by changing the size of the system, as opposed to the usual means of inducing a phase transition by changing the temperature of the system. Second, in the derivation shown above, the phase transition is due to changing an extensive variable, as opposed to the usual description of a phase transition due to changing an intensive variable. To be clear: we envision the (somewhat artificial) construct of a pair of (approximately infinitely large) parallel plates of fixed separation length L. The space between the plates is filled with a noninteracting, massless scalar field, and one measures the pressure on the plates. Since the separation length is fixed, the system cannot actually collapse. However, one is tempted to interpret the right plot of figure (10) as follows. Consider a system of two (approximately infinitely large) parallel plates filled with a noninteracting, massless scalar field. One plate is fixed, but the other plate is freely allowed to move and is exposed to a constant external pressure p; at equilibrium, the parallel plates settle down to an average separation distance L set by p. As one slowly dials up the external pressure p, L decreases but the system continues to find a new, smaller L at which the plates are in equilibrium with the external pressure p 7 . At a certain large enough critical pressure p c , the scalar field inside the plates can no longer resist the external pressure, and the system collapses. One would very much like to confirm this extrapolated interpretation from the canonical ensemble; to do so would require going to the higher ensemble in which the system is in contact with a hypothetical thermal pressure bath. It is perhaps not so easy to see from figure (10) or figure (11), but we will show below that the phase transition occurs at a length of order the thermal de Broglie wavelength, L c ∼ 1/T . Since the phase transition seen from the canonical ensemble calculation occurs at a length of order the thermal de Broglie wavelength, it is possible that fluctuations in the separation distance L from the mean distance L that one would necessarily observe in a system exposed to a thermal pressure bath spoil the observation of a phase transition. A quantitative derivation of the properties of a massless, noninteracting scalar field between two parallel plates of variable separation distance and exposed to a thermal pressure bath via a higher ensemble is left for future work.
Since we observe this phase transition for a massless, noninteracting scalar field theory in a small enough geometric confinement, a natural question to ask is: is this transition one of Bose-Einstein condensation? One can see an indication of an answer in the negative from figure (12). On the left, we plot the entropy S as a function of the area of the parallel plates, V 2 , and the separation of the plates, L, for our isolated scalar field theory system with constant total energy E. For small enough L one can explicitly see how the entropy transitions from a convex to a concave function. We further quantify this phase transition by examining the most positive eigenvalue of the Hessian of S(V 2 , L), i.e. with E fixed. For a concave function, all eigenvalues of the Hessian are negative. The plot on the right of figure (12) shows in gray the region of (V 2 , L) space for which the most positive eigenvalue of the Hessian is greater than zero, which is to say the region for which the entropy is no longer a convex function of V 2 and L 8 . The critical length that forms the phase transition boundary of this region is a function of the parameters V 2 , E, and T × L. In the right hand plot of figure (12), the edge of the region for which the entropy is no longer concave is given by a function L c (E, V 2 , T × L) where T × L ≈ 1; i.e. the phase transition occurs when L ≈ 1/T . Therefore this geometrical phase transition can occur for a Bose system at arbitrarily large temperature. We can conclude then that the phase transition is not one of Bose-Einstein condensation.
One may then naturally next ask if the phase transition can be found in a massless, noninteracting Fermi field. We will show below that, yes, the phase transition persists for a massless, noninteracting Fermi field. In order to repeat the above analysis for a Fermi field, we must slightly alter our quantization condition [50,51]. In order to prevent any Dirac current from leaking through  the plates confining our system, we require that In order to satisfy eq. (5.4), one must take the momentum modes of our Dirac field to satisfy Following the methods of section (A.1) the Fermion partition function for a massless, noninteracting Dirac field constrained between two parallel plates of area V 2 and kept a distance L apart evaluates to whereβ ≡ π/(T L). Note that the overall factor of 2 in eq. (5.6) is due to the spin-1/2 nature of the Dirac field. The results of repeating the scalar field theory analysis are shown in figures (13), (14), and (15). Qualitatively the picture is the same for the Dirac theory as for the scalar theory; quantitatively, one finds that the critical length separating the concave from convex region of the entropy S as a function of plate area V 2 and separation length L is related to a slightly higher temperature of the gas, L c ≈ 0.8/T . One can even go so far as to show that the phase transition exists for massive, nonrelativistic Bosons and Fermions, too, although we leave the details for a future publication.

Summary and prospects
In this work we implemented the first step of a strategy designed to investigate the thermodynamics of small QGP systems, focusing on the finite size corrections to the usual Stefan-Boltzmann thermodynamic properties computed in thermal field theory and also demonstrating the emergence of a new geometric phase transition. In particular, we set up a framework for probing relevant finite size effects by means of an actual spatial boundary for a geometric confinement, imposing Dirichlet boundary conditions. We argued in section (1.2) that Dirichlet BCs are the most appropriate for capturing the relevant finite size effects for heavy ion collisions and quark-gluon plasma phenomenology because these BCs prevent the weakly coupled quantum fields, supposed to account for the relevant degrees of freedom inside a QGP, from propagating outside of the geometric region defining this QGP system. For the sake of simplicity, we considered D − 1 spatial dimensions with c ≤ D − 1 dimensions of finite length L i arranged in standard rectilinear coordinates. We then filled the space inside our geometrically confined region with a massless, noninteracting scalar field. Such a model encodes the main aspects of a massless, noninteracting gas of gluons. In order to best make contact with lattice QCD results and also for simplicity we computed the thermodynamic quantities in the canonical ensemble, in the process extending standard thermal field theory techniques.
We presented the main results of the detailed analytic calculations in the Appendices in section (2). The results were computed by two independent methods, resulting in two different infinite sums for each geometrically confined system. The two different sums have very different numerical convergence properties: one converges exponentially fast for T × L i < 1 while the other converges exponentially fast for T × L i > 1.
In section (3) we discussed some important qualitative consequences of our derivation of the statistical mechanics of our system. In particular, we saw that the First Law of Thermodynamics is generalized by the presence of different pressures: instead of a single scalar pressure as in the case of all spatial dimensions having infinite extent, the pressure in the i th direction may be different from the pressure in the j th direction depending on the various lengths L k . On the other hand, we found that the limited spatial extent of our system did not affect the Third Law of Thermodynamics. Of critical importance, we computed the size of the fluctuations in energy away from the mean energy dictated by the contact of our system with a thermal heat bath. Even though the systems of parallel plates or a tube have some but not all direction(s) of finite size, the volume of these systems is infinite. Therefore the relative size of the energy fluctuations compared to the mean energy is 0. As a result, the canonical ensemble can be used to compute the thermodynamic properties for isolated, constant energy parallel plates or infinite tube systems. At the same time, the energy fluctuations in a finite volume box are not zero, and there must be corrections in the application of the canonical ensemble to an isolated, finite volume box system.
In section (4) we presented a number of quantitative plots of the free energy, entropy, energy, pressure, and heat capacity for a massless, noninteracting scalar field contained inside infinite parallel plates, an infinite rectangular tube, and a finite-sized box. Of particular note was the surprisingly large finite size corrections for the pressure. For a system of the approximate size and temperature of a high multiplicity proton-proton collision at LHC, which has been argued to exhibit hydrodynamic behavior based on thermodynamic quantities computed in the Stefan-Boltzmann limit [44], we found ∼ 20% corrections for an infinite tube and ∼ 40% corrections for a symmetric, finite box. Even for systems of the size of mid-central nucleus-nucleus collisions at LHC the corrections are of order ∼ 10%. Since the size of the azimuthal anisotropy measured in such collisions at LHC are of order ∼ 10% [14], these corrections to the equation of state due to the finite size of the system may have important implications for heavy ion phenomenology.
In section (5) we discovered that an isolated system of noninteracting particles confined within parallel plates at temperature T undergoes a second order phase transition at a critical length L c ∼ 1/T : for L > L c the system is stable and resists compression; for L < L c the system collapses. The phase transition is not one of Bose-Einstein condensation as we showed that a noninteracting Dirac field also experiences such a phase transition. Nor is the transition a relativistic effect as it is also experienced by massive, nonrelativistic Bosons and Fermions. It is tempting to interpret our results as follows: a system of fixed energy can only resist so much external pressure before collapsing. One can see that figures (10) and (13) support this interpretation: the system energy density is on the same order as the pressure, with p E/V , when the system undergoes the phase transition at the critical length. Further clarity on this interpretation requires the use of a higher order ensemble in which the system is put in contact with a thermal pressure bath. Such an investigation would also provide insight into the importance of the fluctuations in the separation distance between the parallel plates about the average, equilibrium length. Should this higher ensemble show that a pressure driven phase transition exists, one could use the higher ensemble to calculate the critical exponent for this second order transition. Other work that has claimed the observation of a first order phase transition in small systems [69] performed the calculation in a finite volume box in the canonical ensemble; it is not clear that the energy fluctuations inherent in using the canonical ensemble invalidates the conclusions reached in that work. In the only other work that we are aware of that examines finite size driven phase transitions [110], the phase transition is due to self interactions of the system; since our system is noninteracting, the phase transition we see is due purely to geometric confinement effects. Future work includes determining whether or not a similar phase transition exists for an isolated system with periodic boundary conditions in one direction or for a tube system. It is important for us to point out that geometrically confining a thermal quantum field into a certain region of space with Dirichlet boundary conditions (even with an infinite volume) appears to provide a solution to the infrared Linde problem [76,111,112]. The origin of the Linde problem is the existence of the zero mode in systems with, e.g., periodic boundary conditions: a momentum mode with zero energy. Imposing Dirichlet boundary conditions naturally provides an infrared cut-off for the momentum modes p > p min ∼ 1/L dictated by the system size L. We explicitly showed in this work that the Matsubara zero mode disappears from the leading order thermal field theory calculations. It would be interesting to quantitatively check that in fact Dirichlet boundary conditions cure the zero mode Linde problem at higher orders in perturbation theory, too, and see that the perturbative series becomes analytic in the coupling.

A Details of calculations
In this section, we detail two different methods of computing the free energy density for a single neutral noninteracting massless scalar field in c geometrically confined dimensions by deriving the logarithm of the partition function. It is important to note that the usual method for computing such a logarithm cannot be applied in a way that is necessarily straightforward since the compactification of any number of spatial dimensions leads to summations that are formally divergent.
For both methods we employ dimensional regularization [113] within the modified minimal subtraction (MS) scheme [114] in order to regulate the divergences, but we also complement this regularization procedure through the use of zeta type of regularization whenever necessary. In the present massless case of interest, the usual (power like) divergences, if any, are set to zero by Dim. Reg. so that there are no divergences that would need to be cured by means of renormalizing the vacuum [115].
In subsection (A.1) we follow the usual method as described by [97,99]. In subsection (A.2) we present an alternative calculation. Naturally, the two methods yield analytically equivalent results which then of course agree numerically-even though their different convergence properties, when the summations are truncated at some finite orders, make them conveniently suitable in different regimes.
It is important to note that we will make significant use of Epstein zeta functions and their various analytic representations. The derivation of these representations often exploits the Poisson summation formula, We will not go into the detail of these derivations, but rather point the interested reader to the extensive literature [100,[116][117][118][119][120][121][122][123] that is the result of three decades of intensive work on these special functions.

A.1 Usual computation
In this section, we compute the canonical representations for the free energy density in different scenarios, following and adapting the general methods of [97,99]. For pedagogical reasons, we present the case of c = 1 in detail, but the result will then be easy to extend to arbitrary c ≤ D − 1.
Our starting point is then either eq. (2.8) or eq. (2.9), leading to We then consider the following two identities and after substituting eq. (A.3) into eq. (A.2), defining for brevity ω 2 ≡ ω 2 k + m 2 , dropping an infinite T -and L i -independent term, and employing eq. (A.4), we obtain where we again dropped infinite T -and L i -independent terms. Assessing the case c = 1, it is time to release the momentary (periodic) compactifications of the D − 2 other spatial dimensions, sending each of the R compactification lengths to its appropriate value (presently L 2 or L 3 ) each being asymptotically large, and get where ω 2 → k 2 + m 2 , D ≡ 4 − 2 being now a parameter ( will be set to zero at the end), and we introduced an arbitrary regularization scale (Λ 2 ≡Λ 2 e γ E /4π) in order to keep the proper dimensionality when D = 4, since we employ Dim. Reg. within the MS scheme. We now focus on the computation of the second term of eq. (A.7). To do so, recall that the volume of an N -dimensional unit sphere is given by Ω N = 2π N/2 /Γ (N/2), and consider the integral Making the substitutions p → T x and a = T b in the above integral, Taylor expanding the logarithm, and rewriting the integration measure, we obtain Then, the substitution z = √ x 2 + b 2 leads to an integral representation of the modified Bessel function of the second kind, namely K, and we get Therefore, in order to compute the second term of eq. (A.7), we need the following result We further notice that the above result is indeed finite, as the second term of eq. (A.7) contains neither Ultra-Violet (UV) nor Infra-Red (IR) divergences.
We focus now on the computation of the first term of eq. (A.7), which obviously contains a UV divergence. There are a number of different ways of regularizing this divergence, and we follow more specifically [99], employing again Dim. Reg. within the MS scheme. We then have We notice an important point related to the above calculation: the integral does not feature any IR divergences in the massless limit as long as L 1 is finite, because 1 never vanishes. Combining eq. (A.7) with eqs. (A.12) and (A.13), we may write the free energy density, f (c=1) (T, L 1 ), of a system with one geometrically confined dimension as where, of course, we took the massless limit. It is then straightforward to extend the above procedure to an arbitrary number c of compactified dimensions, and obtain the result for which we notice that Dim. Reg. took care of the non logarithmic UV divergence, and that the above result will not need to be renormalized 9 .
Computing the summand in the first term of eq. (A.15) will then require some Epstein-Zeta regularization. Note the present Epstein-Zeta function of interest is defined [119] as We reproduce now the result for the analytic continuation of the above function (see [119] for more details), in which the following recursion relation is useful Note that for certain values of D − c the gamma function in the numerator of eq. (A.15) will have a pole. For precisely those values of D − c the zeta function will multiply these infinities by trivial zeros thus rendering the divergences harmless.
(A.18) which reduces to E 0 1 (s; a) = a −s ζ(2s) in the limit m → 0, providing Re (s) < 1/2. We may now use the expressions (A.17) and (A.18) in order to rewrite eq. (A.15) as which is our canonical result for the free energy density with c geometrically compactified dimensions. Notice that the first line in eq. (A. 19) may now be explicitly computed using eqs. (A.17) and (A.18), for any number of compactified dimensions, and we refer the reader to subsection (2.2) for the corresponding final, and furthermore refined, results.

A.2 Alternative computation
In this section, we compute the alternative representations for the free energy density in different scenarios. We generally follow the methods of [97,99] but take great care with the regularization procedure; given the spatial compactification(s), we avoid the formal manipulation of divergences. Moreover, we notice that certain massive results will be given as byproducts, even though these are of no interest for the present work. In fact, straightforwardly taking the asymptotically large R limit in eq. (2.9) could appear problematic for some values of c. Indeed, the subsequent integral and the remaining infinite sums may not have a common strip of convergence. To show that this would not be a problem, we could make use of the monotone convergence (or Beppo-Levi) theorem. Given the asymptotically large R limit and since the infinite summation kernel of eq. (2.9) would then still be positive definite and monotonically increasing 10 , one could take the infinite series as a set of partial summations under a limit and exchange the order between the limit and the integral. This could allow for an analytic continuation without spoiling the convergence of the expression, thereby asserting the existence of the dimensionally regularized free energy density for all c. However, we could, instead, proceed with a convenient trick, making use of the following identity where A is positive definite and α is real. 1/A α can serve the purpose of defining a master expression that can be analytically continued as a function of the dimension D, before applying eq. (A.20) to find the log. Doing so and momentarily assuming a convenient range for α, allows us to work with convergent expressions. Using this identity, and sending all R to their asymptotically large values, we define the following master sum-integral 11 21) related to the free energy density in D dimensions, after analytic continuation in terms of the complex D variable, via the following identity We will then first focus on the evaluation of eq. (A.21). Prior to setting the number of compactified spatial dimensions c to some value in order to do so, and thus to compute the free energy density, let us perform the continuous momentum integration. By doing so, we will be able to analytically continue the expression as a function of the dimensional parameter D, and evaluate the infinite summations outside of their original radii of convergence when needed. For now, we shall assume that there exists some strip (to be specified below) for the complex variable D, in which eq. (A.21) is convergent. We then proceed to rescale the continuous momentum k, in order to factorize all the frequencies and the mass. We then obtain the following expression In order to assure the existence of the strip we take α > 0. Then the integral above is clearly convergent in the D-strip such that Re (D) − 1 − c ∈ (0, 2α). The set of infinite summations, being some Epstein zeta type of function, naturally converges in a D-strip such that 2α+1+c−Re (D) > c.
Consequently, all we need to assume for now is that α and D satisfy the following conditions: 2α > Re (D) − 1 ≥ Re (D) − 1 − c > 0, the middle one being automatically satisfied. We shall then work under this assumption, until we are able to perform an analytic continuation in the dimensional parameter D. Then, the constraint on α will be released, and we will be able to compute the free energy density according to eq. (A. 22), in all the possible spatial compactification cases. Next, we use standard techniques [97,99] to perform the continuous momentum integration in eq. (A.23), and we get being left with the evaluation of the following (1 + c)-dimensional Epstein zeta function for which we need to set c to some particular value prior to any evaluation.

Two infinite parallel plates
In this subsection, we shall analytically continue eq. (A.25) in order to analytically continue eq. (A.24). Taking c = 1 We then split the Matsubara sum into a zero and a nonzero mode contribution, following and make use of the two-dimensional Epstein zeta function representation [118] +∞ n, =1 for the second term, the first one being the usual Riemann zeta function. Analytic continuation of the above formula readily follows from the straightforward continuation of the Euler gamma function, the Riemann zeta and modified Bessel functions, or simply using the functional equation for the corresponding Epstein zeta function [118].
As it turns out, given a certain choice for a and b, the limit of vanishing T or asymptotically high L 1 appears explicitly. Note also the change in the argument of the Bessel function, upon an exchange between a and b. This gives the possibility to slightly enhance the convergence of the remaining two-fold sum, as depending on the choice for a and b, a factor of b/a or a/b (giving either a 2T L 1 or its inverse) appears in the argument of the Bessel function. And given that the larger the argument the better the convergence, we choose the former possibility in order to easily reach the thermodynamic, Stefan-Boltzmann limit. Therefore, the asymptotically high L 1 limit (that is the usual Stefan-Boltzmann finite temperature result) will appear explicitly in the analytically continued result. Using then eq. (A.28) with a ≡ (π/L 1 ) 2 and b ≡ (2πT ) 2 , we further analytically continue eq. (A.26) and can finally extend it together with eq. (A.24) for c = 1, as functions of D from the original D-strip to the whole complex plane. Gathering all together, we get the following result for the corresponding eq. (A.24) with c = 1 and m = 0, valid on the whole D-complex plane. This means that we can actually probe the above result for any value of α. We then apply the relation (A. 22), and since there is no pole around D = 4 (i.e., no logarithmic UV divergence), we readily set the dimension to four and get the renormalized 12 result for the corresponding free energy density. We notice indeed that the first term is the usual Stefan-Boltzmann finite temperature result, relevant to a massless scalar field in a noncompactified space.
In addition, we see that the above result reduces to the Stefan-Boltzmann one when sending L 1 to infinity. Note however, that the zero temperature result relevant to the usual Casimir effect, and equal to −π 2 /1440L 4 1 [118], appears only explicitly when using the same analytical representation for the two-dimensional Epstein zeta function, but done with a = (2πT ) 2 and b = (π/L 1 ) 2 instead of the present choice. Our result contains a two-fold infinite summation which is of course convergent. And in fact, one can also get a closed form when performing one (or the other) summation, with once again two different looking yet equivalent representations, depending on which summation is explicitly done. We choose to explicitly perform the summation over the 1 variable, for convenience when further enhancing the convergence as we are going to explain below. 13 Doing so, we end up with the following one-fold representation Li 3 e −4πT L1n , (A. 31) involving, again, two polylogarithm functions. With the above representation, one would need to change the remaining infinite summations by the corresponding integrals, in order to probe the asymptotic T = 0 limit where the n variable becomes then continuous. Note that this representation, when truncated to some finite order, is relatively well behaved in terms of convergence. However, the situation can easily be further improved. Indeed, making use of a contour integral representation for the polylogarithm functions, one can replace the summations above with another set of summations with much better convergence properties. This new one-fold representation happens to be exponentially enhanced, as can be easily checked. Relabeling the dummy variables, our final result for the free energy density of one neutral massless noninteracting scalar field is (A.32) As a matter of fact, this new representation allows for both the T = 0 and L 1 = ∞ limits to be carried out, which is relevant to the corresponding quantum field when coupled to a heat bath at temperature T , and geometrically confined in between two infinite parallel plates separated by a distance L 1 . We notice that the term proportional to T L −3 1 actually cancels the s −3 term in the second sum, but it is preferable to keep the result as such. The above second sum indeed converges more rapidly.

Infinite rectangular tube
We shall now analytically continue eq. (A.25) with c = 2, that is No counter term was needed, yet Dim. Reg. set the non logarithmic type of divergence to zero which makes the result formally renormalized. 13 Notice that both the Matsubara and the spatial summations have been done. Even though we keep the same dummy variables n and 1 , they do not correspond anymore to the original modes which have been summed over.
which will allow us to investigate the free energy density of a neutral noninteracting scalar field coupled to a heat bath at temperature T , and geometrically confined inside an infinite tube of rectangular section with two sides of respective finite lengths L 1 and L 2 . For computational convenience, let us keep the mass nonzero for a little longer.
We then notice that we can use the following identity noticing the difference of definition for the a i coefficients, with respect to eq. (A.17). Thus, with the help of this above representation, together with the identity (A.34), we can obtain an analytic and symmetric representation for the corresponding master sum-integral I (2) α , and further apply the relation (A.22) for obtaining the free energy density. Finally, given that D = 4 − 2 , we Laurent expand the resulting expression around = 0. After a few more steps, which we will avoid here for the sake of brevity, we finally obtain Let us now simply comment on the massless limit, given the unusual UV structure of the above expression whose renormalization techniques, if necessary, can be found at [49,124] by means of background field methods. However, we notice that under the limit m → 0, the above result becomes convergent. Applying this limit, and performing the summations whenever it is possible to get a closed form, we obtain the following massless finite result In addition, the first set of summations can be conveniently rewritten in a symmetric fashion with (A.38) While the above expressions are finite, as we just mentioned, their convergence is rather slow, but can be enhanced using some Poisson resummation formulas such as in eq. (A.1). Thus, let us first use eq. (A.28) in order to improve the convergence properties of the two-fold summations, once properly re-expressed using eq. (A.34). Then, we can proceed in implementing a specific resummation to the remaining three-fold summations, following where from the second to the third line, we used eq. (A.1) for the n-and j -summations, and singled out the j = n = 0 term. From the third to the last line, we used the integral representation of the modified Bessel function of the second kind. Bringing all the results together, performing some of the summations to get closed forms when it is possible, and relabeling some of the variables, we finally arrive to the following which is the symmetrized and renormalized result for the free energy density of a noninteracting massless neutral scalar field coupled to a heat bath at temperature T , and geometrically confined inside an infinite tube of rectangular section with two sides of respective lengths L 1 and L 2 .

Finite volume box
Let us finally take care of the last case, analytically continuing eq. (A.25) with c = 3, that is 41) in order to investigate the free energy density of a neutral noninteracting scalar field coupled to a heat bath at temperature T , and geometrically confined inside a finite volume box with sides of respective finite lengths L 1 , L 2 , and L 3 . As previously, notice again that we can use a convenient identity, which in this case reads g n 2 , 0, 2 2 , 2 3 ; m + in order to rewrite the multifold summation in eq. (A.41), and make use of an appropriate representation for computing the free energy density with the help of eq. (A.35). Doing so allows us to obtain an analytic and symmetric representation for corresponding master sum-integral I (3) α , which we do not display here for the sake of brevity. Further applying the relation (A. 22), in order to get the free energy density, we Laurent expand the resulting expression around = 0. Again, skipping some minor technical details, we finally obtain Unlike in the case of the tube geometry, we cannot here simply come back to an earlier version of the massive result such as eq. (A.43). Indeed, even though there is no infrared (or even UV) divergence, and the application of the massless limit upon eq. (A.43) would leave an expression which is overall finite, each of the summations therein produces, in fact, a divergence. And only the sum of all of these divergences actually vanishes. The reason is that when m = 0, each of the sums in eq. (A.43) can be represented by a certain zeta function, and we hit the corresponding pole in the → 0 limit. Therefore, we will implement the limit m → 0, keeping nonzero for the sake of regularizing intermediate stage divergences. Doing so, we find where the structure of the intermediate divergences, which we mentioned previously, is now obvious. Notice further that as the sum of all the poles in identically vanishes, the finite part which is by constructionΛ-dependent will also vanish, leaving the overall result not only finite but also thankfully renormalization scale independent.
Following the massless result of the infinite rectangular tube case, we can again use the Poisson resummation formula (A.1), together with expressions such as eq. (A.28) and eq. (A.34) in order to compute the above two-fold summations while singling out the intermediate stage divergences. In addition, we can use eq. (A.38) in order to deal with the three-fold summations, including those which will appear in the decomposition of the four-fold summation. In the end, we are only left with the computation of the last term in the above expression, and to do so we define And while the last three terms can be handled via the procedure which we used for the three-fold summations in the previous subsection, the first three terms need to be computed separately. We then define another set of summations, namely S 1+3 , which denotes the first of the above terms. Analytically continuing this set of sums, using for example eq. (A.35), leads to the following result Finally bringing all the results together, and expanding around = 0, we see that all theintermediate stage divergences indeed cancel. In addition, the renormalization scale dependence disappears as a consequence of the fact that there is no overall logarithmic UV divergence in the present case. We further perform some of the summations, in order to get closed forms whenever it is possible. Then, relabeling some of the variables, we finally obtain which is the symmetrized and renormalized free energy density of a noninteracting massless neutral scalar field coupled to a heat bath at temperature T , and geometrically confined inside a box with sides of respective finite lengths L 1 , L 2 , and L 3 . We refer to eq. (2.17) for a conveniently much more compact form of this result.

B.1 Fourier decomposition
In this section, we first describe the derivation of the form and normalization of the Fourier decomposition of a massless scalar field (adding a mass does not modify the argument), along the direction of a compactified dimension with DBCs. To this end, we apply the method of separation of variables to the scalar field, and consider only, for the sake of argument, the part relevant to one of the DBCs which must also obey the Klein-Gordon equation We now give some useful identities: It should be noted, indeed, that in order to go from eq. (2.4) to eq. (2.5), the following integral identities, defining the usual Kronecker delta function, are needed

B.2 Length derivatives
We present a short derivation of an important result, i.e. the length derivatives of the internal and free energies, which allows us to compute the pressures directly from the free energy. Consider eq. (3.1) which we differentiate following (B.13)