Classical Casimir free energy for two Drude spheres of arbitrary radii: A plane-wave approach

We derive an exact analytic expression for the high-temperature limit of the Casimir interaction between two Drude spheres of arbitrary radii. Specifically, we determine the Casimir free energy by using the scattering approach in the plane-wave basis. Within a round-trip expansion, we are led to consider the combinatorics of certain partitions of the round trips. The relation between the Casimir free energy and the capacitance matrix of two spheres is discussed. Previously known results for the special cases of a sphere-plane geometry as well as two spheres of equal radii are recovered. An asymptotic expansion for small distances between the two spheres is determined and analytical expressions for the coefficients are given.


Introduction
The Casimir effect is often seen as a quantum effect arising from the vacuum fluctuations of the electromagnetic field between two objects. However, for non-zero temperature T also thermal photons with wavelength λ T =hc/k B T contribute to the Casimir force. In fact, for distances larger than the wavelength λ T , the main contribution to the Casimir force is due to thermal fluctuations. This leads to a finite force even in the classical limit ofh → 0 which in view of the definition of the thermal wavelength is equivalent to the high-temperature limit T → ∞. The Casimir free energy, which then no longer depends on Planck's constant, is found to be linear in temperature. Consequently, the Casimir entropy becomes constant, thereby revealing the entropic origin of the Casimir effect in the classical limit [1].
Within the scattering approach to the Casimir effect [2], the high-temperature limit amounts to taking the zero-frequency term of the Matsubara sum. The associated simplification of the problem has allowed to obtain analytical solutions not only for the archetypal plane-plane geometry [3,4] but also for a scalar field with Dirichlet boundary conditions in the sphere-plane and sphere-sphere geometry as well as for the electromagnetic field in the sphere-plane geometry for boundary conditions corresponding to a Drude metal [5]. Even though it was suspected that the extension of the latter to two spheres of different radii might not be possible [6] we will see in the following that an analytical expression for the Casimir free energy in the general setup of two Drude spheres can be obtained within the scattering approach.
Besides the general theoretical interest in analytical solutions, there is also practical interest in such an expression. While most Casimir experiments so far have been carried out using the sphere-plane geometry, the sphere-sphere geometry has received more attention lately in experiments measuring Casimir forces [7,8] or addressing colloidal systems [9,10]. Carrying out the experiment in an aqueous salt solution offers the opportunity to study the zero-frequency contribution even outside the high-temperature limit by changing the salt concentration [7].
Furthermore, theoretical results in the high-temperature limit can provide a crucial ingredient to a semi-analytical approach [6] useful in the analysis of experimental data. There, the terms for non-zero Matsubara frequencies are treated within the derivative expansion. For the zero-frequency contribution, it is found to be advantageous to employ known exact results available for the sphere-plane geometry [5] and two spheres with equal radii [11]. An exact analytical expression for the setup of two spheres with arbitrary radii will thus be valuable.
It is common to treat geometries involving one or more spheres within a spherical or bispherical multipole expansion. With such approaches the high-temperature limit of two spheres with different radii has not been explicitly derived so far. However, it could have been obtained by combining the results of [5] and [12] together with the capacitance matrix discussed in Section 4.4. In [5], bispherical coordinates where employed to determine the free energy for two Dirichlet spheres while [12] applied field theoretical methods to calculate the difference in free energy of two spheres for Dirichlet and Drude boundary conditions.
Here, we will take a different approach by working in the plane-wave basis which has been shown to allow for interesting physical insights [13,14] as well as an efficient numerical method [15]. Our derivation of the Casimir free energy of two Drude spheres with arbitrary radii will entirely be based on the plane-wave basis. A round-trip expansion of the scattering of electromagnetic waves between the two spheres leads to an interesting combinatorial problem which can be solved. Furthermore, our calculation sheds light on a relation between the scattering approach to the Casimir effect and a problem of electrostatics.
The paper is organized as follows. Section 2 introduces the scattering approach within the plane-wave basis, where we express the Casimir free energy as a sum over round trips between the two spheres. In Section 3 we illustrate the basic idea of our approach by deriving an exact expression for the Casimir free energy of a scalar field. Our result is found to be dual to the known result [5] in the sense that the free energy is obtained by summing over round trips instead of bispherical multipoles. By evaluating the spherical monopole contributions in Section 4 and subtracting them from the free energy of the scalar field, we obtain as our main result an exact expression for the Casimir free energy for an electromagnetic field in the presence of two Drude spheres. It turns out that the monopole contributions can be related to the capacitance matrix of the sphere-sphere geometry [12,16]. Furthermore, we show that our result for the Casimir free energy agrees with previously obtained expressions for the sphere-plane geometry [5] and two spheres of equal radii [11]. Finally, in Section 5, the short-distance expansion for the general sphere-sphere geometry is derived with some technical details relegated to the appendix.

Classical Casimir free energy within the plane-wave basis
We start by compiling all ingredients required to evaluate the Casimir free energy in the hightemperature limit within the scattering approach. The geometry of our sphere-sphere setup is shown in Fig. 1 where the two spheres have generally different radii R 1 and R 2 and are placed at a centre-to-centre distance L = R 1 + R 2 + L. L denotes the smallest distance between the two sphere surfaces. The z-axis is chosen to go through the spheres' centres. Furthermore, the spheres are assumed to be made of a Drude-type metal with a dielectric function for imaginary frequencies ω = iξ with the plasma frequency ω P and the relaxation frequency γ. The dielectric function (1) implies a finite dc conductivity ω 2 P /γ. As a consequence, only the electric modes contribute to the Casimir energy in the high-temperature limit as we will see below.
Within the scattering approach to the Casimir effect [2], the free energy is obtained by summation over terms containing the round-trip operator M at the Matsubara frequencies ξ n = 2πn/k B T . Here, k B and T are the Boltzmann constant and the temperature, respectively. In the high-temperature limit L/λ T 1, only the zero-frequency term is relevant and the Casimir free energy becomes The round-trip operator M describes one complete round trip of the electromagnetic waves between the two spheres and is defined as Figure 1: Representation of the geometry with two spheres of radii R 1 and R 2 . The distance between the spheres is given by L and L defines the separation between the sphere centres.
R 1 and R 2 are the reflection operators for the two spheres while the operators T 12 and T 21 describe the translation between the centres of the spheres. In the following, we omit the argument of the round-trip operator because we will exclusively be concerned with the zerofrequency case. For our purpose, it is convenient to expand the logarithm appearing in (2) into a Mercator series. The Casimir free energy then reads which in physical terms amounts to an expansion in the number r of round trips. In order to evaluate the trace in (4), we have to choose a basis. While it may appear as natural to use spherical [17][18][19] or bispherical [5] multipoles, we found it convenient to make use of a plane-wave basis which has been proven useful lately in the study of the sphere-sphere geometry [13,15].
Specifically, we use the angular spectral representation [20] consisting of plane waves denoted by |k, p, φ . Here, k refers to the projection of the wave vector onto the plane perpendicular to the z-axis. The polarization p can be transverse electric (TE) or transverse magnetic (TM) with respect to the Fresnel plane spanned by the z-axis and the incoming wave vector. Introducing the Wick rotated z-component of the wave vector κ, we obtain from the dispersion relation Since the imaginary frequency ξ is preserved during a round trip, we do not include it in the parameters characterizing the plane-wave basis. Furthermore, ξ = 0 in the high-temperature limit considered here, so that κ = |k|. Finally, φ = ± specifies the direction along the z-axis in which the plane wave decays. φ changes its sign at each reflection.
In the angular spectral representation, the trace of the r-th power of the round-trip operator in the plane-wave basis can now be expressed as trM r = p 1 ,...,p 2r where the indices 2r + 1 and 1 are identified to account for the trace. The exponential factors represent the diagonal matrix elements of the two translation operators covering the distance L between the centres of the spheres. This latter choice allows us to make use of the standard reflection operators with the origin of the reference frame at the spheres' centres. The expression (6) requires the knowledge of the matrix elements of the reflection operator. We concentrate on the results found in the limit of vanishing imaginary frequency ξ and refer the reader to [13] for more details. The matrix elements are obtained from the Mie scattering amplitudes by transforming from the polarization basis referring to the Fresnel plane to the polarization basis referring to the scattering plane. The Mie scattering amplitudes can be expressed in terms of a sum over multipoles and consist of the angle functions τ (cos(Θ)) and π (cos(Θ)) accounting for the scattering geometry and the material-dependent Mie coefficients a and b [21]. For imaginary frequencies, the scattering angle Θ is defined through cos The low-frequency behavior of the electric Mie coefficient for spheres made of a Drude metal is given by a ∼ ξ 2 +1 . The magnetic Mie coefficient b contains an additional power of ξ and can thus be neglected with respect to a . The low-frequency behavior of the two angle functions appearing in the Mie scattering amplitudes is found as τ (cos(Θ)) ∼ ξ −2 and π (cos(Θ)) ∼ ξ −2 +2 . Therefore, in the limit of vanishing ξ, only the combination a τ and thus only the Mie scattering amplitude for waves with polarization lying in the scattering plane contributes.
In the polarization basis taken with respect to the Fresnel plane, it follows that in the zero-frequency limit only the matrix element differs from zero. Here, we have expressed the transverse wave vector k i in polar coordinates through the modulus k i and the angle ϕ i . The sum over the multipoles can be carried out and the non-vanishing reflection matrix becomes Note the subtraction of 1 because of the missing monopole term = 0 in (7) which distinguishes the electromagnetic from the scalar case. After inserting the reflection matrix element (8) into the expression (6) for the trace, it is convenient to switch to Cartesian coordinates Here, ρ n = R n /L denotes the dimensionless radius of spheres n = 1, 2 and the argument of the hyperbolic cosines is abbreviated by χ The trace over the r-th power of the round-trip operator is now given by a sum over 2r-dimensional Gaussian integrals. After having determined the matrices associated with the bilinear forms in the exponentials, our main task will be to evaluate the corresponding determinants.
As already remarked above, the subtraction of 1 in the last two factors in (9) arises because the monopole term does not contribute in the case of electromagnetic waves. Including the monopole term amounts to considering the case of a scalar field with Dirichlet boundary conditions on the spheres. In the literature [5], it has been found useful to first evaluate the scalar case and then to determine the correction corresponding to the monopole contributions. In the next section, we will thus consider the scalar case. The plane-wave approach will lead us to an expression for the Casimir free energy which is equivalent to the known result [5].

Scalar field with Dirichlet boundary conditions
According to the discussion in the previous section, the trace over the r-th power of the round-trip operator for a scalar field and two spheres with Dirichlet (D) boundary conditions can be expressed in the plane-wave basis as Expanding the product, one obtains a sum over 2 2r Gaussian integrals where the bilinear form in the exponent can be written with the help of the 2r-dimensional symmetric matrix where all combinations of signs appear in the expansion of the product in (10). As far as the determinant is concerned, it only matters whether the number of minus signs in the upper or lower half of the matrix is even or odd as indicated by the superscript + or − of M ± r , respectively.
For a single round trip, r = 1, the determinant reads where characterizes the geometry of the sphere-sphere arrangement with the effective radius R eff = . For a general number of round trips, it can be useful to view (11) as the Hamiltonian matrix of a periodic tight-binding model and to reexpress the problem in terms of transfer matrices [22,23]. One then finds with µ = arcosh(y). (15) We are now in a position to evaluate the Gaussian integrals in (10). Noting that the bilinear form in the exponent is given by M + r and M − r in half of the terms each, we find with (4) the Casimir free energy for a scalar field and Dirichlet boundary conditions in the high-temperature limit as a sum over round trips Our result (16) can be viewed as a dual representation of the earlier result presented in [5]. Following the notation introduced there, we define and write the Casimir free energy as Expanding the last factor in a geometric series, the sum over r can be evaluated. With the help of the Mercator series, one finally obtains in agreement with the result derived by means of bispherical coordinates [5]. Particularly for small distances, the round-trip representation (16) may be numerically advantageous as compared to the expansion (19) because it possesses better convergence properties. It is also straightforward to read off the Casimir free energy within the proximity force approximation by simply retaining the leading order of the hyperbolic functions. Details of the asymptotic expansion in µ of the Casimir free energy will be discussed in section 5.
4 Electromagnetic case for two general Drude spheres 4

.1 Monopole contributions in the scalar case
The trace over the r-th power of the round-trip operator in the electromagnetic case differs from the scalar case only by the monopole term = 0 as one can see by comparing the corresponding expressions (9) and (10). As already discussed at the end of section 2, does the −1 in the brackets of (9) account for the subtraction of the monopole term. Hence, when expanding the product in (9), all summands containing at least one of these terms −1 yield the negative monopole contributions present in the scalar case. In the following, we will thus focus on the difference Already in previous works it was found convenient to study the difference between the scalar case with Dirichlet boundary conditions and the electromagnetic case for Drude-type objects [5,12]. The difference ∆ r consists of a sum over Gaussian-type integrals where the bilinear form in the exponent is represented by a tridiagonal matrix with the off-diagonal matrix elements arising from the hyperbolic cosines. Whenever in the expansion of the product in (9) a factor −1 appears, the corresponding pair of off-diagonal matrix elements vanishes. In contrast to the matrix (11) in the scalar case, the matrix representing the bilinear form in the exponent of the integrand in (9) is now block-diagonal and can be written as where w denotes an element of a set Π 2r,k containing a multiset of tuples {(n 1 , t 1 ), (n 2 , t 2 ), (n 3 , t 3 ), . . . , (n k , t k )} with i n i = 2r for r round trips. Each block is a symmetric tridiagonal 2-Toeplitz matrix [24] of the form where pairs of off-diagonal matrix elements alternate between ρ 1 and ρ 2 and each pair can come with an arbitrary sign. Each block is characterized by its size n and the index of the first off-diagonal entry, indicated by the superscript 1 or 2 and accounted for by t i in the multiset w. As we will see in more detail later, we cannot set t i to 1 or 2 freely. Rather, its value needs to be compatible with the values of n i−1 and t i−1 . The result of the Gaussian integration will involve the determinant of M w which equals the product of the determinants of the individual blocks. For odd dimension, the determinant of the blocks (22) is given by [25] det m while for even dimension one has to distinguish between blocks starting with ρ 1 or ρ 2 on the off-diagonal det m Here, U k denotes the Chebyshev polynomial of the second kind and order k while y has been introduced in (13) and characterizes the geometry of the sphere-sphere setup. We note that the determinants do not depend on the choice of signs in (22) so that in view of the Gaussian integration we end up with 2 n−1 equivalent blocks m (1/2) n . For the monopole contributions (20) we now obtain together with (9) and (10) where the number of blocks in M w is given by the summation index k and x t denotes the transpose of x. In the derivation, we have taken into account that a factor −2 is associated with each block as one can see by evaluating the product. Furthermore, we have accounted for the multiplicity related to the signs in the block matrices mentioned above. Evaluating the Gaussian integrals, we arrive at

Combinatorics of blocks and diagrammatic representation
The decomposition of the block matrix M w into blocks can be conveniently analyzed in terms of a diagrammatic representation. For r round trips, we consider a graph consisting of a chain of 2r + 1 nodes where the last node should be identified with the first one. These nodes represent the two spheres and are depicted successively in black and white corresponding to spheres 1 and 2, respectively. Each block m (1/2) n is represented by a black or white line where the colors refer to the superscripts 1 and 2, respectively. Therefore, the color of a line equals the color of the node from which the line starts, reading the diagram from left to right. The length of the line is given by the dimension n of the corresponding block.
The connection between the block matrix and its diagrammatic representation is illustrated in Fig. 2. In this example, the round trip starts on sphere 1 and the block m represents half a round trip ending on sphere 2. Given the odd dimension of the first block, the color switches to white. The next block, m 3 , has again an odd dimension and corresponds to one and a half round trips. In general, it is not required that matrices of odd dimension follow each other directly. We are now back to sphere 1 and it follows a black line symbolizing the block m Disregarding the color for a second, there exists an obvious connection to the problem of integer partition. It is well known that so-called ordinary Bell polynomials defined through provide such a partition [26]. We will have the opportunity to employ this relation later when deriving the short-distance behavior of the Casimir free energy. However, for the following discussion we will need to keep the color because according to (24) the determinant of our blocks for even dimension depends on the superscript.
To the best of our knowledge, a generalization of the Bell polynomials to our situation with color does not exist. Fortunately, the dependence of the color of a line on the length and color of the previous line allows us to express the partitions in a recursive way instead.
Let us consider r round trips, i.e. a chain of length 2r. We introduce functions h 2r and h (2) 2r as sums over the inverse determinants of all possible block matrices for r round trips starting on spheres 1 and 2, respectively. These functions will allow us later to express the monopole contributions ∆ r as given by (26). For convenience, we define abbreviations for the inverse of the determinants of the blocks with n ≥ 1 a n = 1 det m which will occur in ∆ r . In our diagrams, the coefficients a n are thus associated with black lines and the coefficients b n with white ones. We can now express the first of the two required functions, which starts on sphere 1, recursively as h (1) where the variable t will later serve to determine the number of blocks. The first term on the right-hand side of (29) accounts for a single block of maximal size while the other two terms correspond to matrices consisting of more than one block. The second term arises from a single block of even dimension followed by a block matrix starting again from sphere 1. The sum runs over all possible sizes of the first block. The relative minus sign between the first and the second term is due to the fact that each new block contributes a minus sign leading to the factor (−1) k for k blocks in (26). The third term differs from the second one in so far as the first block has an odd dimension so that the remaining part starts on sphere 2. Instead of h (1) in the second term, we thus have h (2) in the third term. To close the system of recursive equations, we similarly derive three more equations and for graphs starting on sphere 2 h 2(r−n)+1 (t) and h (2) We note in passing that these recursion relations can be interpreted as the Laplace expansion of appropriately chosen Hessenberg matrices but we will not make use of this fact in the following.
The sum h 2r accounts for all different kinds of block matrices occurring for r round trips. However, we have not yet properly accounted for the multiplicity of the block matrices. So far, we have considered open-chain diagrams with the starting point chosen at a specific node for which a pair of off-diagonal elements vanishes. Since the trace at the origin of (26) implies a closed chain, we can have several starting points. As Fig. 3 demonstrates for r = 3, a configuration can start at three different nodes. Specifically, the white line can start on one of the three white nodes. Since, in general, there are r nodes of one color, we have to multiply the contribution for r round trips by a factor of r.
By proceeding as just described, we include all circular permutations of the blocks, but the functions h (1,2) 2r already include all non-equal circular permutations. Hence, to avoid double counting, we have to remove the k cyclic permutations of a partition in k blocks by dividing the contribution arising from k blocks by k. Fig. 3 illustrates a non-trivial example. On the left, we present the diagrams of the partition diag(m 1 ) and all its circular permutations. As demonstrated by the graphs on the right, there are r = 3 possible ways of choosing a starting point for each partition. However, as one can see, each diagram appears four times. Hence, by dividing by four, the correct number of block matrices is obtained.
We now make use of the parameter t introduced in equations (29)-(32) which ensures that contributions arising from k blocks come with a factor t k . The monopole contributions for a given number of round trips (26) can thus be written as The negative sign arises due to our definition of h (1,2) 2r , where all partitions in odd numbers of blocks occur with a positive sign compared to those with even numbers.
So far, we have considered the monopole contributions for a given number of round trips. The full correction of the Casimir free energy with respect to the result (16) or (19) for the scalar case is obtained by summation over all numbers r of round trips, which we will carry out in the following section.

Monopole correction to the classical Casimir free energy
According to the round-trip expansion (4) of the Casimir free energy, the monopole contributions are determined by Strictly speaking, ∆ equals the negative monopole contributions, but for simplicity we will continue to refer to this quantity as monopole contributions.  (2) 1 ), we reproduce all four configurations. On a closed chain the four permutations appear naturally so that each diagram occurs four times on the right-hand side.
After inserting (33) and interchanging summation and integration, it is convenient to introduce the generating functions for the inverse block-matrix determinants In the second equality, we decompose the sum into contributions from even (e) and odd (o) powers of x. Hence, the monopole term yields Correspondingly, we introduce the generating functions for the inverse determinants of individual blocks (28) and where we apply the same decomposition of the functions as in (35). The generating functions H (1,2) e can be determined by summing over the recurrence relations (29)-(32) and we find with an analogous expression for H Noting that the numerator equals the derivative with respect to t of the denominator, it is straightforward to evaluate the integral in (36) and we find where α = R 2 /R 1 and β = R 1 /R 2 take the ratios of the sphere radii into account and y is defined in (13). The expressions in (42)-(44) are obtained by inserting (28) together with (23) and (24) into (37) and (38). It is instructive to convince oneself that indeed all partitions of round trips are contained in (41) by expanding the logarithm as The first sum accounts for arbitrary repetitions of full round trips starting either on sphere 1 or on sphere 2 as represented by A e or B e , respectively. Expressions containing both A e and B e can only arise if half round trips represented by A o and B o occur as is the case in the second term. Reading this term from left to right, it can clearly be seen that half a round trip induces a change between full round trips starting on sphere 1 and on sphere 2. The number of factors (−1) correctly reflects the number of blocks in the matrices M w .

Relation to the capacitance matrix
It appears that the result (41) was so far not known in the Casimir community. Nevertheless, it can be obtained by combining results from the literature, a fact which we only became aware of after the work presented here had been carried out. As was shown by Fosco et al. [12], the difference between the Casimir free energy of objects made of Drude metals and the Casimir free energy for a scalar field with Dirichlet boundary conditions is related to the capacitance matrix C of the arrangement of conductors. For the special case of two conductors, [12] Even though this was not mentioned in [12], the capacitance matrix elements of two conducting spheres of arbitrary radii were already known to Maxwell [27]. Following the more modern notation in [16], the capacitance coefficients can be expressed as Comparing these coefficients and (46) with our result (41) connects the capacitance coefficients to the scattering of electromagnetic waves in the static limit. It thus highlights the relation between our round-trip description and the method of image charges used by [27] to obtain the capacitance coefficients. We remark that the general result (46) and our result (41) differ by a factor R 1 R 2 T 2 in the logarithm. While this factor would be irrelevant for the Casimir force, it makes a difference for the Casimir entropy. Its origin can be traced back to the different handling of the Casimir free energy of the individual objects [1,28]. While the scattering approach does not contain the free energy of the spheres at an infinite distance, this contribution is present in [12]. For (41), the entropy in the high-temperature limit becomes a constant as expected [1].

Casimir free energy for two Drude spheres of general radii and limiting cases
According to (34), the sum of the expressions (16) and (41) gives the Casimir free energy for two Drude spheres of arbitrary radii and thus constitutes the main result of this paper. Instead of reproducing the two expressions here, it is useful to resum the result as we did in Section 3 for the scalar case and to express it in terms of the variable Z introduced in (17).
Noting that the Chebyshev polynomials of the second kind appearing in equations (42)-(44) can be written as we obtain the classical Casimir free energy for two Drude spheres as where we have introduced the function and correspondingly for g β (Z), where α and β = 1/α are the ratios of sphere radii as defined below (44). We obtain the limit of a sphere of radius R in front of a plane by setting R 2 = R and letting R 1 go to infinity. Then, g β = 1 and B e and B o vanish because there is no second sphere were the electromagnetic waves could be scattered. Since the functional dependence of the scalar part of (51), i.e. the first sum, is not affected, we focus on the monopole contributions for which we obtain where Z depends only on the aspect ratio = L/R through By some minor transformations, one can convince oneself, that (53) agrees with the result found earlier by Bimonte and Emig [5].
Similarly, we obtain the Casimir free energy for two Drude spheres of equal radii by setting R 1 = R 2 = R so that g α = g β = Z 1/2 = Y . In this case, the scattering at the two spheres cannot be distinguished and we have A e = B e and A o = B o . The monopole contributions then read where the parameter Y is a function of the aspect ratio δ = L/2R The result (55) leads to the same Casimir free energy as obtained earlier by using the transformation optics approach [11].

Short-distance expansion
In experiments, the closest distance L between the two spheres is typically small compared to the radii R 1 and R 2 . Therefore, we will now determine a short-distance expansion of the Casimir free energy (51) by separately considering the scalar contribution F (D) and the monopole contributions ∆. The leading-order term will correspond to the proximity-force approximation whose validity can be assessed by the higher-order terms.
In the following, we make use of the fact that F (D) and ∆ can be expressed in terms of the F-series introduced by Garvin [29] as a generalization of the Lambert series. With a choice of coefficients appropriate for our situation, we introduce Here, we follow the notation used by Banerjee and Wilkerson who provide an asymptotic expansion of this series around q = 1 [30]. We start by expanding the Casimir free energy for two Dirichlet spheres. Its representation (18) can be expressed in terms of the series (57) as For small distances L R 1 , R 2 , the variable Z = exp(−µ) is close to unity and µ defined in (15) is small. In the following, we will use µ as our expansion variable.
Making use of the asymptotic expansion of the generalized Lambert series around q = 1 stated in theorem 2.2 of [30], we obtain from (58) for the Casimir free energy in the scalar case with Dirichlet boundary conditions with Glaisher's constant A = 1.28242 . . . and the Bernoulli numbers B k [31]. The first term corresponds to the high-temperature result of the proximity-force approximation with the specific value of the Riemann zeta function ζ(3) = 1.20205 . . . The terms up to order µ 4 were already given in [5] and are consistent with our result. Now we turn to the short-distance expansion of the monopole term as given by the second term of (51) where the argument of the logarithm can again be expressed in terms of a generalized Lambert series (57). Readers not interested in the technical details of the derivation will find the final result in (77).
We bring the sums in the monopole contributions ∆ into the form of a generalized Lambert series by writing the function (52) as where Replacing α by β simply changes the sign of this function, so that g β = Z 1/2−v . For our purpose, we need to expand v(µ) into a Taylor series The coefficients for n > 0 can be expressed in terms of where S(n, k) are the Stirling numbers of the second kind. Before making use of the generalized Lambert series, it is convenient to introduce a notation for the prefactors and sums appearing in the argument of the second logarithmic term in (51). By defining and the monopole contributions (41) can be brought into the form where we introduced the abbreviations J (±) = J(3/2 ± v) and I (±) = I(3/2 ± v). For the further analysis, we now express the series (66) in terms of the generalized Lambert series (57) as Applying the asymptotic expansion of the generalized Lambert series [30], we obtain with the digamma function ψ(c) and the Bernoulli polynomial B 2n (c) [31]. For c = 1, the functions ψ(c) and B 2n (c) are given by the negative Euler-Mascheroni constant, −γ = −0.57721 . . . and the Bernoulli numbers B 2n , respectively. Before proceeding with our calculation, we note that a short-distance expansion of the capacitance coefficients for two general spheres has already been carried out in [32] and [33], where the latter one also applied the asymptotic expansion of the generalized Lambert series. Besides using a different definition of the dimensionless capacitance coefficients and geometric parameters, we determine, in contrast to previous work, a complete expansion of the functions I(c) as well as J(c) in powers of µ. In the definition of our geometric parameters, we follow the notation common in the Casimir community which also simplifies to obtain the limits of equal spheres and of the sphere-plane geometry.
In order to obtain a complete expansion of the argument of the logarithm in (67) in powers of µ, we need to account for the fact that I (±) and J (±) depend on µ through c = 3/2 ± v(µ). By making use of the Taylor series for v given in (62) with the coefficients (63) and (64), one immediately obtains a corresponding Taylor series for c(µ) which is required to determine the Taylor series in powers of µ for the digamma function ψ(c(µ)) and the Bernoulli polynomial B 2n (c(µ)) appearing in (69).
The expansion coefficients n (u) and δ n (u) are defined in (104)-(109) in Appendix A and depend only on the geometric parameter This parameter can take arbitrary values between 0 and 1/4, corresponding to the sphere-plane geometry and equally sized spheres, respectively. Table 1 gives the values of the expansion coefficients for the two limiting cases. The sphere-plane limit, u = 0, is consistent with the results given in [5]. Their numerical constants γ i , i = 1, 2, 3, 4 can now be expressed analytically as  Table 1: Expansion coefficients n (u) and δ n (u) appearing in (77) in the limits of the sphereplane geometry (u = 0) and of equal spheres (u = 1/4).

Conclusions
We have for the first time derived an exact analytical expression for the Casimir free energy of two Drude spheres of arbitrary radii completely within the scattering approach common in Casimir physics. In contrast to previous work on the sphere-plane geometry and two spheres of equal radii, the plane-wave basis was used, which led to a connection with a combinatorial problem. The structure of this combinatorial problem highlights the difference between the general two-sphere case and the corresponding limiting cases. The scattering approach also provides an intuitive interpretation of the structure of the result for the Casimir free energy. Earlier work by Fosco et al. [12] has pointed out the relevance of the capacitance matrix for the Casimir free energy in the high-temperature limit. However, the fact that an analytical expression for the capacitance matrix exists even for two spheres of different radii seems to have largely escaped the attention of the Casimir community. Our work thus provides an interesting connection between Casimir physics and electrostatics. This is in particular the case for the short-distance expansion where by profiting from results obtained within the electrostatics community, we derived a systematic expansion in powers of µ which might be useful in that community as well.