PaperlayerView exact revision on arXiv

Source-linked AI summary

Numerical methods for computing Casimir interactions

Steven G. Johnson

arXiv:1007.0966v2quant-phphysics.comp-ph

TL;DR

The paper addresses how to compute Casimir interactions for arbitrary geometries and materials when practical general methods were historically unavailable. It reviews classical-electromagnetism-based formulations and organizes them by problem, basis, and solution technique, emphasizing complex-frequency adaptations. The review concludes that several methods have been demonstrated, but their differing strengths and weaknesses preclude absolute rankings.

  • Problem

    Practical numerical methods for Casimir forces in arbitrary geometries and materials were historically unavailable, motivating systematic computational approaches.

  • Method

    The review surveys Casimir computational formulations through their relationships to classical electromagnetism and classifies them by problem, basis, and solution technique.

  • Results

    Several effective numerical methods have been proposed and demonstrated, while different approaches retain distinct strengths and weaknesses.

  • Takeaways & Limitations

    Classical electromagnetic methods can be adapted to Casimir physics, with complex- or imaginary-frequency and wide-bandwidth requirements receiving particular emphasis.

  • Takeaways & Limitations

    Sharp corners can limit finite-difference accuracy to O(∆x^p) for some 1 < p < 2.

Abstract

from arXiv · show

We review several different approaches for computing Casimir forces and related fluctuation-induced interactions between bodies of arbitrary shapes and materials. The relationships between this problem and well known computational techniques from classical electromagnetism are emphasized. We also review the basic principles of standard computational methods, categorizing them according to three criteria---choice of problem, basis, and solution technique---that can be used to classify proposals for the Casimir problem as well. In this way, mature classical methods can be exploited to model Casimir physics, with a few important modifications.

7.1 Introduction

Casimir computations historically lacked practical general methods, but recent approaches adapt classical electromagnetism to arbitrary geometries and materials. This review compares methods without seeking a universally best choice.

  • Until recently, practical numerical methods for arbitrary geometries and materials were unavailable, limiting solutions to special geometries and often perfect metals.
  • The review surveys methods for Casimir forces, energies, torques, and related interactions through their relationships to classical electromagnetic techniques.
  • It compares approaches along quasi-independent axes to clarify their strengths, weaknesses, and applicability rather than identify a single best method.
  • Adapting classical solvers requires complex- or imaginary-frequency Maxwell solutions and wide-bandwidth analyses.

7.2 Characterization of numerical methods: Three axes

Casimir numerical methods can be organized by the problem posed, the basis used to discretize unknowns, and the technique used to solve the resulting equations. These choices connect Casimir computations to established classical-EM methods while imposing frequency-integration requirements.

  • Methods differ by the problem solved, the basis expressing unknowns, and the solution technique determining their coefficients.
  • Casimir formulations include eigenfrequency summation, frequency-domain field correlations from Green’s functions, and path-integral interaction energies.
  • All formulations require integration or summation over frequencies, with causal analyticity enabling contour choices in the upper-half complex-frequency plane.
  • PDE approaches use volumetric unknowns, integral equations use surface unknowns, and hybrid eigenmode methods propagate fields between piecewise-constant cross-sections.
  • Finite differences use grid values, FEM/BEM use localized element bases, and spectral methods use truncated non-localized complete bases.
  • Finite differences favor simplicity but uniform resolution and lower-order convergence, whereas FEM and spectral methods offer adaptive or high-order convergence at greater complexity.
  • For large or nonsparse systems, iterative methods can use O(N) storage; conjugate gradient is especially suitable when the matrix is Hermitian positive-definite.

7.3 The impracticality of eigenmode summations

Although Casimir energy can be expressed as a sum over electromagnetic zero-point modes, direct eigenmode summation is numerically impractical for general systems because oscillatory contributions cancel across frequencies.

  • At zero temperature, the Casimir energy can be represented as the sum of zero-point energies of all electromagnetic modes.
  • The computational strategy is to calculate eigenfrequencies numerically, sum them for the energy, and differentiate the result to obtain forces.
  • Discretization makes the number of eigenfrequencies finite, with the grid’s Nyquist-like frequency providing an upper bound.
  • Even after differentiation, the mode summand is wildly oscillatory and receives substantial contributions from nearly every frequency that mostly cancel.
  • Direct dense eigensolvers require O(N^3) time and O(N^2) storage, making eigenmode summation impractical beyond simple one-dimensional problems.
  • Contour-integral reformulations connect mode summation to Green’s functions, Wick rotation, and imaginary-frequency scattering-matrix techniques.

7.4 The complex-frequency plane and contour choices

Casimir frequency integrals can be deformed into the upper complex-frequency half-plane because of causality. Wick rotation is computationally essential: real-frequency contributions are oscillatory and broadband, while imaginary-frequency contributions decay smoothly.

  • For perfect-metal plates, the force contribution f(ω) is integrated over frequency, while finite-permittivity generalizations require numerical integration.
  • Wick rotation converts the frequency integral to imaginary frequencies ω = iξ, a form used in practical Casimir computations.
  • Causality makes the relevant Green’s functions analytic in the upper-half complex-frequency plane, allowing the integration contour to be deformed without changing the result.
  • Near the real-frequency axis, contributions are singular or badly behaved, with strong oscillations, slow decay, and difficult broadband cancellation.
  • In the complex-frequency plot, contours parallel to the real axis are oscillatory and non-decaying, whereas contours parallel to the imaginary axis are non-oscillatory and exponentially decaying.
  • Because the imaginary-frequency integrand is smooth and exponentially decaying, exponentially convergent quadrature can achieve high accuracy using relatively few frequency points.
  • Imaginary-frequency solvers require validity across frequency bands rather than only a narrow real-frequency range.

7.5 Mean energy/force densities and the fluctuation–dissipation theorem

Casimir interactions can be viewed as geometry-dependent electromagnetic-field fluctuations whose average correlations are obtained from classical Green’s functions. Standard classical electromagnetic techniques can therefore compute Casimir quantities after modification for imaginary or complex frequencies.

  • Field fluctuations produce average energy densities and forces that can be integrated to obtain Casimir energies, forces, and related quantities.

7.5.1 Background

The review formulates Casimir energies and forces through Green’s functions, field correlations, and imaginary-frequency integrations. Imaginary frequencies improve numerical behavior, while zero frequency and coincident-point divergences require specific treatment.

  • Green’s functions: The electric and magnetic field correlations are expressed through classical Green’s functions, with the magnetic Green’s function obtainable from derivatives of the electric one.An alternative magnetic-dipole-current definition can be more convenient numerically, and a diagonal delta term does not affect energy differences or forces.
  • Energy density: Casimir energy is obtained by integrating the classical electromagnetic energy density over space and imaginary frequency.The imaginary-frequency formulation remains valid for arbitrary dissipative and dispersive media through its path-integral equivalence.
  • Regularization: The coincident-point Green’s function is formally infinite, but vacuum or self-field contributions can be subtracted analytically or regularized by discretization.Finite-difference and finite-element grids provide numerical regularization, while boundary-element methods explicitly separate vacuum and scattered fields.
  • Imaginary-frequency formulation: Wick rotation converts real-frequency fields into imaginary-frequency Green’s functions that are real, non-oscillatory, and exponentially decaying.At imaginary frequencies, the governing operator is real-symmetric positive-definite for ξ > 0, supporting effective numerical solvers.
  • Zero frequency: At nonzero temperature, the zero-frequency contribution is explicit and may require special low-frequency treatment, with controversy remaining for perfect dissipationless metals.At zero temperature, quadrature can avoid the measure-zero ξ = 0 point.
  • Force and torque: Casimir forces are computed by integrating the mean stress tensor over a surface enclosing the object, while torques use the moment arm crossed with the stress tensor.The stress-tensor field correlations are evaluated at imaginary frequencies.

7.5.2 Finite-difference frequency-domain (FDFD)

FDFD discretizes the Casimir Green’s-function problem on a spatial grid and solves the resulting sparse linear systems at imaginary frequencies. Its simplicity is offset by interface and corner singularities that can reduce convergence, while iterative and multigrid methods improve large-problem costs.

  • Discretization: FDFD divides space into a uniform grid, replaces derivatives with differences, and solves the discretized equations for each desired source.This is the finite-difference frequency-domain method widely used in classical electromagnetism.
  • Discretization: A Yee staggered grid places different electromagnetic components at offset positions to improve three-dimensional curl discretization.The one-dimensional discretization produces a tridiagonal system, while the full problem has multiple field components at each grid point.
  • Domain truncation: The computational domain must be truncated, with rapidly decaying Casimir interactions and perfectly matched layers providing practical boundary treatments.Coordinate transformations can represent PMLs as changes to ε and µ while preserving the governing PDE and discretization.
  • Accuracy: O(∆x^2) center-difference accuracy can fall to O(∆x) across discontinuous material interfaces unless specialized interface treatments or carefully smoothed materials are used.Appropriate interface modifications can restore second-order or higher accuracy.
  • Solution techniques: Sparse FDFD systems have real-symmetric positive-definite matrices at imaginary frequencies, enabling sparse-direct Cholesky factorization and iterative methods.The matrix has only O(N) nonzero entries, where N is the number of degrees of freedom.
  • Accuracy: Sharp geometric corners can limit convergence to O(∆x^p) for 1 < p < 2 even after material discontinuities are treated.The strongest uncompensated solution singularity determines the observed convergence rate.
  • Computational cost: An unpreconditioned iterative stress-integration calculation has O(N^2) time complexity, while an ideal multigrid preconditioner can reduce this to O(N^(2−1/d)).Further savings come from not summing over every point on the stress-integration surface.

7.5.3 Boundary-element methods (BEMs)

BEM discretizes surface currents on interfaces between homogeneous materials and solves the resulting linear system to obtain fields and Casimir observables. Imaginary-frequency evaluation improves matrix properties and enables efficient iterative solvers, while stress-tensor integration can achieve O(N log N) complexity under fixed-geometry refinement.

  • Surface formulation: BEM replaces volumetric unknowns with surface currents on interfaces between homogeneous materials, enforcing boundary conditions through surface integral equations.For perfect conductors, the unknowns can be physical electric surface currents; for permeable media, fictitious electric and magnetic surface currents provide the needed degrees of freedom.
  • Discretization: Triangular surface meshes use localized RWG basis functions associated with mesh edges, each supported only on its two adjacent triangles.Galerkin testing produces one equation per basis coefficient and reduces the field problem to a finite linear system.
  • Imaginary-frequency solution: The BEM matrix becomes real-symmetric and positive-definite at imaginary frequencies, while the Green’s function gains exponential spatial decay.These properties support more effective fast-matrix methods and linear solvers without changing the short-distance singularity requiring specialized integration treatment.
  • Force evaluation: A BEM solution can provide the mean stress tensor on an enclosing surface, whose surface integral yields the Casimir force.With dense-direct solvers, the same matrix factorization can be reused across evaluations; fast matrix-vector products and conjugate gradients can reduce single-point costs.
  • Force evaluation: O(N log N) overall complexity is possible when geometry is fixed and stress-tensor quadrature uses a number of points independent of N.For increasingly complex geometries or smaller feature sizes, the required number of stress-quadrature points may grow with N, and detailed scaling studies were unavailable.
  • Relation to other formulations: The BEM and path-integral approaches solve linear systems with the same matrix and degrees of freedom, while the path-integral formulation also computes interaction energy.The review describes comparable computational cost for dense solvers.

7.5.4 Other possibilities: FEM and spectral methods

FEM can apply classical electromagnetic techniques to the PDE formulation of Casimir Green’s functions, but nonuniform meshes require explicit vacuum regularization. Spectral integral-equation methods are another possibility, though their advantage over spectral path-integral methods is unclear.

  • Spectral methods: Spectral integral-equation methods, including multipole expansions for spheres and cylinders, can compute electromagnetic responses to currents.The review states that their advantages relative to spectral path-integral methods are unclear.
  • Quadrature: Stress-tensor integration uses multidimensional quadrature over a two-dimensional surface, with exponentially decreasing error possible for analytic integrands.Quadrature approximates integrals with weighted samples, and low-dimensional cubature schemes can be used for the surface integral.
  • Finite-element methods: FEM techniques can solve the PDE formulation of the Green’s function using a general mesh and are, in principle, directly applicable from classical electromagnetism.FEM divides space into elements and represents unknowns with localized basis functions, whereas BEM applies the analogous idea to surface unknowns.
  • Finite-element methods: Nonuniform FEM meshes require manual vacuum regularization because the divergent vacuum contribution varies with spatial resolution.A possible procedure is subtracting a calculation on the same mesh with the objects removed; these possibilities remained to be explored for Casimir physics.

7.5.5 Finite-difference time-domain (FDTD) methods

FDTD adapts classical time-domain simulation to Casimir calculations by converting the frequency integral into a time integral and using an artificial conductivity to control decay. The resulting method accumulates pulse responses with a precomputed, geometry-independent kernel.

  • Time-domain formulation: Casimir calculations are broad-bandwidth, so FDTD can replace separate frequency computations with a time-domain response to short pulses.The approach converts the frequency integral to a time integral and provides a time-domain counterpart to complex-frequency evaluation.
  • Time-domain formulation: Field correlations can be obtained from delta-time current responses, Fourier transformed, and combined with frequency-weighting factors to compute force or energy.Linearity makes the Fourier transform of the pulse response equivalent to the frequency-domain field-correlation calculation.
  • Time-domain formulation: The force can be computed directly as Γ(t)g(−t)dt, avoiding an explicit Fourier transform of the accumulated response.The kernel is precomputed, geometry-independent, and may depend on temperature when the thermal coth factor is included.
  • Complex-frequency transformation: Adding artificial conductivity σ everywhere produces dissipation that makes the time-domain response decay exponentially, allowing the simulation to be truncated after a short time.The conductivity arises from a contour choice that avoids exponentially growing negative-frequency solutions.
  • Implementation: The practical method uses an off-the-shelf real-time solver, artificial conductivity, short-pulse responses Γ(t), and a precomputed geometry-independent kernel g(t).FDTD is the most common classical electromagnetic time-domain implementation identified in the review.

7.5.6 Accelerating FD convergence

Finite-difference calculations can be accelerated through extrapolation, subtraction of isolated-object error forces, and reduced-order stress-tensor integration. These techniques target discretization error and the number of required linear solves rather than convergence itself.

  • Acceleration techniques: Richardson extrapolation fits results at multiple spatial resolutions to estimate the limit Δx → 0 and can improve O(Δx) convergence to O(Δx^2).The method increases convergence order cheaply when the resolution-dependent error rate is known.
  • Acceleration techniques: Subtracting the residual force computed for isolated objects reduces discretization-asymmetry error in the interaction-force calculation.Computing the A–B force and subtracting the forces for A alone and B alone triples the number of computations but can reduce the resolution needed for accuracy.
  • Acceleration techniques: Direct stress-tensor summation over every grid point is wasteful because each point requires a separate linear solve or time-domain simulation.The stress tensor is conceptually smooth in space, motivating methods that approximate its surface integral with fewer responses.
  • Acceleration techniques: Discrete cosine transforms expand fields on grid-aligned integration-surface faces so that only a small number of terms are needed for accurate stress integration.This reduces the number of right-hand sides that must be solved compared with evaluating the stress tensor independently at every point.

7.6 Path integrals and scattering matrices

Path-integral formulations express Casimir interactions through classical Green’s functions and scattering matrices after imposing boundary conditions with auxiliary interface fields. Their computational trade-offs depend on basis choice, geometry, solver scaling, and method-specific convergence or memory limits.

  • Path-integral formulation: Path integrals convert the Casimir interaction into a functional integral over vector potentials, with auxiliary interface fields enforcing material boundary conditions.For perfect conductors, the auxiliary fields can be interpreted as surface currents and the Green’s function is evaluated in the exterior medium.
  • Path-integral formulation: After choosing a finite current basis, the functional integral becomes a Gaussian integral whose evaluation produces a matrix log determinant.The resulting matrix is formed from Green’s-function matrix elements in the chosen basis.
  • Computational simplifications: Only the inter-object coupling block depends on relative object positions, allowing computations to be reused across multiple separations.For two objects, the self-coupling blocks remain unchanged as separation varies, while the cross-coupling block changes.
  • Basis and geometry: Spectral multipole bases converge exponentially when unaccounted singularities are absent, but convergence slows for corners, cusps, or near-touching surfaces.Spherical harmonics provide uniform angular resolution, whereas close surfaces require greater local resolution near the narrow gap.
  • Computational trade-offs: Dense-matrix path-integral calculations can avoid stress-integration surfaces and reduce system size through block decompositions, but their cost may grow as O(N^2 log N).The stated scaling assumes N linear systems solved with an O(N log N) fast solver and O(N) storage; low-rank trace or log-determinant methods may improve efficiency.
  • Computational trade-offs: Path-integral force expressions are exactly equivalent to stress-tensor surface integration in the continuum limit, linking the two computational formulations.The force expression becomes a surface integral over Green’s-function diagonals as the surface discretization is refined.
  • Basis and geometry: For corrugated periodic surfaces, eigenmode expansion computes reflection matrices by expanding fields in cross-sectional eigenmodes and matching boundary conditions where the cross-section changes.The method reduces the problem to scattering of incident planewaves from each surface.
  • Basis and geometry: RCWA handles uniform-cross-section regions using at most a two-dimensional cross-sectional discretization, but dense reflection matrices require O(N^2) storage and O(N^3) computation.Periodic structures with simple unit cells are favorable; nonperiodic structures can use PML boundaries at greater cross-sectional cost.

7.7 Periodicity and other symmetries

Periodicity and symmetry reduce Casimir calculations by decomposing the problem into smaller unit-cell or symmetry sectors. The resulting contributions are integrated or summed over the corresponding wavevectors or representation labels.

  • Translational periodicity: A structure periodic with period Λ can be reduced to one unit cell using Bloch-periodic boundary conditions.The fields at x = Λ equal those at x = 0 multiplied by the phase factor e^{ikxΛ}.
  • Translational periodicity: The Casimir energy or force is computed separately for each Bloch wavevector and then integrated over the Brillouin zone.For general periodic unit cells, the integration can be restricted to the irreducible Brillouin zone.
  • Continuous symmetry: Continuous translational symmetry corresponds to Λ → 0 and requires integration over all kx.Perfect-metal structures with continuous translational symmetry can further decompose into two polarizations.
  • Rotational and group symmetry: Rotational symmetry decomposes fields into angular sectors with e^{imφ} dependence, and the total interaction is the sum over integer m.More generally, Casimir contributions can be separated into irreducible representations of the structure’s symmetry group.

7.8 Nonzero-temperature corrections

At nonzero temperature, the zero-temperature imaginary-frequency integral becomes a Matsubara-frequency sum. The resulting corrections are usually small at room temperature for submicron separations, though geometry and materials can increase them.

  • Finite-temperature formulation: At T > 0, the zero-temperature imaginary-frequency integral is replaced by a sum over Matsubara frequencies.The n = 0 term carries weight 1/2, and ξn = 2πkBTn/ℏ.
  • Finite-temperature formulation: The finite-temperature sum is mathematically equivalent to a trapezoidal-rule approximation of the zero-temperature integral with spacing ∆ξ = 2π/λT.Its convergence is O(∆ξ^2), so the two results are close when the integrand varies slowly on that frequency scale.
  • Magnitude of corrections: At room temperature, λT ≈ 7.6 µm, making temperature corrections typically negligible for submicron Casimir separations.Careful choices of materials and geometry can nevertheless produce larger temperature effects.
  • Scope: Nonequilibrium configurations with objects at different temperatures are outside the scope of the review.The review discusses finite-temperature equilibrium corrections but does not cover these nonequilibrium situations.

7.9 Concluding remarks

Numerical Casimir computations remain open to many geometry and material combinations, while classical-electromagnetism parallels suggest continued method development. The review cautions against absolute rankings because methods have distinct strengths, weaknesses, and implementation costs.

  • Future directions: Many geometry and material combinations remain unexplored, leaving newly answerable questions about modifying Casimir phenomena through nanofabrication.The review identifies the computational area as rich with opportunities.
  • Future directions: Parallels with computational electromagnetism suggest ongoing improvements and developments in numerical Casimir methods.The review connects future progress to established classical-EM techniques.
  • Method selection: Different numerical techniques should not receive absolute rankings because they exhibit unique theoretical and practical strengths and weaknesses.The review notes that implementation simplicity and general availability can outweigh theoretical suboptimality in practice.
Loading 1007.0966v2…