Use S4 to calculate optical force

In this tutorial, we will simulate the behaviors of lateral and normal optical forces between two coupled photonic crystal slabs, in an exemplary system as shown in Fig.1. Each slab consists of equally sized alternating regions of high-index material (ɛ=12) and air, with lattice constant a the thickness of each slab is 0.5a. The two slabs are separated from each other by a distance d. Moreover, the two slabs are possibly laterally shifted by Δx along the x-direction. The excitation is assumed to be a normally incident plane wave with electric field in the y-direction. Figure 1.Schematic of the double slab system. The red arrows indicate the directions of the incident light. Each slab consists of an array of high index rods surrounded by air. The empty rectangle indicates the surface over which the integration of the stress tensor is performed. In the right panel, the bottom slab is shifted relative to the first slab by a distance ∆x. [See Liu, et. al. Optics Express 17, 21897 (2009)]

The simulated lateral and normal forces as functions of relative shift between the slabs is shown below. Figure 2.Lateral and normal forces as functions of relative shift in the horizontal direction between the slabs. The vertical spacing between the slabs is d=0.5a. The forces are periodic with respect to the displacement. Only one period in ∆x is shown. [See Liu, et. al. Optics Express 17, 21897 (2009)]

The Lua script for this case is shown below and can be downloaded here.

```S = S4.NewSimulation()
S:SetLattice({1,0}, {0,0}) -- 1D lattice
S:SetNumG(27)
-- Material definition
S:AddMaterial("Silicon", {12,0}) -- real and imag parts
'AirAbove', --name
0, --thickness
'Vacuum') --background material
S:SetLayerPatternRectangle('Slab', -- which layer to alter
'Silicon', -- material in rectangle
{0,0}, -- center
0, -- tilt angle (degrees)
{0.25, 0.5}) -- half-widths
S:SetLayerPatternRectangle('Slab2', -- which layer to alter
'Silicon', -- material in rectangle
{0.15,0}, -- center
0, -- tilt angle (degrees)
{0.25, 0.5}) -- half-widths
-- E polarized along the grating "rods"
S:SetExcitationPlanewave(
{0,0}, -- incidence angles (spherical coordinates: phi in [0,180], theta in [0,360])
{1,0}, -- s-polarization amplitude and phase (in degrees)
{0,0}) -- p-polarization amplitude and phase
S:SetFrequency(0.57)
for dx=-0.5,0.5,0.01 do
S:SetLayer('Slab2', 0.5)
S:SetLayerPatternRectangle('Slab2',
'Silicon', -- material in rectangle
{dx,0}, -- center
0, -- tilt angle (degrees)
{0.25, 0.5}) -- half-widths
T1x,T1y,T1z = S:GetStressTensorIntegral('Spacer', 0.5)
-- Returns the integral of the electromagnetic stress tensor over a unit cell surface normal to the z-direction
--0.5 is surface offset from the beginning of the layer
T2x,T2y,T2z = S:GetStressTensorIntegral('AirBelow', 0)
print(dx..'\t'..T2x-T1x..'\t'..T2z-T1z)
-- A factor of 0.5 due to time averaging is taken into account here
end
```