An Accelerated Method for Investigating Spectral Properties of Dynamically Evolving Nanostructures

The discrete-dipole approximation (DDA) is widely applied to study the spectral properties of plasmonic nanostructures. However, the high computational cost limits the application of DDA in static geometries, making it impractical for investigating spectral properties during structural transformations. Here we developed an efficient method to simulate spectra of dynamically evolving structures by formulating an iterative calculation process based on the rank-one decomposition of matrices and DDA. By representing structural transformation as the change of dipoles and their properties, the updated polarizations can be computed efficiently. The improvement in computational efficiency was benchmarked, demonstrating up to several hundred times acceleration for a system comprising ca. 4000 dipoles. The rank-one decomposition accelerated DDA method (RD-DDA) can be used directly to investigate the optical properties of nanostructural transformations defined by atomic- or continuum-scale processes, which is essential for understanding the growth mechanisms of nanoparticles and algorithm-driven structural optimization toward enhanced optical properties.


.1. Method
In the discrete-dipole approximation (DDA) [1][2][3] , the geometry is discretized into polarizable cubic lattices which represent the point dipoles. Every point dipole's polarizability is associated with its local dielectric constant. The polarization of the ℎ dipole is induced by the electric field from the incident beam and also the rest ( − 1) dipoles (Eq. (1)-(2)).
= (1) where is the polarization of the ℎ dipole, is the ℎ dipole's polarizability and denotes the electric field at the ℎ dipole's position which is defined in Eq. (2): where , is the electric field of the incident beam at the location and is the contribution from the ℎ dipole located at . They can be formulated as below: where is the imaginary number, E is the amplitude of the incident beam, is the wavenumber with the same direction of the wave propagation and its amplitude is defined as = 2 , where is the wavelength of the incident beam, , ≡ − , , ≡ | , |, and is the polarization of the ℎ dipole.
To calculate the polarizations in a self-consistent manner, the system can be described by simplifying the Maxwell equations into a set of linear equations (Eq. (5)).
where is a 3 vector describing the local electric field of the incident wave in every dipole position, is a 3 × 3 symmetric matrix depending on the geometry and materials of the geometry. is a For any matrix , at the diagonal position, it is defined by the inverse of the polarizabilities (α , , α , , α , ) of the ℎ dipole in X, Y and Z directions (Eq. (7)). For simplicity, it is generally assumed that the dipole is isotropic in the X, Y and Z directions, so that α , = α , = α , = α .
For the matrix , ≠ that is off-diagonal, it is defined by the relative distance between the ℎ and ℎ dipoles (Eq. (8)).
where ̂, = − | , | and is the 3 × 3 identity matrix. It should be noted these off-diagonal matrices are independent of the polarizability of the dipoles, and thus are kept constant when the properties of the dipoles are changed.
After calculating , we can further evaluate the extinction and absorption cross-sections via the following Eq. (9) and Eq. (10), respectively. C = 4 |E | 2 ∑ Im( , * • ) (9) where C is the extinction cross-section, Im( ) denotes the imaginary part of and * is the conjugate of .
where C is the absorption cross-section and α is the polarizability of the ℎ dipole.
where α , is the Clausius-Mossotti polarizability 2 defined as below: α , = 3 3 4 ( 2 − 1) where is the complex refractive index of the ℎ dipole and is the dipole length. The term in Eq. (11) is defined as: where is the dipole length. If not mentioned, the refractive indexes used in this work were all from reference 6 .
Once C and C are calculated, the extinction and absorption efficiency factors can be calculated from the corresponding cross-sections via Eq. (14).
where is the effective radius of a sphere with the same volume as the geometry, which can be defined by During the study of the spectral properties of nanostructures, the extinction efficiency factors by uniformly sampling 1000 possible orientations of the simulated nanoparticle with respect to the incident beam are calculated and averaged to give the final UV-Vis spectrum. To enable accuracy, the validation criteria (Eq. (15)) should always be satisfied.
where is the complex refractive index of the material, is the wavenumber of the incident beam and is the dipole length.
Based on the DDA described above, how the polarizations will change when the replacement/addition/removal of the dipoles in the geometry happens will be discussed below, as indicated in Figure S1. Figure S1. The scheme for the discrete-dipole approximation. The original geometry is represented by a series of discrete dipoles. The replacement (indicated by changing the colour of a dipole from yellow to blue), addition or removal (indicated by making one dipole transparent) of the dipoles from the system will change the polarizations of the dipoles.

Formulating new solutions when the polarizability of a dipole is changed
A replacement reaction within the nanostructure can be represented by the change in the properties of the dipoles. Here, the relationship between the new solution of dipole polarizations with the original solution before the replacement is derived. Eq. (5) can be reformulated as below to calculate the polarizations: where − is the inverse of . Changing the property of a dipole is equivalent to updating its polarizability tensor , where the diagonal elements of are updated. The off-diagonal elements are kept the same as indicated by Eq. (8). The diagonal part of that is only dependent on the polarizability and regardless of the geometry can be expressed as follows: The off-diagonal blocks are independent of polarizability tensors, thus they are not shown here. Also, S7 for simplicity, it is assumed the dipoles are isotropic in the following derivation. However, the formulation is general and similar expressions can be derived for anisotropic dipoles where only the diagonal elements exist but their values are different.
Next, assume the property of dipole is changed so that it has a new polarizability α ′ , then the new matrix of ′ will be: where α ′ is the new polarizability of dipole .
Calculating the inverse matrix of ′ directly is straightforward but computationally inefficient.
Instead, a method which can analytically describe the relationship between the original solution and the updated solution should be derived. The new inverse matrix ′ − would be written as a function of the original inverse matrix − . To do so, Lemma 1 7 was introduced: Lemma 1. If A and A+B are invertible, and B has rank one, then let g = tr( − ). If g ≠ −1, we By comparing Eq. (17) and (18), we can write: By writing so, we can define ′ as the summation of and three rank-one matrices: where , and correspond to the sparse matrices with only one non-zero element of α ′ −1 − α −1 in Eq. (19). For anisotropic dipoles where only the diagonal elements exist but with different values, the elements of , and should be calculated in every corresponding direction of X, Y and Z respectively.
With Lemma 1, the next step is to add the first rank-one matrix ( ), whose physical significance is to change the polarizability of the ℎ dipole from α to α ′ in the X direction. Its corresponding inverse matrix ( ′ ) − can be calculated as: where, ( − ) , is the element of − at the ℎ row and the ℎ column. Then for the new solution to the polarizations of dipoles ( ′ ) , it can be written as: where is a vector of 3N elements that defines the polarizations of the dipoles in the original system S10 and 3 −2 indicates the (3 − 2) ℎ element in the vector .
The change in polarizability tensor along Y and Z directions can be introduced similarly with varied indexes, and solutions for the system with the change of polarizabilities in multiple dipoles can be obtained by applying the same procedure as above iteratively.
By adding a rank-one matrix and calculating the new inverse matrix and dipole solution according to the procedure described above, we can solve the system where the property of either a single dipole or a set of dipoles are changed. Here, we only need to calculate the initial inverse matrix and the rest is only matrix multiplication, which is less computationally expensive compared to solving a linear system and can be accelerated using GPU/multiple CPUs. Furthermore, the approach described above not only can find a new solution to the system when the polarizability of dipole is changed, but also can be used as a valid approximation to calculate the new solution of the system when the addition or removal of dipoles happens, which will be discussed later.

Problem
The growth of a nanostructure can be represented by the addition of new dipoles to the system. Here, the relationship between the new solution to the polarizations and the original solution is derived using the rank-one decomposition when multiple dipoles are added to the system. Again, it is assumed the original linear system was solved to obtain the corresponding initial − and . When a new dipole is added, the new inverse matrix ( ′ − ) can be calculated from them and used to evaluate the new polarizations of the dipoles ( ′ ).
First, when a dipole is added to the system, the corresponding matrix for the linear coefficients in the new system was denoted as ′ and can be written as: where + is the polarizability tensor of the newly added dipole and is the 3 × 3 matrix which is constructed with the original dipoles. ′ +1, is the 3 × 3 matrix that describes the dipole interaction between the newly added ( + 1) ℎ dipole and the ℎ dipole (see the general description of DDA and Eq. (7)). Note here the property that both and ′ are symmetric are used.
The influence of the new dipole on the original system is added sequentially for its polarizability in the X, Y and Z directions. The procedure to evaluate the solution of the new system after adding a dipole goes as follows: 1. A (3 + 1) × (3 + 1) matrix ′ is defined, which takes the top-left (3 + 1) × (3 + 1) block of ′ , and its inverse matrix ( ′ − ) will be evaluated based on − .
3. Similar to step 2, once ′ , − is solved, the inverse of ′ will be evaluated based on it.
This procedure expands the original matrix to ′ by adding a single row and column every time and calculating their corresponding inverse matrices through matrix rank-one decomposition. The corresponding ′ , ′ , , and ′ are defined below: By adding an extra column and row to the original matrix, the relation between the original inverse matrix and the current inverse matrix will be derived below. Once the relation is derived, it can be applied repeatedly to sequentially solve ′ − , ′ , − , and ′ − starting from − .

Method
We will begin to derive ′ − based on − , but it should be noted the same procedure can be used to calculate the rest inverse matrices ( ′ , − and ′ − ). ′ is defined in Eq. (26), which can be decomposed into the summation of ′ , and two rank-one matrices defined below: where its last row and column are filled with 0s except the diagonal one. Its inverse matrix can be written explicitly based on − : The corresponding solution of dipoles, in this case, can also be written as: With this auxiliary matrix, ′ can be decomposed to the auxiliary matrix and two rank-one matrices: where ( × ) indicates a × block matrix composed of only 0s.

S14
For simplicity and generality, we denote the extra row and column except for the last element as , so that: Then after applying Lemma 1 twice, we will be able to evaluate ′ − based on ′ , , which can be explicitly expressed by − (Eq. (30)) below.
First, the inverse matrix of ′ ] is calculated: where C k is the ℎ element of vector , which is defined as below: Note we use the property that − is symmetric so that − = ( − ) .
The solution of the dipoles ( ′ , ) corresponding to ′ , can be calculated via: where D is defined as D = • and α +1 E +1, − α +1 • is the last element of ′ , .
Thus, the inverse matrix and its corresponding dipole distribution of ′ , which has an extra column and row compared to , can be solved. The relation defined by Eq. (34)-(38) is general for any symmetric matrices of and ′ , where ′ has an extra column and row compared to .
By applying this relation three times, the inverse matrix of ′ and its corresponding dipole distribution ′ can be obtained. In the case of adding more dipoles, the linear system can be solved by applying the relation multiple times iteratively.

Problem
In the previous section, we formulate the change in the polarizations when dipoles are added to the system, which represents nanostructure growth. Similarly, with the etching of the nanostructure, the equivalent process can be represented by removing the old dipoles from the system. So, through the matrix rank-one decomposition, the relationship between the new solution of polarizations and the original solution in an etching process can be derived. Assuming the original linear system including the old dipoles is solved to obtain the corresponding − and . Here we assume there is ( + 1) dipoles in the original system, with one dipole to be removed. When one dipole is removed, the new inverse matrix ( ′ − ) can be calculated from − and used to evaluate the new solution of dipoles ( ′ ). Thus, the size of is 3( + 1) × 3( + 1), while the size of ′ is 3 × 3 .
3. Similar to step 2, once ′ , − is solved, the inverse of ′ will be evaluated based on it.
This procedure shrinks the original matrix to ′ by deleting a single row and column every time and calculating their corresponding inverse matrices through matrix rank-one decomposition. The corresponding ′ , ′ , , and ′ are defined below: By removing a column and row from the original matrix, the current inverse matrix can be expressed as a function of the original inverse matrix and will be derived. Once the relation is derived, we apply it repeatedly to sequentially find ′ − , ′ , − , and ′ − starting from − .

Method
We start to derive ′ − based on − . The relation between ′ and is shown in Eq. (42), where ′ is a (3 + 2) × (3 + 2) matrix. The inverse matrix ′ − will be calculated through an auxiliary matrix, which does not consider the interaction between the ( + 1) ℎ dipole and the rest of the dipoles: Its inverse matrix can be evaluated directly as: can be expressed as the summation of and two rank-one matrices, so that: where is the negative of the last row without the last element in . Following the same procedure described in the last section, the ′ , can be evaluated as follows: ) is defined as: As indicated by Eq. (34), the inverse matrix of ′ , has the form as follows: where X is unknown constants that will be solved but not be used further. So, we only need to S19 calculate ′ , − and take the top-left (3 + 2) × (3 + 2) block to get ′ − , instead of actually calculating ′ , − further by adding a second rank-one matrix.
Again, ′ , − can be evaluated via Lemma 1: As an example, the equivalence between the direct implementation of DDA and iterative rank-one decomposition accelerated DDA method (RD-DDA) has been shown for a limited number of dipoles, which is available at https://github.com/croningp/RD-DDA.

The validation and numerical approximation for replacement, growth, and etching
As described above, we can solve the exact solution when a replacement, addition or removal of dipoles happens through the rank-one decomposition method. However, considering the growth/etching process, the size of the coefficient matrix (A) will increase/decrease at each step, this makes the algorithmic implementation computationally less efficient. Instead, the approach described in Eq.(21)-(22) offers a computationally efficient method for replacement, growth and etching which will be discussed below.
First, by defining the outer bounds, all the possible positions where the dipoles can exist are considered to create the initial matrix. The absence of the ℎ dipole in a specific position can be approximated by setting its polarizability as , where ≈ 0. During the growth process, when the nanostructure growth front reaches the position of the ℎ dipole, its polarizability can be changed from to α . Similarly, during the etching process, when the nanostructure etching front reaches the ℎ dipole happens, its polarizability can be changed from α to . In both cases, the system can be solved iteratively according to Eq. (21)-(22).
Here, for good numerical stability in estimating the polarizations of dipoles, the selection of is crucial. A series of s in the case of the growth of Au cubes (10 × 10 × 10 nm 3 ) to Au nanorods (10 × 10 × 30 nm 3 ) in the water medium, with a dipole length of 1 nm, were benchmarked. The UV-Vis spectra from the exact solution were calculated from the method described in Section 1.3 as S20 standard data, where the initial and final dipole numbers were 1000 and 3000 respectively. The complex refractive index data for Au was used based on reference 6 . The calculated extinction efficiency factors (Q ) with layer wise growth (sequential dipole addition to complete one layer) are shown in Figure S2.

Figure S2. The simulated UV-Vis (extinction efficiency factors) in a growth trajectory.
The initial geometry is a cube with a 10 nm length, and the final geometry is a rod with the size of 10 × 10 × 30 nm 3 . The growth happens by adding dipoles of 1 nm in length to the initial cube. The blue-to-red transition indicates the growth of the formation of nanorods. The results shown in the figure are sampled after adding a complete layer composed of 100 dipoles.
Then, the same well-defined growth trajectory can be simulated via two possible strategies: Strategy 1: Select Au cube (10 × 10 × 10 nm 3 ) as the initial geometry with outer bounds defined by nanorod geometry (10 × 10 × 30 nm 3 ), and the UV-Vis spectra from the same growth trajectory were simulated using a small enough for the non-Au dipoles.

Strategy 2:
Reverse the growth trajectory to create an "etching trajectory". The corresponding UV-Vis spectra from this etching trajectory with a small enough value of for non-Au dipoles can also be simulated. Reverse the sequence of the UV-Vis spectra will give the UV-Vis change in the growth process.
By changing the polarizability of the dipoles from to α, the UV-Vis spectra of such growth or etching trajectory can be created. Imaging there is one dipole composed of the medium, whose polarizabilities is ≈ 0, we can write matrix as below: When the polarizability of the ℎ dipole is changed from to α, there would be a factor defined as (21) and (22), which can be calculated as: Given → 0, we should estimate the value of this term to check if it is divergent or convergent to evaluate the numerical stability. We need to estimate the value of ( − ) 3 −2,3 −2 , which is the It can be shown ( − ) 3 −2,3 −2 can be approximated as when → 0 as follows: First, a matrix labelled as was constructed by sorting the existing dipoles so that the dipoles purely composed of Au are indexed lower and those of the medium are indexed higher. We assume there are dipoles in total and of them are composed of pure Au and ( − ) composed of the medium, then the matrix can be written as: where , is the 3 × 3 matrix that describes the interaction between the ℎ and ℎ dipole as indicated in Eq. (7). We define the top-left 3 × 3 block of as ′ , which consists of × block matrices. The shape of these small block matrices is 3 × 3. ′ is identical to an matrix constructed from only considering Au dipoles listed in the same order.
For any matrix of ( ) , where > , the diagonal element is −1 , with approximating 0. Later, the inverse matrix of can be estimated as follows: 1. We rewrite as: where ′′ is a 3( − ) × 3( − ) matrix identical to an matrix constructed from medium dipoles listed in the same order, whose diagonal elements are all equal to −1 . and consider the interaction terms between the Au dipole set and the medium dipole set.
2. According to the inverse of the block matrix, − can be written as: here, we need to calculate the diagonal elements of this matrix. Since approximates 0, the off-diagonal elements of ( ′′ − ′ − ) are infinitely small compared to −1 , while the diagonal elements can be approximated by −1 . Thus, we can write: where is the 3 × 3 identity matrix. Based on it, we can have: It indicates that any diagonal element that corresponds to the dipole composed of medium and whose index is larger than 3 in − , can be approximated by .
The original matrix can be constructed by swapping the columns and rows of . During this transformation, the inverse matrix − can also be constructed by swapping the corresponding rows and columns in − . If the ℎ dipole that constructed matrix is composed of the medium, we have: Thus, ( − ) 3 −2,3 −2 that correspond to a dipole composed of the medium will be , which indicates Eq. (52) can be approximated as: which approximates infinity intrinsically and gives a large numerical error when is too small and the numeric precision is not high enough. It should be noted that it is just an estimation of when the polarizability of the ℎ dipole in the X direction is changed from to α. Later, it can be observed that −1 − α −2 is good to estimate the order of magnitude of in the subsequent steps: their orders go up similarly when decreases.
Additionally, the magnitude of increases when smaller is used.
However, for strategy 2 where a dipole is removed, we have α = α and α ′ = . The value of ( − ) 3 −2,3 −2 is not as small as , which means can be approximated as: where 1 ( −1 −α −1 ) is neglected for → 0, which ensures numerical stability during the etching process.
In the benchmark, we found out that the inherent numerical errors can negatively impact the first strategy, as suggested by Eq. (60) and will be discussed later in detail. UV-Vis spectra were generated through strategy 2, which means during the calculation, the actual initial geometry is a 10 × 10 × 30 nm 3 rod with 3000 dipoles. The calculations were implemented with 64-digit numeric precision.
In strategy 2, we created a system composed of all 3000 dipoles representing the nanorod. As part of the etching process, we sequentially remove dipoles from the system by changing its polarizability α to , which should be close to 0 compared to the original polarizability (α) for Au dipole. The polarizations of dipoles where different was used to represent the polarizability of the medium ( Figure S3a) were solved via Eq. (21)-(22), which give the extinction efficiency factors (Q ) again.
Then we calculated the absolute difference between the efficiency factors solved from applying different s and the standard data. The relative absolute difference is further calculated by dividing absolute difference with the standard data. The average of the relative absolute difference among multiple wavelengths from the trajectory is shown in Figure S3b. When is equal or smaller than 10 −4 α, the calculated Q converged to the standard data, which proves that Eq. (21)-(22) can be used as an efficient numerical approximation given that is smaller as compared to α (~10 −4 α or less).
Strategy 1 was implemented by initialising the system using various s as the polarizability of the medium. Sequentially, after changing the polarizability from a small value for to a constant value α for Au in the growth process, the numerical instability of the UV-Vis spectra was observed for a given small , which is shown in Figure S4a. For ≤ 10 −4 α, such numerical instability leads to large deviations from the expected values and even unphysical NaN (infinity) values due to the presence of large values during the calculation. To dampen such problems, we reduce the order of S25 magnitude of by multiplying it by a scale factor, so that the magnitude of the elements of is closer to 1. Here the scale factor ( ) is set as: where is the dipole length in the unit of meter, and is a constant which we set as 10 19 / 3 considering the dipole length is in nanoscale to make as close to 1 as possible. By introducing such a scale factor, the UV-Vis spectra from the same trajectory with various s were simulated again, as shown in Figure S4b. For comparison, the errors before and after rescaling is shown in Figure   S4c-h. Since there is no numerical error for ≥ 10 −2 α, rescaling the system does not influence the results (Figure S4c-d). By rescaling the system, the numerical instability problem can be damped.
Alternatively, we can increase the numeric precision during the calculation. By using 128-digit precision, the UV-Vis spectra in the growth process were simulated again. The numerical instability problem was not observed, as seen in Figure S5a. By applying a scale factor ( ), the results did not show obvious difference, as shown in Figure S5b-h. Only when was set as 10 −10 α , small numerical instability is observed and applying the scale factor can slightly influence the results Figure S5h.  In the etching case, = 10 −8 α and 64-digit numeric precision was used in the etching process to enable numerical stability. matrix was not rescaled during the calculation. The calculated factor , with only neglectable differences as shown in Figure   S6a. The index corresponds to the matrix indexes relevant to the 6000 updates needed to etch away the 2000 dipoles in the X, Y and Z directions.
For the growth process, we compared the order of magnitude between and −1 − α −2 when the calculation is numerically stable for a series of s (i.e., using 128-digit numeric precision). Again the index corresponds to the matrix indexes relevant to the 6000 updates needed to grow the 2000 dipoles. was varied from 10 −1 α to 10 −10 α, and its corresponding estimation value ( −1 − α −2 ) as well as the actual value of during the growth were recorded.
It is observed the order of magnitude of increases when smaller is applied and shows similar order of ( −1 − α −2 ), as shown in Figure S6b. (labelled in different colours) and the estimated −1 − α −2 (black lines) to show the similarity of the orders. In both plots, the value is averaged among 51 wavelengths from 400 nm to 900 nm with an interval of 10 nm.
In summary, we have discussed RD-DDA from various perspectives. For a replacement process, the method described in Section 1.2 gives the exact and computationally efficient strategy to obtain the new solutions. For a growth process, Section 1.3 gives the exact solution to track the dipole change but the change of matrix size may dampen computation efficiency. Implementing the method described in Section 1.2 by setting a small polarizability of the medium can avoid the change of the matrix size, but can also cause numerical instability, so higher (128-digit) precision during the calculation is recommended. For an etching process, the exact solution to track the removal of a dipole is discussed in Section 1.4, with the varied matrix size. The method described in Section 1.2 with a small polarizability of the medium can be a good approximation to the exact solutions and ensures numerical stability in the etching process but not the growth process. However, for any growth process, an equivalent etching process can be created by reversing the trajectory, indicating that the system for the growth process can be solved by solving it as an etching process.
By far, we have discussed both general analytical solutions to track the dipole change under replacement, addition, and removal, and their computationally efficient implementation. The computational time from RD-DDA and from directly solving a new system in every step will be estimated, compared, and discussed below. Then we implemented RD-DDA in two different ways with both 64-digit and 128-digit precision: 1. Solving the system analytically by expanding the system as described in Section 1.3, which is labelled as rank-one decomposition iterative solution 1(RS 1).
Create an etching process trajectory by reversing the growth trajectory and removing dipoles by changing the corresponding polarizability from α to (see Section 1.5), which is labelled as rank-one decomposition iterative solution 2 (RS 2). This setting is due to the numerical instability problem described above.
We defined an acceleration factor (F , Eq. (63)) to quantify the computational efficiency of RD-DDA as compared to the direct solution method (implementing the DDA directly):  An additional example to quantify the computational efficiency consists of the growth process of nanospheres from a radius of 10 nm to 11 nm, which was benchmarked similarly as described above.
The sequential addition of 1 nm dipoles occurred in the ascending order of the distance between the  For the case of coating Au octahedra to form Au@Ag nanostructures as discussed above, the comparison of time cost to simulate intermediates using the direct solution method and RD-DDA is shown in Figure S9, with different numerical precisions. When both methods are set in 64-digit precision, the acceleration factor (F , Eq. (63)) is larger than one when the number of intermediates is larger than 18 and 9 for RS 1 and RS 2 respectively. When the methods are set to 128-digit precision, F is larger than one if the number of intermediates exceeds 5 and 4 for RS 1 and RS 2 respectively.
When RD-DDA is in 128-digit precision and the direct solution method is in 64-digit precision, more intermediates are required to make F larger than one, which is 33 and 24 for RS 1 and RS 2 respectively. It should be noted that in this comparison, the total number of the solved systems will be equal to the number of intermediates plus one (which corresponds to the initially solved nanostructure in RD-DDA). For the second case of growing Au nanospheres to larger nanospheres, a similar comparison can be made ( Figure S10). Here, F is larger than one if the number of intermediates exceeds 27 and 14 for RS 1 and RS 2 respectively.

Estimating the time cost of the direct solution method and RD-DDA
In this section, we will discuss a general strategy to estimate the time cost for both methods including the direct solution method and RD-DDA. Since the computation is mainly matrix manipulations, we can assume the time cost for every update in RD-DDA depends mainly on the matrix size (or equivalently, the dipole number). In the meantime, with the benchmark cases using RD-DDA simulation discussed in the previous section, we observed the time cost to obtain one solution recursively after one dipole change is linearly correlated to the number of existing dipoles in the system (equivalently, the matrix size), as shown in Figure S11. The time cost and its linear fitting to perform one update in the growth of Au nanospheres with RS 1.
The same raw data is used in Figure S8.
Thus, by assuming the time cost for one update (after a dipole change) towards a different nanostructure is linearly correlated to the existing dipole number (thus the matrix size), one possible strategy to estimate the time cost for RD-DDA is proposed as follows: 1. If we know the trajectory, we can estimate the lowest and highest time cost for one update by performing the update on the smallest system (with 1 dipoles) and the largest system (with 2. Then we can create a linear curve using ( 1 , 1 )and ( 2 , 2 ). The linear curve will be used to Step (2). With this method, is precise and ∑ is estimated.
We will validate if the linear interpolation strategy is good to estimate the RD-DDA time cost in the two benchmark cases described above. Specifically, we will estimate the time cost (∑ ) for all the possible small trajectories from the whole growth trajectories for validation. These small trajectories include at least two updates to create the linear curve.
On the one hand, we have recorded the actual update time cost for these small trajectories (i.e., the S35 time cost for RS 1 in Figure 1b and Figure S8). On the other hand, we have also recorded the time cost for one update for the smallest and largest structures (two points), which can be used to create a linear curve to estimate the ∑ term as described above. By comparing the actual and estimated time cost for all the small trajectories, the relative error due to the estimation can be calculated (see Figure S12). The maximal relative error from the estimation is 3.77% for the case of the growth of Ag on the surface of Au octahedra and 2.00% for the case of the growth of nanospheres, respectively.

Efficient tracking of UV-Vis spectra using the RD-DDA
In the last section, we formulated a rank-one decomposition method to accelerate the discrete-dipole approximation for dynamically evolving nanostructures. Using this approach, it is possible to efficiently track the spectral properties like UV-Vis spectra of nanostructures for any given morphological change. In this section, we create trajectories defining morphological changes in three different cases: 1. Custom-built trajectories corresponding to different growth modes.

2.
A trajectory generated from an empirical crystallographic surface growth model.

3.
A trajectory generated from an atomic model using the kinetic Monte Carlo simulation.

Custom-built trajectories
In this case, we studied the growth of a thin layer of Ag on the surface of Au octahedra. It was achieved by adding a layer of 1 nm Ag dipoles to the surface of the Au core. The Au core is created by a series of dipoles composed of pure Au atoms with a dipole length of 1 nm. The edge length of the final octahedral Au@Ag core-shell nanostructure is 20 nm. By varying the growth strategies, three types of trajectories were generated.
1. Adding the layer of Ag dipoles randomly (labelled as Random).
2. Adding the layer of Ag dipoles in ascending order of the distance from the centre of the surface (labelled as Centre).
3. Adding the layer of Ag dipoles in descending order of the distance from the centre of the surface (labelled as Tip).
The growth trajectories with random dipole addition were generated 10 times to create statistically significant results. During the growth process, the number of dipoles increased from 3303 to 4089.
An etching trajectory generated by reversing the growth trajectory was used to calculate the polarizations with a 64-digit precision. The polarizability of the medium is set as = 10 −10 α. It should be noted the changed dipoles were purely composed of Ag, so α specifies the polarizability of Ag dipoles in this case.
The corresponding UV-Vis spectra can thus be simulated using RD-DDA (See Figure 2 in the manuscript). The full details of the spectra are available in Supplementary Information (SI) Video S1-S3.

An empirical crystallographic surface growth model
Crystallographic surface growth models are usually empirical but powerful to elucidate the morphological changes during crystal growth. Previously, using Transmission Electron Microscopy S37 (TEM) the morphological transformation of Au arrow-headed nanorods into Au octahedra 9 has been investigated. Inspired by this experimental observation, here we constructed a trajectory using an empirical model by defining the growth of multiple crystal surfaces to describe this transformation.
Once the trajectory was generated, the UV-Vis spectra of the intermediates were calculated.
The initial arrow-headed nanorods were enclosed within two types of crystallographic surfaces: (110) and (111) 9 . We set constant growth rates for these surfaces as 1 and 2 respectively. The region of the nanostructure is defined by the following enclosed planes in the space: where 1 and 2 are the growth rate for (110) and (111) surfaces, 1 and 2 are the initial conditions. is a non-negative integer defining the growth step. We set 1 = 0.0205 / , 2 = 0.01 1 to ensure that we can capture the atomic layer growth precisely and faster growth rate of (110) surface, and 1 = 2.87 and 2 = 6.56 to generate an initial arrow-headed rod shape. The fine sampling in the growth model equivalently represents the atomic layer growth which is shown in Figure S13a-c, where the intermediates after every 80 growth steps together with the initial Au arrow-headed rods and the final octahedra are shown. However, at the atomic scale, simultaneous growth over the complete crystallographic plane represents an ideal process. In general, at the atomic scale, the crystal growth is a stochastic process, where the order of sequential addition of atoms towards a single atomic layer growth to the crystallographic surface of the nanostructure could potentially create different trajectories. Here, we generate various growth trajectories by shuffling the sequence of the addition of atoms for the same layer to investigate the bounds on the spectroscopic deviation due to stochasticity in the growth process at the atomic scale. The way of generating these trajectories is described below: 1. Since we implemented a continuum crystal surface growth model with atomic layer precision, an extra layer of atoms can be added after one step. We labelled the structure before and after adding this layer as and −1 .
2. The atomic growth trajectory was then generated by randomly adding the extra atoms one by one sequentially, which enables the transformation from −1 to . However, the sequence of adding these atoms can be different, which generates different growth trajectories from −1 to . 3. After the addition of a single atom, the influence on the structure's dipole representation should be checked, which can cause the addition of one dipole or dipoles. If a face-centred lattice is fully occupied after adding the atom, a dipole representing this lattice should be added. If the addition of one atom causes the addition of multiple dipoles, they will be added sequentially with a shuffled sequence. Thus, multiple trajectories composed of the sequential addition of dipoles for the transformation from −1 to can be generated, which S39 approximates the atomic growth process.
4. The overall trajectory from 0 to can be generated by combining the trajectories to transform 0 to 1 , then 1 to 2 ,…, and eventually −1 to .
For the RD-DDA calculations, the polarizability of the medium was set as = 10 −10 α with 64-digit numeric precision, where α is the polarizability of a 0.41 nm dipole purely composed of Au. The wavelength was sampled from 400 nm to 900 nm with an interval of 10 nm. We reversed the growth trajectory to avoid the numeric instability problem during the calculation.
In the simulated spectra, the transverse and longitudinal modes of the original Au arrow-headed rods were observed which are typical spectral features for anisotropic nanoparticles. When the arrowheaded rods transformed into the Au octahedra leading to an isotropic shape, these modes merge into a single peak. Their peak prominence also varies over the growth process, as shown in Figure S14.
It is observed that the complete coating of a compact layer of dipoles on the surface can enhance the extinction of the longitudinal mode. The full details of the spectra from one trajectory are available in SI Video S4-S5.

. Au nanostructures
The Monte Carlo simulation to investigate nanostructural transformations is based on the previous study 10 , and only its implementation is been discussed here.
The Monte Carlo method includes two types of possible events: 1. The addition of a single Au atom on the surface vacant sites.
2. The removal of a single Au atom from the surface of the nanostructure.
The probability of these two events depends on the difference between the chemical potentials at the surface and the environment (e.g., the solution phase). Each step in the Monte Carlo Simulation consists of the following sub-steps: 1. Sample the event type with equal probability: addition or removal 2. If an addition event is sampled, select a vacant site from all the available vacant sites. The probability of addition of the atom on the selected site is calculated as: where and are the numbers of vacant/surface sites of the current structure, is the chemical potential difference of the atom between the solution and the surface of the nanostructure, is the binding energy from the neighbouring atom, is the coordination number of the site after adding one atom, and = 1 with is the Boltzmann constant and T is the temperature.
3. If a removal event is sampled, select a surface site from all the available surface atoms. The probability of removal of the atom on the selected surface site is calculated as: 4. Steps 1-3 are repeated until the simulation is finished.
It should be noted that compared to the original method, we used the ratio of the current surface/vacant S41 sites to approximate the detailed balance to accelerate the simulation. In the original work, they used the sites before/after the Monte Carlo steps for the detailed balance. Due to the relatively large numbers of the sites, this approximation holds. The surface sites were found by selecting atoms with coordination numbers smaller than 12, while the surface vacant sites were in the surroundings of existing atoms. The etching process of Au octahedra with an edge length of ~ 9.3 nm was selected as an example. In the simulation, we set = × 300 = 0.0259 , = 0.3275 and = −6 = −1.965 . Approximately two million Monte Carlo steps are required to etch the octahedra into spheres.
Starting with well-defined Au octahedra, the kinetic Monte Carlo simulation was performed three times with the same parameters but with varied random seeds for sampling. An example trajectory, with atomic-scale variations during the Monte Carlo simulation, an equivalent trajectory at dipole scale, and their corresponding UV-Vis change are shown in SI Video S6-S7.
To find the equivalent dipole transformation of the atomic transformation, we first generated a set of dipoles representing the initial Au octahedron. After the Monte Carlo simulation, it was observed that these initial dipoles included all the dipoles involved in the structural transformation. Thus, these initial dipoles were used to calculate the A matrix. Then, we recorded the addition/removal of atoms, and if the dipole was not fully occupied, it was regarded to represent the medium, which generate the trajectory of dipole transformation. During the dipole transformation, when a dipole is removed or added, its polarizability is changed between α and accordingly.
It should be noted here we directly removed or added a dipole, thus the polarizability is changed between α and . However, if the dipole is not fully filled with atoms, a series of intermediate polarizabilities (α ) can be estimated based on the number of atoms within the dipole. The scattering signals when the polarizability of the dipole is changed among α , α and can also be tracked with the proposed method above. This strategy will be discussed in detail when the replacement event happens in the example below (Section 2.3.2).
During the simulation, 128-digit precision is used. = 10 −6 α is used to define the polarizability of the medium. A cubical volume of face-centred cubic (FCC) lattice with a length of 0.41 nm was set as a single dipole, and when the lattice is fully occupied, the polarizability will be set as α, otherwise . Using the RD-DDA simulation, we sampled at the wavelength from 450 nm to 700 nm with an interval of 10 nm, with additional points between 520 nm and 570 nm with an interval of 2 nm to precisely capture the peaks. The simulated UV-Vis spectra are shown in Figure S15. Figure S15. The UV-Vis change from three different trajectories labelled from 1 to 3. The spectra from the trajectory 1, 2 and 3 are shown in (a), (b) and (c) respectively. The three repeats all showed a consistent tendency of the UV-Vis change for the transformation.

Au-M nanostructures
In this section, we extended the kinetic Monte Carlo combined with the RD-DDA approach to multimetallic systems. In the case of a bimetallic system with gold as one metal, we labelled the first atom b. For the replacement event, we assumed the attachment of a specific type of atom to the surface triggered the replacement process, thus , and , are proportional to selecting the type of atom that is originally in the solution phase and will be inserted into the nanostructure, e.g., a higher , will cause more trials to replace the surface site M with Au.

S43
We denote the selected atom type that will be added, removed and replaced as and the unselected atom type as for the following descriptions.
3. If an addition event with atom type is tried, select a vacant site on the surface with uniform distribution. The probability of accepting this event is calculated as: The corresponding parameters are described above.
5. If a replacement event is tried, select one surface site of atom type from the nanostructure to be replaced by atom type . The probability of accepting this event is calculated as: where , and , are the numbers of surface sites with atom type and respectively. The additional parameters are defined so that is the chemical potential difference in the solution and the surface of nanostructure for atom type , − and − are the binding energy of atom type with Au and M, while − and − are the