# Use of Photonic

# Introduction

Photonic is a perl package for calculations on photonics and metamaterials available from CPAN and from Github. It is a modular package using an object system to perform calculations about optical properties of systems with arbitrary geometry and composition. I prepared a couple of examples which may serve as an introduction. A more comprehensive manual is in my plans.

# Simple example on the use of the package.

As every perl program, it must begin by loading the modules to use:

```
use warnings;
use strict;
use feature qw(say);
use PDL;
use PDL::NiceSlice;
use PDL::Complex;
use PDL::Graphics::Gnuplot;
use Photonic::Geometry::FromB;
use Photonic::LE::NR2::AllH;
use Photonic::LE::NR2::EpsL;
```

The package `Photonic::Geometry::FromB`

defines all the geometric
attributes of the system given a characteristic funcion (B). In the
code below, we take a circle of radius (R), taken as a fraction of the lattice
parameter (a) of a square lattice.

```
my $N=50; # We take 2*$N+1 as the number of pixels
my $a=2*$N+1; # on each side
my $R=0.3*$a; # Define the radius of the circle
my $B=zeroes($a,$a) # a square array of (*$2N+1)x(2*$N+1)
->rvals # compute the distance to the center
# # of the array
<$R; # compare to $R
```

Thus, our characteristic function takes the value 1 within a circle of a given radius and 0 outside. We can visualize the geometry by plotting it.

```
my $win=gpwin('png', output=>'../../../../assets/image/20210220ejemplo0.png'); # initialize plotter
$win->plot(
{title=>"Characteriztic function", view=>['equal','xy']},
with=>'image', $B
);
```

The result is in Fig. 6.

Now we construct a geometry object.

```
my $k=pdl(1,0); #Direction of the wavevector k+0 for G=0
my $geometry=Photonic::Geometry::FromB->new(B=>$B, Direction0=>$k);
```

The geometry object can produce lots of information about the
object. For example, it may calculate the reciprocal lattice
**G**, or a set
of unit vectors in the direction of the reciprocal vectors. However,
the case of the vector **G** =0 is special, as it has no intrinsic
direction. Thus, we supplement the constructor above with the
direction we require for this null vector. This corresponds to the
longitudinal direction for the macroscopic field, i.e., it is the
direction of the wavevector **k** of the macroscopic field in the
nonretarded limit **k→0**.

Given a geometry, we can initialize a Haydock calculator object.

```
my $haydock_calculator=Photonic::LE::NR2::AllH->new(
geometry=>$geometry,
nh=>4*$N
);
```

The notation above is the following: `LE`

means we are using a set of
modules based on the longitudinal dielectric function, `NR2`

means to
use the non retarded methods for binary systems, and `AllH`

means this
module obtains all Haydock coefficients corresponding to a given
geometry, stopping after the number of coefficients indicated by the
parameter `nh`

or before, if there is convergence.

With a calculator of Haydock coefficients, we can proceed to construct an object to calculate the dielectric function.

```
my $epsilon_calculator=Photonic::LE::NR2::EpsL->new(
nr=>$haydock_calculator,
nh=>4*$N,
);
```

This object may be used to obtain the complex dielectric function for any pair of values of the dielectric functions of the components (A) (the host) and (B) (the particles). We may use it as follows:

```
my $epsA=pdl(1)->r2C; # vacuum
say "epsA epsB epsM";
foreach(2..5){
my $epsB=pdl($_)->r2C; # some dielectric
my $epsM=$epsilon_calculator->evaluate($epsA, $epsB);
say sprintf "%.3f %.3f %.3f", $epsA->re, $epsB->re, $epsM->re;
}
```

Or else, we may obtain full spectra as follows:

```
my $w=zeroes(100)->xlinvals(.1,1); # frequencies
my $tau=40; #lifetime
my $epsB=drude($w, $tau); # metallic particle
my $epsM=epsilonM($epsA, $epsB, $epsilon_calculator);
```

I have to define the function that computes the dielectric function of the circular inclusions. I use the Drude model for a given frequency and lifetime.

```
sub drude {
my ($w, $tau)=@_;
my $eps=defined $tau? 1-1/($w*($w+i/$tau)) : (1-1/$w**2)->r2C;
return $eps;
}
```

I also need to define a function to calculate the dielectric function
for all the frequencies. For a single frequency I could call
`$epsilon_calculator->evaluate($eps_A,$eps_B)`

, but I want to
calculate for all frequencies in a pdl-like manner.

```
sub epsilonM {
my($epsA, $epsB, $eps, $res)=@_;
$res//=null; #initialize result
epsilonMaux($epsA,$epsB, $res, $eps);
return $res->complex;
}
```

This function uses an auxiliary function that threads over the additional functions.

```
BEGIN {
thread_define 'epsilonMaux(ea(2);eb(2);[o]em(2)), NOtherPars=>1',
over {
my $ea=shift;
my $eb=shift;
my $er=shift;
my $eps=shift;
$er.=$eps->evaluate($ea->cplx, $eb->cplx)
};
}
```

This is a strange but useful way to define threadable functions in
pdl. The idea is to define the function as if it would act on single
objects. In this case, on the single complex dielectric function ε_{a}
and the single complex dielectric dielectric function
ε_{b}. Nevertheless, if supplied with a large pdl, with lots of
dimensions and dielectric functions, then it iterates automatically
over all extra dimensions.

As a comparison, we can calculate the response according to the 2D Maxwell-Garnett approximation.

```
my $f=$geometry->f; # filling fraction
my $beta=$f*($epsB-$epsA)/($epsA+$epsB); #~polarizability
my $epsMG=$epsA*(1+$beta)/(1-$beta); #Maxwell Garnett
```

Now we can plot the results using gnuplot.

```
$win=gpwin('png', output=>'../../../../assets/image/20210220ejemplo1.png'); # initialize plotter
$win->plot(
{title=>
"Dielectric response of square array of Drude cylinders",
xlabel=>"Frequency ω/ω_p",
ylabel=>"ε",
},
with=>"lines", legend=>"Ph Re", $w, $epsM->re,
with=>"lines", legend=>"Ph Im", $w, $epsM->im,
with=>"lines", legend=>"MG Re", $w, $epsMG->re,
with=>"lines", legend=>"MG Im", $w, $epsMG->im,
);
```

The code may be run as

```
./ejemplo.pl
```

The table above results from the `say`

statements in the code
7. The spectrum is plotted in the next figure.

As can be seen, for this filling fraction Maxwell-Garnett is very
close to the numerical calculation. For larger radii, the numerical
calculation shows a dipolar peak that
is further red-shifted from the surface plasmon of a single cylinder
at ω_{sp}=ω_{p}/√2, and from the dipolar resonance of MG
theory, while additional peaks appear closer to $ω_{sp}
corresponding to the excitations of multipoles higher than the
dipole.

# Example of a field calculation

I can also calculate the field. To that end I initialize a field calculator and a Haydock calculator. First I add similr boilerplate code as in the previous example.

```
use warnings;
use strict;
use feature qw(say);
use PDL;
use PDL::NiceSlice;
use PDL::Complex;
use PDL::Graphics::Gnuplot;
use Photonic::Geometry::FromB;
use Photonic::LE::NR2::AllH;
use Photonic::LE::NR2::Field;
use Photonic::Utils qw(tile vectors2Dlist);
```

I use the same geometry as above.

```
my $N=100; # We take 2*$N+1 as the number of pixels
my $a=2*$N+1; # on each side
my $R=0.3*$a; # Define the radius of the circle
my $B=zeroes($a,$a) # a square array of (*$2N+1)x(2*$N+1)
->rvals # compute the distance to the center
# # of the array
<$R; # compare to $R
my $k=pdl(1,0); #Direction of the wavevector k+0 for G=0
my $geometry=Photonic::Geometry::FromB->new(B=>$B, Direction0=>$k);
```

The Haydock calculator has to save the states, so I add the keepStates flag

```
my $haydock_calculator=Photonic::LE::NR2::AllH->new(
geometry=>$geometry,
keepStates=>1,
nh=>4*$N
);
```

I initialize the field constructor using this Haydock calculator.

```
my $field_calculator=Photonic::LE::NR2::Field->new(
nr=>$haydock_calculator,
nh=>4*$N);
```

I compute the field for given dielectric functions. I choose a dielectric function close to -1, the condition for surface plasmon resonance of a single cylinder.

```
my $eps_A_f=1+0*i; #vacuum
my $eps_B_f=-1.5+.1*i; #metallic close to resonance
my $field=$field_calculator->evaluate($eps_A_f, $eps_B_f);
```

To plot the fields I use some utility routines. I use `tile`

to repeat
the unit cell a couple of times and I use `vectors2Dlist`

to build a
vector field ready for plotting with gnuplot. As `tile`

repeats a real
scalar field, I have to move the complex and cartesians indices before
and after tiling.

```
my $repetitions_x=2;
my $repetitions_y=2;
my $repeated_field=tile($field->mv(0,-1)->mv(0,-1), $repetitions_x, $repetitions_y)
->mv(-1,0)->mv(-1,0);
my $repeated_field_abs=$repeated_field->Cabs2->sumover->sqrt;
my $repeated_field_dir=$repeated_field->re->norm;
my $scale=4; #size of arrows
my $decimate=10; #how many pixels between arrows
my @vector_list=vectors2Dlist($repeated_field_dir, $scale, $decimate);
my $win=gpwin('png', output=>'../../../../assets/image/20210220ejemplof.png'); #initialize plotter
$win->plot(
{title=>"Field magnitud and direction", view=>['equal', 'xy'], , cbrange=>[0,25]},
with=>'image', $repeated_field_abs,
{lc=>3, lt=>1}, with=>'vectors', @vector_list);
```

For heavy calculations requiring many states and/or with a very large grid, maybe in 3D, it may be useful to save the states in a file instead of saving them in memory. This is easily done when initializing the Haydock calculator, by adding a filename, as in

```
my $stateFN="states.dat"; # File to save states
my $haydock_calculator_file=Photonic::LE::NR2::AllH->new(
geometry=>$geometry,
stateFN=>$stateFN, keepStates=>1,
nh=>4*$N
);
```

Then the rest of the program is identical but for this haydock calculator.

```
my $field_calculator_file=Photonic::LE::NR2::Field->new(
nr=>$haydock_calculator_file,
nh=>4*$N);
my $field_file=$field_calculator_file->evaluate($eps_A_f, $eps_B_f);
my $repeated_field_file=tile($field_file->mv(0,-1)->mv(0,-1),
$repetitions_x, $repetitions_y)->mv(-1,0)->mv(-1,0);
my $repeated_field_file_abs=$repeated_field->Cabs2->sumover->sqrt;
my $repeated_field_file_dir=$repeated_field->re->norm;
my @vector_list_file=vectors2Dlist($repeated_field_file_dir, $scale, $decimate);
my $win1=gpwin('png', output=>'../../../../assets/image/20210220ejemplof1.png'); #initialize plotter
$win1->plot(
{title=>"Field magnitud and direction", view=>['equal', 'xy'], , cbrange=>[0,25]},
with=>'image', $repeated_field_file_abs,
{lc=>3, lt=>1}, with=>'vectors', @vector_list_file);
```

I run the program.

```
./ejemplof.pl
```

The result is in the following figures.

As expected, both images are identical. Furthermore, we have created a file with the states which may be used for other calculations.

```
ls -l states.dat
file states.dat
states.dat: cannot open `states.dat' (No such file or directory)
```