# Light scattering by large particles – a refinement to the Mean-Field T-matrix method –

Robert Botet

► To cite this version:

Robert Botet. Light scattering by large particles – a refinement to the Mean-Field T-matrix method – Journal of Quantitative Spectroscopy and Radiative Transfer, 2025, 347, pp.109589. 10.1016/j.jqsrt.2025.109589. hal-05247390

HAL Id: hal-05247390

https://hal.science/hal-05247390v1

Submitted on 9 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Light Scattering by Large Particles
– a refinement to the Mean-Field T-matrix method –

Robert Botetª

ªUniversité Paris-Saclay, Laboratoire de Physique des Solides, CNRS UMR8502, 91405 Orsay, France

## Abstract

In the realm of light scattering computational techniques, the T-matrix method serves as a pivotal tool for analysing dust particles composed of aggregated spherical monomers. However, an inherent and significant limitation restricts its application to small aggregates due to the computational demands of processing the translation matrix, whose size becomes prohibitive for aggregates with a large number of spheres. The Mean-Field T-matrix method (MFTM) was proposed as an approximation to address this issue, enabling the treatment of arbitrarily large particles. In the present study, the MFTM method is extended to encompass dust particles that are large, irregular aggregates of polydisperse spherical monomers made from various materials. Additionally, a recently proposed accurate pair correlation function is integrated to describe the spatial distribution of monomers within the aggregate. These advancements render the refined MFTM method both efficient and capable of processing realistic, large dust aggregates. This paper presents a comprehensive application of the enhanced MFTM method to fractal aggregates of spheres with radii following a Schulz distribution. Notably, the refined method is expected to facilitate the derivation of optical properties of atmospheric aerosols and the interpretation of observations from protoplanetary disks.

## Keywords:

fractal; light scattering; polydispersity; T-matrix; cross-sections; pair correlation function

## 1. Introduction

Dust particles in natural environments – ranging from Earth's atmosphere to protoplanetary disks – exhibit intricate shapes, diverse compositions, and highly variable sizes. Rather than existing as simple grains, these particles often form large, porous aggregates made of smaller monomers, displaying fractal or irregular structures. Such aggregates play a crucial role in numerous physical processes, influencing light scattering, thermal emission, and radiative transfer in both atmospheric and astrophysical settings.

In the Earth's atmosphere, the growth of such aggregates is primarily limited by aerodynamic forces and gravity, with carbonaceous aerosols and mineral dust particles reaching sizes up to a few hundred micrometers [1, 2, 3]. These constraints are supported by laboratory studies showing aggregate bouncing or compaction beyond certain size thresholds [4, 5, 6]. Aggregates larger than a few hundred micrometers either disintegrate [7, 8, 9] or undergo compression [10, 11]. In contrast, the low-pressure conditions of astronomical environments allow the formation of even larger and more tenuous aggregates. For instance, particles collected by the ROSETTA mission exceeded 100 micrometers in diameter [12, 13], while low-density aggregates capable of lifting from comet nuclei were estimated to reach up to
730 micrometers [14]. Theoretical models in protoplanetary contexts predict aggregates spanning up to millimeters or more [15, 16, 17, 18, 19, 20]. Such large particles, with low porosity or fractal structures, can explain polarization and light flux observations in protoplanetary disks [17, 21], with sizes ranging from 1 to 100 micrometers commonly acknowledged [22, 23].

These large aggregates, though few in number, can dominate the optical response of a dust population due to their large scattering cross-sections. Moreover, the situation is further complicated by polydispersity in monomer sizes [22, 23, 24? , 26], the presence of multiple grain types, and subtle structural features such as necking or overlapping between monomers [27]. A few works on the optics of small aggregates with polydisperse monomers have been published up to now [28, 29, 30, 31]. The general consensus is that scattering cross-sections increase with monomer polydispersity, but no reliable result exists yet for large aggregates. Despite their relevance, these ‘minor structures’ are often overlooked in optical modeling [32, 33], leaving a gap in our understanding of how dust clouds interact with light across different environments.

Several rigorous numerical methods have been developed to model light scattering by complex aggregates, such as the Discrete Dipole Approximation (DDA) [34] and the T-matrix method [35]. The latter, in particular, has proven effective for aggregates of spherical monomers and is widely used for interpreting observational data. However, its computational cost becomes prohibitive for large aggregates due to the complexity of matrix inversion operations [36].

To overcome these limitations, a Mean-Field T-matrix (MFTM) method was proposed [37, 38] as an asymptotic, statistical approach rooted in mean-field theory. The mean-field framework was first developed in statistical physics to address the problem of a large number of objects interacting together. As Kadanoff wrote in a famous review on the general mean-field approach [39]: “More is the same”, or, less rhetorically: “use of effective fields can be quite accu- rate, at least in the classical applications involving large pieces of matter”. In the context of light scattering by large fractal aggregates, this approach allows the derivation of analytical or semi-analytical expressions for optical cross- sections. While previous MFTM implementations addressed aggregates of identical monomers, realistic applications require extending the method to polydisperse and compositionally heterogeneous aggregates.

This paper presents a refined version of the MFTM method, which accommodates spherical monomers of varying sizes and materials, and incorporates realistic spatial distributions via a pair correlation function introduced in [40, 41]. These enhancements enable efficient yet accurate modeling of large, irregular aggregates beyond the practical limits of traditional T-matrix computations, making the refined MFTM a valuable complementary tool.

The remainder of this article is structured as follows.

Section 2 presents the governing equations of the refined MFTM, with some technical derivations deferred to Appendix A for clarity.

Section 3 applies the *refined* method to fractal aggregates with monomer sizes distributed according to a Schulz distribution. Indeed, this realistic radius distribution allows for analytical simplifications.

Section 4 outlines the general procedure for computing optical properties using the refined MFTM.

Section 5 illustrates its application through example calculations and comparisons with standard T-matrix results.

Conclusions and outlook are summarized in Section 6.

Additional mathematical details are provided in the appendices.

## 2. MFTM equations for aggregates of polydisperse spherical monomers

### 2.1. *Description of the system*

The particle is an aggregate of *N* spherical homogeneous monomers. It is large (i.e. *N* ≫ 1), rigid and statistically isotropic (orientation-averaged).

Monomers are labelled according to their distance to the centre of mass of the particle (sphere 1 is closest to the centre of mass, while sphere *N* is furthest away). Each monomer, say monomer *j*, is characterized by its radius, *s*<sub>j</sub>, and the refractive index, *m*<sub>j</sub>, of the material of which this monomer is made.

For the sake of simplicity, the centre of the sphere 1 is selected to be the origin of the Cartesian coordinate frame of the aggregate. The Euclidean vector distance between centre of sphere *l* and centre of sphere *j* is noted: $\mathbf{r}_{1l} \equiv \mathbf{r}_1 - \mathbf{r}_j$, where $\mathbf{r}_1$ is the vector distance between centre of sphere *l* and the origin of the frame (i.e. $\mathbf{r}_1$ is short notation for $\mathbf{r}_{1l}$).
The external incident field onto the aggregate is assumed to be a monochromatic plane wave of wavelength $\lambda$. The amplitude of the electromagnetic field is $E_0 = 1$. The time factor $\exp(-i\omega t)$ is assumed and is suppressed throughout the paper [42].

The incident wave vector defines the z-axis of the coordinates frame, that is: $\mathbf{k}_{\text{inc}} = k \mathbf{u}_z$, with $k = 2\pi/\lambda$ the free-space wavenumber. In the same way polarization direction of the incident wave can be selected to define the direction $\mathbf{u}_x$ since the particle is statistically isotropic. This choice is convenient for it makes the expansion of the incident field onto the vector spherical harmonics simple. A linear polarization angle $\gamma \neq 0$ of the incident plane wave with respect to the x-z plane would lead to slightly more complicated incident field coefficients (4) (see for example [43]) resulting in additive $\gamma$-dependent correction to the phase $\phi$ in the amplitude matrix elements (14)-(15) of the wave. After orientation-averaging over the azimuthal angle $\phi$ – as required for statistically isotropic aggregates –, the final MFTM equations do not depend on the value of $\gamma$, then the choice $\gamma = 0$ is sufficient to obtain the general results for isotropic particles.

In this aggregate reference frame, the scattering wave vector writes: $\mathbf{k}_{\text{sca}} = k(\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)$.

## 2.2. Basics of T-matrix method for an aggregate of N monomers

Within the Transition Matrix (T-matrix) method, the incident electromagnetic field and the field scattered by sphere j are written on the basis of vector spherical harmonics, $\mathbf{N}_{mn}^{(3)}(\mathbf{r}_j)$, $\mathbf{M}_{mn}^{(3)}(\mathbf{r}_j)$ evaluated in the spherical coordinate system referenced to sphere j (see for example [44] and references therein):

$$ \mathbf{E}_{\text{sca}}^{j} = \sum_{n=1}^{\infty} \sum_{m=-n}^{n} [a_{mn}^{j} \mathbf{N}_{mn}^{(3)}(\mathbf{r}_j) + b_{mn}^{j} \mathbf{M}_{mn}^{(3)}(\mathbf{r}_j)] . \quad (1) $$

As usual, upper index (3) (resp. (1)) means that spherical Hankel function $h_n^{(1)}$ of the first kind (resp. Bessel function $j_n$) is used in the spherical harmonics). In (1), $a_{mn}^j$ and $b_{mn}^j$ denote the multipole expansion coefficients of order n and degree m of monomer j. Because of the eventual multiple scattering between monomers in the aggregate, these coefficients are different from the corresponding multipole Mie scattering coefficients $\bar{a}_n^j$, $\bar{b}_n^j$ of the isolated monomer j [42]. Note that diverging Mie scattering coefficients $\bar{a}_n^j$ or $\bar{b}_n^j$ are not physical, even close to a Fröhlich resonance [45].

The T-matrix relates the expansion coefficients of the incident field to the expansion coefficients of the scattered field. This is realized through boundary conditions. Binding the incident and scattered electromagnetic fields at the sphere surfaces give a set of inhomogeneous linear equation involving all the $a_{mn}^j$ and $b_{mn}^j$ coefficients, namely:

$$ a_{mn}^{j} + \bar{a}_{n}^{j} \sum_{\substack{l=1 \\ l \neq j}}^{N} \sum_{\substack{\gamma=1 \\ \mu=-\gamma}}^{\infty} \sum_{\substack{l'=1 \\ l \neq j}}^{+\gamma} (A_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj}) a_{\mu\gamma}^{l} + B_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj}) b_{\mu\gamma}^{l}) = \bar{a}_{n}^{j} p_{mn}^{j}, \quad (2) $$

$$ b_{mn}^{j} + \bar{b}_{n}^{j} \sum_{\substack{l=1 \\ l \neq j}}^{N} \sum_{\substack{\gamma=1 \\ \mu=-\gamma}}^{\infty} \sum_{\substack{l'=1 \\ l \neq j}}^{+\gamma} (B_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj}) a_{\mu\gamma}^{l} + A_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj}) b_{\mu\gamma}^{l}) = \bar{b}_{n}^{j} q_{mn}^{j}, \quad (3) $$

where $p_{mn}^j, q_{mn}^j$ are the multipole expansion coefficients of the incident electromagnetic field expressed in the coordinate system centred at the sphere j. By definition of the global coordinate system, the field $\mathbf{E}_{\text{inc}}$ incident on the aggregate is propagating in the z direction, and electric polarization is aligned along the x direction. Then, the multipole coefficients of the incident plane wave are [46, 47, 43, 48]:

$$ p_{1n}^{j} = i^{n+1} \frac{2n+1}{2n(n+1)} e^{i \mathbf{k}_{\text{inc}} \cdot \mathbf{r}_j} ; \quad p_{-1n}^{j} = -i^{n+1} \frac{2n+1}{2} e^{i \mathbf{k}_{\text{inc}} \cdot \mathbf{r}_j}, \quad (4) $$

$$ q_{1n}^{j} = p_{1n}^{j} ; \quad q_{-1n}^{j} = -p_{-1n}^{j}, \quad (5) $$

$$ p_{mn}^{j} = q_{mn}^{j} = 0 \quad \text{if } m \neq \pm 1. \quad (6) $$

In (2)-(3), the addition coefficients $A_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj})$ and $B_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj})$ (discussed in Appendix A) come from the expansion of the vector spherical harmonics $\mathbf{N}_{mn}^{(3)}(\mathbf{r}_j)$ and $\mathbf{M}_{mn}^{(3)}(\mathbf{r}_j)$ onto the basis of the vector spherical harmonics $\mathbf{N}_{\mu\nu}^{(1)}(\mathbf{r}_l)$ and $\mathbf{M}_{\mu\nu}^{(1)}(\mathbf{r}_l)$ evaluated in the spherical coordinate system referenced to sphere l [49].
### 2.2.1. General scattered field

The total electric field scattered by the aggregate is obtained after summing up the electric fields $E_{sca}^j$ scattered by all the monomers:

$$E_{sca} = \sum_{j=1}^{N} E_{sca}^{j}. \qquad (7)$$

Translating the coordinate system referenced to sphere j into the one referenced to sphere 1, it writes [48, 50]:

$$E_{sca} = \sum_{n=1}^{\infty} \sum_{m=-n}^{n} [a_{mn}N_{mn}^{(3)}(\mathbf{r}_1) + b_{mn}M_{mn}^{(3)}(\mathbf{r}_1)], \qquad (8)$$

with

$$a_{mn} = a_{mn}^{1} + \sum_{j=2}^{N} \sum_{v=1}^{\infty} \sum_{\mu=-v}^{v} [A_{\mu\nu,mn}^{(1)}(\mathbf{r}_{j1}) a_{\mu\nu}^{j} + B_{\mu\nu,mn}^{(1)}(\mathbf{r}_{j1}) b_{\mu\nu}^{j}], \qquad (9)$$

$$b_{mn} = b_{mn}^{1} + \sum_{j=2}^{N} \sum_{v=1}^{\infty} \sum_{\mu=-v}^{v} [B_{\mu\nu,mn}^{(1)}(\mathbf{r}_{j1}) a_{\mu\nu}^{j} + A_{\mu\nu,mn}^{(1)}(\mathbf{r}_{j1}) b_{\mu\nu}^{j}]. \qquad (10)$$

### 2.2.2. Far-field scattered field

Let us consider a point S in space at a distance $r_S = r(\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)$, far from the aggregate. Vector $r_S = r\mathbf{u}_{sca}$ defines a scattering direction, $\mathbf{u}_{sca}$, of the electromagnetic wave. Using asymptotic behaviour of the spherical Hankel function $h_n^{(1)}$ at large arguments, the far-field total electric field (7) scattered by the entire aggregate at the location S, is transverse and the $\theta$ and $\phi$ components are:

$$E_{sca,\theta} \sim \sum_{j=1}^{N} \frac{e^{ikr_{JS}}}{-ikr_{JS}} \sum_{n=1}^{\infty} \sum_{m=-n}^{n} (-i)^{n+1} [a_{mn}^{j} \tau_{mn}(\theta_{JS}) + b_{mn}^{j} \pi_{mn}(\theta_{JS})] e^{im\phi_{JS}}, \qquad (11)$$

$$E_{sca,\phi} \sim \sum_{j=1}^{N} \frac{e^{ikr_{JS}}}{-ikr_{JS}} \sum_{n=1}^{\infty} \sum_{m=-n}^{n} (-i)^n [a_{mn}^{j} \pi_{mn}(\theta_{JS}) + b_{mn}^{j} \tau_{mn}(\theta_{JS})] e^{im\phi_{JS}}, \qquad (12)$$

where the angular functions $\tau_{mn}(\theta)$ and $\phi_{mn}(\theta)$ are:

$$\pi_{mn}(\theta) = \frac{m P_{n}^{m}(\cos\theta)}{\sin\theta}; \quad \tau_{mn}(\theta) = \frac{d}{d\theta} P_{n}^{m}(\cos\theta), \qquad (13)$$

and the angular variables are defined by: $\mathbf{r}_{JS} = r_{JS}(\sin\theta_{JS}\cos\phi_{JS}, \sin\theta_{JS}\sin\phi_{JS}, \cos\theta_{JS})$.

Writing $\mathbf{r}_{JS} = \mathbf{r}_S - \mathbf{r}_j$, and $r \equiv r_S \gg r_j$, one has: $r_{JS} = r - \mathbf{u}_{sca} \cdot \mathbf{r}_j + O(r_j^2/r)$ and $\theta_{JS} = \theta_S + O(r_j/r)$, $\phi_{JS} = \phi_S + O(r_j/r)$. In the far-field situation, the terms $O(r_j/r)$ can be neglected in front of $O(1)$, then the amplitude matrix elements of the wave are:

$$E_{sca,\theta} \sim \frac{e^{ikr}}{-ikr} \sum_{n=1}^{\infty} \sum_{m=-n}^{n} (-i)^{n+1} [a_{mn}\tau_{mn}(\theta) + b_{mn}\pi_{mn}(\theta)] e^{im\phi}, \qquad (14)$$

$$E_{sca,\phi} \sim \frac{e^{ikr}}{-ikr} \sum_{n=1}^{\infty} \sum_{m=-n}^{n} (-i)^n [a_{mn}\pi_{mn}(\theta) + b_{mn}\tau_{mn}(\theta)] e^{im\phi}. \qquad (15)$$

where the total scattering coefficients, $a_{mn}$, $b_{mn}$, are [47, 51, 52]:

$$a_{mn} = \sum_{j=1}^{N} a_{mn}^{j} e^{-i \mathbf{k}_{sca} \cdot \mathbf{r}_{j}}, \quad b_{mn} = \sum_{j=1}^{N} b_{mn}^{j} e^{-i \mathbf{k}_{sca} \cdot \mathbf{r}_{j}}. \qquad (16)$$
### 2.2.3. Computational limitations

T-matrix method is limited to relatively small aggregate sizes. Indeed, the size of the matrix (2)-(3) to be inverted to obtain the multipole scattering coefficients is $4Nn_c^2 \times 4Nn_c^2$, in which the maximum relevant order is [53, 54, 55]:

$$n_c = \lfloor ks_{\max} + 4.05(ks_{\max})^{1/3} \rfloor + 8, \quad (17)$$

where $s_{\max}$ is the largest monomer radius in the aggregate, and $\lfloor \cdots \rfloor$ denotes the integer part. It can result in a huge matrix when the number, $N$, of monomers is large. Even using accelerated algorithm combined with large-scale parallel processing platforms, simulated aggregates of a million monomers, constitute presently the extreme limit for large-scale computers [36].

### 2.3. Mean-Field T-matrix equations

The Mean-Field T-matrix (MFTM) method is an approximation to the T-matrix method capable of rapidly calculating estimates of the optical quantities of a large, orientation-averaged, aggregate of spherical grains when multiple scattering is present [38]. Basically, the MFTM method supposes that every monomer experiences a similar incident electromagnetic field (coming from outside and from all the other monomers of the aggregate). It means in particular that the number of monomers on the boundary of the particle is much smaller than the number of monomers in the volume, i.e. the particle is expected to be large.

It is convenient to define two auxiliary sets of coefficients, namely:

$$\bar{q}_{mn}^j = i^{n+1} \frac{2n+1}{2} \frac{(n-m)!}{(n+1)!} e^{i \mathbf{k}_{\text{inc}} \cdot \mathbf{r}_j} ; \quad \bar{p}_{mn}^j = \zeta_m \bar{q}_{mn}^j, \quad (18)$$

where $\zeta_m$ denotes the sign of $m$ (that is: $\zeta_m = +1$ if $m \ge 0$ and $= -1$ for negative values of $m$). These coefficients generalize the current multipole coefficients (4)-(6) of the incident wave, that is:

$$\frac{p_{mn}^j}{\bar{p}_{mn}^j} = \frac{p_{mn}^j}{\bar{p}_{mn}^j} = 1 \quad \text{for } m = \pm 1, \quad (19)$$

$$= 0 \quad \text{for } m \neq \pm 1. \quad (20)$$

The multipole coefficients are scaled with the coefficients $\bar{p}_{mn}^j$ and $\bar{q}_{mn}^j$, namely:

$$a_{mn}^j = \bar{a}_n^j \bar{p}_{mn}^j \delta_{mn}^{(a)} ; \quad b_{mn}^j = \bar{b}_n^j \bar{q}_{mn}^j \delta_{mn}^{(b)}, \quad (21)$$

defining the coefficients $\delta_{mn}^{(a,b)}$ where upper indices (a) and (b) correspond to the TM and TE modes respectively.

Note that the Rayleigh-Gans-Debye (RGD) theory [56] corresponds to the special case of single scattering: $\delta_{mn}^{(a,b)} = 1$. The MFTM method goes beyond this approximation by calculating the values of the scattering coefficients $\delta_{mn}^{(a,b)}$ from multiple scattering acting within the aggregate. The derivation of these scattering coefficients is described below.

The mean-field approximation in this context is to assert that the coefficients $\delta_{m,n}^{(a,b)}$ are monomer independent [37, 38]. This means that the MFTM method assumes that each monomer is exposed to the same environment.

Expressing (2)-(3) in terms of $\delta_{m,n}^{(a,b)}$, using (18)-(19), and replacing sums $\sum_{l \neq j} \cdots$ over the monomers $l$ with averages $\langle \cdots \rangle$ over the spatial distribution of the monomers and of their characteristics, the set of equations defining the values of $\delta_{m,n}^{(a,b)}$ write [38]:

$$\delta_{m,n}^{(a)} + (N-1) \sum_{v=1}^{\infty} i^{v-n} \frac{2v+1}{2n+1} \frac{(n+1)!}{(v+1)!} \frac{(v-m)!}{(n-m)!} \left( \langle \langle A_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{jl}} \bar{a}_v^l \rangle \rangle \delta_{m,v}^{(a)} + \zeta_m \langle \langle B_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{jl}} \bar{b}_v^l \rangle \rangle \delta_{m,v}^{(b)} \right) = 1, \quad (22)$$

$$\delta_{m,n}^{(b)} + (N-1) \sum_{v=1}^{\infty} i^{v-n} \frac{2v+1}{2n+1} \frac{(n+1)!}{(v+1)!} \frac{(v-m)!}{(n-m)!} \left( \zeta_m \langle \langle B_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{jl}} \bar{a}_v^l \rangle \rangle \delta_{m,v}^{(a)} + \langle \langle A_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{jl}} \bar{b}_v^l \rangle \rangle \delta_{m,v}^{(b)} \right) = 1, \quad (23)$$
where $m$ can take the values $\pm 1$ only, because $\delta_{m,n}^{(a,b)} = 0$ when $m \neq \pm 1$ (this result comes from (102) in Appendix A).

Applying the symmetry relations (106)-(107) (derived in Appendix A) to the equations (22)-(23), the following result ensues:

$$ \delta_{1n}^{(a,b)} = \delta_{-1n}^{(a,b)} \qquad (24) $$

for any value of order $n$. Let us write $\delta_n^{(a,b)}$ the common values of $\delta_{\pm 1n}^{(a,b)}$, that is:

$$ \delta_n^{(a)} = \delta_{1n}^{(a)} ; \quad \delta_n^{(b)} = \delta_{1n}^{(b)} . \qquad (25) $$

These multiple-scattering coefficients are solutions of the set of linear coupled equations:

$$ \delta_n^{(a)} + (N-1) \sum_{v=1}^{\infty} [A_n^v \langle \bar{a}_v \rangle \delta_v^{(a)} + B_n^v \langle \bar{b}_v \rangle \delta_v^{(b)}] = 1, \qquad (26) $$

$$ \delta_n^{(b)} + (N-1) \sum_{v=1}^{\infty} [B_n^v \langle \bar{a}_v \rangle \delta_v^{(a)} + A_n^v \langle \bar{b}_v \rangle \delta_v^{(b)}] = 1. \qquad (27) $$

In (26)-(27), $\langle \bar{a}_v \rangle$ and $\langle \bar{b}_v \rangle$ are the mean values of the Mie coefficients $\bar{a}_v^j$ and $\bar{b}_v^j$ over monomer sizes and refractive indices.

After some algebra on the relations (104)-(105), using (108)-(110) and known identities on the associated Legendre polynomials, the complex-valued coefficients $A_n^v, B_n^v$ appearing in (26)-(27) can be written under the forms:

$$ A_n^v = \frac{2v+1}{2} \sum_{\substack{p=|n-v| \\ p+n+v \text{ even}}}^{n+v} \alpha(n, v; p) S_p ; \quad B_n^v = \frac{2v+1}{2} \sum_{\substack{p=|n-v+1| \\ p+n+v \text{ odd}}}^{n+v-1} \beta(n, v; p) S_p, \qquad (28) $$

in which the positive rational numbers $\alpha(n, v; p)$ and $\beta(n, v; p)$ do not depend on the aggregate, while the complex-valued coefficients $S_p$ depend on the wavelength and the geometrical structure of the aggregate. They are respectively given by:

$$ \alpha(n, v; p) = (2p+1) \frac{n(n+1) + v(v+1) - p(p+1)}{2n(n+1)v(v+1)} \int_{-1}^{+1} P_v^1(x) P_n^1(x) P_p(x) dx, \qquad (29) $$

$$ \beta(n, v; p) = (2p+1) \frac{p+1}{n(n+1)v(v+1)} \int_{-1}^{+1} P_v^1(x) P_n^1(x) \frac{x P_p(x) - P_{p+1}(x)}{1-x^2} dx. \qquad (30) $$

$$ S_p = \langle j_p(k r_{lj}) h_p^{(1)}(k r_{lj}) \rangle, \qquad (31) $$

where, in (31), symbol $\langle \cdots \rangle$ means average over the spatial distribution of the monomers, whatever their size and refractive index. The coefficients $\alpha(n, v; p), \beta(n, v; p)$ can also be expressed in terms of 3-j symbols, for example [49, 47]:

$$ \alpha(n, v; p) = -(2p+1) \frac{n(n+1) + v(v+1) - p(p+1)}{\sqrt{n(n+1)v(v+1)}} \begin{pmatrix} v & n & p \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} v & n & p \\ 1 & -1 & 0 \end{pmatrix}, \qquad (32) $$

and it is possible – though more complicated – to write the coefficient $\beta(n, v; p)$ in a similar form [57]. The identities:

$$ \sum_{p=|n-v|}^{n+v} \alpha(n, v; p) = 1 ; \quad \sum_{p=|n-v|}^{n+v} \beta(n, v; p) = 1, \qquad (33) $$

and: $\alpha(n, v; p) = 0$ when $n + v + p$ is an even number; $\beta(n, v; p) = 0$ when $n + v + p$ is an odd number, can be used as checks. There are also symmetries: $\alpha(n, v; p) = \alpha(v, n; p)$ and $\beta(n, v; p) = \beta(v, n; p)$.
### 2.3.1. Three notes about the matrix MFTM equation

An important fact should be mentioned here: the MFTM equations (26)-(27) constitute in principle an infinite set of coupled linear equations. However, the orders $n$ and $v$ are in practice restricted to the range $1 \le n, v \le n_c$, with $n_c$ given in (17), because the contribution of all scattering coefficients $\bar{a}_n$ and $\bar{b}_n$ with $n > n_c$ is negligibly small. Thus, the effective mean-field matrix (26)-(27) to be inverted is concretely of size $2n_c \times 2n_c$ that is correlated to the size parameter of the largest monomer in the aggregate. This is a dramatic simplification over the regular T-matrix method which requires inverting a $4n_c^2 N \times 4n_c^2 N$ matrix to find the relevant $a_{m,n}^j, b_{m,n}^j$ coefficients.

General mathematical conditions for the infinite set of linear equations (26)-(27) to be invertible are not known. Only a few partial results are available. An example is discussed below in Section 5.1: for large fractal aggregates of fractal dimension $d_f > 2$, a sufficient condition for the existence of the inverse of the MFTM equations (26)-(27) is written.

It is worth noting the special decoupling of the average quantities in (26)-(27) due to the no size- or material-segregation assumptions (this point is discussed in [58], and shown on Fig. 1 below): the average addition coefficients $(N-1)A_n^v$ and $(N-1)B_n^v$ depend on the morphology of the aggregate only, while $\langle \bar{a}_v \rangle$ and $\langle \bar{b}_v \rangle$ depend on the monomer-radius and refractive-index distributions.

### 2.3.2. $S_p$ coefficients

When the aggregate is large, the quantities to be averaged over all the couples $(l, j)$ of monomers whatever their respective radius can be expressed in terms of the pair correlation function $g(r)$. This radial function is such that: $4\pi r^2 g(r)$ is the probability to find a pair of monomers (no matter the sizes) at distance $r$. It is normalized according to:

$$ \int_0^\infty 4\pi r^2 g(r) dr = 1. \qquad (34) $$

The general expression of $S_p$, as written in (31), is then given by:

$$ S_p = \int_0^\infty j_p(k r) h_p^{(1)}(k r) 4\pi r^2 g(r) dr. \qquad (35) $$

The equations (24)-(30) and (35) solve completely the problem of light scattering by an ensemble of orientation-averaged aggregates of $N$ monomers of random radii and refractive indices, within the MFTM method, when the pair-correlation function is known.

## 2.4. Mean-Field optical quantities

From the far-field amplitudes (11)-(12), the mean-field amplitude scattering matrix $\{S_1, S_2, S_3, S_4\}$ takes the simple form:

$$ S_1 = \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} \left( \alpha_n \delta_n^{(a)} \pi_{1n} + \beta_n \delta_n^{(b)} \tau_{1n} \right), \quad S_2 = \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} \left( \alpha_n \delta_n^{(a)} \tau_{1n} + \beta_n \delta_n^{(b)} \pi_{1n} \right), \qquad (36) $$

$$ S_3 = S_4 = 0, \qquad (37) $$

with the sums:

$$ \alpha_n = \sum_{j=1}^{N} \bar{a}_n^j e^{i(k_{\text{inc}} - k_{\text{sca}}) \cdot \mathbf{r}_j} ; \quad \beta_n = \sum_{j=1}^{N} \bar{b}_n^j e^{i(k_{\text{inc}} - k_{\text{sca}}) \cdot \mathbf{r}_j}. \qquad (38) $$
### 2.4.1. Linear polarization profile

Putting aside the $N$ terms $j = l$ out of the $N(N-1)$ terms $j \neq l$, and using the no size-segregation assumption, the average values (over positions and characteristics of the monomers) of $|S_1|^2$ and $|S_2|^2$ are:

$$ \langle |S_1|^2 \rangle = N \langle |\bar{S}_1(s, \theta)|^2 \rangle + (N-1) S(\theta) \langle |\bar{S}_1(s, \theta)|^2 \rangle, \quad (39) $$

$$ \langle |S_2|^2 \rangle = N \langle |\bar{S}_2(s, \theta)|^2 \rangle + (N-1) S(\theta) \langle |\bar{S}_2(s, \theta)|^2 \rangle. \quad (40) $$

In the relations (39)-(40), the functions $\bar{S}_1(s, \theta)$ and $\bar{S}_2(s, \theta)$ are the amplitude scattering elements of a single spherical particle of radius $s$ and corrected Mie scattering coefficients $\{\delta_n^{(a)} \bar{a}_n(s), \delta_n^{(b)} \bar{b}_n(s)\}$, that is:

$$ \bar{S}_1(s, \theta) = \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} (\bar{a}_n(s) \delta_n^{(a)} \pi_{1n} + \bar{b}_n(s) \delta_n^{(b)} \tau_{1n}), \quad \bar{S}_2(s, \theta) = \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} (\bar{a}_n(s) \delta_n^{(a)} \tau_{1n} + \bar{b}_n(s) \delta_n^{(b)} \pi_{1n}). \quad (41) $$

Also, the real-valued $\theta$-depending function $S(\theta) = \langle e^{i(k_{\text{inc}}-k_{\text{sca}}) \cdot r_{ij}} \rangle$ – average over the distances $r_{ij}$ and $l \neq j$ – is essentially the Fourier-transform of the pair-correlation function. After structure and orientation averaging, it is written [38]:

$$ S(\theta) = \int_{0}^{\infty} j_{0}(2kr \sin(\theta/2)) 4\pi r^{2} g(r) dr. \quad (42) $$

Note the inequality:

$$ |S(\theta)| \le 1. \quad (43) $$

The size distribution average quantities appearing in (39)-(40) can be written explicitly as:

$$ \langle |\bar{S}_1(s, \theta)|^2 \rangle = \sum_{n=1}^{\infty} \sum_{n'=1}^{\infty} \frac{2n+1}{n(n+1)} \frac{2n'+1}{n'(n'+1)} \left[ \langle \bar{a}_n \bar{a}_{n'}^* \rangle \delta_n^{(a)} \delta_{n'}^{(a)*} \pi_{1n} \pi_{1n'} + \langle \bar{b}_n \bar{b}_{n'}^* \rangle \delta_n^{(b)} \delta_{n'}^{(b)*} \tau_{1n} \tau_{1n'} + \right. \\ \left. + \langle \bar{a}_n \bar{b}_{n'}^* \rangle \delta_n^{(a)} \delta_{n'}^{(b)*} + \langle \bar{a}_{n'} \bar{b}_n \rangle \delta_{n'}^{(a)*} \delta_n^{(b)} \pi_{1n} \tau_{1n'} \right], \quad (44) $$

$$ |\langle \bar{S}_1(s, \theta) \rangle|^2 = \left| \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} \left( \langle \bar{a}_n \rangle \delta_n^{(a)} \pi_{1n} + \langle \bar{b}_n \rangle \delta_n^{(b)} \tau_{1n} \right) \right|^2, \quad (45) $$

and similar expressions for $S_2$ after replacement $\pi_{1n} \leftrightarrow \tau_{1n}$. For $\theta = 0$, the results: $\pi_{1n} = \tau_{1n} = -n(n+1)/2$ and $S(0) = 1$ can be used.

Note that prior calculation of the average quantities $\langle \bar{a}_n \bar{a}_{n'}^* \rangle$, $\langle \bar{b}_n \bar{b}_{n'}^* \rangle$, $\langle \bar{a}_n \bar{b}_{n'}^* \rangle$ (for $1 \le n, n' \le n_c$) are required before estimating values of (44).

From expressions (39)-(40) and (44)-(45), the average degree of linear polarization, $P$, of the scattered wave is:

$$ P = \frac{\langle |S_1|^2 \rangle - \langle |S_2|^2 \rangle}{\langle |S_1|^2 \rangle + \langle |S_2|^2 \rangle}. \quad (46) $$

### 2.4.2. Scattering cross-section, asymmetry parameter and phase function

Using (from the definitions (13)) the expansions of $\pi_{1n}$ and $\tau_{1n}$ in terms of Legendre polynomials:

$$ \pi_{1n}(\theta) = \frac{n(n+1)}{2n+1} \frac{1}{\sin^2 \theta} (P_{n+1}(\cos \theta) - P_{n-1}(\cos \theta)), \quad (47) $$

$$ \tau_{1n}(\theta) = - \frac{n(n+1)}{2n+1} \frac{1}{\sin^2 \theta} \left( (P_{n+1}(\cos \theta) - P_{n-1}(\cos \theta)) \cos \theta + (2n+1)P_n(\cos \theta) \sin^2 \theta \right), \quad (48) $$
At last, the scattering phase function (with normalization: $\int I(\theta) d\Omega = 1$) is:

$$I(\theta) = \frac{S_{11}(\theta)}{k^2 C_{\text{sca,N}}}. \qquad (59)$$

When the complete phase function $I(\theta)$ is required for a study, expression (59) can be used to calculate by integration the scattering cross-section $C_{\text{sca,N}}$ – instead of performing (or for checking) the calculations (53)-(54) – and the asymmetry parameter $g$.

### 2.4.3. Extinction cross-sections

The extinction cross-section is obtained from the sum of the extinction cross-sections of the individual monomers. Alternatively, the mean-field solution (36)-(38) corresponds to the field scattered by a sphere with modified Mie coefficients $\{\delta_n^{(a)} \alpha_n, \delta_n^{(b)} \beta_n\}$. The aggregate cross-section extinction for aggregate of size $N$ is then given by the optical theorem (setting $\mathbf{k}_{\text{sca}} = \mathbf{k}_{\text{inc}}$ in (38)):

$$C_{\text{ext,N}} = \frac{2\pi N}{k^2} \sum_{n=1}^{\infty} (2n + 1) \text{Re}\{\delta_n^{(a)} \langle a_n \rangle + \delta_n^{(b)} \langle b_n \rangle\}. \qquad (60)$$

### 2.4.4. Absorption cross-sections

The absorption cross-section of a given monomer is obtained from calculating the net rate at which electromagnetic energy crosses the surface of the monomer [42]. The T-matrix calculation was done in [48, 43]. Replacing the derivatives of Ricatti-Bessel functions with spherical Bessel and Hankel functions in the formulae given in these papers – and after a bit of algebra –, the MFTM approximation of the T-matrix result gives the absorption cross-section of an aggregate of size $N$:

$$C_{\text{abs,N}} = \frac{2\pi N}{k^2} \sum_{n=1}^{\infty} (2n + 1) \left( c_n^{(1)} |\delta_n^{(a)}|^2 + c_n^{(2)} |\delta_n^{(b)}|^2 \right), \qquad (61)$$

$$c_n^{(1)} = \left\langle \frac{\text{Im}\{m^2\} (n+1) |j_n(mx)|^2 - |m|^2 x \text{Im}\{m j_n(mx) j_{n+1}^*(mx)\}}{x |(n+1)(m^2-1) j_n(mx) h_n^{(1)}(x) + m x (h_n^{(1)}(x) j_{n+1}(mx) - m h_{n+1}^{(1)}(x) j_n(mx))|^2} \right\rangle, \qquad (62)$$

$$c_n^{(2)} = \left\langle \frac{-\text{Im}\{m^* j_n(mx) j_{n+1}^*(mx)\}}{x^2 |m h_n^{(1)}(x) j_{n+1}(mx) - h_{n+1}^{(1)}(x) j_n(mx)|^2} \right\rangle. \qquad (63)$$

The parameter $x \equiv k_s$ in (62)-(63) denotes the size parameter of the monomer of radius $s$, and $m$ is the refractive index of the monomer. Symbol $\langle \cdots \rangle$ denotes average over the monomer size and composition distributions. When $m$ is a pure real number (i.e. non-absorbing material), it is clear from (62)-(63) that $c_n^{(1)} = c_n^{(2)} = 0$, hence $C_{\text{abs,N}} = 0$, as it must be.

Note that the expected relation: $C_{\text{ext,N}} = C_{\text{sca,N}} + C_{\text{abs,N}}$ does not hold exactly from the results (60), (53) and (61) because MFTM is an approximation of the rigorous theory. In this situation, (61) is more accurate for estimating $C_{\text{abs,N}}$ than using the MFTM value of $C_{\text{ext,N}} - C_{\text{sca,N}}$.

## 3. MFTM equations for large fractal aggregates made of Schulz-distributed monomer radii and uniform material

In this Section, the values of the monomer radii, $s$, are distributed according to the normalized distribution function $f(s)$, and all the monomers are made of a same material of refractive index $m$. The mean monomer radius $s_0$ and the
polydispersity, $\sigma$, of the monomer radius are defined through the moments of the radius distribution:

$$ \langle s^n \rangle = \int_{0}^{\infty} s^n f(s) ds \qquad (64) $$

$$ s_0 = \langle s \rangle ; \quad \sigma = \frac{\sqrt{\langle s^2 \rangle - \langle s \rangle^2}}{\langle s \rangle} . \qquad (65) $$

Let us consider the particular case of the two-parameters ($z > 0, s_0 > 0$) Schulz distribution [59]:

$$ f(s) = \frac{1}{\Gamma(z)} \left(\frac{z}{s_0}\right)^z s^{z-1} e^{-zs/s_0}, \qquad (66) $$

defined over the range $s > 0$. Polydispersity of that distribution is: $\sigma = 1/\sqrt{z}$. The curve shape representing (66) is unimodal when $z > 1$. The larger $z$, the closer the Schulz distribution is to the Gaussian distribution $N(s_0, s_0^2/z)$.

An example of fractal cluster generated by the regular Ballistic Cluster-Cluster Aggregation algorithm, is shown on Fig. 1. The fractal dimension of the aggregate is $d_f = 1.9$. The spherical monomers are distributed following the Schulz distribution (66) with polydispersity $\sigma = 0.5$ (that is: $z = 4$) and mean radius $s_0 = 0.5$. No size segregation (i.e. no blob of monomers of approximately the same radii) is visible in the figure. This point is important to check or to argue before using the MFTM method, because it is a basic condition for the method be reliable.

In the simulated aggregates, simple check of size segregation can be done partitioning the actual monomer radius distribution into bins, and calculating the mean number, $n_i(j)$ of monomers belonging to bin $j$ and tangent to a bin-$i$ monomer. Segregation of the monomers from (e.g.) the bin $i$ results in a strong increase in the value of $n_i(i)$. For regular BCCA aggregates, such a check is shown in Fig. 2. The function $n_i(j)$ appears to be almost constant, in accordance to the no-segregation assumption. This insures that the MFTM method can be used dependably in this case. Limited bias in the data of $n_i(j)$ between small and large values of $(i, j)$ in Fig. 2 is irrelevant and due to the differences in the respective collision cross sections of the monomers [60].

Note that possible monomer-size segregation in Cluster-Cluster Aggregation processes has not been studied in detail so far. This is because the polydispersity of monomers is not expected to radically alter the inherent randomness of fractal aggregates under usual physical conditions. Similar problems arise in closely related contexts, such as the evidence that the polydispersity of spheres does not prevent hyperuniform organization of sphere packings [61].

For comparison and illustration, a modified BCCA process in which size segregation is forced for monomers belonging to a particular bin, is discussed in Appendix D.

### 3.1. Largest monomer in an aggregate

Let $s_{max}$ be the ensemble average radius of the largest monomer in an aggregate of size $N$. It is given by the condition [62]:

$$ N \int_{s_{\max}}^{\infty} f(s) ds = 1. \qquad (67) $$

Using the size-distribution (66), the integral in (67) is the S-shaped function: $\Gamma(z, zs_{max}/s_0)/\Gamma(z)$, with $\Gamma$ the (incomplete or complete) Gamma function.

When $z$ and $N$ are both large, and $\log N \gg z$, the approximation [63] $\Gamma(z, zx)/\Gamma(z) \approx (x e^{-x} z)^z/(z(x-1))$ leads to: $s_{max} \approx s_0(1 + \log N/z)$. Hence the estimation of the value of the maximum multipole order:

$$ n_c \approx \frac{ks_0}{z} \log N, \qquad (68) $$

increasing logarithmically with the aggregate size $N$.

Relation (68) gives information on how the mean-field translation matrix increases with the aggregate size, but in the applications, it is more precise to calculate the positive solution of:

$$ \left( z \frac{s_{\max}}{s_0} e^{-s_{\max}/s_0} \right)^z = \frac{z}{N} \left( \frac{s_{\max}}{s_0} - 1 \right), \qquad (69) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. Example of 3D BCCA fractal aggregate ($d_f = 1.9$) made of 16 384 Schulz-distributed monomers (polydispersity: 0.5). Colours correspond to values of the monomer radius (16 colours are used; from blue = small monomers to red = large monomers). The different monomer sizes look randomly distributed throughout the aggregate, that is no segregation of the monomer sizes is apparent. This result is proved quantitatively on Fig. 2.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Mean number, $n_i(j)$, of monomers belonging to bin $j$ and tangent to bin-$i$ monomers. Average over 100 independent 3D BCCA aggregates of size $N = 16\,384$. Monomer radii are Schulz-distributed (polydispersity: 0.5). The complete monomer radii range is divided into 16 bins, each bin consisting of 1024 monomers of about the same radius. Any spatial segregation of the monomers belonging to bin $i$ would result in a peak in $n_i(i)$ characteristic of an excess of monomers of bin $i$ tangent to other monomers in the same bin $i$. In the figure, there is no appearance of irregularity of the $n_i(j)$ function, i.e. no monomer-size segregation occurs in the regular BCCA process. This figure can be compared to Fig. 8 of Appendix D (using the same range of values) that corresponds to an aggregation process with forced size segregation.</figcaption></figure>
in the variable $\langle s_{\text{max}} \rangle$, and use (17) to write:

$$n_c = [ks_{\text{max}} + 4.05(ks_{\text{max}})^{1/3}] + 8, \quad (70)$$

to obtain the value of the largest order $n$ relevant to the current situation.

## 3.2. Pair correlation function of the fractal aggregate

For aggregates of fractal dimension $d_f < 3$ made of monomers following the Schulz radius distribution (66), the global pair correlation function writes (see [41] and details in Appendix B):

$$4\pi r^2 g(r) = \frac{1}{N} (4\pi r^2 g_L(r) + 4\pi r^2 g_R(r)), \quad (71)$$

$$4\pi r^2 g_L(r) = \frac{z^{2z}}{\Gamma(2z)} \frac{2}{s_0} \left(\frac{r}{s_0}\right)^{2z-1} e^{-zr/s_0}, \quad (72)$$

$$4\pi r^2 g_R(r) = \frac{\kappa d_f k_f}{s_0} \left(\frac{r}{s_0}\right)^{d_f-1} e^{-\kappa(r/R_g)^{d_f}} \left(1 - \frac{\Gamma(2z, zr/s_0)}{\Gamma(2z)}\right). \quad (73)$$

In (73), the coefficient $\kappa$ is a function of the only fractal dimension: $\kappa = (\Gamma(2/d_f)/d_f)^{d_f/2}$ for normalization of $g(r)$. The mean gyration radius, $R_g$, of the fractal aggregate is given by the relation:

$$N = k_f \left(\frac{R_g}{s_0}\right)^{d_f}, \quad (74)$$

with $k_f$ the fractal prefactor. In most of the applications, the value of $k_f$ is close to 1 [38]. Note however that the value of the prefactor is known to decrease slowly with monomer polydispersity [64].

In (72), the function $g_L(r)$ results from the local arrangement of the monomers, that is $g_L(r)$ is independent of the fractal dimension $d_f$. On the other hand, the regular part $g_R(r)$ (in (73)) depends essentially on the global (fractal) morphology of the particle.

## 3.3. Coefficients $S_p$ for aggregates of Schulz-polydispersed monomers

The coefficients $S_p$ are given in (31) and, more explicitly, in (35). Using the analytical form (71)-(73) of the pair correlation function for the Schulz-distributed monomer radii, the coefficients $S_p$ write:

$$S_p = \frac{1}{N} \left( S_p^{(L)}(ks_0) + S_p^{(R)}(ks_0, R_g/s_0) \right), \quad (75)$$

with $p$ a non-negative integer number. In (75), the first component, $S_p^{(L)}$, is:

$$S_p^{(L)}(ks_0) = 2 \frac{z^{2z}}{\Gamma(2z)} \int_0^\infty j_p(ks_0 u) h_p^{(1)}(ks_0 u) u^{2z-1} e^{-zu} du. \quad (76)$$

The integral in (76) can be expressed as the sum of two hypergeometric functions ${}_3F_2$:

$$S_p^{(L)}(x) = \frac{\Gamma(2z + 2p)}{\Gamma(2z)} \frac{p!^2}{(2p + 1)!^2} \left(\frac{2x}{z}\right)^{2p} {}_3F_2\left[p+z, p+z+1/2, p+1; -\left(\frac{2x}{z}\right)^2\right] - \\ - \frac{i}{(2z-1)(2p+1)} \frac{z}{x} {}_3F_2\left[z-1/2, z, 1/2; p+3/2, 1/2-p; -\left(\frac{2x}{z}\right)^2\right]. \quad (77)$$
The second component in (75) is:

$$S_p^{(R)}(ks_0, R_g/s_0) = \kappa d_f k_f \int_0^\infty j_p(ks_0 u) h_p^{(1)}(ks_0 u) u^{d_f-1} e^{-\kappa(\kappa u_{s_0}/R_g)^{d_f}} \left(1 - \frac{\Gamma(2z, zu)}{\Gamma(2z)}\right) du, \quad (78)$$

which can be written – using the Heaviside step function $1 - H(u/2)$ as an approximation of $\Gamma(2z, zu)/\Gamma(2z)$ for large $z$ – under the simple form:

$$S_p^{(R)}(ks_0, R_g/s_0) \approx \kappa d_f k_f \int_2^\infty j_p(ks_0 u) h_p^{(1)}(ks_0 u) u^{d_f-1} e^{-\kappa(\kappa u_{s_0}/R_g)^{d_f}} du, \quad (79)$$

convenient for numerical integration (accuracy of the Heaviside function approximation is shown on Fig. 5 of Appendix B).

From inequality (C.2) (in Appendix C) on the product function $j_p(x) h_p^{(1)}(x)$, one has:

$$|S_p^{(L)}(ks_0)| < \frac{z}{z-1/2} \sqrt{\frac{\pi}{2p+1}} \frac{1}{ks_0}, \quad (80)$$

$$|S_p^{(R)}(ks_0, R_g/s_0)| < (\kappa k_f)^{1/d_f} \Gamma(1 - 1/d_f) \sqrt{\frac{\pi}{2p+1}} \frac{N^{1-1/d_f}}{ks_0}. \quad (81)$$

These two inequalities can be used as consistency checks of numerical values obtained for $S_p^{(L)}$ and $S_p^{(R)}$.

## 4. General MFTM protocol for large aggregates of polydisperse, non-uniform monomers

The main intention of this paper is to provide a practical tool for estimating the optical properties of large dust particles in the form of aggregates of disparate spherical monomers. The general MFTM protocol is described below for the Schulz-distribution case, with recommendations for practical computation. Application to other non-Schulz distribution of the monomer sizes is straightforward in principle (when the size distribution is not singular). The aggregates are assumed to be large – typically larger than a few thousand monomers –, and can be fractal or porous non-fractal. The calculation sequence is as follows (checks are optional):

* 1) estimation of the largest monomer size, $s_{max}$, using equation (69). The Newton-Raphson method is efficient to solve that equation (check: $s_{max} \ge s_0$).

* 2) estimation of the largest multipole order, $n_c$, using (17) with the value of $s_{max}$ found in the previous step.

* 3) calculation of the $S_p$ coefficients for $p = 0, \dots, 2n_c$ using (75), (77), and (78) or (79). The pair-distribution function $g(r)$ decreasing faster than exponentially at large distance, Gauss-Laguerre quadrature is numerically stable to calculate these coefficients (check: (80)-(81)).

* 4) writing the $A_n^v$, $B_n^v$ coefficients for $n, v = 1, \dots, n_c$ using (28)-(30) and the values of $S_p$ found in the previous step (check: (33) on the coefficients $\alpha(n, v; p)$ and $\beta(n, v; p)$).

* 5) calculation of the coefficients $\langle a_n \rangle$, $\langle b_n \rangle$, $\langle a_n a_{n'}^* \rangle$, $\langle a_n b_{n'}^* \rangle$, $\langle b_n b_{n'}^* \rangle$ for $n, n' = 1, \dots, n_c$. These quantities require integration over the monomer radius distribution $f(s)$. Since the Schulz distribution has no singularity at zero radius and decreases exponentially at large radius, Gauss-Laguerre quadrature is efficient (check Cauchy-Schwarz inequalities: $|\langle a_n \rangle|^2 \le |\langle a_n \rangle|^2$ and $|\langle a_n b_{n'}^* \rangle|^2 \le |\langle a_n a_{n'}^* \rangle| |\langle b_{n'} b_{n}^* \rangle|$).

* 6) the multiple-scattering correction factors $\delta_n^{(a)}$, $\delta_n^{(b)}$ for $n = 1, \dots, n_c$ are found after computing the inverse of the $2n_c \times 2n_c$ matrix (26)-(27). Since $n_c$ is generally small and the matrix not singular, any stable method of matrix inversion can be used (check: $\delta_n^{(a)}$, $\delta_n^{(b)}$ verifies (26)-(27)).

* 7) once $\delta_n^{(a)}$, $\delta_n^{(b)}$ known, the optical cross-sections are given using (60), (53), (61); asymmetry parameter is defined in (56) after calculating (58), (57); linear polarization is (39), (40), (42), (44), (45) then (46); scattering phase function is calculated using (59) and (49)-(51).

That way, the MFTM values of the optical quantities for a fractal aggregate of any number of monomers are typically obtained within a few minutes on a laptop (the cpu time depends only logarithmically on the aggregate size).
## 5. Application of the MFTM method and comparison with T-Matrix and RGD calculations

Extensive comparison between MFTM and T-Matrix results will be the topic of a subsequent article [65]. However, scattering and absorption cross-sections in cases relevant for soot [66] or for tholins [67] are discussed here to give a few brief illustrations of the MFTM method.

MFTM approach can essentially be used in two ways, either as approximation (i.e. to calculate approximate optical quantities for large, $N \gg 1$, aggregates) or as extrapolation (i.e. to derive analytical behaviours of optical quantities for large aggregates). Both ways are illustrated below.

### 5.1. An example with $d_f > 2$

When the fractal dimension is larger than 2, particle boundaries are expected to play a leading role in light scattering processes [68]. This can be seen when generally discussing the MFTM approach as an extrapolation scheme for large aggregates, as follows.

For any fractal dimension $d_f$, relation (79) writes:

$$S_p^{(R)} \approx \kappa d_f k_f \left(\frac{R_g}{s_0}\right)^{d_f} \int_{2s_0/R_g}^{\infty} j_p(k R_g v) h_p^{(1)}(k R_g v) v^{d_f-1} e^{-\kappa v^{d_f}} dv. \quad (82)$$

The function $v^{d_f-1} \exp(-\kappa v^{d_f})$ reaches maximum for the finite positive value $v^* = ((d_f - 1)/d_f)^{1/d_f} \sqrt{d_f/\Gamma(2/d_f)}$ which is an increasing function of the fractal dimension, going from $v^* = 0$ for $d_f = 1$ to $v^* = 1.30$ for $d_f = 3$. Around $v^*$, the value of the argument $k R_g v$ is large when $N \gg 1$, so the large-$x$ approximation $2x^2 j_p(x) h_p^{(1)}(x) \approx 1 + (-1)^{p+1} \exp(2ix)$ can be used in (82). The fractal dimension being here larger than 2, the leading behaviour of $S_p^{(R)}$ is $O(N^{1-2/d_f})$, then the $N$-independent $S_p^{(L)}$ is insignificant in front of $S_p^{(R)}$. Thus (using (75)):

$$S_p \approx \frac{\pi}{2d_f \sin 2\pi/d_f} \frac{1}{(k R_g)^2} \quad (83)$$

is independent of the index $p$. The asymptotic behaviours of the coefficients $A_n^v, B_n^v$ result from (83), the definitions (28) and identities (33):

$$(N-1)A_n^v \approx (N-1)B_n^v \approx \frac{1}{c_f (k R_g)^2} (2v+1), \quad (84)$$

with the constant $c_f = 8 j_0(2\pi/d_f)$.

In a like manner [69], the coefficients $K_{n,n'}^{(1,2)}$ defined in (54) – and required to calculate the scattering cross-section – can be proven to be $O(1/(k R_g)^2)$ when $N \gg 1$.

Using (84), the MFTM equations (26)-(27) write:

$$\delta_n^{(a)} = \delta_n^{(b)} \qquad , \quad (85)$$

$$\delta_n^{(a)} = 1 - \frac{N}{c_f (k R_g)^2} \sum_{v=1}^{\infty} (2v+1)(\langle \bar{a}_v \rangle + \langle \bar{b}_v \rangle) \delta_v^{(a)}. \quad (86)$$

The self-consistent equation (86) is simple enough to discuss briefly the existence of a stable solution $\delta_n^{(a)}$, i.e. the question of whether the Mean-Field T-matrix equations (26)-(27) are invertible.

Let us define the mean size parameter $x_0 \equiv k s_0$, and the amplitude $\Sigma_0 = (1/2) \sum_v (2v+1)(\langle \bar{a}_v \rangle + \langle \bar{b}_v \rangle)$ of the wave scattered by the aggregate in the forward direction. Let us also define the auxiliary parameter: $\Delta = k_f^{2/d_f} N^{1-2/d_f} \Sigma_0 / (c_f x_0^2)$. Solving iteratively (86) by successive scattering orders, $l$, starting from the single scattering case ($l = 0$; $\delta_n^{(a)} = 1$), one finds, at the order $l$:

$$\delta_n^{(a)} = 1 - \Delta + \Delta^2 - \dots + (-1)^l \Delta^l. \quad (87)$$
It is clear from (87) that scattering order iterations can be continued ad infinitum when $|\Delta| < 1$. This is a sufficient condition for (86) to have a finite solution in the scattering coefficients $\delta_n^{(a)}$ [70]. This condition can be written more explicitly as the upper bound:

$$N < \left| \frac{c_f}{k_f^{2/d_f}} \frac{\lambda_0^2}{\Sigma_0} \right|^{d_f/(d_f-2)} \quad (88)$$

When (88) is realized, the scattering coefficient $\delta_n^{(a)}$ is certainly defined and is independent of the multipole index $n$, hence the general asymptotic solution (generalizing a result known for the monodisperse monomers case [38, 71]):

$$\delta_n^{(a)} = \delta_n^{(b)} \approx \frac{1}{1 + \frac{\Sigma_0}{c_f} \frac{N}{(k R_g)^2}} \quad (89)$$

Note however that $\delta_n^{(a)} = 1/(1 + \Delta)$ is always solution of (86), even if condition (88) is not fulfilled. This raises the unresolved question of the relevance of the mathematical solution when it is not the consequence of the iterative method in scattering order (this point is beyond the topic of this work).

When solution of (86) exists, the scattering coefficients (89) – which stem for multiple scattering contribution within the aggregate – behave like $\sim N^{2/d_f}/N$, that is like the ratio between an effective surface area ($\sim R_g^2$) of the particle and its effective volume ($\sim a_{\text{eff}}^3$).

Using (89), the general formulae (61) and (53) lead to the large-$N$ behaviours of the absorption and scattering cross-sections in the case $d_f > 2$:

$$C_{\text{abs},N} \approx \frac{c_f^2}{|\Sigma_0|^2} \frac{(k R_g)^4}{N} C_{\text{abs},1}, \quad (90)$$

$$C_{\text{sca},N} \approx 2\pi c_f R_g^2, \quad (91)$$

with $C_{\text{abs},1} = 2\pi \sum_n (2n + 1)(c_n^{(1)} + c_n^{(2)})/k^2$, according to (61).

Since the geometric cross-section, $G$, of a fractal aggregate of fractal dimension $d_f > 2$ scales as [71]: $G \sim R_g^2$, it results from (91) that the scattering efficiency defined as $Q_{\text{sca}} = C_{\text{sca}}/G$ tends to a constant value when $N \to \infty$. This resembles the case of geometrical optics where scattering, absorption and extinction efficiencies tend to constant values as the size parameter tends to infinity. However, the absorption efficiency resulting from (90) tends to 0 when $N \to \infty$, meaning that absorption by a fractal aggregate is less efficient than for a compact particle of same size.

MFTM asymptotic behaviours (90)-(91) are general when $d_f > 2$ provided the linear MFTM equations (26)-(27) can be solved. Let us apply these results to the case studied in [72], that is: soot refractive index $m = 1.6 + 0.6i$, green-laser wavelength $\lambda = 532$ nm, monomers of mean radius $s_0 = 15$ nm, fractal dimension $d_f = 2.1$, and prefactor $k_f = 2.3$. From (90) and $4/d_f - 2 = -0.095$, the behaviour:

$$C_{\text{abs},N}/(N C_{\text{abs},1}) \sim 1/N^{0.095} \quad (92)$$

is expected when $N \gg 1$ (sufficient condition (88) leads to $N \le 3693019$ though no failure of matrix inversion was observed in this case up to $N = 10^9$). MFTM approach clearly departs from the regular RGD method, as the latter predicts: $C_{\text{abs},N}/(N C_{\text{abs},1}) = 1$ for any size $N$ of the aggregate. In [73], the authors argued that RGD approach has to be corrected to take into account multiple scattering within the aggregate, and that these corrections should lead to smooth decrease of the ratio $C_{\text{abs},N}/(N C_{\text{abs},1})$. This conclusion is in qualitative agreement with the analytical result (92).

Figure 3 shows comparison of T-matrix results of $C_{\text{abs},N}/(N C_{\text{abs},1})$ for various values of $N$ (black squares, $N = 15, 50, 100, 200, 400, 600$ from [72] using GMM01f code [74]) with the MFTM extrapolation result (90) as the dashed line and direct MFTM results (continuous red line) using (26)-(27), (28), (30), (35), (71), (75) and (61) with polydispersity $\sigma = 0.005$ – that is close to the monodisperse case –, and aggregate size up to $N = 1000000$. Even if the mean-field approach is supposed not to hold for small systems, one can see on the figure that MFTM realizes a good approximation even for moderate aggregate sizes.
<figure><img src="image_5.png" /><figcaption>Figure 3. Double logarithmic plot of fractal ($d_f = 2.1$) aggregate normalized absorption cross-sections, $C_{abs,N}/(NC_{abs,1})$, vs. aggregate size $N$. Monomers are all of radius $s_0 = 15$ nm and material with refractive index $m = 1.6 + 0.6i$. Black squares are T-Matrix results reproduced from [72], for $N = 15, 50, 100, 200, 400, 600$. Red continuous curve represents the MFTM data between $N = 15$ and $N = 1000000$. Dashed line is the extrapolation formula (90) ~ $N^{-0.095}$. The line $C_{abs,N}/(NC_{abs,1}) = 1$ is the Rayleigh-Gans-Debye (RGD) approximation.</figcaption></figure>

## 5.2. Examples with $d_f < 2$

In [75], numerical soot aggregates similar to those above are considered, but with fractal dimension $d_f = 1.78$ (see also [43] for simulated aggregates of soot with fractal dimension $d_f = 1.9$). In Table 1, T-matrix values of the absorption and scattering cross-sections are compared with the corresponding values using the MFTM method. The error bars on the T-matrix cross-sections result from different realizations of the fractal aggregates of given $N$.

For $N \ge 200$, the differences between exact and approximate values are found to be less than 7%. Even in the cases of smaller aggregates, values obtained from both methods agree well.

<table><caption>Table 1. Absorption and scattering cross-sections ($d_f = 1.78, m = 1.6 + 0.6i, s_0 = 15$ nm). T-matrix results vs MFTM results</caption><thead><tr><th>N</th><th>C<sub>abs,N</sub> (nm<sup>2</sup>) T-matrix</th><th>C<sub>abs,N</sub> (nm<sup>2</sup>) MFTM</th><th>C<sub>sca,N</sub> (nm<sup>2</sup>) T-matrix</th><th>C<sub>sca,N</sub> (nm<sup>2</sup>) MFTM</th></tr></thead><tbody><tr><td>10</td><td>1 508.9 ± 9.5</td><td>1 476.4</td><td>43.0 ± 0.1</td><td>27.8</td></tr><tr><td>20</td><td>3 046.5 ± 40.0</td><td>3 066.4</td><td>154.2 ± 1.6</td><td>120.3</td></tr><tr><td>50</td><td>7 396.3 ± 28.1</td><td>7 936.7</td><td>675.9 ± 13.9</td><td>626.9</td></tr><tr><td>100</td><td>15 012.0 ± 94.5</td><td>15 939.2</td><td>1 769.2 ± 79.8</td><td>1 781.1</td></tr><tr><td>200</td><td>29 551.9 ± 171.4</td><td>31 584.0</td><td>4 312.4 ± 188.4</td><td>4 512.2</td></tr><tr><td>400</td><td>58 097.1 ± 290.5</td><td>62 000.4</td><td>10 004.6 ± 268.0</td><td>10 622.3</td></tr><tr><td>600</td><td>86 387.7 ± 457.9</td><td>91 697.0</td><td>16 182.5 ± 856.0</td><td>17 088.8</td></tr></tbody></table>

Another example is worth giving here to briefly discuss the role of monomer polydispersity. Let us consider the case of fractal aggregates of monomers with radius polydispersity $\sigma$. Fractal dimension is $d_f = 1.9$ (they are simulated using the Ballistic Cluster-Cluster Aggregation model [76]) and the monomers, of mean size parameter 0.048, are made of a material of refractive index $m = 1.2052+0.2568i$ (this is a possible value of $m$ for Titan tholins in the 3.4 µm absorption band [67]). T-matrix calculations are done using MSTM-v4.0 code [74]. For monodisperse (i.e. $\sigma = 0$) monomers, the T-matrix values of the absorption and scattering efficiencies are: $Q_{abs} = 1.5195$ and $Q_{sca} = 0.02676$ respectively (average over three independent aggregates of size $N = 131072$ and 490 orientations). These values are very close to the MFTM results for $\sigma = 0.01$, namely: $Q_{abs} = 1.5200$ and $Q_{sca} = 0.02670$ respectively.
<figure><img src="image_6.png" /><figcaption>Figure 4. MFTM results of absorption (blue squares) and scattering (red circles) cross-sections vs. polydispersity $\bar{\sigma}$ for fractal aggregates $d_f = 1.9$, and refractive index $m = 1.2052 + 0.2568i$. Mean size parameter $k s_0 = 0.048$. The monomer radii are distributed following the Schulz distribution (66). For comparison, the two black circles represent the T-matrix values of the respective cross-sections for the monodisperse case ($\bar{\sigma} = 0$).</figcaption></figure>

Evolution of these efficiencies with increasing polydispersity is shown on Fig. 4. For this very absorbing material, absorption and scattering efficiencies increase significantly with the polydispersity of the monomers. It means that one should use numerical codes including polydispersity (such as the present version of MFTM), since real aggregates in nature are generally made of polydisperse grains.

## 6. Conclusion

The refined Mean-Field T-Matrix (MFTM) method developed in this work represents a substantial and original advancement in the modelling of light scattering by large, irregular dust aggregates. While it builds on the foundational MFTM framework introduced in earlier work, the present version significantly expands its applicability by incorporating polydisperse monomers, heterogeneous compositions, and a physically realistic pair correlation function, (71)-(73), to describe aggregate structure. To date, no other asymptotic method in the field offers this level of generality and efficiency.

This approach is particularly well-suited for statistically isotropic fractal aggregates composed of thousands of monomers or more—a regime where exact numerical methods become computationally prohibitive, and where the mean-field approximation is not only valid but highly effective. Under these conditions, the total electromagnetic field experienced by each monomer can be treated as position-independent, enabling a compact and accurate description of the aggregate's optical behavior.

Although the current formulation is optimized for fractal geometries, it can be naturally extended to porous non-fractal aggregates by taking a fractal dimension $d_f = 3$. In this case, porosity emerges as the key structural parameter, and the cutoff behavior of the pair correlation function becomes critical for distinguishing between aggregates with similar porosity. This direction will be the focus of future developments.

The original MFTM method [38] has already been successfully applied to a variety of problems – ranging from latex [77, 78] or polystyrene [79] aggregates, from colloidal nanoparticles [80] to planetary[81, 82], astrophysical environments [83] and brown dwarf atmospheres [84], as well as astrophysical dust. The present refinement opens the way to a new class of applications, particularly in atmospheric and astrophysical science, where large, irregular, and polydisperse aggregates are not only expected but known to influence observational signatures. By offering an accurate, versatile, and computationally efficient solution in a regime inaccessible to exact methods, the refined MFTM stands out as a uniquely capable tool in the field of light scattering by complex particulate media.
# 7. Acknowledgements

The author would like to thank Pascal Rannou for helpful discussions on this work and for providing the T-matrix results discussed in Fig. 4.

# Appendix A Mean-field expressions of the coefficients $A_{\mu\nu,mn}^{(J)}(\mathbf{r}_{lj})$ and $B_{\mu\nu,mn}^{(J)}(\mathbf{r}_{lj})$

## A.1 T-matrix

The addition coefficients $A_{\mu\nu,mn}^{(J)}(\mathbf{r}_{lj})$ and $B_{\mu\nu,mn}^{(J)}(\mathbf{r}_{lj})$ can be written under the forms [38]:

$$A_{\mu\nu,mn}^{(J)}(\mathbf{r}_{lj}) = (-1)^m i^{n-v} \frac{2n+1}{2n(n+1)} \sum_{\substack{p=[n-v] \\ p+n+v \text{ even}}}^{n+v} i^p [n(n+1) + v(v+1) - p(p+1)] a(\mu, v; -m, n; p) \times \Phi_{\mu-m,p}^{(J)}(\mathbf{r}_{lj}), \quad (93)$$

$$B_{\mu\nu,mn}^{(J)}(\mathbf{r}_{lj}) = (-1)^{m+1} i^{n-v} \frac{2n+1}{2n(n+1)} \sum_{\substack{p=[n-v]+1 \\ p+n+v \text{ odd}}}^{n+v-1} i^p \frac{2p+1}{2p-1} b(\mu, v; -m, n; p) \times \Phi_{\mu-m,p}^{(J)}(\mathbf{r}_{lj}), \quad (94)$$

$$b(\mu, v; -m, n; p) = (n-m)(n+m+1) a(\mu, v; -m-1, n; p-1) + 2m(p-\mu+m+1) a(\mu, v; -m, n; p-1) - (p-\mu+m+1)(p-\mu+m) a(\mu, v; -m+1, n; p-1), \quad (95)$$

where the auxiliary function $\Phi_{\mu-m,p}^{(J)}$ is defined by:

$$\Phi_{\mu-m,p}^{(J)}(\mathbf{r}_{lj}) = z_p^{(J)}(kr_{lj}) P_p^{\mu-m}(\cos \theta_{lj}) e^{i(\mu-m)\phi_{ij}}, \quad (96)$$

and the function $z_p^{(J)}$ is related to either the spherical Bessel or the spherical Hankel function of the first kind according to the value of the upper index (J), namely: $z_p^{(1)}(x) = j_p(x)$ and $z_p^{(3)}(x) = h_p^{(1)}(x)$. The vector linking monomer j and monomer l is: $\mathbf{r}_{lj} = \mathbf{r}_j - \mathbf{r}_l = r_{lj}$ ($\sin \theta_{lj} \cos \phi_{lj}, \sin \theta_{lj} \sin \phi_{lj}, \cos \theta_{lj}$). The numerical coefficients $a(\mu, v; m, n; p)$ in (93) and (95) are defined through the formal expansion of associated Legendre polynomials:

$$P_{\nu}^{\mu}(x) P_{n}^{m} = \sum_{p=[n-v]}^{n+v} a(\mu, v; m, n; p) P_{p}^{\mu+m}(x). \quad (97)$$

From the identity $(v-\mu)!P_{\nu}^{\mu}(x) = (-1)^{\mu}(v+\mu)!P_{\nu}^{-\mu}(x)$ applied to (97), the following symmetry relations are found:

$$\frac{(p - m + \mu)!}{(\nu + \mu)!(n - m)!} a(\mu, \nu; -m, n; p) = \frac{(p + m - \mu)!}{(\nu - \mu)!(n + m)!} a(-\mu, \nu; m, n; p), \quad (98)$$

$$\frac{(p - m + \mu)!}{(\nu + \mu)!(n - m)!} b(\mu, \nu; -m, n; p) = \frac{(p + m - \mu)!}{(\nu - \mu)!(n + m)!} b(-\mu, \nu; m, n; p). \quad (99)$$

## A.2 Mean-Field addition coefficients (general degree)

The mean-field equations (22)-(23) use the average addition coefficients:

$$\langle\langle A_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj}) e^{-i k r_{lj} \cos \theta_{lj}} \bar{a}_{\nu}^{l'} \rangle\rangle; \quad \langle\langle B_{\mu\nu,mn}^{(3)}(\mathbf{r}_{lj}) e^{i k r_{lj} \cos \theta_{lj}} \bar{b}_{\nu}^{l'} \rangle\rangle, \quad (100)$$

and similar forms where Mie coefficients $\bar{a}$ are replaced with $\bar{b}$. Symbols $\langle\langle \cdots \rangle\rangle$ denote averages over orientation, the pair distribution function $g(r)$ and the size- and refractive index-distribution. In both expressions (100), monomer locations appear in the only combination $\Phi_{\mu-m,p}^{(J)}(\mathbf{r}_{lj})$.

Using the identity:

$$i^p \int_0^\pi P_p(\cos \theta) e^{-i k r \cos \theta} \sin \theta d\theta = 2 j_p(kr), \quad (101)$$
and the assumption that no size segregation is present in the aggregate, the average value of $\Phi_{\mu-m,p}^{(J)}$ writes:

$$ \langle \Phi_{\mu-m,p}^{(J)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle = (-i)^p \langle j_p(kr) z_p^{(J)}(kr) \rangle \langle \bar{a}_{\nu} \rangle \delta_{\mu-m,0}, \quad (102) $$

where $\langle \cdots \rangle$ denotes average over the orientation and pair-correlation function. In particular, note that all the quantities (100) vanish for $\mu \neq m$ because of the angular averaging: $\langle \exp(i(\mu - m)\phi_{lj}) \rangle = 0$ in this case.

After introducing the quantity:

$$ S_p = \langle j_p(kr) h_p^{(1)}(kr) \rangle, \quad (103) $$

where average occurs over the distribution of the monomer-monomer distances $r_{lj}$, relations (100) write in the case $\mu = m$:

$$ \langle A_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle = (-1)^m i^{n-\nu} \frac{2n+1}{2n(n+1)} \langle \bar{a}_{\nu} \rangle \sum_{\substack{p=|n-\nu| \\ p+n+\nu \text{ even}}}^{n+\nu} [n(n+1) + \nu(\nu+1) - p(p+1)] a(m, \nu; -m, n; p) S_p, \quad (104) $$

$$ \langle B_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{b}_{\nu}^{l} \rangle = (-1)^{m+1} i^{n-\nu} \frac{2n+1}{2n(n+1)} \langle \bar{b}_{\nu} \rangle \sum_{\substack{p=|n-\nu|+1 \\ p+n+\nu \text{ odd}}}^{n+\nu-1} \frac{2p+1}{2p-1} b(m, \nu; -m, n; p) S_p. \quad (105) $$

with the symmetry relations resulting from (98)-(99):

$$ \langle A_{-mv,-mn}^{(3)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle = \frac{(v-m)!}{(v+m)!} \frac{(n+m)!}{(n-m)!} \langle A_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle, \quad (106) $$

$$ \langle B_{-mv,-mn}^{(3)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle = - \frac{(v-m)!}{(v+m)!} \frac{(n+m)!}{(n-m)!} \langle B_{mv,mn}^{(3)}(\mathbf{r}_{lj}) e^{-i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle, \quad (107) $$

and the expressions for the numerical coefficients $a(\cdots)$ [85]:

$$ a(m, \nu; -m-1, n; p) = - \frac{2p+1}{2} \int_{-1}^{+1} P_{\nu}^{m}(x) P_{n}^{-m-1}(x) P_{p}^{1}(x) dx, \quad (108) $$

$$ a(m, \nu; -m, n; p) = \frac{2p+1}{2} \int_{-1}^{+1} P_{\nu}^{m}(x) P_{n}^{-m}(x) P_{p}(x) dx, \quad (109) $$

$$ a(m, \nu; -m+1, n; p) = - \frac{2p+1}{2} \int_{-1}^{+1} P_{\nu}^{m}(x) P_{n}^{-m+1}(x) P_{p}^{-1}(x) dx. \quad (110) $$

### A.3 Mean-Field addition coefficients (degree $m = \pm 1$)

The following symmetry relations are found from (106)-(107):

$$ \langle A_{-1\nu,-1n}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle = \frac{n(n+1)}{\nu(\nu+1)} \langle A_{1\nu,1n}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle \quad (111) $$

$$ \langle B_{-1\nu,-1n}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle = - \frac{n(n+1)}{\nu(\nu+1)} \langle B_{1\nu,1n}^{(3)}(\mathbf{r}_{lj}) e^{i\mathbf{k}_{\text{inc}} \cdot \mathbf{r}_{lj}} \bar{a}_{\nu}^{l} \rangle \quad (112) $$

## Appendix B Pair correlation function of a fractal particle made of aggregated polydisperse monomers (Schulz-distribution of monomer radius)

### Partial pair correlation functions

Assuming homogeneous spatial distribution of the monomer-radius values inside the aggregate [58], the partial pair correlation functions write [40]:

$$ 4\pi r^2 g_{s,s'}(r) = \frac{2}{N} \delta(r - (s + s')) , \quad \text{if } 0 < r \le s + s', \quad (B.1) $$

$$ = \left(1 - \frac{2}{N}\right) c_{s+s'} 4\pi r^2 g_W(r) , \quad \text{if } r > s + s', \quad (B.2) $$
where the regular function $g_W(r)$ – characterizing the fractal distribution of the monomers inside the aggregate – is the normalized Weibull distribution [38]:

$$4\pi r^2 g_W(r) = \frac{\kappa d_f}{R_g^{d_f}} r^{d_f-1} e^{-\kappa(r/R_g)^{d_f}}, \qquad (B.3)$$

with $d_f$ the fractal dimension of the aggregate. The parameter $\kappa$ is [40]: $\kappa = (\Gamma(2/d_f)/d_f)^{d_f/2}$, and the coefficient $c_{s+s'}$ in (B.2) is calculated from the normalization of $4\pi r^2 g_{s,s'}(r)$:

$$c_{s+s'} = e^{\kappa((s+s')/R_g)^{d_f}}. \qquad (B.4)$$

In principle, Weibull distribution in (B.2) is correct for $s + s' < r < \gamma R_g$, with $\gamma \approx 2.5$ for all known regular aggregate models. Beyond $\gamma R_g$, the pair correlation function vanishes much faster than Weibull function. For the sake of simplicity, the approximation $\gamma R_g = \infty$ is used in (B.2) since $r > \gamma R_g$ corresponds to a very small amount of very distant monomer pairs when the aggregate is large.

### Total pair correlation function (Schulz distribution of the monomer sizes)

Using the Schulz distribution (66), the total pair correlation function is:

$$g(r) = \frac{1}{\Gamma^2(z)} \left(\frac{z}{s_0}\right)^{2z} \int_0^\infty \int_0^\infty s^{z-1} s'^{z-1} e^{-z(s+s')/s_0} g_{s,s'}(r) ds ds', \qquad (B.5)$$

that is, using (B.1)-(B.2) [41]:

$$4\pi r^2 g(r) = \frac{2}{N} 4\pi r^2 g_L(r) + \left(1 - \frac{2}{N}\right) 4\pi r^2 g_R(r), \qquad (B.6)$$

$$4\pi r^2 g_L(r) = \frac{1}{\Gamma(2z)} \left(\frac{z}{s_0}\right)^{2z} r^{2z-1} e^{-zr/s_0}, \qquad (B.7)$$

$$4\pi r^2 g_R(r) = \frac{1}{\Gamma(2z)} 4\pi r^2 g_W(r) \int_0^{zr/s_0} u^{2z-1} e^{-u} e^{\kappa(s_0 u/z R_g)^{d_f}} du, \qquad (B.8)$$

with the functions $g_L(r)$ and $g_R(r)$ defined over the interval $r \in (0, \infty)$. Normalization condition $\int_0^\infty 4\pi r^2 g(r) dr = 1$ results from integration by parts.

The function (B.7) corresponds to the local singular part of the pair correlation distribution function [86], while the function (B.8) is the regular part. Comparing the magnitudes of these two terms is insightful: on Fig. 5, it is clear that the local part of $g(r)$ is dominant in the pair-distance domain around $r \sim 2s_0$.

When the aggregate is large (more precisely, when $zR_g/s_0 \gg 1$), the approximation: $e^{\kappa(s_0 u/z R_g)^{d_f}} \approx 1$ used in (B.8) leads to:

$$4\pi r^2 g_R(r) \approx 4\pi r^2 g_W(r) \left(1 - \frac{\Gamma(2z, zr/s_0)}{\Gamma(2z)}\right). \qquad (B.9)$$

An example is shown in Fig. 6 to compare the effects of the two options (B.8) or (B.9). The conclusion is that (B.9) is accurate for usual applications (even for moderate sizes, such as here: $N = 100$).

## Appendix C Some inequalities

### C.1 Inequalities for the coefficients $K_{n,n'}^{(1,2)}$

The auxiliary functions $\Phi_{n,n'}$ and $\Psi_{n,n'}$ defined in (50)-(51) are written under the alternative forms:

$$\Phi_{n,n'}(\theta) = \frac{(2n+1)(2n'+1)}{2} \left[ \frac{P_n' P_{n'}'}{n(n+1)n'(n'+1)} (1+\cos^2\theta) + P_n P_{n'} + \left( \frac{P_n' P_{n'}}{n(n+1)} + \frac{P_n P_{n'}}{n'(n'+1)} \right) \cos\theta \right] \qquad (C.1)$$

$$\Psi_{n,n'}(\theta) = \frac{(2n+1)(2n'+1)}{2} \left[ - \frac{2P_n' P_{n'}'}{n(n+1)n'(n'+1)} \cos\theta + \frac{P_n' P_{n'}}{n(n+1)} + \frac{P_n P_{n'}}{n'(n'+1)} \right] \qquad (C.2)$$
<figure><img src="image_7.png" /><figcaption>Figure 5. Plots of the local distribution ($4\pi r^2 g_L(r)$, red continuous line) and of the limited Weibull distribution ($4\pi r^2 g_R(r)$, black thick continuous line). Parameters are: $s_0 = 1/2$, $d_f = 1.9$, $N = 100$, and Schulz parameter: $z = 100$. The local part is prevailing around $r \approx 2s_0$. In the case of monodisperse monomers (i.e. $z \to \infty$), the local singular part of the pair correlation function is the Dirac delta function: $4\pi r^2 g_L(r) \to \delta(r - 2s_0)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Comparison between the approximated pair correlation function (B.6) using (B.9) (black dashed curve) and the Schulz pair correlation function using (B.8) (red continuous curve). Parameters are: $d_f = 1.9$, $N = 100$ (hence: $R_g/s_0 = 11.3$). Schulz parameter is $z = 100$. Both functions are almost identical.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7. Blow-up of a 3D modified BCCA fractal aggregate (see text) made of 16 384 Schulz-distributed monomers (polydispersity: 0.5), and exhibiting radius-segregation for monomers belonging to bin 3 (blue-coloured monomers). Colours are the same as in Fig. 1. The corresponding monomer-size segregation is studied quantitatively on Fig. 8.</figcaption></figure>

in which $P_n \equiv P_n(\cos \theta)$ and $P'_n \equiv dP_n/dx$ at $x = \cos \theta$. Using the known inequalities [87]: $|P_n| \le 1$ and $|P'_n| \le n(n+1)/2$, one finds: $|\Phi_{n,n'}(\theta)| \le (2n+1)(2n'+1)(5+|\cos \theta|+\cos^2 \theta)/8$ and $|\Psi_{n,n'}(\theta)| \le (2n+1)(2n'+1)(2+|\cos \theta|)/4$. Since $|S(\theta)| \le 1$ (see (43)), the following inequalities are found after integration over $\theta$:

$$ \mathcal{K}_{n,n'}^{(1)} \le \frac{35}{24}(2n+1)(2n'+1); \quad \mathcal{K}_{n,n'}^{(2)} \le \frac{5}{4}(2n+1)(2n'+1) \qquad (C.3) $$

## C.2 Inequality for the $j_p(x) h_p^{(1)}(x)$ functions

Starting from the representation (after [87] (18.17.20 and 18.17.21)):

$$ j_p(x) h_p^{(1)}(x) = \frac{1}{ix} \int_0^1 e^{2ixu} P_p(1-2u^2) du, \qquad (C.1) $$

and using the improved Bernstein's inequality [88] $(1-x^2)^{1/4}|P_p(x)| < 2/\sqrt{\pi(2p+1)}$ for Legendre polynomials, expression (C.1) leads to the inequality:

$$ |j_p(x) h_p^{(1)}(x)| < \sqrt{\frac{\pi}{2p+1}} \frac{1}{x}, \qquad (C.2) $$

for $x > 0$ and $p \ge 0$.

# Appendix D A model of Cluster-Cluster aggregation with forced monomer-size segregation

Regular Cluster-Cluster aggregation models can be modified to force monomers segregation by their sizes. A simple way is to change the sticking rules on collisions. In the regular models, sticking probability (i.e. the probability to stick irreversibly during a collision event) is = 1 whatever the colliding monomers.
<figure><img src="image_10.png" /><figcaption>Figure 8. Mean number, $n_i(j)$, of monomers belonging to bin $j$ and tangent to bin-$i$ monomers for modified-BCCA aggregates (see text). Average over 100 independent 3D aggregates of size $N = 16\,384$. Monomer radii are Schulz-distributed (polydispersity: 0.5). The complete monomer radii range is divided into 16 bins, each bin consisting of 1024 monomers of about the same size. In the figure, the peak in the $n_i(j)$ function characterizes the monomer-size segregation. This figure can be compared to the Fig. 2 (using the same range of values) corresponding to the regular BCCA model without forced monomer-size segregation.</figcaption></figure>

Let us consider the following example of modified BCCA model. The set of sizes of the ensemble of monomers to be aggregated is split into 16 bins. Each bin corresponds to monomers of about the same radius. In the modified BCCA model, the probability for two monomers to stick during a collision = 1 except in the case where a monomer belonging (e.g.) to the bin 3 collides a monomer from another bin. In the latter case, the sticking probability is taken to be very small (say $10^{-3}$, for example). This way, it is very unlikely that sticking occurs between a bin 3 monomer and a monomer from a bin other than 3. Consequently, there is segregation of the monomers of bin 3, forming a fractal blob of monomers of about the same size inside the complete aggregate. An example of such aggregate of $N = 16\,384$ polydisperse monomers is shown on Fig. 7. The corresponding function number $n_i(j)$ of monomers of bin $j$ tangent to a monomer of bin $i$ is shown on Fig. 8. Size segregation results in the high peak $n_3(3)$, that is a monomer of the bin 3 is in average connected to about two other monomers belonging to the same bin. One can also note the decreases of $n_3(j)$ and of $n_j(3)$ when $j \neq 3$.

This kind of two-step aggregation process (first step is creation of a BCCA aggregate, $B_s$, made of segregated monomers, followed by the second step : BCCA aggregation of the aggregate $B_s$ and all the other polydisperse monomers) may be realized in some conditions, but this topic is still overlooked (see [89] for an example).
## References

[1] Long, CM, Nascarella, MA, Valberg, PA. Carbon black vs. black carbon and other airborne materials containing elemental carbon: physical and chemical distinctions. *Environmental Pollution* 2013; 181:271-286. https://doi.org/10.1016/j.envpol.2013.06.009

[2] International Carbon Black Association (ICBA), 2004. *Carbon Black User's Guide: Safety, Health, & Environmental Information*. International Carbon Black Association (ICBA).

[3] Kok, JF. A scaling theory for the size distribution of emitted dust aerosols suggests climate models underestimate the size of the global dust cycle. *Proceedings of the National Academy of Sciences of the United States of America (PNAS)* 2010;108:1016-1021. https://doi.org/10.1073/pnas.101479810

[4] Weidling, R, Güttler, C, Blum, J. Free collisions in a microgravity many-particle experiment. I. Dust aggregate sticking at low velocities. *Icarus* 2012;218:688-700. https://doi.org/10.1016/j.icarus.2011.10.002

[5] Brisset, J, Heißelmann, D, Kothe, S, Weidling, R, Blum, J. Low-velocity collision behaviour of clusters composed of sub-millimetre sized dust aggregates. *A&A* 2017;603:A66-1-11. https://doi.org/10.1051/0004-6361/201630345

[6] Langkowski, D, Teiser, J Blum, J. The Physics of Protoplanetary Dust Agglomerates. II. Low-Velocity Collision Properties 2008). *The Astrophysical Journal* 2008;675:764-776. https://doi.org/10.1086/525841

[7] Blum, J, Wurm, G. The Growth Mechanisms of Macroscopic Bodies in Protoplanetary Disks. *Annu. Rev. Astron. Astrophys.* 2008;46:21-56. https://doi.org/10.1146/annurev.astro.46.060407.145152

[8] Testi, L, Birnstiel, T, Ricci, L, Andrews, S, Blum, J, Carpenter, J, Dominik, C, Isella, A, Natta, A, Williams, JP Wilner, DJ. Dust evolution in protoplanetary disks. in Protostars and Planets VI, University of Arizona Press (2014), eds. H. Beuther, R. Klessen, C. Dullemond, Th. Henning; 339-361. https://doi.org/10.2458/azu_uapress_9780816531240-ch015

[9] Zsom, A, Dullemond, CP. A representative particle approach to coagulation and fragmentation of dust aggregates and fluid droplets. *A&A* 2008;489:931-941. https://doi.org/10.1051/0004-6361:200809921

[10] Dominik C, Tielens, AGGM. The Physics of Dust Coagulation and the Structure of Dust Aggregates in Space. *The Astrophysical Journal* 1997;480:647-673. https://doi.org/10.1086/303996

[11] Suyama, T, Wada, K, Tanaka, H. Numerical Simulation of Density Evolution of Dust Aggregates in Protoplanetary Disks. I. Head-on Collisions. *The Astrophysical Journal* 2008;684:1310-1322. https://doi.org/10.1086/590143

[12] Güttler, C, Mannel, T, Rotundi, A, Merouane, S, Fulle, M, Bockelee-Morvan, D, Lasue, J, Levasseur-Regourd, AC, Blum, J, Naletto, G, Sierks, H, Hilchenbach, M, Tubiana, C, Capaccioni, F, Paquette, JA, Flandes, A, Moreno, F, Agarwal, J, Bodewits, D, Bertini, I, Tozzi, GP, Hornung, K, Langevin, Y, Kruger, H, Longobardo, A, V Corte, D, Tóth, I, Filacchione, G, Ivanovski, SL, Mottola, S, Rinaldi, G. Synthesis of the morphological description of cometary dust at comet 67P/Churyumov-Gerasimenko. *A&A* 2019;630:A24-1-15. https://doi.org/10.1051/0004-6361/201834751

[13] Lasue, J, Maroger, I. Botet, R, Garnier, Ph, Merouane, S, Mannel, Th, Levasseur-Regourd, AC Bentley, MS. Flattened loose particles from numerical simulations compared to particles collected by Rosetta. *A&A* 2019;630:A28-1-9. https://doi.org/10.1051/0004-6361/201834766

[14] Podolak, M, Flandes, A, Della Corte, V, Krüger, H. Simple Model for Understanding the DIM Dust Measurement at Comet 67P/Churyumov-Gerasimenko. *Planetary and Space Science* 2016;133:85-89. https://doi.org/10.1016/j.pss.2016.04.006

[15] Dominik, C, Blum, J, Cuzzi, JN Wurm, G. Growth of dust as the initial step toward planet formation. in *Protostars and Planets V* 2007, ed. B. Reipurth, D. Jewitt and K. Keil, 783-800. ISBN: 9780816526543

[16] Natta, A, Testi, L, Calvet, N, Henning, Th, Waters, R Wilner, D. Dust in proto-planetary disks: properties and evolution. in *Protostars and Planets V* 2007, ed. B. Reipurth, D. Jewitt and K. Keil, 767-781. ISBN: 9780816526543

[17] Tazaki, R, Tanaka, H, Kataoka, A, Okuzumi, S, Muto T. Unveiling Dust Aggregate Structure in Protoplanetary Disks by Millimeter-wave Scattering Polarization. *The Astrophysical Journal* 2019;885:52-1-13. https://doi.org/10.3847/1538-4357/ab45f0

[18] Duchêne, G, Ménard, F, Stapelfeldt, KR, Villenave, L, Wolff, SG, Perrin, MD, Pinte, C, Tazaki, R, Padgett, DL. JWST Imaging of Edge-on Protoplanetary Disks. I. Fully Vertically Mixed 10 µm Grains in the Outer Regions of a 1000 au Disk. *The Astronomical Journal* 2024;167:77-1-17. https://doi.org/10.3847/1538-3881/acf9a7

[19] Tazaki, R, Lazarian, A, Nomura H. Radiative Grain Alignment in Protoplanetary Disks: Implications for Polarimetric Observations. *The Astrophysical Journal* 2017;839:56-1-17. https://doi.org/10.3847/1538-4357/839/1/56

[20] Tazaki, R, Tanaka, H, Muto, T, Kataoka, A, Okuzumi, S. Effect of dust size and structure on scattered-light images of protoplanetary disks. *Monthly Notices of the Royal Astronomical Society* 2019;485:4951-4966. https://doi.org/10.1093/mnras/stz662

[21] Tazaki, R, Ménard, F, Duchêne, G, Villenave, M, Ribas, Á, Stapelfeldt, KR, Perrin, MD, Pinte, C, Wolff, SG, Padgett, DL, Ma, J, Martinien, L, Roumesy M. JWST Imaging of Edge-on Protoplanetary Disks. IV. Mid-infrared Dust Scattering in the HH 30 disk. *arXiv:2412.07523v2* (2025)

[22] Désert, FX, Boulanger, F, Puget, JL. Interstellar Dust Models for Extinction and Emission. *Astronomy and Astrophysics* 1990;237:215-236.

[23] Tielens, AGGM. The large and the small of interstellar dust. *A&A* 2009;500:335-336. https://doi.org/10.1051/0004-6361/200912167 II. The opacity of the dust aggregate size distribution. *A & A* 2011;532:A43-1-10. https://doi.org/10.1051/0004-6361/201117058

[24] Ormel, CW, Peszun, D, Dominik, C, Tielens, AGGM. Dust coagulation and fragmentation in molecular clouds. I. How collisions between dust aggregates alter the dust size distribution. *A&A* 2009;502:845-869. https://doi.org/10.1051/0004-6361/200811158

[25] Liu, C, Yin, Y, Hu, F, Jin, H, Sorensen, C M. The Effects of Monomer Size Distribution on the Radiative Properties of Black Carbon Aggregates. *Aerosol Science and Technology* 2015;49:928-940. https://doi.org/10.1080/02786826.2015.1085953

[26] Tazaki, R, Ginski, C, Dominik, C. Fractal Aggregates of Submicron-sized Grains in the Young Planet-forming Disk around IM Lup. *The Astrophysical Journal Letters* 2023;944:L43-1-14. https://doi.org/10.3847/2041-8213/acb824

[27] Teng, S, Liu, C, Schnaiter, M, Chakrabarty, RK, Liu, F. Accounting for the effects of nonideal minor structures on the optical properties of black carbon aerosols. *Atmos. Chem. Phys. Discuss.* 2019;19:2917–2931. https://doi.org/10.5194/acp-19-2917-2019

[28] Farias, TL, Köylü, ÜÖ, Carvalo, MG. Effects of polydispersity of aggregates and primary particles on radiative properties of simulated soot. *Journal of Quantitative Spectroscopy & Radiative Transfer* 1996;54:357-371. https://doi.org/10.1016/0022-4073(95)00166-2
[29] Liu, C, Yin, Y, Hu, F, Jin, H, Sorensen, CM. The Effects of Monomer Size Distribution on the Radiative Properties of Black Carbon Aggregates. *Aerosol Science and Technology* 2015;49:928-940. https://doi.org/10.1080/02786826.2015.1085953

[30] Doner, N, Liu, F. Impact of morphology on the radiative properties of fractal soot aggregates. *Journal of Quantitative Spectroscopy & Radiative Transfer* 2017;187:10-19. https://doi.org/10.1016/j.jqsrt.2016.09.005

[31] Yon, J, Liu, F, Morán Cofre, JC, Fuentes, A. Impact of the primary particle polydispersity on the radiative properties of soot aggregates. *Proceedings of the Combustion Institute* 2019;37:1151-1159. https://doi.org/10.1016/j.proci.2018.07.065

[32] Muinonen, K, Zubko, E. Optimizing the discrete-dipole approximation for sequences of scatterers with identical shapes but differing sizes or refractive indices. *Journal of Quantitative Spectroscopy & Radiative Transfer* 2006;100:288-294. https://doi.org/10.1016/j.jqsrt.2005.11.045

[33] Perry, JD, Gostomski, E, Matthews, LS, Hyde, TW. The influence of monomer shape on aggregate morphologies. *A&A* 2012;539:A99-1-7. http://dx.doi.org/10.1051/0004-6361/201117940

[34] Draine BT. The Discrete Dipole Approximation for Light Scattering by Irregular Targets. in: *Light Scattering by Nonspherical Particles, Theory, Measurements, and Applications* (2000), Mishchenko M. I., Hovenier, J. W., Travis, L. D., editors. New York: Academic Press, 131-145. https://doi.org/10.1088/0957-0233/11/12/705

[35] Mishchenko MI, Travis LD, Mackowski DW. T-matrix computations of light scattering by nonspherical particles: a review. *Journal of Quantitative Spectroscopy & Radiative Transfer* 1996;55:535-575. https://doi.org/10.1016/0022-4073(96)00002-7

[36] Mackowski, DW, Kolokolova, L. Application of the multiple sphere superposition solution to large-scale systems of spheres via an accelerated algorithm. *Journal of Quantitative Spectroscopy & Radiative Transfer* 2022;287:108221. https://doi.org/10.1016/j.jqsrt.2022.108221

[37] Berry, M, Percival, IC. Optics of Fractal Clusters Such as Smoke *Optica Acta* 1986;33:577-591. https://doi.org/10.1080/713821987

[38] Botet, R, Rannou, P, Cabane, M. Mean-field approximation of Mie scattering by fractal aggregates of identical spheres. *Applied Optics* 1997;36:8791-8797. https://doi.org/10.1364/AO.36.008791

[39] Kadanoff, LP. More is the Same: Phase Transitions and Mean Field Theories. *Journal of Statistical Physics* 2009;137:777-797. https://doi.org/10.1007/s10955-009-9814-1

[40] Botet, R, Rannou, P, Tazaki, R. The importance of the local structure of fractal aggregates. *Journal of Physica A: Math. and General* 2024;57:115001. https://doi.org/10.1088/1751-8121/ad2c82

[41] Rannou, P, Botet, R, Tazaki, R. Improvement of the Mean Field T-matrix method for scattering by fractal aggregates of identical spheres in astrophysical environments. *Icarus* 2024;424:116247. https://doi.org/10.1016/j.icarus.2024.116247

[42] Bohren, CF, Huffman, DR. *Absorption and Scattering of Light by Small Particles* (1983). Wiley Interscience, New York. https://doi.org/10.1002/9783527618156

[43] Mackowski, DW. Calculation of total cross-sections of. multiple sphere clusters. *Journal of the Optical Society of America A* 1994;11:2851-2861. https://doi.org/10.1364/JOSAA.11.002851

[44] Mishchenko, MI, Hovenier JW, Travis, LD. *Light Scattering by Nonspherical Particles – Theory, Measurements and Applications* (2000). Academic Press, San Diego. https://doi.org/10.1088/0957-0233/11/12/705

[45] Arnoldus, HF. Mie scattering near the Fröhlich mode *Journal of the Optical Society of America A* 2025;42:580-586. https://doi.org/10.1364/JOSAA.558101

[46] Saxon, DS. 1955 Lectures on the scattering of light. University California at Los Angeles, Department of Meteorology, Scientific Report no. 9.

[47] Bruning, JH, Yo, YT. Multiple Scattering of EM Waves by Spheres Part I-Multipole Expansion and Ray-Optical Solutions. *IEEE Transactions on Antennas and Propagation* 1971;19:378-390. https://doi.org/10.1109/TAP.1971.1139944

[48] Mackowski, DW. Analysis of radiative scattering for multiple sphere configurations. *Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences* 1991;433:599. https://doi.org/10.1098/rspa.1991.0066

[49] Cruzan, OR. Translational addition theorems for spherical vector wave functions. *Quarterly of Applied Mathematics* 1962;20:33-40. https://doi.org/10.1090/qam/132851

[50] Mackowski, D, Kolokolova, L, Sparks, W. T-matrix approach to calculating circular polarization of aggregates made of optically active materials. *JQSRT* 2011;112:1726-1732. https://doi.org/10.1016/j.jqsrt.2011.02.003

[51] Xu, YL. Electromagnetic scattering by an aggregate of spheres: asymmetry parameter. *Physics Letters* 1998;A 249:30-36. https://doi.org/10.1016/S0375-9601(98)00708-7

[52] Nguyen, T D, Tanaka M, Tanaka, K. Analysis and computation of electromagnetic scattering by multiple spheres configuration. in *International Conference on Communications and Electronics* 2010, Nha Trang, Vietnam, 2010, 335-340. https://doi.org/10.1109/ICCE.2010.5670658.

[53] Wiscombe, WJ. Improved Mie scattering algorithms. *Applied Optics* 1980;19:1505-1509. https://doi.org/10.1364/AO.19.001505

[54] Mishchenko, MI, Travis, LD, Lacis, AA. *Scattering, absorption, and emission of light by small particles* (2002), Cambridge University Press. ISBN 0 521 78252

[55] Zhang, Y, Ding, J, Yang, P, Panetta, RL. Vector spherical wave function truncation in the invariant imbedding T-matrix method. *Optics Express* 2022;17:30020-30037. https://doi.org/10.1364/OE.459648

[56] van Hulle, P, Weille, M-E, Talbaut, M, Coppalle, A. Comparison of Numerical Studies Characterizing Optical Properties of Soot Aggregates for Improved EXSCA Measurements. *Part. Part. Syst. Charact.* 2002;19:47-57. https://doi.org/10.1002/1521-4117(200204)19:1<47::AID-PPSC47>3.0.CO;2-W

[57] Xu, YL. Efficient Evaluation of Vector Translation Coefficients in Multiparticle Light-Scattering Theories. *Journal of Computational Physics* 1998;139:137-165. https://doi.org/10.1006/jcph.1997.5867

[58] Bushell, G, Amal, R. Fractal Aggregates of Polydisperse Particles. *Journal of Colloid and Interface Science* 1998;205:459-469. https://doi.org/10.1006/jcis.1998.5667

[59] Kotlarchyk, M, Stephens, RB, Huang, JS. Study of Schultz distribution to model polydispersity of microemulsion droplets. *The Journal of Physical Chemistry* 1988;92:1533-1538. https://doi.org/10.1021/j100317a032

[60] Wu, MK, Friedlander, SK. Enhanced power law agglomerate growth in the free molecule regime. *J. Aerosol Sci.* 1993;24:173-182. https://doi.org/10.1016/0021-8502(93)90002-Q

[61] Berthier, B, Chaudhuri, P, Coulais, C, Dauchot, O, Sollich, P. Suppressed Compressibility at Large Scale in Jammed Packings of Size-Disperse
Spheres. Phys. Rev. Lett. 2011;106:120601, 1-4. https://doi.org/10.1103/PhysRevLett.106.120601

[62] Haan, L., Ferreira, A. *Extreme Value Theory. An Introduction* (2006). Springer Series in Operations Research and Financial Engineering. Springer-Verlag New York. ISBN 144192020X

[63] Tricomi, FG. Asymptotische Eigenschaften der unvollständigen Gammafunktion. *Mathematische Zeitschrift* 1950;53:136-148.
https://doi.org/10.1007/BF01162409

[64] Eggersdorfer, ML, Pratsinis, SE. *Aerosol Sci. Technol.* 2012;46:347–353. https://doi.org/10.1080/02786826.2011.631956

[65] Botet, R, Rannou, P. Light Scattering by Large Particles – Mean-Field T-matrix vs T-matrix method – in preparation

[66] Charalampapotdos, TT. Morphology and Dynamics of Agglomerated Particulate in Combustion Systems Using Light Scattering Techniques. *Prog. Energy Combus. Sci.* 1992;18:13–45. https://doi.org/10.1016/0360-1285(92)90031-U

[67] Brassé, C, Muñoz, O, Coll, P, Raulin, F. Optical constants of Titan aerosols and their tholins analogs: Experimental results and modeling/observational data. *Planetary and Space Science* (2015);109:110-159-174. http://dx.doi.org/10.1016/j.pss.2015.02.012

[68] Sorensen, CM. Light Scattering by Fractal Aggregates: A Review. *Aerosol Science & Technology* 2001;35:648-687.
https://doi.org/10.1080/02786820117868

[69] the asymptotic behaviours of the coefficients $K_{n,n'}^{(1,2)}$ result from the results (beyond the scope of this paper):
$$ \int_0^\pi \sin(2kr \sin(\theta/2)) \Phi_{n,n'}(\theta) \cos(\theta/2) d\theta \sim (2n+1)(2n'+1)(1+(-1)^{\delta_{n,n'}} \cos(2\theta))/(2kr)^2 \text{ and } \int_0^\pi \sin(2kr \sin(\theta/2)) \Psi_{n,n'}(\theta) \cos(\theta/2) d\theta \sim (2n+1)(2n'+1)(1-(-1)^{\delta_{n,n'}} \cos(2\theta))/(2kr)^2 \text{ when } r \to \infty, \text{ with } \delta_{n,n'} \text{ the Kronecker symbol.} $$

[70] the relation $|Δ| < 1$ is a sufficient but not necessary condition for (86) to have a finite solution. See textbooks about summing diverging series, e.g. Bender, CM, Orszag, SA. *Advanced Mathematical Methods for Scientists and Engineers* (1978). McGraw-Hill Kogakusha, Ltd

[71] Tazaki, R. Analytic expressions for geometric cross-sections of fractal dust aggregates. *Monthly Notices of the Royal Astronomical Society* 2021;504:2811–2821. https://doi.org/10.1093/mnras/stab1069

[72] Soewono, A, Rogak, SN. Morphology and Optical Properties of Numerically Simulated Soot Aggregates. *Aerosol Science and Technology* 2013;47:267–274. https://doi.org/10.1080/02786826.2012.749972

[73] Yon, J, Rozé, C, Girasole, T, Coppalle, A, Mées, L. Extension of RDG-FA for Scattering Prediction of Aggregates of Soot Taking into Account Interactions of Large Monomers. *Part. Part. Syst. Charact.* 2008;25:54-67. https://doi.org/10.1002/ppsc.200700011

[74] Multiple-sphere numerical codes are currently available at https://scattport.org/

[75] Liu, F, Smallwood, GJ. Radiative Properties of Numerically Generated Fractal Soot Aggregates: The Importance of Configuration Averaging. *Journal of Heat Transfer* 2010;132:023308-1-6. https://doi.org/10.1115/1.4000245

[76] Meakin, P. Fractal aggregates in Geophysics. *Reviews of Geophysics* 1991;29:317-354. https://doi.org/10.1029/91RG00688

[77] Lambert, S, Thill, A, Ginestet, P, Audic, JM, Bottero, JY. Structural Interpretations of Static Light Scattering Patterns of Fractal Aggregates: I. Introduction of a Mean Optical Index: Numerical Simulations. *Journal of Colloid and Interface Science* (2000);228: 379-385.
https://doi.org/10.1006/jcis.2000.6965

[78] Thill, A, Lambert, S, Moustier, S, Ginestet, P, Audic, JM, Bottero, JY. Structural Interpretations of Static Light Scattering Patterns of Fractal Aggregates: II. Experimental Study. *Journal of Colloid and Interface Science* 2000;228:386-392. https://doi.org/10.1006/jcis.2000.6836

[79] Soos, M, Lattuada, M, Sefcik, J. Interpretation of Light Scattering and Turbidity Measurements in Aggregated Systems: Effect of Intra-Cluster Multiple-Light Scattering. *J. Phys. Chem. B.* 2009;113:14962-14970. https://doi.org/10.1021/jp907284t

[80] Lattuada, M, Ehrl, L. Scattering properties of dense clusters of colloidal nanoparticles. *J. Phys. Chem. B.* 2009;113:5938-5850.
https://doi.org/10.1021/jp810976z

[81] Rannou, P, McKay, CP, Botet, R, Cabane, M. Semi-empirical model of absorption and scattering by isotropic fractal aggregates of spheres. *Planetary and Space Science* 1999;47:385-396. https://doi.org/10.1016/S0032-0633(99)00007-0

[82] Rannou, P, Cabane, M, Botet, R, Chassefière, E. A new interpretation of scattered light measurements at Titan's limb. *Journal of Geophysical Research: Planets* 1997;102:10875-11026. https://doi.org/10.1029/97JE00719

[83] Tazaki, T, Tanaka, H. Light Scattering by Fractal Dust Aggregates. II. Opacity and Asymmetry Parameter. *The Astrophysical Journal* 2018;860:79. https://doi.org/10.3847/1538-4357/aac32d

[84] Lodge, MG, Wakeford, HR, Leinhardt, ZM. Aerosols are not spherical cows: using discrete dipole approximation to model the properties of fractal particles. *Monthly Notices of the Royal Astronomical Society* 2024;527:11113–11137. https://doi.org/10.1093/mnras/stad3743

[85] Xu, YL. Electromagnetic scattering by an aggregate of spheres. *Applied Optics* 1995;34:4573-4588. https://doi.org/10.1364/AO.34.004573

[86] Tomchuk, OV, Avdeev, MV, Bulavin, LA. Modeling fractal aggregates of polydisperse particles with tunable dimension. *Colloids and Surfaces A* 2020;605:125331. https://doi.org/10.1016/j.colsurfa.2020.125331

[87] Abramowitz, M, Stegun, IA. *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables* (1972). US Government Printing Office, Washington, 10th printing. ISBN 0486612724

[88] Lorch, L. Alternative proof of a sharpened form of Bernstein's inequality for Legendre polynomials. *Applicable Analysis* 1983;14: 237-240.
https://doi.org/10.1080/00036818308839426

[89] Heinson, WR, Sorensen, CM, Chakrabarti, A. Computer Simulation of Aggregation with Consecutive Coalescence and Non-Coalescence Stages in Aerosols. *Aerosol Science and Technology* 2010;44:380-387. https://doi.org/10.1080/02786821003662942

Declaration of interests

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.