Interpretation of Mie scattering

Question

I’m usually on Physics SE, but this question seems to fit better into Earth Sciences:

I’m interested in radiation transfer and in particular Mie scattering in the atmosphere. So I did some calculation with Miepython (https://pypi.org/project/miepython/). Below is an example result for extinction, scattering and backscattering on small water droplets of 1000 nm diameter.

Unfortunately, I’m a bit puzzled about the results.

First, is there any common definition of the angular range in which scattering and backscattering is? The Miepython manual doesn’t state anything about. E.g., is scattering maybe in 0$^\circ$-180$^\circ$ and backscattering in 90$^\circ$-180$^\circ$?

Second, the backscattering shows a lot of variations with wavelength. If I increase the particle size, these get even more. What is the physical explanation for this zig-zag behaviour of alternating maxima and minima?

Update

After receiving helpful comments I calculated a polar plot, following instructions from https://miepython.readthedocs.io/en/latest/03_angular_scattering.html (incl. axes labeling...) for unpolarized radiation.

The result is:

However, I tried to integrate (sum up) over different ranges and multiplied by the geometric cross section (what miepython requires you to do, and which has been done already for the plots above), but couldn't reproduce numbers of the above plot at 500 nm.

Answer 1

Based on the comments and after a lot of trial and error plus reading carefully the Miepython documentation (!), I managed to answer my question:

For the problem that I wasn't able to reproduce scattering and backscatter efficiencies $qsca$ and $qback$ (or corresponding cross sections $\sigma$ in the first plot from values of the phase function in the second plot (update) in my question: I had two problems, or better say there were two pitfalls:

• I made an equidistant $\vartheta$-vector and converted to $mu = cos(\vartheta)$, calculated the phase function $p(\mu)$ and plotted it as polar plot and histogram, which I tried then to integrate (sum up). That was of course not working as $\mu = cos(\vartheta)$ is no longer equidistant. Now I'm doing an equidistant $\mu=cos(\vartheta)$-vector for calculation and integration. For plotting against $\vartheta$ I just calculate $\vartheta$ from $\mu$ using $arccos()$. The result is non-equidistant datapoints in the histogram then, but when the resolution is high enough, you don't see it..
• There is a normalisation of the phase function in miepython, which is by default set to the single scattering albedo $a = qsca/qext$. So after summing up (integrating), to get $qsca$ you need to multiply by $qext$. (But there are other options you can pass to Miepython for normalisation, the default is just $qsca/qext$... see https://miepython.readthedocs.io/en/latest/03a_normalization.html)

With this, I could reproduce values for $qsca$ and $qback$ at 500 nm. Here is the code:

import numpy as np
import miepython
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 10})

###############################################################################
# Parameters
###############################################################################
m = 1.335 - 1E-9j  # complex refractive index for water droplet at 500 nm from https://refractiveindex.info/?shelf=main&book=H2O&page=Daimon - 24.0C
lambda0 = 0.5 # in mu, i.e. 500 nm
diameter = 1 # in mu
r = 0.5*diameter  # radius of aerosol
x = 2*np.pi*r/lambda0 # size parameter

###############################################################################
# Miepython calculation
# - Calculate scattering efficiencies
# - Calculate phase function and plot as polar plot and histogram
###############################################################################
mu = np.linspace(1,-1,100001) # make mu eqidistant and calculate (for the plots) theta  from it, not vice versa (important as we integrate later the histogram, which is a function of mu, so mu needs to be equidistant for easy sum up -> integrating)
theta = np.arccos(mu)*180/np.pi

# Calculate extinction, scattering and backscatter efficiency
qext500, qsca500, qback500, g2_500 = miepython.mie(m,x)

# Now calculate the phase function
normalisation='albedo' # default is albedo a = qsca/qext, see https://miepython.readthedocs.io/en/latest/03a_normalization.html
p_mu = miepython.i_unpolarized(m, x, mu, normalisation) # p_mu is now normalised to a = qsca/qext
# Test alternative way of calculating the phase function:
s1, s2 = miepython.mie_S1_S2(m, x, mu)
phase = (np.abs(s1)**2 + np.abs(s2)**2 )/2
if np.allclose(p_mu, phase):
    print('p_mu and phase are the same vectors')

# plot phase function as polar plot and histogram
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12,5))
ax1=plt.subplot(121, projection='polar')
ax1.plot(theta/180*np.pi, p_mu, color='blue')
ax1.set_title(r'm = ' + str(np.around(m,9)) + ', d = '+str(diameter) + r' $\mu$m, $\lambda$ = ' + str(np.around(1000*lambda0,0)) + r' nm')

ax2 = plt.subplot(122)
ax2.plot(theta, p_mu)
ax2.set_xlabel(r'$\vartheta$ ($^\circ$)')
ax2.set_ylabel('Unpolarized Scattered light [1/sr]')
ax2.grid(which = 'minor', color='grey', linestyle=':')
ax2.grid(which = 'major', color='grey', linestyle=':')
ax2.set_title(r'm = ' + str(np.around(m,9)) + ', d = '+str(diameter) + r' $\mu$m, $\lambda$ = ' + str(np.around(1000*lambda0,0)) + r' nm')

plt.tight_layout()
plt.show()

fig.savefig('phase function.png', dpi=fig.dpi, bbox_inches='tight')
  
###############################################################################
# Now integrate histogram to reproduce scattering and backscatter effiencies 
###############################################################################
# print values for comparison
print('qback500 = ' + str(qback500))
print('qsca500 = ' + str(qsca500))
print('qext500 = ' + str(qext500))
print('qsca and qext is almost the same, because of almost no absorption (imaginary part of complex refractive index, so extinction (loss) is mainly caused by scattering...)')

print('\nNow calculating the scattering efficency from the phase function:')

# spacing dmu for integration
oben = np.abs(mu[1:])
unten = np.abs(mu[:-1])
dmu = np.mean(np.abs(oben - unten))
print('\ndmu = ' + str(dmu))
print('dmu is the spacing on the mu = cos(theta) axis')

# qsca from histogram:
print('Integration of phase function over full solid angle:')
print('- There is no dependency on the azimuth, so integration over azimuth gives 2*pi and the remaining integration is over mu only.')
print('- The integration over mu is from 1 to -1, i.e. from theta=0° to 180°')
integral_p_mu = 2*np.pi*dmu*np.sum(p_mu)
print('Integral over p_mu histogram = ' + str(integral_p_mu))
# To get qsca from this, remind that p_mu was normalised to a = qsca/qext. Thus, to get qsca, one needs to multiply the obtained result here with qext:
print('\nqsca calculated from phase function = ' + str(integral_p_mu * qext500))

print('\nBackscattering efficency calculated from physe function')
idx180 = np.where(np.abs(180 - theta) < 0.01)[0][0] # index of backward direction
# Again, we need to multiply by qext because of the normalisation 
# We use only the 180° datapoint of the phase function, but still need to integrate over the solid angle, which just gives 4*pi')
print('p_mu[idx180]*qext500*4*np.pi = ' + str(p_mu[idx180]*qext500*4*np.pi))

And the answer to my original question is:

• $qsca$ is the total scattering efficiency in all directions (integration over a full sphere of $4\pi$), i.e. not only forward half-dome or so..
• $qback$ is scattering in the backward direction of $180^\circ$ only (one single direction, not the entire rear half-dome or so...)