{'fraction_non_alphanumeric': 0.0474335790791487, 'fraction_numerical': 0.0162748643761302, 'mean_word_length': 4.584077669902912, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of nonzero-temperature Casimir forces in the time domain

We show how to compute Casimir forces at nonzero temperatures with time-domain electromagnetic simulations, for example using a finite-difference time-domain (FDTD) method. Compared to our previous zero-temperature time-domain method, only a small modification is required, but we explain that some care is required to properly capture the zero-frequency contribution. We validate the method against analytical and numerical frequency-domain calculations, and show a surprising high-temperature disappearance of a non-monotonic behavior previously demonstrated in a piston-like geometry.

We show how to compute Casimir forces at nonzero temperatures with time-domain electromagnetic simulations, for example using a finite-difference time-domain (FDTD) method. Compared to our previous zero-temperature time-domain method, only a small modification is required, but we explain that some care is required to properly capture the zero-frequency contribution. We validate the method against analytical and numerical frequency-domain calculations, and show a surprising high-temperature disappearance of a non-monotonic behavior previously demonstrated in a piston-like geometry.
In this paper, we show how to compute nonzerotemperature (T > 0) corrections to Casimir forces via time-domain calculations, generalizing a computational approach based on the finite-difference time-domain (FDTD) method that we previously demonstrated for T = 0 [1,2]. New computational methods for Casimir interactions have become important in order to model nonplanar micromechanical systems where unusual Casimir effects have been predicted [1][2][3][4][5][6][7][8][9][10][11][12], and there has been increasing interest in T > 0 corrections [13][14][15][16][17][18][19], especially in recently identified systems where these effects are non-negligible [16]. Although T > 0 effects are easy to incorporate in the imaginary frequency domain, where they merely turn an integral into a sum over Matsubara frequencies [20], they turn out to be nontrivial to handle in the time-domain because of the singularity of the zero-frequency contribution, and we show that a naive approach leads to incorrect results. We validate our approach both with a one-dimensional system where analytical solutions are available, and also in a two-dimensional (2D) piston-like geometry [4,8,21,22] where we compare to a frequency-domain numerical method. In the 2D piston geometry, we observe an interesting effect in which a non-monotonic phenomenon previously identified at T = 0 disappears for a sufficiently large T .
The Casimir force is a combination of fluctuations at all frequencies ω, and the T = 0 force can be expressed as an integral F (0) = ∞ 0 f (ξ)dξ over Wick-rotated imaginary frequencies ω = iξ [20]. At a nonzero T , this integral is replaced by a sum over "Matsubara frequencies" ξ n = nπω T for integers n, where ω T = 2k B T / and k B is Boltzmann's constant [20]: The transformation from the T = 0 integral to a summation can be derived directly by considering thermodynamics in the Matsubara formalism. Eq. (1) corresponds to a trapezoidal-rule approximation of the T = 0 integral [23]. At room temperature, the Matsubara frequency corresponds to a wavelength 2π/ξ = 7 µm, much larger than separations where the Casimir effect is typically observed, so usually T > 0 corrections are negligible [13,14]. However, experiments are pushing towards > 1µm separations [24][25][26] in an attempt to observe these corrections. Also, a recent theoretical prediction shows much larger T corrections with appropriate material and geometry choices [16].
To compute Casimir forces in arbitrary geometries, it is desirable to exploit mature methods from computational classical electromagnetism (EM), and a number of approaches have been suggested [27][28][29][30][31]. One technique is to use the fluctuation-dissipation theorem: the mean-square electric and magnetic fields E 2 and H 2 can thereby be computed from classical Green's functions [20], and the mean stress tensor can be computed and integrated to obtain the force [1,2,8]. In particular, at each ω, the correlation function of the fields E 2 is given by: where G E jk = (G E k ) j is the classical dyadic "photon" Green's function, proportional to the electric field in the j direction at x due to an electric-dipole current in the k direction at x ′ , and solves where ε is the electric permittivity tensor, µ is the magnetic permeability tensor, andê k is a unit vector in direction k. The magnetic-field correlation H 2 has a similar form [1,2]. Note that the temperature dependence appears as a coth factor (from a Bose-Einstein distribution). If this is Wick-rotated to imaginary frequency ω = iξ, the poles in the coth function gives the sum (1) over Matsubara frequencies [32]. In our EM simulation, what is actually computed is the electric or magnetic field in response to an electric or magnetic dipole current, respectively. This is related to G ij by where E jk (ω, x, x ′ ) denotes the electric field response in the jth direction due to a dipole current source This equation can be solved for each point on a surface to integrate the stress tensor, and for each frequency to integrate the contributions of fluctuations at all frequencies. Instead of computing each ω separately, one can use a pulse source in time, whose Fourier transform contains all frequencies. As derived in detail elsewhere [1,2], it turns out that this corresponds to a sequence of time-domain simulations, where pulses of current are injected and some function Γ(t) of the resulting fields (corresponding to the stress tensor) is integrated in time, multiplied by an appropriate weighting factor g(t). We perform these simulations by using the standard FDTD technique [27], which discretizes space and time on a uniform grid. In frequency domain, Wick rotation to complex ω(ξ) is crucial for numerical computations in order to obtain a tractable frequency integrand [8,23], and the analogue in time domain is equally important to obtain rapidly decaying fields (and hence short simulations) [1,2]. In time domain, one must implement complex ω indirectly: because ω only appears explicitly with ε in Eq. (3), converting ω to the complex contour ω(ξ) ≡ ξ 1 + iσ ξ is equivalent to operating at a real frequency ξ with an artificial conductivity ε(r) → ε(r)(1 + iσ ξ ) [1,2]. One cannot use purely imaginary frequencies ω = iξ in the time domain, because the corresponding material ε → −ε has exponentially growing solutions in time [1]. Thus, by adding an artificial conductivity everywhere, and including a corresponding Jacobian factor in g(t), one obtains the same (physical) force result in a much shorter time (with the fields decaying exponentially due to the conductivity). Now, we introduce the basic idea of how T > 0 is incorporated in the time domain, and explain where the difficulty arises. The standard T > 0 analysis of Eq. (1) is expressed in the frequency domain, so we start there by exploiting the fact that the time-domain approach is derived from a Fourier transform of the frequency-domain approach. In particular, g(t) is the Fourier transform of a weighting factor g[ω(ξ)] in the Fourier domain [1,2]. At real frequency, the effect of T > 0 is simply to include an additional weighting factor coth ω(ξ) ωT in the ω(ξ) integral from Eq. (2). So, a straightforward, but naive, approach is to replace g[ω(ξ)] with: using the T = 0 g [ω(ξ)] expression from Ref. 1, and then Fourier transform this to yield g(t). However, there is an obvious problem with this approach: the 1/ω singularity in coth ω(ξ) ωT means that Eq. (5) is not locally integrable around ξ = 0, and therefore its Fourier transform is not well-defined. If we naively ignore this problem, and compute the Fourier transform via a discrete Fourier transform as in [1,2], simply assigning an arbitrary finite value for the ξ = 0 term, this unsurprisingly gives an incorrect force for T > 0 compared to the analytical Lifshitz formula for the case of parallel perfect-metal plates in 1D [33], as shown in Fig. 1 (green dashed line).
Instead, a natural solution is to handle ω = 0 by the coth factor as in Eq. (5), but to subtract the ω = 0 pole and handle this contribution separately. As explained below, we will extract the correct ω = 0 contribution from the frequency-domain expression Eq. (1), convert it to time domain, and add it back in as a manual correction to g(t). In particular , the coth ω(ξ) ωT function has poles at ω = inπω T for integers n. When the frequency integral is Wick-rotated to imaginary frequency, the residues of these poles give the Matsubara sum Eq. (1) via contour integration [32]. If we subtract the n = 0 pole from the coth, obtaining the result of the time-domain integration of g n>0 (t)Γ(t) will therefore correspond to all of the n > 0 terms in Eq. (1), nor is there any problem with the Fourier transformation to g n>0 (t). Precisely this result is shown for the 1D parallel plates in Fig. 1, and we see that it indeed matches the n > 0 terms from the analytical expression.
We check Eq. (11) for the 1D parallel plate case in Fig. 1 against the analytical Lifshitz formula [33]. As noted above, the g n>0 term (6) correctly gives the n > 0 terms, and we also see that the σk B T term gives the correct n = 0 contribution, and hence the total force is correct.
As another check, we consider a more complicated geometry: a piston-like configuration from Ref. 4, shown schematically in the inset of Fig. 2. This system consists of two square rods adjacent between two sidewalls, which we solve here for the 2D case of z-invariant fluctuations. At T = 0, such geometries were shown to exhibit an interesting non-monotonic variation of the force between the two blocks as a function of sidewall separation d [4,8,21,22], which does not arise in the simple pairwiseinteraction heuristic picture of the Casimir force. This can be seen in the solid lines of Fig. 2, where the nonmonotonicity arises from a competition between forces from transverse-electric (TE) and transverse-magnetic (TM) field polarizations [4], which in turn can be explained by a method-of-images argument [21]. In Fig. 2, the solid lines are computed by a T = 0 frequencydomain boundary element method (BEM) evaluating a path-integral expression [9], whereas the circles are computed by the T = 0 FDTD method [1,2], and both methods agree. We also compute the force at T = 1×πc /k B a where the ξ = 0 + term dominates. We see that the FDTD method with the T > 0 modification Eq. (11) (diamonds) agrees with the frequency-domain BEM results (dashed lines), where the latter simply use the Matsubara sum (1) to handle T > 0.
Interestingly, Fig. 2 shows that the non-monotonic effect disappears for T = 1×πc /k B a, despite the fact that the method-of-images argument of Ref. 21 ostensibly applies to the ξ = 0 + quasi-static limit (which dominates at this large T ) as well as to ξ > 0. The argument used the fact that TM fluctuations can be described by a scalar field with Dirichlet boundary conditions (vanishing at the metal), and in this case the sidewalls introduce oppositesign mirror sources that reduce the interaction as d decreases; in contrast, TE corresponds to a Neumann scalar field (vanishing slope), which requires same-sign mirror sources that increase the interaction as d decreases [21]. In Fig. 2, however, while the T = 1 × πc /k B a TM force still decreases as d decreases, the TE force no longer increases for decreasing d at T = 1×πc /k B a. The problem is that the image-source argument most directly applies to z-directed dipole sources in the scalar-field pictureelectric J E z currents for TM and magnetic J H z currents for TE-while the situation for in-plane sources (corresponding to derivative of the scalar field from dipole-like sources) is more complicated [34]. For a sufficiently large T dominated by the ξ = 0 + contribution (as is the case here), we find numerically that the J H z sources as ξ → 0 + no longer contribute to the force. Intuitively, as ξ → 0 + a magnetic dipole source produces a more and more constant (long wavelength) field, which automatically satisfies the Neumann boundary conditions and hence is not affected by the geometry. Instead, numerical calculations show that the TE ξ = 0 + contribution is dominated by J E x sources and the corresponding electric stress-tensor terms, which turn out to slightly decrease in strength as d decreases. (A related effect is that, for small d, it can be observed in Fig. 2 that the T = 1 force is actually smaller than the T = 0 force, again due to the suppression of the TE contribution. Since the force diverges as T → ∞, this means that the force changes non-monotonically with T at small d; a similar non-monotonic temperature depen-dence was previously observed for Dirichlet scalar-field fluctuations in a sphere-plate geometry [19].) In contrast, if we consider the 3D constant crosssection problem with z-dependent fluctuations, corresponding to integrating e ikz z fluctuations over k z [20], then we find that the non-monotonic effect is preserved at all T . This is easily explained by the fact that, for perfect metals, k z = 0 is mathematically equivalent to a problem at k z = 0 and ξ → ξ 2 + k 2 z [8,35], and so the n = 0 Matsubara term still contains contributions equivalent to ξ > 0 in which the J H z mirror argument applies and the situation is similar to T = 0. In any case, this 2D disappearance of non-monotonicity seems unlikely to be experimentally relevant, because we find that it only occurs for T 0.7 × πc /k B a , which for a = 1µm corresponds to T 5000K.
The main point of this paper is that a simple (but not too simple) modification to our previous time-domain method allows off-the-shelf FDTD software to easily calculate Casimir forces at nonzero temperatures. Although the disappearance of non-monotonicity employed here as a test case appears unrealistic, recent predictions of other realistic geometry/material effects [16], combined with the fact that temperature effects in complex geometries are almost unexplored at present, lead us to hope that future work will reveal further surprising temperature effects that are observable in micromechanical systems.