PDF

Ideal Cloak1
Introduction
Electromagnetic or optical invisibility can be achieved by coating an object with a transparent gradient-index structure that bends the rays of light around the concealed object (Ref. 1). The structure has to be able to do so with radiation incident from any direction, which can be achieved by making it rotationally invariant. Inside this transparent shell, situated is a nontransparent object whose scattering properties in free space are inessential as long as the cloak operates perfectly by preventing any light rays from hitting the object. The concept of invisibility based on omnidirectional cloaking was introduced by Sir John Pendry (Imperial College, UK) and his collaborators in 2006 (Ref. 1).
The cloak of invisibility modeled here is a concentric spherical shell, whose interior surface represents the concealed cavity. Omnidirectional cloaks require anisotropic material properties, which can be calculated using the transformation optics theory (Ref. 1). Although ray and beam bending can be achieved with a gradient of isotropic refractive index, index gradient alone is not sufficient for omnidirectional invisibility. This can be shown by means of the uniqueness theorem that applies to scattering problems involving bodies composed of isotropic materials (Ref. 2).
The refractive index in the azimuthal directions (normal to the radial direction) experiences a gradual change from unity on the exterior surface of the cloak, where it matches free space, down to zero on the interior surface. With a proper choice of index distribution, you can ensure that any ray hitting the cloak never reaches the interior surface, and thus never probes the object. The refractive index in the radial direction is not continuous in this particular cloak design. The resulting index discontinuity at the exterior surface does not lead to reflections because only the index in the tangential direction affects reflectivity.
This model demonstrates the use of optical tracing for studying optically large gradient-index structures with anisotropic optical properties. Additionally, the model introduces a smoothing technique for handling discontinuities of refractive index on curved surfaces, which are typical in conventional optical devices such as lenses.
Model Definition
There is no explicit support for modeling geometrical optics in the Particle Tracing Module, but an analogy between the Hamilton equations and the equations for rays in the zero wavelength limit allows us to solve the problem. The analogy is as follows:
The wave vector, k (SI unit: 1/m) plays the same role in geometrical optics as the momentum, p, of particles in classical mechanics.
The angular frequency, ω (SI unit: 1/s) plays the role of the Hamiltonian, H.
For a classical particle, Hamilton’s equations are:
and using the analogy above:
For geometrical optics, the angular frequency is given by
where n (dimensionless) is the refractive index of the material and c (SI unit: m/s) is the speed of light in a vacuum. For vacuum, the refractive index is simply 1. Inside the cloak, the refractive index is anisotropic, so it is more convenient to express the wave vector using spherical coordinates:
The angular frequency is hence given by:
Results and Discussion
The ray trajectories are plotted in Figure 1. The rays reach the cloak and bend around the inner sphere, which would appear invisible to an observer.
Figure 1: Plot of the light rays traveling through the cloak.
A better way of determining whether the incoming rays are returned to their original trajectory is by using a Poincaré Map or Phase Portrait plot. Figure 2 shows a Poincaré Map in the yz-plane at the initial time step (red dots) and at the final time step (blue dots). In this Poincaré map, the horizontal position indicates the particle’s y-coordinate and the vertical position indicates the particle’s z-coordinate. The particle position after traveling through the cloak is almost exactly the same as it was initially.
Figure 3 shows the change in the particles’ position in the yz-plane after traveling through the cloaking device. The particles at the maximum and minimum y-coordinates have greater absolute error in their final positions despite being deflected at lower angles compared to the particles passing through the middle of the cloak. The higher absolute error may be due to these particles entering the anisotropic domain at a very oblique angle of incidence.
Figure 2: Poincaré map in the yz-plane at for Poincaré sections at x = 1 (red) and x = 1 (blue).
Figure 3: Change in the y-component of the particle position after traveling through the cloaking device.
References
1. J. Pendry, D. Schurig, and D.R. Smith, “Controlling Electromagnetic Fields,” Science, vol. 312, no. 5781, pp. 1780–1782, 2006.
2. A.I. Nachman, “Reconstruction from Boundary Measurements,” Ann. Math., vol. 128, pp. 531–576, 1988.
Application Library path: Particle_Tracing_Module/Tutorials/ideal_cloak
Modeling Instructions
This model comes courtesy of Yaroslav Urzhumov, Center for Metamaterials and Integrated Plasmonics, Duke University.
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Mathematics>Mathematical Particle Tracing (pt).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Time Dependent.
6
Global Definitions
Specify the dimensions of the air domain and the cloak.
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Geometry 1
Block 1 (blk1)
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type 2*L.
4
In the Depth text field, type 2*L.
5
In the Height text field, type 2*L.
6
Locate the Position section. From the Base list, choose Center.
Sphere 1 (sph1)
1
In the Geometry toolbar, click  Sphere.
2
In the Settings window for Sphere, locate the Size section.
3
In the Radius text field, type a.
Sphere 2 (sph2)
1
In the Geometry toolbar, click  Sphere.
2
In the Settings window for Sphere, locate the Size section.
3
In the Radius text field, type b.
4
Click  Build All Objects.
5
Click the  Go to Default View button in the Graphics toolbar.
Definitions
Now add the expressions which transform the refractive index of the cloak from Cartesian to spherical coordinates. The wave vector must also be transformed.
Variables 1
1
In the Home toolbar, click  Variables and choose Local Variables.
Load the variable definitions from a file.
2
In the Settings window for Variables, locate the Variables section.
3
Click  Load from File.
4
Mathematical Particle Tracing (pt)
Using the analogy presented in the introduction section above, enter an expression for the angular frequency, which is called H_photon in this case.
1
In the Model Builder window, under Component 1 (comp1) click Mathematical Particle Tracing (pt).
2
In the Settings window for Mathematical Particle Tracing, locate the Particle Release and Propagation section.
3
From the Formulation list, choose Hamiltonian.
Particle Properties 1
1
In the Model Builder window, under Component 1 (comp1)>Mathematical Particle Tracing (pt) click Particle Properties 1.
2
In the Settings window for Particle Properties, locate the Hamiltonian section.
3
In the H text field, type H_photon.
4
Locate the Particle Mass section. In the mp text field, type 1.
Release from Grid 1
Next, release the particles in the y direction for fixed x- and z-coordinates.
1
In the Physics toolbar, click  Global and choose Release from Grid.
2
In the Settings window for Release from Grid, locate the Initial Coordinates section.
3
In the qx,0 text field, type -1.
4
In the qy,0 text field, type range(-0.38,0.01,-0.02) range(0.02,0.01,0.38).
5
Locate the Initial Velocity section. Specify the v0 vector as
Because the Hamiltonian formulation is being used to model rays, the settings entered for the Initial velocity determine the initial wave vector direction, not the velocity of the model particles.
Mesh 1
The mesh needs to be fine in the cloak region so that the particle trajectories can be computed to a high degree of accuracy.
Free Tetrahedral 1
In the Mesh toolbar, click  Free Tetrahedral.
Size 1
1
Right-click Free Tetrahedral 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
5
Click to expand the Element Size Parameters section. Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section. Select the Maximum element size check box.
7
8
Select the Minimum element size check box.
9
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section. In the Maximum element size text field, type hmax.
5
In the Minimum element size text field, type hmax/2.
6
In the Model Builder window, right-click Mesh 1 and choose Build All.
Study 1
To accurately compute the particle trajectories in an anisotropic medium, the default solver tolerances need to be made more strict. Do this by first showing the default solver, and then reducing the relative and absolute tolerances.
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
Click  Range.
4
In the Range dialog box, type 2[m]/c_const in the Stop text field.
5
From the Entry method list, choose Number of values.
6
In the Number of values text field, type 301.
7
Click Replace.
8
In the Settings window for Time Dependent, locate the Study Settings section.
9
From the Tolerance list, choose User controlled.
10
In the Relative tolerance text field, type 1e-6.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node.
3
In the Model Builder window, click Time-Dependent Solver 1.
4
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
5
From the Method list, choose Runge-Kutta.
6
From the Runge-Kutta method list, choose Dormand-Prince 5.
7
In the Model Builder window, click Fully Coupled 1.
8
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
9
In the Study toolbar, click  Compute.
Results
Particle Trajectories (pt)
The path of the rays is best visualized by adding selections for the cloak inner and outer surfaces.
1
In the Settings window for 3D Plot Group, locate the Plot Settings section.
2
Clear the Plot dataset edges check box.
3
Click to expand the Title section. From the Title type list, choose None.
4
In the Model Builder window, expand the Particle Trajectories (pt) node.
Particle Trajectories 1
1
In the Model Builder window, expand the Results>Particle Trajectories (pt)>Particle Trajectories 1 node, then click Particle Trajectories 1.
2
In the Settings window for Particle Trajectories, locate the Coloring and Style section.
3
Find the Line style subsection. From the Type list, choose Line.
4
Find the Point style subsection. From the Type list, choose None.
Color Expression 1
1
In the Model Builder window, click Color Expression 1.
2
In the Settings window for Color Expression, locate the Coloring and Style section.
3
Clear the Color legend check box.
Surface 1
1
In the Model Builder window, right-click Particle Trajectories (pt) and choose Surface.
2
In the Settings window for Surface, locate the Data section.
3
From the Dataset list, choose Study 1/Solution 1 (sol1).
4
Locate the Expression section. In the Expression text field, type cos_phi.
5
Locate the Coloring and Style section. From the Color table list, choose WaveLight.
6
Clear the Color legend check box.
Selection 1
1
Right-click Surface 1 and choose Selection.
2
Surface 2
1
In the Model Builder window, right-click Particle Trajectories (pt) and choose Surface.
2
In the Settings window for Surface, locate the Data section.
3
From the Dataset list, choose Study 1/Solution 1 (sol1).
4
Locate the Coloring and Style section. From the Coloring list, choose Uniform.
5
From the Color list, choose Gray.
Selection 1
1
Right-click Surface 2 and choose Selection.
2
Definitions
View 1
1
In the Model Builder window, under Component 1 (comp1)>Definitions click View 1.
2
In the Settings window for View, locate the View section.
3
Clear the Show grid check box.
Results
Particle Trajectories 1
1
In the Model Builder window, under Results>Particle Trajectories (pt) click Particle Trajectories 1.
2
In the Particle Trajectories (pt) toolbar, click  Plot.
3
Click the  Go to XY View button in the Graphics toolbar. The plot should look like Figure 1.
To see how well the cloak performs, look at the rays in phase space before and after they pass through the cloak. You can do this in two different ways.
The first method is to define a pair of Cut Plane datasets on the incoming and outgoing sides of the cloak, then plot a Poincaré Map of the ray positions as they intersect each plane.
Cut Plane 1
1
In the Model Builder window, expand the Results>Datasets node.
2
Right-click Datasets and choose Cut Plane.
3
In the Settings window for Cut Plane, locate the Data section.
4
From the Dataset list, choose Particle 1.
5
Locate the Plane Data section. In the x-coordinate text field, type -0.99.
Cut Plane 2
1
Right-click Cut Plane 1 and choose Duplicate.
2
In the Settings window for Cut Plane, locate the Plane Data section.
3
In the x-coordinate text field, type 0.99.
Ray Position Relative to Initial Position
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Ray Position Relative to Initial Position in the Label text field.
3
Locate the Plot Settings section. Select the x-axis label check box.
4
5
Select the y-axis label check box.
6
Poincaré Map 1
1
In the Ray Position Relative to Initial Position toolbar, click  More Plots and choose Poincaré Map.
2
In the Settings window for Poincaré Map, locate the Data section.
3
From the Cut plane list, choose Cut Plane 1.
4
In the Ray Position Relative to Initial Position toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
Poincaré Map 2
1
Right-click Poincaré Map 1 and choose Duplicate.
2
In the Settings window for Poincaré Map, locate the Data section.
3
From the Cut plane list, choose Cut Plane 2.
4
Locate the Coloring and Style section. From the Color list, choose Blue.
5
Select the Radius scale factor check box.
6
7
Click to expand the Title section. From the Title type list, choose None.
8
In the Ray Position Relative to Initial Position toolbar, click  Plot.
9
Click the  Zoom Extents button in the Graphics toolbar. The plot should look like Figure 2.
The second method is to construct a Phase Portrait of the rays and verify that their positions and velocities are the same before and after passing through the cloak.
Change in Lateral Position
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Change in Lateral Position in the Label text field.
3
Locate the Data section. From the Dataset list, choose Particle 1.
4
From the Time selection list, choose Last.
5
Locate the Plot Settings section. Select the x-axis label check box.
6
7
Select the y-axis label check box.
8
In the associated text field, type Change in position (mm).
Particle 1
1
In the Change in Lateral Position toolbar, click  More Plots and choose Particle.
2
In the Settings window for Particle, locate the y-Axis Data section.
3
In the Expression text field, type qy-at(0,qy).
4
From the Unit list, choose mm.
5
Click to expand the Title section. From the Title type list, choose None.
6
Locate the x-Axis Data section. From the Parameter list, choose Expression.
7
In the Expression text field, type at(0,qy).
8
From the Unit list, choose mm.
9
Click to expand the Coloring and Style section. Find the Line markers subsection. From the Marker list, choose Point.
10
From the Positioning list, choose In data points.
11
In the Change in Lateral Position toolbar, click  Plot. The plot should look like Figure 3. You can see that the rays do indeed return to their original position after passing through the cloak.
 

1
This model is courtesy of Yaroslav Urzhumov, Center for Metamaterials and Integrated Plasmonics, Duke University Durham, NC.