Numerical studies of Casimir interactions

0.13 MB·English

by  S. Pasquali

Numerical studies ofCasimirinteractions
S. Pasquali, A. C. Maggs
LaboratoiredePhysico-ChimeThe´orique,Gulliver,CNRS-ESPCI,10rueVauquelin,75231ParisCedex05,France.
We study numerically the Casimir interaction between dielectrics in both two and three dimensions. We
demonstrate how sparse matrix factorizations enable one to study torsional interactions in three dimensions.
Intwo dimensions westudy the full cross-over between non-retarded and retarded interactions asafunction
ofseparation. Weuseconstrainedfactorizationsinordertomeasuretheinteractionofaparticlewitharough
dielectricsurfaceandcomparewithascalingargument.
Dispersionforceshaveastheiroriginthefluctuationsofpo- InthreedimensionswediscretizeMaxwell’sequationstoa
8
larization in materials coupled to the long-rangedelectrody- cubicYee-lattice[14],latticeconstanta = 1,associatingthe
0
namicinteractiondescribedbyMaxwell’sequations.Thefirst electricdegreesoffreedomto thelinks, magneticdegreesof
0
2 convincingdemonstrationof theirimportancewasthe calcu- freedomare localized onthe facesof the lattice. We remind
lationbyKeesomoftheinteractionbetweenfluctuatingclas- thereaderthatthefinitedifferenceapproximationtothe∇×
n
sicaldipoles[1]. Theintroductionofquantumfluctuationsby operator,heredesignatedCurl,mapstheelectricfieldonfour
a
J London [2] accounted for the long-ranged, 1/r6, part of the links surroundingthe face of the cube to the magnetic field.

8 van der Waals interaction in most materials. Later, Casimir CurlisneededintheMaxwellequation
2 and Polder [3] showed that retardation modifies the interac-
tions in an important manner– leading to a decay in the in- ∂H =−cCurlE
h] teraction which is asymptotically 1/r7 at zero temperature. ∂t
p Further advances were made by the Russian school [4] who
The adjoint operatormaps fields from the faces to the links.
- showed how to formulatethe interactionsin terms of the di- ∗
t WewilldenoteitCurl . Itintervenesinthesecondtimede-
n electric response of materials. Overviews with many refer-
pendentMaxwellequation.
a
encestotheoreticalandexperimentaldevelopmentsaretobe
u
q foundin[5,6,7]. RetardedCasimirinteractionsarethedom- ∂D =cCurl∗H
[ inant interaction between neutral surfaces at the submicron ∂t

scale.
1 The importance of clearly distinguishing the two operators
v Whilst the analytic basis of the theory is largely estab-
will become apparentwhen we discuss the two-dimensional
5 lished its application is difficult in experimentally interest-
case below. We use Heaviside-Lorentz units in which
8 ing geometries. One is constrained to work with perturba-
3 Maxwell’sequationsare directlyparameterizedbythespeed
tive expansion about exactly solvable geometries [8], or use
4 oflightinvacuum,c.
ad hoc schemes such as the proximity force approximation.
. From these two equations Lifshitz theory [15] shows that
1
Onlyafewgeometrieshavebeenattackedwithexactanalytic
0 the free energy of interaction between dielectric bodies is
techniques [9]. Recently several attempts have been made
8 found from from the imaginary time wave equation for the
0 to study numericallythe interactionsby using methodsfrom vectorpotentialinthetemporalgaugewhereE=−A˙/cand
: modern computational science– including fast multigrid lat-
v φ=0
tice solvers [10] in order to calculate Green functions and
i
arX feonrecregsi,eosr[1th1e].useofdiscretizeddeterminantstodeterminefree (cid:26)ǫ(r~,2ωc)2ω2 +Curl∗Curl(cid:27)A=DAA=0
Inthis Letterwe will presenta series oftechniqueswhich
enableonetoevaluatetheinteractionbetweendielectricbod- Alternatively one introduces a magnetic formulation and
iesinfullvectorialelectrodynamics. Firstly,wecalculatethe workswithapotentialsuchthatH = G˙ /candconsidersthe
torsional potential between two three-dimensional bodies in waveequation
the retarded regime, using a full discretization of Maxwell’s
equations, we note that the Casimir torque has recently re- ω2 +Curl 1 Curl∗ G=D G=0
(cid:26)~2c2 ǫ(r,ω) (cid:27) G
ceivedtheattentionofexperimentalists[12,13].Formorede-
tailedstudieswepresentresultsfortwo-dimensionalsystems.
Inourworkwealwaysconsiderthedifferencesinfreeen-
Thisallowsusto studythecross-overbetweenthe near-and
ergybetweenpairsofconfigurations;wethusavoidafullac-
far-fieldregimesandalsotomeasuretheinteractionbetween
count of the self-energy variations of dielectric media [11].
a particle and a rough surface. With these two-dimensional
Thefreeenergydifferencebetweentwoconfigurations1,2is
systemsweimplementgeneralstrategieswhichsubstantially
foundfrom
increase the efficiency of simulations, at the same time de-
∞
creasing the sensitivity of the results to numerical round-off dω
U1,2 = lndet D1(ω)−lndet D2(ω) (1)
errors. Z0 2π
(cid:8) (cid:9)
2
1.2
1
0.8
Max 0.6
U
U/ 0.4
FIG. 1: A Pair of structured dielectric rings. Each quadrant has 0.2
differentdielectricproperties.
0
−0.2
0 1/8 1/4 3/8 1/2 5/8 3/4 7/8 1
for either choice of wave operator, DA or DG; while self- angle/π
energycontributionsaredifferentinthetwoformulationswe
haveverifiedwithourcodesthatbothgivethesameresultfor
thelong-rangedpartoftheinteractionsthatweareinterested
FIG.2: Interactionenergyasafunctionofangleasaringoffigure1
in.
rotates. Inthelinearpartsofthecurvethetorqueisalmostindepen-
We performthefrequencyintegrationineq.(1) bychang-
dentoftheangle. Ringdiameter36a, separationandthickness2a.
ing integration variables to z, where ω = αz/(1−z) with Roundingisdeterminedbytheratioofseparationtodiameterofthe
0 < z < 1. The parameter α is chosen so that the ma- rings.10daysofcomputationwithNg =8.Choleskyfactor9GB.
jor features in the integrand occur for values of z near 1/2.
We thenuse N -pointLegendre-Gaussquadratureto replace
g
the integral by a weighted sum over discrete frequencies.
angleisnoticeablytriangularinshapebetweenπ/8to3π/8.
We evaluatedeterminantsby findingthe Choleskyfactoriza-
This is understood by the fact that the interaction energy is
tion L of D(ω) such that L is lower triangular [16] and
D D dominated by the interactions directly across the gap. The
L LT =D(ω). ThedeterminantofDisthengivenby
D D fluctuations in the curve about the expected linear behavior,
together with its slight asymmetry give an idea of the noise
lndetD(ω)=2 ln(L )
D,i,i comingfromirregularitiesoftheinterpolationofthedisksto
Xi
thelattice. Thisirregularityisparticularlyclearoncomparing
When we examine the detailed structure of Maxwell’s thepointsforπ/4and3π/4.
equations discretized to V = L3 sites in three dimensions
Wenowturntotwo-dimensionalelectrodynamicswherewe
we discover that the Curl operator is a matrix of dimension
can study systems with larger linear dimensions. Such large
3V × 3V and has 12V non-zero elements. The operator
systemsizesareneededinordertofollowthecross-oversbe-
∗
(Curl Curl) has 39V non-zero elements. The major tech-
tweendifferentregimesintheinteractionofparticlesorifone
nicaldifficultycomesfromthefactthatthematriceswework
wishestosimulatestructuredordisorderedmaterialsinorder
with have dimensions which are very large, ∼ 106 × 106.
understandtheefficiencyofanalyticapproximations.
All numericalwork was performedwith an IntelXeon-5140
In three dimensions the two formulation in terms of D
workstation. A
andD arelargelyequivalent.Intwo-dimensionalelectrody-
WenowcalculatetheCasimirtorquebetweentwoparallel G
namicsthisisnolongerthecase. Consideranelectrodynamic
ringscenteredon a commonaxis, figure 1. Each ring is di-
systeminwhichtherearetwocomponentsoftheelectricfield
videdintoquadrantswithalternatingdielectricproperties.We
in the x−y plane; the magnetic field then has just a single
takepermittivitieswhichareindependentoffrequency,corre-
componentin the z direction. The Curl operatorbecomesa
sponding to the full retarded regime [15] with ǫ1(ω) = 5,
rectangularmatrixofdimensionsV×2V wherenowV =L2.
ǫ2(ω) = 10; the space around the rings is a vacuum with
Thestandardformulationintermsofthevectorpotentialgives
ǫ = 1. Wemeasuretheenergyofinteractionasthetopring
r
toanoperatorD ofdimensions2V ×2V with14V non-zero
is rotated with respect to the lower. The zero of the inter- A
elements;thealternativeformulationintermsofD leadsto
action correspondsto aligned rings. As the ringsare rotated G
determinants of dimensions V × V involving just 5V non-
theinterfacebetweenthedielectricmaterials, asinterpolated
zeroelements;thesizeofthematrixthatwemustworkwith
to the lattice, undergoes some re-arrangement changing the
issmallerintheD formulation. WeusedD inthefollow-
selfenergyoftherings. Wethusperformtworuns. Thefirst G G
ingnumericalwork,havingcheckedthatweobtainequivalent
run of a single rotating ring determines this variation in the
results.
self-energy. Thesecondrunwiththebothringsallowsoneto
measuretheinteractionenergybysubtraction. Westartedbymeasuringthecross-overbetweentheshort-
WeworkedwithasystemofdimensionsV =55×55×55, ranged non-retarded interaction to the long-ranged Casimir
figure 2. Thegraphoftheinteractionenergyasafunctionof force. We studied a pair of dielectric particles described by
3
thesinglepoleapproximationtothedielectricconstant
100
χ
ǫ(ω)=1+
1+ω2/ω2~2
0
whereχ isthezerofrequencyelectricsusceptibility. Thein- 5U r
teractionisretardedforseparationsD ≫ c/ω0,non-retarded −10−1
forD ≪c/ω0.
We measured the interaction between two dielectric par-
ticles in a square, periodic cell of dimensions L × L using
SuiteSparse[17]toperformboththeorderingandthefactor-
10−2
ization of the matrices. We placed a first particle at the ori- 101 102 103
gin, andconsideredtwopossiblepositionsofasecondparti- r
cletocalculateafreeenergydifferenceusingeq.(1).Thefirst
FIG.3: Scaledinteractionfreeenergy,−Ur5forapairofdielectric
results were disappointing– rather small systems (L = 50)
particles(ǫ(0)=8)inaboxofdimensions2000×2000asafunction
were sensitive to numerical round-off errors. The origin of
of separation. Curves from top to bottom correspond to ω0/c =
thisproblemwasquiteclear. Inalargesystemthereisanex- 10, 0.3, 0.1, 0.03, 0.010.003. Forlargeω0/c, Ur5 isconstant, (cid:3).
tensive self-energy∼ L2. Pair interactionscalculated as the Forsmallerω0/cweseebothretardedandnon-retardedinteractions.
differencebetweentwolargenumbersareunreliable. SolidlinecorrespondstoU ∼1/r4. 10GBforCholeskyfactor.
vdw
Weavoidedthisproblembyseparatingthefreeenergycon- Sixhoursofcalculation.Ng =25.
tributionsfromtheneighborhoodofthethreeinterestingsites
and the rest of the system. We did this by introducing a
block-wise factorization of D that enabled us to both solve decayingasUvdw =1/r4. Asinthreedimensionsretardation
the round-offproblemwhile re-usingmuchofthe numerical leads to an accelerated decay so that the Casimir interaction
effort need to generate the Cholesky factors thus improving varies as Uc ∼ 1/r5. In our simulations we used values of
theefficiencyofthecode. ω0/cvaryingfrom0.003to10, figure3. We determinedthe
We now write the symmetric matrix from the wave equa- energyof interaction of particles U, as a functionof separa-
X Y tionrwhilemovingthesecondparticleinthesimulationcell
tion in block form, D = . Its determinant
(cid:18)YT Z(cid:19) out to (L/5,L/5); the zero of energy is calculated for two
is det(D) = det(X)det(S) where the Schur complement particlesseparatedby(L/2,L/2). We scaleouttheretarded
S = Z −YTX−1Y [18]. We group sites so that the great behavior,plotting−U(r)r5. Weseethatforthelargestω0/c
majorityiswithintheblockX andsitesthatweareinterested the interactions are retarded for all separations, (cid:3). For the
in are in the block Z. It is the term in det(X) the gives the smaller values of ω0/c the interactionvaries as 1/r4. In the
largeextensivefreeenergywhichcausedournumericalprob- scaledcurvethisgivesthelinearriseclearlyvisibleinthefig-
lems. Itisindependentofthepropertiesofourtestparticles. ure, ⋄. For 0.1 < ω0/c < 0.01 we see both the near- and
Alltheinterestinginformationonenergydifferencesisinthe far-field behaviorsclearly displayed within a single sample–
Schurcomplement,S. permitting the detailed study of cross-over phenomena with
We start by finding the Cholesky factorization of X, L . frequencydependentdielectric behavior. No assumptionsof
x
TheSchurcomplementiscalculatedbysolvingthetriangular symmetryaremadeinthecalculation;themethodcanbeused
equationsL U =Y byforwardsubstitution,thencalculating withbodiesofarbitrarygeometry.
x
S = Z −UTU. Ourseparationofenergiesintoanextensive Wenowturntoaproblemwhereanalyticresultsaremuch
constantandasmallsetofinteractingsitesallowsustostudy more difficultto find: The interaction of a dielectric particle
the interaction of systems of sizes up to L = 2000 before witharoughsurface,figure4. Wegeneratedroughsurfacesas
round-offbecomesaproblem. realizationsofsolid-on-solidrandomwalksona lattice. Ap-
In order to generate data we generalized the method to a proximatelyhalfofthesimulationboxcontainsdielectricma-
threelevelscheme–firstlycollectthesetofsites(here∼100) terial with ǫ = 8,ω0 = ∞; the rest of the box has ǫ = 1.
of all the separations required to generate a curve into the We measure the interaction with a test particle as a function
blockZ,andformtheSchurcomplementformingasmallef- ofthedistancefromtheroughsurfaceusingtheabovemethod
fectivetheoryforalltheseremainingsites. Withinthesmaller ofblockSchurcomplementstoperformasinglelargefactor-
matrix that has been generated we again re-order to succes- izationperfrequencyforeachrealizationofthedisorder. We
sivelyputeachinterestingsetsofvariablesinthebottom-right generated1000roughsurfacesandmeasuredtheaveragein-
corneroftheeffectivetheoryandfindtheSchurcomplement teractionwiththesurfacehUi,asafunctionofseparation,as
of these remainingvariables. We can then calculate interac- wellasthevarianceinthepotentials.
tionsbetweentheparticleswhileminimizinground-offerrors. We understand the results, figure 5, with a scaling argu-
We remind the reader that in two dimensions the electro- ment. When the particle is a distance r from the surface the
static potential is logarithmicbetween two charges, and that interactionisdominatedbyafrontoflengthr alongthesur-
dipole-dipole fluctuations lead to van der Waals interactions face. Since the surface is a random walk its average posi-
4
600 blewith[19].
We have demonstrated the power of direct methods from
550
linearalgebrawhenappliedtothestudyofdispersionforces.
Inthreedimensionswehavemeasuredinteractionsinexperi-
500 mentallyrealizablegeometries–thoughsystem sizesare still
too small to accurately measure cross-overs between differ-
entscaling regimes. In two dimensionswe have shownhow
450 to measure the cross-over between London dispersion and
100 200 300 400 500 600 700 800 900 1000
Casimirinteractions,andhavedeterminedcorrectiontoscal-
FIG.4: Realizationofroughinterfaceandsetofmeasurementposi- ingexponentsfortheinteractionsofadisorderedsystems.
tions,×,fortheinteractionenergywhichwillbeseparatedintothe WorkfinancedinpartbyVolkswagenstiftung.
blockZ.Anisotropichorizontalandverticalscales.
[1] W.H.Keesom,Physik.Zeits.22,129(1921).
100
[2] F.London,Trans.FaradaySoc.33,8(1937).
[3] H. B. G. Casimir and D. Polder, Physical Review 73, 360
3F r [4] (I1.D94.8D).zyaloshinskii,E.M.Lifshitz,andL.P.Pitaevskii,Soviet
10−1 Phys.Usp.4(1961).
[5] J. Mahanty and B. Ninham, Dispersion Forces (Academic
Press,1976).
[6] M.Bordag,U.Mohideen,andV.M.Mostepanenko,Phys.Rep.
10−2 353,1(2001).
101 102 [7] K.A.Milton,JournalofPhysicsA37,R209(2004).
r
[8] T.Emig,A.Hanke,R.Golestanian,andM.Kardar,Phys.Rev.
A67,022114(2003).
FIG.5: (1)◦,−hUir3,averagedinteractionbetweendielectricpar-
[9] T.Emig,A.Hanke,R.Golestanian,andM.Kardar,Phys.Rev.
ticleandroughdielectricsurface.(2)⋄,−Usr3,interactionbetween Lett.87,260402(2001).
particle and flat surface. (3) △, σur3, variance of interaction for [10] A. Rodriguez, M. Ibanescu, D. Iannuzzi, J. D. Joannopoulos,
rough surfaces. (4) (cid:3), δUr3, difference in mean interaction en-
andS.G.Johnson,Phys.Rev.A76,032106(2007).
ergybetweenaflatandaroughsurface. Solidlines: r−3.5andr−4.
[11] S.Pasquali,F.Nitti,andA.C.Maggs,Phys.Rev.E77,016705
L=1000. Twoweeksofsimulationtime. Choleskyfactor2.5GB.
(pages11)(2008).
Ng =20. [12] F.Capasso,J.N.Munday, D.Iannuzzi,andH.B.Chan,IEEE
J.SelectedTopicsinQuant.Elec.13,400(2007).
[13] C.-G.Shao,A.-H.Tong,andJ.Luo,Phys.Rev.A72,022102
tionisdisplacedbyδr ∼ ±r1/2 comparedtotheflatsurface. (2005).
Theinteractionbetweenasmoothsurfaceandaparticlevaries [14] K.S.Yee,IEEETrans.AntennasandPropag.14,302(1966).
as U ∼ 1/r3 in the Casimir regime. The interactionof the [15] E.M.LifshitzandL.P.Pitaevskii,StatisticalPhysics,Part2:
s
particle should thus be U ∼ 1/(r +δr)3. If we expand to Volume9(PergamonPress,1980).
[16] D.Irony,G.Shklarski,andS.Toledo,FutureGenerationCom-
first orderwe find that the varianceof the interactionshould
scale as, (△) σ ∼ r−3.5 while the second order expansion puterSystems20,425(2004).
u [17] T.A.Davis,DirectMethodsforSparseLinearSystems(SIAM,
givesashiftinthemeanpotential,hUi,whichvariesas,((cid:3)), Philadelphia,2006).
δU ∼1/r4.Thenumericaldataarecompatiblewiththisscal- [18] G. H. Golub and C. F. V. Loan, Matrix Computations (Johns
ing.Theargumentiseasilygeneralizedtoaffinesurfaceswith HopkinsUniversitypress,1983).
other,lesstrivialroughnessexponentsgivingresultscompati- [19] H.LiandM.Kardar,Phys.Rev.Lett.67,3275(1991).