About the scattering cross section as function of the angle - optical theorem not satisfied

[Mattia Andrini]

Hello,

I am using ADDA to compute the scattering, absorption and extinction efficiencies of spheres and other shapes (usually aggregation of speheres), in presence of a substrate (non-absorbing) and for a spectrum of wavelength. The absorption and extinction are a direct output, so there is no problem, while for the scattering efficiency I need to separate the forward component from the backward. For this reason, I am using the amplitude matrix and integrating in the proper angular subset. I have already double checked with isolated sphere without substrate to see the accuracy with respect to the Mie results, and here it is consistent as expected.

In presence of the substrate and with more complex shapes, though, I don’t get consistent results, and in particular the optical theorem is not satisfied for certain wavelengths, and Csca may be even bigger that Cext. I have seen this was already recognized as a problem (https://github.com/adda-team/adda/issues/195#event-4053496459), and seems to be related to a overestimation of the amplitude matrix component in the angular range ±30°, but I have not fully understood the issue. Does anyone has a possible solution or a way to devise this problem?

Thank you very much.

[Maxim Yurkin]

Dear Mattia,

First, I assume that you are using the Mueller scattering matrix (in file `mueller` or similar), which describes the scattering intensity. By contrast, the amplitude scattering matrix (file `ampl`) describes the fields, which are generally not additive (to produce any meaningful quantity).

Second, you are mentioning complex shapes, have you tried them in the free space? This can be a nice next test to verify your post-processing. In this case, ADDA can do the same with `-Csca`, which should produce results close to Cext-Cabs (within errors of integration and convergence of the iterative solver). See Section 11.4 "Integral scattering quantities" of the manual. Note that for non-symmetric shapes, you need to integrate the corresponding elements of the Mueller matrix over the whole solid angle (`-store_scat_grid`) and the proper integrand is given by Eq.(66) (in manual v.1.4.0).

For instance, you can try the command line:

```

adda -shape ellipsoid 2 3 -m 1.5 0.2 -grid 16 -prop 1 2 -3 -Csca -store_scat_grid

```

with default files `alldir_params.dat` and `scat_params.dat` from `input/` and compare the results of post-processing `mueller_scatgrid` to Csca found in `CrossSec-X/Y` files. The optical theorem for this run is satisfied with relative accuracy within 0.01%, well below the integration uncertainty in the end of files `log_int_Csca-X/Y`.

 

Third, concerning the substrate, if it is absorbing, then you cannot hope for satisfaction of optical theorem at all. This, and other potential ambiguities, has been recently discussed in Section 2.1 of 
Virkki A.K. and Yurkin M.A. Microwave scattering by rough polyhedral particles on a surface, J. Quant. Spectrosc. Radiat. Transfer 345, 109547 (2025).

For non-absorbing substrate, the optical theorem should hold, but the proper Csca is defined by Eq.(58) of the manual. For the normal incidence, the integrand is similar to that in free space, but with extra factor for theta > 90 deg (scattering into the substrate). Note, that for the surface mode, the angles are always defined in the laboratory reference frame (See Section 11.1 "Definition of scattering plane and angles") of the manual. The main problem here is for defining the proper integrand for non-normal incidence, since it may involve some non-trivial rotation from the scattering plane (through incidence and scattered direction) to the laboratory reference frame. ADDA does something similar internally, so something can be deduced by looking at this code.

Similarly to the above, the command line:

```

adda -shape ellipsoid 2 3 -m 1.5 0.2 -grid 16 -prop 1 2 -3 -Csca -store_scat_grid -surf 10 2 0

```

can be used for testing. Here the optical theorem is satisfied only within 5%, in accordance with much larger integration uncertainty (probably due to abrupt change of scattering intensity at theta=90 and 150 degrees (the latter due to the total internal reflection).  However, the accuracy can be improved to better than 0.4% by incrementing both Jmax in `alldir_params.dat` by 2.

To conclude the above, all the formulaes to compute Csca (equal to Cext-Cabs) are known and implemented in ADDA for free-spaced or non-absorbing substrate. Still, it would be great (helpful for ADDA users) to implement a (Python) script to compute the same from `mueller_scatgrid`. Then, it can be easily adapted to compute the scattering integral in any solid angle. Comparison with built-in Csca from ADDA is an easy test for it. If you manage to do it, at least for some of the cases mentioned above, please consider contributing it, e.g., to `examples/` folder. In some simplified form it can also qualify for `tests/equiv`, but this can be discussed later.

Finally, there is also a workaround, if you are interested only in some ranges of theta (with full range of phi, or at least with the one independent of theta) - relative to the incident direction in the free space and relative to the surface normal near the substrate. Then you can use built-in ADDA calculation of Csca with preliminary modification of ranges in `alldir_params.dat`. The only trick here is to rescale the result according to the resulting solid angle. The fact that such rescaling is necessary - is exactly the topic of the issue that you referred to.

Let me know, if you need more details on any of the above,

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.  

[Mattia Andrini]

Dear Yurkin,

Thank you very much for all the detailed explanation. I am sorry that the answer took more time than I thought, but I had to check with different simulation.

I was not using Mueller matrix, but I have now reverted to it as it is indeed much more simple and straightforward to work with. My substrate is not absorbing (k=0), thus the optical theorem should hold, and the light is normal with respect to the substrate surface. 

Now that I have defined the scattering cross section in the proper way (eq. 66 of the manual), I confirm that the Csca output from adda (using option -Csca) is consistent with the Cext-Cabs (taken from CrossSec-Y/X), at least up to some 0.something %. 

I have also made a python script that takes the S11 component as function of theta and phi (using the -store_scatgrid option) and integrates it in the backward and forward hemisphere with respect to the impinging light direction, taking into account the eventual substrate presence and the proper reference frame. As expected, I have noticed that the simulation time can skyrocket with a too complex theta-phi grid. I was wondering your opinion about this, now I am working with 721 values for theta and 361 for phi. Probably a bit excessive... I will try to run some test to determine the accuracy and avoid wasting too much time.

It may also be of interest (even if expected) that considering the phi dependence is fundamental for very asymmetric particles (in the plane perpendicular to the light direction). If the option -store_scatgrid is not included and the Csca integration is thus performed only considering the theta dependence, the results are not consistent, but this still will need further investigation.

I have double checked that instead by considering both theta and phi, the integrated Csca is consistent with the Csca calculated internally by ADDA and with Cext-Cabs computation, and it does take into account that the cross section coming from the light propagating into the substrate should be divided by the substrate refractive index squared. I have also checked with asymmetric and not standard particles. This script is meant to perform simulation for a spectrum of wavelength in parallel (user can decide how many at the same time), but of course can be used for a single wavelength, and complex shapes can be taken care of (as long as user has the .geom file - but I am also working in this direction). It is basically a standalone script that runs simulation in parallel, reads all the folder outputs, do the integration and save the desired quantities in simple .txt files as function of the wavelength. If you think it can be useful to the community, I will gladly share it as soon as I write some basic documentation (at this stage it is still not really user friendly).

 If something is not clear, please let me know

Again, I thank you for your kind help.

Cheers,

Mattia

[Maxim Yurkin]

Hi Mattia,

> I have also made a python script that takes the S11 component as function of theta and phi (using the -store_scatgrid
> option) and integrates it in the backward and forward hemisphere with respect to the impinging light direction, taking
> into account the eventual substrate presence and the proper reference frame. As expected, I have noticed that the
> simulation time can skyrocket with a too complex theta-phi grid. I was wondering your opinion about this, now I am
> working with 721 values for theta and 361 for phi. Probably a bit excessive... I will try to run some test to
> determine the accuracy and avoid wasting too much time.

Yes, this seems large, but ultimately it depends on your accuracy requirements and the specific algorithm you use for
integration. You can compare them with the default values in file `alldir_params.dat`, which are (if you grind through
the definition of Jmax in the comments) 65 and 32 for theta and phi, respectively. However, as mentioned previously, I
changed it to 257x128 for particles on substrate to get good accuracy. In principle, built-in integration is based on
adaptive (high-order) Romberg algorithm (for theta), but it is expected to be not efficient for discontinuous functions.
When you handle the jumps explicitly (cutting the range at theta=90, but maybe also at total-internal reflection angle),
you can definitely make it more efficient in combination with some high-order quadrature (or at least some integration
rule). For phi integration, the simplest trapezoidal rule on the uniform grid is the most efficient one due to the
periodicity of the integral, so it may be sufficient to use about 32 angles for wavelength-sized particles (but this
needs more testing). For example, I have just tested the grid 257x32 and it resulted in almost the same integration
accuracy as a denser one. So the number of phi values may even be decreased further after some tests.

> It may also be of interest (even if expected) that considering the phi dependence is fundamental for very asymmetric
> particles (in the plane perpendicular to the light direction). If the option -store_scatgrid is not included and the
> Csca integration is thus performed only considering the theta dependence, the results are not consistent, but this
> still will need further investigation.

I guess, it would be interesting to have a separate (faster) simulation option for symmetric particles (based on the
Mueller file), but you do not have as much flexibility here in terms of varying the grid of theta (it must always be
uniform). Then, some simple function to test it against the solid-angle integration is also desirable.

Similarly, you can also implement use of `-Csca` with modified theta ranges in `alldir_params.dat` (that I discussed
previously) and testing against the explicit integration. (based on your text below, I am not sure if you've already
tried it).

> I have double checked that instead by considering both theta and phi, the integrated Csca is consistent with the Csca
> calculated internally by ADDA and with Cext-Cabs computation, and it does take into account that the cross section
> coming from the light propagating into the substrate should be divided by the substrate refractive index squared. I
> have also checked with asymmetric and not standard particles. This script is meant to perform simulation for a
> spectrum of wavelength in parallel (user can decide how many at the same time), but of course can be used for a single
> wavelength, and complex shapes can be taken care of (as long as user has the .geom file - but I am also working in
> this direction). It is basically a standalone script that runs simulation in parallel, reads all the folder outputs,
> do the integration and save the desired quantities in simple .txt files as function of the wavelength. If you think it
> can be useful to the community, I will gladly share it as soon as I write some basic documentation (at this stage it
> is still not really user friendly).

This sounds great. I would only put a bit more emphasis on core functions, like the ones to (efficiently) calculate
hemi-spherical values for any ADDA command line. The wrappers on top of it (like scanning over wavelengths) can also be
useful for users, depending on their applications, but they should not hide the core functions. Thus, the latter should
be easy to reuse in a wide variety of applications. Once you are ready with the code, take a look at
https://github.com/adda-team/adda/wiki/InstructionCommitters#examples .

Maxim.

[Mattia Andrini]

 "built-in integration is based on adaptive (high-order) Romberg algorithm (for theta), but it is expected to be not efficient for discontinuous functions. When you handle the jumps explicitly (cutting the range at theta=90, but maybe also at total-internal reflection angle), you can definitely make it more efficient in combination with some high-order quadrature (or at least some integration rule)". 

To be fair, I am trying to integrate the scattering cross section in the easiest possible way: I sum (S11*sin(theta)*dtheta*dphi), in this way it is just sufficient to handle with care the S11 indexes to sum everything in the right order. I still have to perform more test, but from the preliminary results I found discrepancies between my calculation, ADDA built-in integration and (Cext-Cabs) to be around 0.01%, at least for asymmetric particles without the substrate and with angle grid of 721x361.

So I guess I have margin for a smaller grid, as soon as I check accuracy including the substrate.

I guess, it would be interesting to have a separate (faster) simulation option for symmetric particles (based on the Mueller file), but you do not have as much flexibility here in terms of varying the grid of theta (it must always be uniform). Then, some simple function to test it against the solid-angle integration is also desirable.

This is a good suggestion, indeed I have also this option in the code. I have verified that the Mueller file with only theta dependence is giving decent accuracy for symmetric or kind-of symmetric particles (for example, three spheres arranged in a triangular fashion), now I need to quantify this accuracy.

Thank you again for your time, I will update this thread after a vacation period with some quantification of the integration accuracy.

Best,

Mattia

[Mattia Andrini]

Dear Maxim,

I am sorry this has taken a while. I have committed a change in the example folder on github. I am not sure if this is the correct way to add the example, if I should do it differently please let me know.

Thank you again for your help

Kind regards,

Mattia

Il giorno martedì 29 luglio 2025 alle 12:22:00 UTC+2 yur...@gmail.com ha scritto:

[Maxim Yurkin]

That's great, Mattia. I see the branch patch-1 in your fork of ADDA. But you need to additionally create pull request.
You probably can press the button "Create pull request" at
https://github.com/adda-team/adda/compare/master...m51vortice:adda:patch-1

Maxim.