Casimir calculations in Meep

From AbInitio

(Difference between revisions)
Jump to: navigation, search
Revision as of 16:34, 23 September 2009 (edit)
Alexrod7 (Talk | contribs)
(Example: dispersive materials)
← Previous diff
Revision as of 16:35, 23 September 2009 (edit)
Alexrod7 (Talk | contribs)
(Parallelization)
Next diff →
Line 203: Line 203:
==Parallelization== ==Parallelization==
-[coming soon]+ 
 +If you look at the sample file rods-plates.ctl
==Example: dispersive materials== ==Example: dispersive materials==

Revision as of 16:35, 23 September 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:

Example structure: A two-dimensional piston-like configuration in which we will compute the Casimir force between two blocks between two parallel plates.
Enlarge
Example structure: A two-dimensional piston-like configuration in which we will compute the Casimir force between two blocks between two parallel plates.

This is described in the sample file rods-plates.ctl ([1]). The dashed red lines indicate the surface S. This system consists of two metal a\times a squares in between metallic sidewalls. To run a simulation in which the blocks are (nondispersive) dielectrics one can simply change their materials in the definitions as in a normal Meep simulation. For dispersive dielectrics a few extra steps are needed, which is discussed in a later section.

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, we can use a very short runtime 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 user does not have to explicitly construct the f_n(\textbf{x}) source basis for Casimir calculations. This is all done inside of Meep (see below). The user only has to specify how many harmonic moments n to use for the simulation. Here we briefly describe what is done inside of Meep.

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:

Schematic of Casimir calculation via a sequence of time-domain sources around the integration surface.
Enlarge
Schematic of Casimir calculation via a sequence of time-domain sources around the integration surface.

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)

Snapshots of computed fields from sources around the stress surface.
Enlarge
Snapshots of computed fields from sources around the stress surface.

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) (list 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.

Important Numerical Considerations

Before going on to discuss more complicated simulations, we pause a moment to discuss some important practical considerations for running a simulation.

Vacuum Subtraction

In the 2d blocks example discussed above, as long as the source surface is square an close to the right block, the force converges rapidly with resolution. However, if the user varies the source surface (for example, using a rectangular source surface or one that is not placed symmetrically around the right block), they may get a drastically different result. This is entirely a numerical artifact, as the correct Casimir force should be independent of the source surface used (as long as that surface does not intersect any objects).

This numerical artifact can be removed by a simple procedure, called vacuum subtraction. This procedure consists of running two additional simulations with certain objects deleted and subtracting the results, illustrated below:

The procedure for the double blocks case is illustrated below: Image:Vac.jpg

In an ideal case of infinite spatial resolution, only the first term is nonzero. However, due to discretization effects, for finite resolution they can be quite large and will depend on the source surface. The total vacuum-subtracted force, however, will still be well-behaved for finite resolution.

For example, in the blocks example previously discussed, if one were to shift the center of the soruce surface so that it doesn't align with the right block (but still doesn't overlap with any other objects), the force will get contributions from both vacuum subtraction terms.

The vacuum terms should go to zero as resolution is increased, however they may do so quite slowly. Also it is important to remember that only the total force (ie summed over all values of n) that should converge to zero, so several values of n must be retained in order to verify this.

Although this procedure adds to the computational cost of the problem, it is very easy to implement (it is demonstrated in the rods-plates.ctl example file) and can oftentimes be essential to get a meaningful result.

Simulation Time

As discussed in Part I, the time-convergence of a simulation depends on the size of the system. This can generally be optimized by varying the value of the conductivity σ. The convergence rate also depends on the dielectric function of the medium. This makes sense because a high-dielectric medium has a longer optical thickness than vacuum. The optimal value of σ depends on the system, but generally decreases if the dielectric of the medium is increased (for example, for simulations where the dielectric of the medium is ethanol, the optimal value is somewhere around σ˜0.2).

Running a long simulation to determine temporal convergence is not a large numerical task, because the temporal convergence is largely independent of the spatial resolution and harmonic moment n, so one can simply run a long low resolution job.

Boundary Effects

The Casimir force between metals can be very strong, while the force between dielectric bodies may be much weaker. It is then important to use the appropriate boundary conditions. For example, the default boundary conditions in Meep are metallic - if the boundaries are close enough, the Casimir force between the bodies and the metal walls will overwhelm their mutual their mutual forces. Usually to avoid this problem (unless the problem specifically requires the objects to be in a metal box), specifying periodic boundary conditions by calling (set! k-point (vector3 0 0 0) works the best, even in the presence of PML.

Parallelization

If you look at the sample file rods-plates.ctl

Example: dispersive materials

As discussed in Part I, the treatment of dispersive materials is very similar to the normal case of metals or non-dispersive dielectrics. Here we show how to add them in to the computation. We treat only the case of lossless dielectrics; the case of loss is more complicated, and the loss terms tend not to greatly affect the Casimir force.

A lossless dielectric material with a single Lorentzian resonance, e.g., Silicon, is defined by a resonant frequency ω0 (this is the angular frequency in radians), an oscillator strength C and a high frequency dielectric εf:

\epsilon(\xi) = \epsilon_f + \frac{C \xi_0^2}{\omega_0^2 - \xi^2}

(here we follow the notations of Parts I and II, rather than the rest of the wiki, in which ξ denotes real frequency in FDTD, and ω denotes the more abstract complex frequency mapping).

The conductivity mapping ω2 = ξ2 + iσξ must be applied to all dispersions. Applying this to ε(ξ) above gives the new dispersion function:

\epsilon (\xi) = \epsilon_f + \frac{C \xi_0^2}{\omega_0^2 - \xi^2 - i\sigma \xi}

So the new dispersion is a Lorentzian, but with an additional loss term. This is the correct material to define in Meep.

It is easy to define a dispersive material in Meep (discussed further in Materials in Meep, with examples in Meep Tutorial/Material dispersion). Here is how we go about it (further material examples are in materials.scm ([2]), and rods-plates.ctl ([3]) demonstrates their use.

 (define length-scale 1e-6) ;length scale - units of 1 micron
 (define w0 (* wconv 6.6e15)) ;convert angular frequency (in radians / second)
 (define eps0 11.87) ;DC dielectric
 (define epsf 1.035) ;high frequency dielectric
 (define Silicon (make medium (epsilon epsf)
                   (E-polarizations
                     (make polarizability
                       (omega (/ w0 (* 2 pi))) (sigma (- eps0 epsf)) 
                       (gamma (/ Sigma (* 2 pi)))))))

There are two important things to note about the definitions above. First, "sigma", the oscillator strength variable in the polarizability class, is different from "Sigma", which is our global conductivity (of course "Sigma" is not built in, so it can be renamed). Also, both w0 and Sigma are divided by a factor of . This is because both the resonant frequency and dissipation passed into Meep must be specified in units of 2pic / a, whereas both w0 and Sigma are given in units of c / a.

Now the Silicon that we have defined can be use in the Casimir calculations like any other material type.

Example: cylindrical symmetry

[coming soon]

Example: three-dimensional computations

[coming soon]

Personal tools