Circular dichroism (CD) spectra simulation

[Hau T]

Hey ADDA team,

I am working on some chiral particles, and would like to know the feasibility of simulating their (far-field) CD spectra using ADDA.

I read some topics on adda discussion group (#30, #192). Knowing I'm not good at the mathematics behind the extinction matrix method, my plan is to construct the left and right circular polarized light (LCP/RCP) using X and Y plane wave with +-90 degree phase difference, then calculate the difference in their extinction cross-section. However, it turned out my extinction cross-section of the particle under LCP and RCP is 0 over all ranges.

I was wondering if my approach of using ADDA to simulate CD is fundamentally wrong or not supported by ADDA yet. I was able to obtain the X and Y dipole polarization of the discretized particle, is there a way from there to calculate the CD spectra?

Some related works I've aware are 

1. Circular Dichroism Simulated Spectra of Chiral Gold Nanoclusters:  A Dipole Approximation. J. Phys. Chem. B 2003, 107, 44, 12035–12038
   
2. pyGDM -- new functionalities and major improvements to the python toolkit for nano-optics full-field simulations.  arXiv:2105.04587 (It include a section discussing the near-field optical chirality)

Thanks a lot,

Hau

[Maxim Yurkin]

Dear Hau,

Indeed, this should be possible to calculate CD with ADDA, at least in fixed orientation, although I am not aware of anybody who has actually done it (if somebody from the group can prove me wrong, please do). However, it was done with DDSCAT - see, e.g., https://doi.org/10.1007/s11468-014-9821-1 .

Both methods are fine (through building incident circular polarization or through the amplitude matrix), although the latter is generally easier. I will add the specific formulae for that to the manual. However, in order to test them, it would be great if you can provide some simple shape file (relatively small) of a chiral particle.

And concerning LCP and RCP please also provide the input beam files (for the same shape, mentioned above) and the command lines that you used for simulation (for a single wavelength), or you may provide corresponding log files.

Maxim.

P.S. This answer has been forwarded to your e-mail address for your convenience. However, if you want to continue the discussion please reply to the group's e-mail. You are also advised to check the corresponding discussion thread at http://groups.google.com/group/adda-discuss and/or subscribe to it to receive automatic notifications, since other answers/comments will not be necessarily forwarded to your e-mail.

--
You received this message because you are subscribed to the Google Groups "ADDA questions and answers" group.
To unsubscribe from this group and stop receiving emails from it, send an email to adda-discuss...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/adda-discuss/CABVeukzf%2BrvzR6XFHUGzEW3UmoK2B3zQTCGqPgcY1i9%2Bt9ysdA%40mail.gmail.com.

[Hau T]

Hi Maxim,

Thank you for your explanation.

Attached is a zip file 'particle.zip' containing 

1. 'shape.dat': the sample I drew using Fusion360 and converted into DDSCAT shape format using pip. It is a 20*(40+20) nm twisted rod, with 10 dipoles along the diameter and 30 dipoles along the length.
2. 'IncBeam-Y' (same for LCP and RCP), 'IncBeam-X-rcp' and 'IncBeam-X-lcp':  -store_beam generated beam file for (0,0,0) orientation (I followed the AddingBeam procedure to add the LCP/RCP beam type in order to try orientation averaging)
   
3. 'log-lcp' and 'log-rcp': log files for LCP and RCP

PyGDM shows this structure has some near-field optical chirality, although I am not sure how much of that can translate to the CD spectra.

Regards,

Hau

[Maxim Yurkin]

Hi Hau,

Here are some intermediate results:

1) your beam files are exactly negative of each other - that is not correct (and explains why you get the same Cext for them). So there seems to be an error in the code that you added.

2) the formula for CD is Cext,L - Cext,R = (4pi/k^2)*Im(S4-S3), if R and L are defined by viewing from the source (along the propagation direction). There seems to be two notation conventions with respect to the latter - https://en.wikipedia.org/wiki/Circular_polarization#Handedness_conventions . If you know, which one is commonly used in CD literature, please let me know.

I have run your shape with "-scat_matr ampl -ntheta 2" and the above formula leads to non-zero (positive) result. Actually, Im(S3) = - Im(S4) to good accuracy in this case. You may test this formula after you fix the beam generation.

3) If orientation-averaged result is needed for such small particles, it can be obtained from two or three orientations (or incident directions), as specified, e.g., in the DDSCAT paper that I mentioned above. However, it can also be obtained from the full polarizability tensor of the particle, which in turn can be obtained from the amplitude scattering matrix. I will provide the corresponding formulae later.

Maxim.

[Maxim Yurkin]

Hau, here is another update. I have just commited a test script `tests/equiv/ext_CD.py` which shows how Cext for both
linear and circular polarizations can be obtained from the forward-direction amplitude matrix. It also includes your
helix shape - hope that you don't mind.

Several notes:
1) The handiness convention there is reverse from my previous message, since it seems to be more common in CD literature
(although I still need a broader overview of this issue). So there the R and L polarizations are defined by viewing from
the detector [ e_R,L = (e_y +- i*e_x)/sqrt(2) ], which corresponds to V parameter of the Stokes vector equal to +1 and
-1, respectively. This definition leads to
CD = Cext,L - Cext,R = (4pi/k^2)*Im(S3-S4)

2) I have noticed that the default iterative solver breaks down for circular incident polarizations. I haven't studied
this issue in details, but it is easily fixed by using BiCGStab instead.

3) As with other particle symmetries and reciprocity relations, the above formula (equivalent ways to obtain Cext) is
correct to the residual of the iterative solver. While the default value (1e-5) of the latter is sufficient for vast
majority of applications, it should be decreased if many-digit agreement is required (or tested). For instance, this
test script tests for relative difference to be less than 1e-8, thus, it uses '-eps 10' in ADDA runs.

I would appreciate your feedback on the script and whether you will be able to replicate any published results on CD (to
verify the handiness convention).

Maxim.

[Hau T]

Hi Maxim,

Sorry for the last reply, it has been busy for me during the past weeks.

Your suggestions are very helpful. It seems for asymmetric shapes under plane waves, X and Y polarization will be applied separately to calculate X and Y cross section. For circular polarized light, after I combined X and Y together in a single Einc, I can now get CD signals same as yours.

There is one problem I am still stuck at -  is it normal that the X and Y extinction cross section under a circular polarized light are always the same for an asymmetric shape (e.g., longitudinal and traverse of the rod)? Under CPL, when I use a rod with its longitudinal direction on z-axis, the XY extinction cross section are same traverse distinction, but when I rotated it with '-orientation 0 90 0', the XY extinction cross section changed to be both longitudinal extinction (which I would expect one to be traverse and the other to be longitudinal). Meanwhile, this is not an issue for plane wave and I can obtain both longitudinal and traverse distinction with such orientation.

Thanks,
Hau

[Maxim Yurkin]

Dear Hau,

Not completely sure that I correctly understood your question. But probably it is about Cext for right- and
left-circular polarizations (Cext_R and Cext_L) for such asymmetric particles that standard x- and y-polarizations lead
to different Cext. Then the answer is that those are two different kind of symmetries: one is the asymmetry between
sizes along different axes, another - existence of the symmetry plane (or chirality).

Generally, if the particle is symmetric with respect to the scattering plane (yz by default in ADDA), then its amplitude
matrix is diagonal, since the amplitude matrix for two symmetric configurations have different signs of off-diagonal
elements - see (a) and (c) on p.49 of (van de Hulst 1981) book. If there is a symmetry plane of the particle containing
the incident direction, then a diagonal amplitude matrix at forward direction for such scattering plane can be rotated
into the other scattering plane using formulae on p.52 of that book. This implies that S3=S4 for any such scattering
plane and, hence, Cext_R = Cext_L. This agrees with physical reasoning that Cext_R and Cext_L are invariant to rotation
around the propagation direction.

This can be tested by
adda -shape box 2 3 -scat_matr ampl -ntheta 2 -eps 10
adda -shape box 2 3 -scat_matr ampl -ntheta 2 -orient 30 0 0 eps 10

first line leads to S3=S4=0, while second - to S3=S4.
shape is 1:2:3 cuboid
'-orient 30 0 0' leads to rotation of the scattering plane (since only alpha angle is involved)
'-eps 10' leads to better convergence of the iterative solver minimizing the "numerical" symmetry breaking; the
relations are then satisfied to relative error of the order 1e-10

A more complicated question is what happens when the particle has symmetry plane(s), but none of them contain the
propagation direction. I have tried:
adda -shape box 2 3 -scat_matr ampl -ntheta 2 -orient 10 20 30 -eps 10
The existing theory (which I am aware of) proves that circular dichroism is absent when some averaging over orientations
is involved (but not in fixed orientation as discussed here). Simulations do show finite difference between Cext_R and
Cext_L (0.05% in this example), but it tends to decrease with refining discretization - 0.015% and 0.004% for twice and
four times smaller dipoles, based on the runs:
adda -shape box 2 3 -scat_matr ampl -ntheta 2 -orient 20 40 60 -eps 10 -dpl 30 -grid 32
adda -shape box 2 3 -scat_matr ampl -ntheta 2 -orient 20 40 60 -eps 10 -dpl 60 -grid 64
So that seems to be a discretization artefact since the set of dipole do not necessarily have the same symmetry as the
original particle.

To obtain a definite answer one need to perform systematic simulations with refining discretization to obtain converged
result for S3 and S4, and then see whether they significantly differ. While for a cuboid, the situation is rather
certain, I have also tried limited simulations of a triangular prism on my laptop - here I was not able to obtain a
definite converged result for S3 and S4 (they jump abruptly from one dipole size to another). So this remains an open
question - what is a minimal set of particle symmetries that guarantee Cext_R = Cext_L for a fixed orientation?

Coming back to your original example of a rod - it has an axis of symmetry, hence the plane through this axis and the
axis of propagation is always a symmetry plane (for any orientation). For some orientations (like in your example) it is
exactly conforming to the discretization (i.e. a set of dipoles is also perfectly symmetrical with respect to this
plane) - then any difference between Cext_R and Cext_L may span only from the residual of the iterative solver, but it
is on the order of 1e-5 by default (hardly noticeable). For other orientations, the difference will exist, but it must
vanish with refining discretization.

Maxim.

[Maxim Yurkin]

One more update about circular polarizations. They can be directly obtained from current version of ADDA (in the master branch) using the built-in Bessel beams. Most straightforward is to add the following to the command line:

```

-scat_matr none -sym enf -iter bicgstab -beam besselM 0 0 0 0 0 0.7071067812 0 0 7071067812 0

-scat_matr none -sym enf -iter bicgstab -beam besselM 0 0 0 0 0 0.7071067812 0 0 -7071067812 0

```

for right and left circular polarizations, respectively. The corresponding Cext, etc. can be found in files CrossSec-Y. 

A few comments:

- the numerical values are 1/sqrt(2) (+ or -)

- I have explained the need for `-iter bicgstab` previously (we actually noticed such behavior for Bessel beams long ago, but I have only recently remembered that)

- `-sym enf` ensures that only one polarization is computed (for each run), since the other one in the framework of Bessel beams will be obtained by rotation (leading to the one proportional to the original one for circular polarization) and is, thus, redundant.

- ` -scat_matr none` turns off the calculation of both Mueller and amplitude matrices. On the one hand, they will be hard to interpret anyway. If one tries to extend the definition of these matrices to some beams other than linearly-polarized plane waves, it is natural to require two incident beams to be related by rotation. However, this makes little sense for circular polarizations (for which rotation is equivalent to multiplication by a scalar). This has been discussed in Section IV.A of Glukhova S.A. and Yurkin M.A. Vector Bessel beams: General classification and scattering simulations, Phys. Rev. A 106, 033508 (2022). (PDF). 

On the other hand, if the scattering matrices are desired, it is best to calculate them for standard linear polarizations. Then they can be used to compute any measureable quantities for circular polarizations, as we discussed previously. 

This way to compute CD values has been implemented in the script `/tests/equiv/ext_CD.py` (I have just committed a new version), where it is successfully tested against the previous options. The updated script also tests several other ways to imitate circular polarization through Bessel beams (see comments inside). Some of them are significantly more complicated.

Note also, that you can add `-store_beam -prognosis` to the above command lines to obtain incident circular polarizations without any Python post-processing.

Maxim.