Casimir calculations in Meep

From AbInitio

(Difference between revisions)
Jump to: navigation, search
Revision as of 23:05, 31 August 2009 (edit)
Stevenj (Talk | contribs)
(Introduction)
← Previous diff
Revision as of 23:07, 31 August 2009 (edit)
Stevenj (Talk | contribs)
(Introduction)
Next diff →
Line 28: Line 28:
The classical force is then: The classical force is then:
-<math> F_i^{Classical} = \int_S \int_0^\infty M_{ij}(\omega,x) d\omega dS_j </math> +:<math> F_i^{Classical} = \int_S \int_0^\infty M_{ij}(\omega,x) d\omega dS_j </math>
where the Maxwell stress tensor is given by: where the Maxwell stress tensor is given by:
-<math> M_{ij}(\omega,x) \equiv E_i(\omega,x)E_j(\omega,x) + H_i(\omega,x)H_j(\omega,x) - \frac{1}{2}\delta_{ij}\sum_{k=1}^3 \left(E_k^2(\omega,x) + H_k^2(\omega,x) \right) </math>+:<math> M_{ij}(\omega,x) \equiv E_i(\omega,x)E_j(\omega,x) + H_i(\omega,x)H_j(\omega,x) - \frac{1}{2}\delta_{ij}\sum_{k=1}^3 \left(E_k^2(\omega,x) + H_k^2(\omega,x) \right) </math>
-[ref] showed the the Casimir force can be obtained by a similar integral, except that the stress tensor is replaced with the expectation value of the corresponding operator in the quantum-mechanical grounds state of the system:+As has been known for several decades (reviewed [http://math.mit.edu/~stevenj/papers/RodriguezIb07-pra.pdf here]), the the Casimir force can be obtained by a similar integral, except that the stress tensor is replaced with the expectation value of the corresponding operator in the quantum-mechanical ground state of the system:
-<math> F_i^{Casimir} = \int_S \int_0^\infty < M_{ij}(\omega,x) > d\omega dS_j </math>+:<math> F_i^{Casimir} = \int_S \int_0^\infty < M_{ij}(\omega,x) > d\omega dS_j </math>
-In Part I, it is shown how this formulation can be turned into an algorithm for computing Casimir forces in the time-domain. The reader is referred there for details. Below we present the basic steps used in our implementation.+In Part I, it is shown how this formulation can be turned into an algorithm for computing Casimir forces in the time domain. The reader is referred there for details. Below we present the basic steps used in our implementation.
==Implementation== ==Implementation==

Revision as of 23:07, 31 August 2009

Meep
Download
Release notes
FAQ
Meep manual
Introduction
Installation
Tutorial
Reference
C++ Tutorial
C++ Reference
Acknowledgements
License and Copyright

Contents

Computing Casimir forces with Meep

It is possible to use the Meep time-domain simulation code in order to calculate Casimir forces (and related quantities), a quantum-mechanical force that can arise even between neutral bodies due to quantum vacuum fluctuations in the electromagnetic field, or equivalently as a result of geometry dependence in the quantum vacuum energy.

Calculating Casimir forces in a classical finite-difference time-domain (FDTD) Maxwell simulation like Meep is possible because of a new algorithm described in:

These papers describe how any time-domain code may be used to efficiently compute Casimir forces without modification of the internal code. The newest version of Meep contains several optimizations of these algorithms, allowing for very rapid computation of Casimir forces (reasonably-sized two-dimensional systems can be solved in a matter of seconds).

This page will provide some tutorial examples showing how these calculations are performed for simple geometries. For a derivation of these methods, the reader is referred to the papers above, which will be referred to as Part I and Part II in this webpage.

Introduction

In this section, we introduce the equations and basic considerations involved in computing the force using the method presented in Rodriguez et. al. (arXiv:0904.0267). Note that we keep the details of the derivation to a minimum and instead focus on the calculational aspects of the resulting algorithm.

The general setup for a Casimir force computation is shown in the following figure:

Schematic of Casimir-force calculation on an object, by integrating the Maxwell stress tensor on a surface S around the object.
Enlarge
Schematic of Casimir-force calculation on an object, by integrating the Maxwell stress tensor on a surface S around the object.

The goal is to determine the Casimir force on one object (shown in red) due to the presence of other objects (blue).

Classically, the force on the red object due to the electromagnetic field can be computed by integrating the Maxwell stress tensor Mij (ref: Griffiths) over frequency and over any surface enclosing only that object as shown above.

The classical force is then:

F_i^{Classical} = \int_S \int_0^\infty M_{ij}(\omega,x) d\omega dS_j

where the Maxwell stress tensor is given by:

M_{ij}(\omega,x) \equiv E_i(\omega,x)E_j(\omega,x) + H_i(\omega,x)H_j(\omega,x) - \frac{1}{2}\delta_{ij}\sum_{k=1}^3 \left(E_k^2(\omega,x) + H_k^2(\omega,x) \right)

As has been known for several decades (reviewed here), the the Casimir force can be obtained by a similar integral, except that the stress tensor is replaced with the expectation value of the corresponding operator in the quantum-mechanical ground state of the system:

F_i^{Casimir} = \int_S \int_0^\infty < M_{ij}(\omega,x) > d\omega dS_j

In Part I, it is shown how this formulation can be turned into an algorithm for computing Casimir forces in the time domain. The reader is referred there for details. Below we present the basic steps used in our implementation.

Implementation

The basic steps involved in computing the Casimir force on an object are:

1. Surround the object for which the force is to be computed with a simple, closed surface S. It is often convenient to make S a rectangle in two dimensions or a rectangular prism in three dimensions.

2. Add a uniform, frequency-independent conductivity σ to the dielectric response of every object (which is easily done Meep). The purpose of this is to rapidly reduces the time required for the simulations below. As discussed in Part I, adding this conductivity in the right way leaves the result for the force unchanged.

3. Determine the Green's function along S. This is done by measureing the electric E and magnetic H fields on S in response to a set of different current distributions on S (more on the form of these currents later).

4. Integrate these fields over the enclosing surface S at each time step, and then integrate this result, multiplied by a known function g( − t), over time t.

The Casimir force is given by an expression of the form:

F_i = \sum_n \mathrm{Im} \int_0^\infty dt \, g(t) \int_S dS_j(\textbf{x}) f_n(\textbf{x}) \Gamma_{ij;n}(t,\textbf{x})

where the Γij;n are the fields in response to sources related to the functions fn (discussed in detail later), S is an arbitrary closed surface enclosing the object for which we want to compute the Casimir force, g(t) is a known function, and the index n ranges over all of the integers.

The functions \Gamma_{ij;n}(t,\textbf{x}) are related to the Maxwell stress tensor introduced in the previous section. Here the frequency integration has been turned into an integration over time, which is especially suited for our purposes.

Note that the precise implementation of step (3) will greatly affect the efficiency of the method. For example, computing the fields due to each source at each point on the surface separately requires a separate Meep calculation for each source (and polarization). This corresponds to taking f_n(\textbf{x}) = \delta(\textbf{x}-\textbf{x}_n) for each point \textbf{x}_n \in S, and the sum over n becomes an integration over S. As described in Part II (arXiv:0906.5170), we are free to take fn to be any basis of orthogonal functions. Picking an extended basis (e.g. cosine functions or complex exponentials eikx) greatly reduces the number of simulations required.

Example: two-dimensional blocks

In this section we calculate the Casimir force in the two-dimensional Casimir piston configuration (Rodriguez et. al, Physical Review Letters, vol 99, p 080401 2007) shown below:

This is described in the sample file rods-plates.ctl included in the examples directory. The dashed red lines indicate the surface S. This system consists of two metal a\times a squares in between metallic sidewalls.

First define the geometry, consisting of the two metal sidewalls (each 2 pixels thick) and the two blocks:

 (set-param! resolution 40) 
 (define a 1)
 (define-param h 0.5)
 (set! geometry
    (list (make block (center 0 (+ (/ a 2) h (/ resolution))) ;upper sidewall 
                                  (size infinity (/ 2 resolution) infinity) (material metal))
            (make block (center 0 (- (+ (/ a 2) h (/ resolution)))) ;lower sidewall
                                  (size infinity (/ 2 resolution) infinity) (material metal))
            (make block (center a 0) (size a a infinity) (material metal)) ;right block
            (make block (center (- a) 0) (size a a infinity) (material metal)))) ;left block

Define an air buffer on either side of the blocks and the PML thickness. Then set the computational cell size, and add in pml on the left/right sides:

 (define buffer 1)
 (define dpml 1)
 (set! geometry-lattice (make lattice (size (+ dpml buffer a a a buffer dpml) (+ (/ 2 resolution) h a h (/ 2 resolution)) no-size))) 
 (set! pml-layers (list (make pml (thickness dpml) (direction X))))

Define the source surface S; here we take it to be a square with edges 0.05 away from the right block on each side:

 (define S (volume (center a 0) (size (+ a 0.1) (+ a 0.1))))

As described in Part II, we add a uniform, frequency-independent D-conductivity σ everywhere:

 (define Sigma 1)

(note that "sigma" is another built-in Meep parameter, so here we use a capital S).

As discussed in Part I, the optimal value of σ depends on the system under consideration. In our case, σ = 1 is optimal or nearly optimal. With this choice value of Sigma<\code>, we can use a very short runtime <code>T for the simulation:

 (define-param T 20)  

The only thing left to define is the function g( − t). This is done with the Meep function (make-casimir-g T dt Sigma ft):

 (define gt (make-casimir-g T (/ Courant resolution) Sigma Ex))

Here we can pass in either field type Ex or Hx. Since the E/H fields are defined for integer/half-integer units of the time step dt, we technically require a different g(t) for both polarizations, and this option is allowed. However, for this example we get sufficient accuracy by using the same function for all polarizations.

Sources

The built-in source basis used in Meep for this type of computation consists of a Cosine basis. For each side of S and each non-negative integer n, this defines a source distribution:

f_n(x) = \sqrt{\frac{c_n}{L}} \cos \left(\frac{n\pi x}{L}\right), ~n = 0,1,\ldots

where cn = 1 if n = 0 and cn = 2 otherwise, L is the side length (if each side has a different length, then the functions fn(x) will differ for each side). An illustration of these functions for the system under consideration, compared to point sources, is shown below:

For the simulation, we must truncate the sum over n to some finite upper limit n-max. Typically, a value of 10 is sufficient to get results to within 1%:

 (define-param n-max 10) 

To illustrate the field profiles, below we show four snapshots at different times for what we term \Gamma^E_{yy;n=2}(\textbf{x},t), the y-component of the electric field response to a y-polarized current source with spatial dependence f2(x)

As the fields continue to propagate, the finite conductivity Sigma causes the fields to decay rapidly. By T\simeq 20 the fields have all but vanished. Note that the user is not confined to use this Cosine basis, but that this type is already built into meep and therefore offers the greatest convenience.

Running the Simulation

Computing the Casimir force involves running several independent Meep simulations for each set of parameters. The fact that they are independent makes the problem very easy to run in parallel, which is why we implemented the calculation as a series of independent Meep simulations. These parameters are grouped into the following lists:

 (define pol-list (list Ex Ey Ez Hx Hy Hz))          ;source polarizations
 (define side-list (list 0 1 2 3))                            ;each side of the square surface S
 (define n-list (if (eq? n-max 0) (list 0) (interpolate (- n-max 1) 0 n-max)))    ;number of terms in source sum

For each value of n, side number, and polarization, we run a short meep simulation. For convenience, the source construction, simulation, and field integration are all taken care of by the Scheme function casimir-force-contrib, defined in "/libctl/casimir.scm" (included in the Meep-1.1 release):

 (casimir-force-contrib force-direction integration-vol N Sigma T source-component gt)

Here the desired component of the Casimir force is the X direction (all others are zero by symmetry), and the source volume is S. The integer N is defined to be N\equiv 4\times n + \mathrm{side}.

To find the casimir force, one simply iterates over the parameter lists constructed above, each time adding the result of casimir-force-contrib to the total force:

 (define total-force 0)
 (do ((n 0 (1+ n))) ((= n (length n-list))) ;source components
     (do ((side 0 (1+ side))) ((= side (length side-list))) ;sides of S
        (do ((p 0 (1+ p))) ((= p (length pol-list))) ;field polarizations
            (let* ((curr-pol (list-ref pol-list p))
                      (N (+ (* 4 n) side)))
                    (set! total-force
                       (+ total-force
                            (casimir-force-contrib X S N Sigma T curr-pol gt)))))))

Results

The result, when sampled over many values of h, is a force curve that varies non monotonically in h:

Here the force is measured relative to the value obtained by the PFA (proximity-force approximation), a simple approximation for the force used in lieu of an exact analytic expression or a numerically accurate value. As is clear from the figure, the behavior of the force deviates significantly from the PFA, indicating the need for accurate algorithms for computing the Casimir force.

Using the built-in functions above, it takes roughly 20 seconds to run simulation at resolution 40 for each value of n, including all sides of S and all field polarizations. Typically, only n\leq 10 are needed to get the force to high precision, so the Casimir force for this system (for a single value of h) can be determined in under two minutes on a single processor.

Parallelization

[coming soon]

Pitfalls to avoid

[coming soon]

Example: dispersive materials

[coming soon]

Example: cylindrical symmetry

[coming soon]

Example: three-dimensional computations

[coming soon]

Personal tools