Study of Morphology and Optical Properties of Gold Nanoparticle

Subscriber access provided by READING UNIV

Interface Components: Nanoparticles, Colloids, Emulsions, Surfactants, Proteins, Polymers

Study of morphology and optical properties of gold nanoparticle aggregates under different pH conditions Dezheng Darson Li, Xi Gu, Victoria Timchenko, Qing Nian Chan, Anthony Chin Yin Yuen, and Guan-Heng Yeoh Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.8b01457 • Publication Date (Web): 04 Aug 2018 Downloaded from http://pubs.acs.org on August 4, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

RTICLE

Study of morphology and optical properties of gold nanoparticle aggregates under different pH conditions D. Darson Li *†, Xi Gu †, Victoria Timchenko †, Qing N. Chan †, Anthony C. Y. Yuen †,and Guan H. Yeoh †‡ †

School of Mechanical and Manufacturing Engineering, the University of New South Wales, UNSW Sydney, NSW 2052, Australia.

‡

Australian Nuclear Science and Technology Organisation (ANSTO), Kirrawee DC, NSW 2232, Australia

* Email: [email protected]. Phone: +61 2 9385 4763. KEYWORDS. Gold nanoparticle, particle aggregation, computational fluid dynamics, discrete phase modelling, discrete dipole approximation.

ABSTRACT. Gold nanoparticle aggregation has a strong influence on the plasmonic resonance and hence the effectiveness in various photothermal applications. In relation to this, a comprehensive numerical model is developed to simulate and characterise the gold nanoparticle aggregation process at various particle volume fractions and base fluid pH levels. Computational fluid dynamics techniques are utilised to model the base fluid, while discrete phase modelling is adopted in determining the nanoparticle trajectories. Two-way coupling is implemented to handle the particle-fluid interactions. Discrete dipole approximation approach is utilised to further examine the absorption and scattering efficiency of various gold nanoparticle aggregate structures. At 1

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 40

lower particle volume fraction, short chain-like structures are formed in the particle aggregation process, with a more complex interconnected 'particle network’ structure formed at higher particle volume fractions. With the three base fluid pH levels investigated, gold nanoparticle aggregates are more compact with larger fractal dimensions and higher mean coordination numbers at pH = 3.5, whereas a more ‘loose’ structure formed at pH = 6.7 and 9.4 due to larger electrostatic repulsive forces as a result of changes in the zeta potential and Debye length of the gold nanoparticles. Among the typical gold nanoparticle aggregate structures characterised in this paper, the chain-like tetramer demonstrates the highest absorption efficiency of 1.83 at 700 nm wavelength – comparable to the plasmonic resonance of a nanorod – which lies in the optical window of biological tissue. Predictions of gold nanoparticle optical properties are found to be in good agreement with published experimental data.

Introduction The addition of solid particles to enhance the base fluid properties is a subject of significant interest for more than a century, ever since the field of colloidal science was introduced and pioneered by Thomas Graham back in the 1860s1. This practice has been reported to be effective in changing the mechanical properties such as viscosity2-3 and density4 of the base fluids, as well as thermal5 and optical6 properties. Recent advancements in manufacturing technologies have enabled the large-scale fabrication of micrometre and nanometre-sized particles, which allowed a large range of applications of such systems, such as magnetorheological (MR) fluid in MR fluid dampers7-8, nanofluid heat sink in cooling9, magnetic particle drag targeting10, and gold nanoparticle (GNP) hyperthermia in cancer treatment11. During these applications, particle aggregation is often observed, which may lead to favourable or undesirable outcomes, depending on the applications. In order to optimise the performance of particle 2


Page 3 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

suspensions, it is of fundamental importance to understand the influencing factors of the particle aggregation characteristics, including the size, shape and material of the suspended particles, the solid volume fraction of the particles, the pH of the base fluid, the flow behaviour of the base fluid, and the application of external magnetic/electric fields12. Nam et al.13 investigated the pH-induced aggregation of GNPs and their effect on photothermal cancer therapies. Albanese and Chan14 examined the effect of GNP aggregation on cell uptake and toxicity. Mocanu et al.15 examined the GNP aggregation characteristics under the influence of cysteine, and its effect on the absorbance and over the optical spectrum. At low pH (acidic conditions), negatively charged gold nanoparticles attract and absorb positively charged H+ ions to form a compact layer, known as the Stern layer, which partially screens the negative charges on the nanoparticle surface. A second, much looser layer of ions is also attracted near the surface of the nanoparticle via the Coulomb force, known as the diffuse layer. The zeta potential, defined as the electric potential in this interfacial double-layer, is a key indicator for the stability of GNPs against aggregation. Since the zeta potential is sensitive to the ionic strength, and hence the pH, of the base fluid, modifying the pH of the base fluid has a strong impact on the magnitude of the electrostatic repulsive forces between GNPs according to the DLVO theory, resulting in a drastic change in the aggregation behaviour of GNPs. Although the influencing factors and effects of the particle aggregation process have been experimentally examined on a macro level, a comprehensive computational model is yet to be developed for better characterisations at the micro-level. It is often very challenging for experimentalists to conduct measurements on or even visualise the particle suspensions due to limitations in instrumentations. Whereas for numerical modelling, quantitative data is readily available for every particle at any given time through the duration of the simulation. The traditional numerical modelling approach for particle-fluid mixtures based on the Eulerian-Eulerian framework solves the continuum equations for both fluid and particle phases; they are usually treated as a 3


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 40

single-phase fluid with adjusted properties. Such an approach assumes the mixture to be homogenous, and hence completely ignores the aggregation behaviour of GNPs. Discrete phase modelling approach based on the Eulerian-Lagrangian framework however solves the trajectory of each GNP through time. Various attempts have been made to study the aggregation behaviour of nanoparticles under this framework16-18. In this paper, a comprehensive computational model, which would be able to resolve all the essential physics of particle and fluid movement including Brownian motion, long-range particle-particle interactions, and particle-fluid interactions, is developed to purposefully understand the aggregation behaviour of GNPs. Owing to the increasing computational power and potential utilisation of irradiated GNPs for effective treatment of melanoma, the aggregation behaviour of GNPs and the plasmonic resonance resulting from the formation of selected structures of GNP aggregates are investigated. Conservation equations based on the continuum assumption are utilised to model the continuous phase (base fluid of the solution) which will be coupled with the discrete phase modelling to predict the translational and rotational movements of GNPs. Both long-range and contact particle-particle interactions are considered. The effects of different particle volume fractions and base fluid pH levels are examined. Discrete Dipole Approximation approach is adopted to examine the absorption and scattering efficiency of various typical GNP aggregation structures.

Mathematical Modelling The model for determining the plasmonic resonance of GNP aggregates in a surrounding medium with varied pH conditions comprises three components. The first component is based on the utilisation of the soft-sphere model to compute the particle trajectories and rotation rates based on the forces and torques acting on the particles. The second component describes the flow field within the physical domain which is achieved via the Navier-Stokes transport equations. A two-way coupling is adopted to couple the particle and fluid phases to predict the structure of GNP aggregates in different pH solutions. Once the formation of aggregates is 4


Page 5 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

established, the optical properties are evaluated using the discrete dipole approximation (DDA) method based on Mie theory, which constitutes the third component of the model. In the current study, the DDSCAT software, which has been developed and validated by Draine and Flatau19, is applied to calculate the absorption and scattering efficiency of GNP aggregates. Discrete Phase Modelling For the soft-sphere model, the motion of each discrete GNP is governed by Newton’s second law as described by:

m p ,i

d 2 x p ,i dt 2

= fvdw,ij + fe,ij + fc,ij + f d ,i + f B ,i + flift ,i + f g ,i (1)

I p ,i

dΩ p ,i dt

= M f .i + M c ,ij

(2)

In eq (1), m p ,i and x p,i are the mass and position of the ith discrete GNP. In eq (2), 𝐼𝑝,𝑖 and Ω p,i are the moment of inertia and angular velocity of the ith discrete GNP. For a homogeneous solid sphere, the moment of inertia I p ,i is given by I p,i = ( 2 / 5) mp,i rp2,i where rp,i is the radius of the discrete GNP and the mass mi is evaluated as mp,i = ( 4 / 3) rp3,i  p,i where  p ,i is the density of the discrete GNP. The right-hand sides of eqs (1) and (2) represent the consideration of forces (van der Waals force f vdw , electrostatic force f e , collision force f c , drag force f d , Brownian force f B , Saffman lift force flift and buoyancy and gravity force f g ) and torques

(fluid-induced torque M f and collision-induced torque M c ) acting on the discrete GNP in this study. Note that in eqs (1) and (2), the subscript ij denotes the sum of forces or torques between the ith discrete GNP and all other surrounding GNPs. Continuous Phase Modelling Following the approach suggested by Marshall20, the Navier-Stokes transport equations governing the base fluid of the solution can be written in the form as

5


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

  f f +  f f vf = 0 t

(

  f  f vf +   t  = − f p +     f 

(

)

)

(

f

(

)

 f vf  vf

(

Page 6 of 40

(3)

)

   f v f + v f 

(

)

T

)

2  −  f   vf    + S 3 

(4)

where  f = 1 −  s is the void fraction,  f is the fluid density, v f is the fluid velocity, p is the pressure,  f is the dynamic viscosity of the fluid, and S is a source term which accounts for the two-way coupling due to the hydrodynamic interactions with the discrete phase. Particle-Fluid Interactions GNPs will interact with the surrounding fluid through various forces and torques. This section describes the forces and torques that are relevant to the model being considered in this study. For sub-micron particles, the Stokes drag force f d ,i acting on the GNP caused by the movement of the base fluid is determined using the Cunningham drag law21 which can be expressed in the form as

f d ,i =

6 f rp ,i Cc

(v

f ,i

− v p ,i

)

(5)

where v f ,i and v p ,i are the velocities of the base fluid and the discrete GNP. The Cunningham correction C c which corrects the Stokes drag force is given by Cc = 1 +

 rp ,i

(1.257 + 0.4

(

−1.1 rp ,i / 

)

)

(6)

where  represents the molecular mean free path. The Brownian force f B,i for sub-micron particles can be modelled as a Gaussian white noise process according to Ounis et al. 21,

6


Page 7 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

f B ,i =  i

So t p

(7)

in which  i are Gaussian random numbers with zero-mean and unit variance, t p is the time-step for the particle simulation and the term So is given as

So =

27  f kBT f

(8)

4 2 rp5,i  p2.i Cc

where k B is the Boltzmann constant which can be taken to be 1 10−23 J/K and T f is the temperature of the solution. Since the Cunningham correction C c appears in eq (8), the magnitude of the Brownian force greatly depends on the molecular mean free path. It is evident that hydrodynamic lubrication forces will suppress the Brownian force at small particle separation distances. In addition, when an aggregate is formed, the entire aggregate behaves like a single particle with a significantly larger diameter, which will also dramatically decrease the magnitude of the Brownian force according to eq (7) and (8). Thus, the Brownian force will be limited for the aggregate within a 1 nm distance from other GNPs. According to Marshall20, the Saffman lift force can be expressed in terms of the lift coefficient K L as

flift ,i = K L

f ( v f ,i − v p ,i )   f , i m p.i  p ,i Re,i

(9)

where K L =2.18, Ω f ,i =  v f ,i is the angular velocity of the fluid at the position of the particle centre and Re,i = 2 f f ,i rp2,i  . This force becomes important especially for low Reynolds number where a small

rotating particle experiences a lift force both due to the pressure difference between the top and bottom of the particle when it rotates with the fluid and local gradients of transitional fluid velocities. The buoyancy and gravity forces f g ,i for sub-micron particles can be expressed by 7


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

 p ,i −  f m p ,i g f

f g ,i =

Page 8 of 40

(10)

The torque on a GNP caused by the local rotation of the solution can be computed as

M f .i

8 f rp3,i  1  =  Ω f ,i − Ω p ,i  Cc  2 

(11)

Special attention is required to determine the particle volume fraction  s , especially where the aggregates could grow to the point where they could occupy the entire grid cell. To circumvent the possible occurrence of numerical instability, the approach by Xiao and Sun22 is adopted herein. The particle volume fraction at each grid cell can be determined as s =

1 cn 3 i ( 4 / 3)  rp,i Vc i =1

(12)

with

i =

K ( x − xc / b )



Nc m =1K

(13)

( x − xc / b )

where Vc is the volume of the grid cell, cn is the total number of GNPs in the grid cell, i is the weight of the contribution from each ith discrete GNP and N c is the total number of cells in the computational domain. The following kernel function by means of a Gaussian function is given by  2  2 K ( ) =   

(

)

−3 2

  2  if   1 exp  −  2  if   1 0

(14)

8


Page 9 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

where  = x − xc / b . In this study, the bandwidth b is chosen to be b = 2rp,i / 2 according to the proposal of Khare et al.23. For the source term S in eq (4), the momentum transfer from the particle phase to the fluid phase can be evaluated using cn

(

)

S =  f  f d ,i + flift ,i + f g ,i / Vc i =1

(15)

Equation (15) represents the sum of forces derived from the drag and gravity exerted on the solution from all the total number GNPs represented by cn that reside within the grid cell. Particle-Particle Interactions The fundamentals of particle-particle interactions can be represented using a combination of particle and continuum mechanics. Using the model with stiff particles and soft contacts as proposed by Marshall20, the influence of elastic repulsion for GNPs in a representative particle contact can be demonstrated via the contact force equilibrium. Also, a sphere-sphere adhesion model without any contact deformation for GNPs can be combined with a plate-plate adhesion model for nano-contact flattening. With the increasing flattening, the contact forces and torques acting on GNPs can be described via the Hertzian theory of elastic force-displacement and viscous damping along the line normal to the particle centres and the resistance due to sliding, twisting and rolling of particles over one another. For long-range particle-particle interactions, namely the van der Waals attraction and electrostatic repulsion forces, fvdw,ij and f e ,ij , are derived from the corresponding interaction potential energy ( U vdw,ij and U e,ij ). As predicted by the DLVO theory, f =−

dU dRij

(16)

9


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

where Rij is the centre distance between the ith and jth GNPs. The van der Waals attractive and electrostatic repulsive potential energies, U vdw,ij and U e,ij , are given in eqs (17) and (18), respectively 24.

U vdw,ij

2  4rp ,i 2   2rp ,i 2 AH  2rp ,i  =− + + ln 1 −   (17) 6  Rij 2 − 4rp ,i 2 Rij 2 Rij 2      

( (

))

U e,ij = 2 rp ,i 0 2 ln 1 + exp − Rij − 2rp ,i  (18)  

where AH is the Hamaker constant,  =  0 rel is the permittivity of the medium,  0 is the surface potential of the GNPs, and  is the reciprocal of the Debye length. DLVO theory predicts the combined effects of the van der Waals attractive and electrostatic repulsive forces. However, this model is only used when GNPs are not in contact with each other. When GNPs collide, the particle contact force f c ,ij becomes effective. The approach to compute the particle contact force is adopted from Marshall20. During contact, the van der Waals attractive force is still active, and hence the contact force needs to be modified accordingly. The particle contact force f c ,ij can be calculated as:

f c ,ij

 a 3  a 3/2  = 12 reff   −     a0   a0    

(19)

where  is adhesion surface energy of the GNPs, reff is the effective particle radius between binary particles, which can be calculated using eq (20), a is the radius of the flattened contact region between the GNPs, and a0 is the value of a in the equilibrium state. The value of a0 can be calculated using eq (21), in which Eeff is

the effective Young’s modulus of the colliding particles, and  p is the Poisson’s ratio of the particles. 1 1 1 = + reff rp ,i rp , j

(20) 10


Page 11 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

1/3

 9 reff 2  a0 =    Eeff   

(21)

1 −  2p,i 1 −  2p , j 1 = + Eeff E p ,i E p, j

(22)

The radius of the flattened contact region a in eq (19) depends on the magnitude of the particle overlap  :   a  4 a   = 61/3  2   −   c   a0  3  a0   2

1/2

   

(23)

Particle-Wall Interactions GNPs in close proximity to wall boundaries will interact with the wall through various long-range and contact mechanisms similar to particle-particle interactions as specified in the previous section. Although the primary focus of this study is to examine the particle-particle aggregation behaviour in the bulk, it is essential to build a robust particle-wall interaction model to ensure numerical stability at nearwall regions. For simplicity, only repulsive interactions between GNPs and walls are considered. Colloidal particles suspended in the base fluid in the near-wall region experience a strong short-range repulsive force, known as the hydration or structural force.25 The mathematical expression to characterise the potential energy U h,iw resulted from this repulsive interaction is adopted from Kroupa et. al.26 U h,iw =  rp,i F0 0 exp(− Hiw /  0 )

(24)

where F0 is the hydration force constant,  0 is the characteristic decay length, and H iw is the surface separation distance between the ith particle and the wall. The values for F0 and  0 are experimentally determined by Thompson and Collins27. Plasmonic Resonance of Particle Aggregates DDA, introduced by DeVoe28, is considered as one of the most flexible and reliable methods in calculating the optical properties of particles with arbitrary geometries. The 11


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 40

target is replaced by an array of point dipoles (polarisable points) and the electromagnetic scattering problem for an incident periodic wave interacting with this array of point dipoles is solved to calculate the absorption and scattering efficiency of the aggregates of GNPs. Through DDA, the size of the target is characterised by the effective radius Reff of any arbitrary structure:

Reff =

3

3V 4

(25)

where V is the equivalent volume. The shape of gold nanorod (GNR) is characterised by the aspect ratio (AR = axis length/diameter). Since the solution is water based, the refractive index of water, 1.33 + 0i , was applied at all wavelengths. McPeak et al.29 measured the values of the complex dielectric function of gold at various wavelengths using a thin film of gold. However, when the target dimension was smaller than the mean free path length of electrons in the bulk material, the limitations to electron collisions within the particle would change the dampening constant of conduction electrons. This would result in size-dependent dielectric constants which are different from the dielectric constants of the bulk material. Kreibig and Vollmer30 developed a model to modify the real and imaginary components of the dielectric constant. This size effect need to be considered, and the complex dielectric function measured by McPeak et al.29 need be modified in the current case. Based on the model developed by Kreibig and Vollmer30, the real ( r ) and imaginary (  i ) components of the bulk dielectric constant  bulk can be calculated using the refractive index of the bulk material:  r = n2 − k 2

(26)

 i = 2nk

(27)

12


Page 13 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

where n is the index of refraction and k is the index of absorption. The modified dielectric constant  m can thus be expressed as  m =  bulk +  p

1 1 − p 2  + i bulk  + i eff 2

(28)

in which  is the frequency of the incident electromagnetic wave. According to Ordal et al.31, the value of bulk plasmonic frequency  p can been taken to be 9.03 eV and the relaxation frequency of the bulk material  bulk as 0.072 eV. The effective relaxation frequency  eff can be written in the following form:

 eff =  bulk +

As u f leff

(29)

where the value of the Fermi velocity u f is 1.4  106 m/s and the value of the surface scattering parameter As is 0.332. The effective mean free path of the photons leff is taken to be equal to the length of GNP aggregates.

Numerical Modelling Computational Domain Although the computer code has been principally developed to analyse any generic three-dimensional problems, two-dimensional simulations have been conducted in this study to better understand the aggregation behaviour of GNPs as represented by the application of a thin-layer of GNP solution on the surface of human skin. This restricts the motion of GNPs and fluid flow to be predicted only on the x-y plane. The resulting particle and fluid velocities, as well as the corresponding particle-particle and particlefluid interaction forces in the z-direction (perpendicular to human skin surface) are set to zero. Fluid Domain The computational domain is a small square window of the applied solution. The dimension of the fluid domain is 1.5μm by 1.5μm by 33.51nm. The computational domain is discretised with 40 divisions 13


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 40

of uniform spacings along the x-direction and y-direction with a division of 1 uniform spacing along the zdirection (Refer to Figure S1 in the supporting information). The mesh quality is demonstrated to be adequate by a mesh convergence study. Water is used as the based fluid with density  f = 998.2kg/m3, dynamic viscosity  f = 1.003cP and relative permittivity  rel of water at 25°C of 80. Discrete Phase The GNPs considered in this study are spheres with a diameter of 20nm. The properties are based on GNPs synthesised and surface treated according to an established method33. While the size of the fluid domain remains unchanged for all the simulation cases conducted in this study, three different particle volume fractions, namely  s = 5%, 10% and 20%, have been considered which corresponds to a total of 900, 1,800 and 3,600 particles respectively. Key modelling parameters for the continuous and discrete phases are summarised in Table 1. Table 1. Key Modelling Parameters Symbol Quantity AH Hamaker constant dp Particle diameter Young’s modulus34 E F0 Hydration force constant27  Particle surface energy 0 Characteristic decay length27 0 Vacuum permittivity  rel Relative permittivity of water at 25°C pH 3.5 35 −1 Debye Length pH 6.7  pH 9.4 f Dynamic viscosity of base fluid f Density of base fluid

Value 2.5×10-19 20nm 1.16GPa 2×108N/m2 1.44mJ/m2 1.0nm 8.85×10-12 F/m 80 15.1nm 60.7nm 18.8nm 1.003cP

p

19,320kg/m3 -19.8mV -169.7mV -61.8mV

0

Density of particles Particle surface potential

35

pH 3.5 pH 6.7 pH 9.4

998.2kg/m3

14


Page 15 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Boundary and Initial Conditions Since the computational domain represents a small square window of the irradiated GNPs in the real-life treatment of melanoma, symmetry boundary conditions are applied on all boundaries of the fluid domain. The initial velocity of the fluid was set to zero at the beginning of the simulation. For the initial placement of GNPs in the small square window, the locations for each GNP are randomly determined within the fluid domain via an in-house computer code. Each GNP is checked against all neighbouring GNPs to avoid possible particle overlaps before the numerical calculation proceeds. Please refer to Figure S2 in the supporting information for the initial particle positions at the start of the computation. Numerical Schemes The Finite Volume Method (FVM) is used to discretise the Navier-Stokes equations, and the pressure-velocity coupling is computed using the PISO scheme36. For spatial discretisation, 2nd order upwind scheme is employed to determine the face values based on cell-centre values, and the least squares cell-based approach is used to discretise the diffusion and advection terms. First order implicit integration is used for the temporal discretisation. A two-way coupling approach is adopted between the discrete phase trajectories and the fluid continuum, as shown in Figure S3 in the supporting information. At the start of the computation for each time step, the flow field of the continuum (base fluid) is solved under the computational fluid dynamics principles using the Algebraic Multi-Grid (AMG) numerical solver. The mass and momentum conservation equations are solved using an iterative method, which provides essential input variables for the computation of the particle-fluid interaction forces. Together with long-range and contact particle-particle interaction forces calculated using the DLVO theory and the soft-sphere approach under discrete phase modelling, the total force acting on each nanoparticle, and hence the acceleration, velocity, and updated position of each nanoparticle are determined. To close the two-way coupling loop, the particle trajectory information is fed back to the computation of the

15


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 40

fluid domain through the source term of the Navier–Stokes momentum equations for each mesh cell. The fluid and discrete equations are solved alternately until the numerical solutions for both phases converge. While the time step adopted for the fluid continuum t f is 1×10-9s, the time step used for the marching of each discrete GNP is taken to be considerably smaller at 1×10-11s. As recommended by Marshall20, the particle time step t p is characterised by the smaller of particle hydrodynamic time and particle convective time scale  tc , as shown in eq (30), where U p ,i is the characteristic macroscopic velocity of the GNPs.   2p ,i tc = 2rp ,i  2  E p ,iU p ,i 

   

(30)

Results and Discussion There are multiple influencing factors that can significantly affect the dynamics of GNP aggregation, including particle properties (e.g. material, size, coating), fluid properties (e.g. viscosity, pH) and external environmental factors (e.g. temperature, magnetic/electric field). The general GNP aggregation characteristics are first presented. Then the effect on GNP aggregation due to varying particle volume fraction  s and pH of the base fluid is further analysed. Lastly, predictive results on the plasmonic resonance resulting from the formation of selected GNP aggregation structures and the comparison with experimental data are discussed. Particle Aggregation Characteristics The results presented in this section is based on 10% particle volume fraction with base fluid pH = 3.5. A set of snapshots of the GNP structure at different stages of the particle aggregation process is included in Figure 1. Initially, the particles are randomly distributed within the computational domain. Initial particle overlaps are avoided during particle injection. Shortly after the injection, random motion of GNPs is observed due to Brownian motion, which leads to the collision of particles. Due to van der Waals attraction forces, a small number of particles remain in contact and form aggregates after the 16


Page 17 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

collision. These small particle cluster then join each other to form larger aggregates. It can be seen that aggregation happened rapidly between t=1×10-6s and t=1×10-5s, after which, an equilibrium state is reached with minimal changes afterwards.

(a)

(b)

(c)

(d) (e) (f) Figure 1. Particle structure with particle volume fraction of 0.10 and base fluid pH of 3.5 at (a) t=1×105

ms, (b) t=1×10-4ms, (c) t=0.001ms, (d) t=0.005s, (e) t=0.01s, (f) t=0.2s.

Figure 2 (a) shows the number of ‘singular’ GNPs within the computational domain as a function of time. It is shown that the GNP aggregation process is very rapid, and within 0.1ms, most particles are contained within an aggregate, with the number of single unaggregated GNPs becoming negligible. Figure 2 (b) demonstrates the change in the average number of particles in an aggregate with respect to time. As small aggregates join each other to form larger aggregates, the average number of particles in each particle cluster increases 17


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 40

throughout the simulation. These particle aggregation characteristics are common among all simulations cases in this study, with a range of different particle volume fractions and base fluid pH levels.

(a)

(b)

Figure 2. (a) Number of single particles not in an aggregate v.s. time, and (b) average number of particles in each particle cluster v.s. time, with particle volume fraction of 0.10 and base fluid pH of 3.5. To understand the particle aggregation process quantitatively, a fractal analysis is conducted. The fractal dimension D f of a structure measures the self-similarity and irregularity of its geometry. A higher D f indicates more complex and interconnected aggregate structure, while lower D f indicates a simpler chain-like aggregate structure. By definition, Al

Df

(31)

And hence

log l A =

log ( A ) log ( l )

 Df

(32)

18


Page 19 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

where A is the two-dimensional projected area of the particle clusters and the longest dimension of a particle cluster l is defined as the largest distance between any two points on the outline of the particle cluster as shown in Figure 3, and a typical plot utilising eq (32) is shown in Figure 4. A and l are determined using an in-house post-processing code37-39.

Figure 3. Image post-processing output showing the longest dimension of each particle cluster with particle volume fraction of 0.10 and base fluid pH of 3.5 at t=0.2ms.

19


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 40

Figure 4. Logarithmic plot between the two-dimensional projected area of the particle clusters and its longest dimensions with particle volume fraction of 0.10 and base fluid pH of 3.5 at t=0.2ms, where the gradient of the line of best fit indicates the fractal dimension. The fractal dimensions at different stages of the aggregation process are then determined by the gradient of the line of best fit using linear regression, and are summarised in Table 2. Table 2. Longest dimension, fractal dimension and regression R2 as a function of time with particle volume fraction of 0.10 and base fluid pH of 3.5. s

0.10

l (μm)

Df

1  10−5

0.095

1.041

R2 0.984

1  10−4

0.120

1.089

0.978

−4

0.126 0.144 0.302 0.304 0.292 0.307 0.307 0.306

1.150 1.168 1.211 1.236 1.281 1.284 1.289 1.283

0.970 0.964 0.947 0.945 0.942 0.945 0.947 0.943

t (ms)

5  10 0.001 0.005 0.01 0.05 0.1 0.15 0.2

pH

3.5

As shown in Table 2, initially, the longest dimension l for the largest particle cluster in the computational domain rapidly increases with respect to time. After 0.005ms, the longest dimension stabilised at around 0.3µm. 20


Page 21 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

The magnitude of the fluctuations in the longest dimension is less than the diameter of the GNPs of 20nm, and therefore, these fluctuations are not caused by the attachment of new primary particles or the breakage of existing particle clusters. They are caused by the slight change of particle orientations within the largest cluster. This can be confirmed by the observations from Figure 1. In order to determine the fractal dimensions at each stage of particle aggregation process, linear regression is used to find the line of best fit within the log( A) v.s. log(l ) plot. The coefficient of determination R 2 is greater than 0.94 for all cases, indicating high reliability of the predicted results. The fractal dimension increases from 1.04 at the start of the simulation until it converges to a value of about 1.28 after 0.05ms, indicating a more complex and more interconnected GNP aggregate is formed.

Effect of Particle Volume Fraction on GNP Aggregation In this study, three different particle volume fractions are considered, namely  s = 0.05, 0.10 and 0.20. A set of snapshots of the particle structure at t = 0.2ms with base fluid pH = 3.5 are included in Figure 5. It can be seen that the structure of GNP aggregates changes dramatically with respect to different particle volume fractions. As shown in Figure 5 (a), with a dilute suspension of  s = 0.05, the GNP clusters are dominated by dimers, trimers and other short chain-like structures. As  s increases to 0.10, much longer and more complex chain structures are formed, as the shorter chains formed at the initial stage are in much closer proximity to each other, resulting in a larger van der Waals attraction between them. As  s further increases to 0.20, an interconnected ‘particle network’ is formed. As shown in Figure 6, only ten particle clusters exist in the computational domain at t=0.2ms with particle volume fraction of 0.20 and base fluid pH of 3.5. One primary particle cluster spanned across the entire computational domain with nine other much smaller particle clusters filling in the space between the branches of the primary particle network.

21


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

Page 22 of 40

(c)

Figure 5. Particle structure with particle volume fraction of (a) 0.05, (b) 0.10, and (c) 0.20 with base fluid pH of 3.5 at t=0.2ms.

Figure 6. Image post-processing output coloured by cluster number showing the longest dimension of each particle cluster with particle volume fraction of 0.20 and base fluid pH of 3.5 at t=0.2ms. In order to quantitatively understand the effect of particle volume fraction on particle aggregation, the mean coordination number Z is introduced. The mean coordination number is an important physical quantity in measuring the mechanical and thermal properties of self-assembly structures, and is calculated as: 22


Page 23 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Np

N c ,i

i =1

Np

Z =

(33)

where N p is the total number of particles in the computational domain, and N c ,i is the number of particles in contact with the ith particle. Thus, dimers have a coordination number of 1. With a long straight chain of particles, Z = 2 , and D f = 1 With more complex particle cluster geometries, the Z generally increases. The theoretical maximum mean coordination number for particles with uniform diameters under a two-dimensional setup is 6, and 12 under a three-dimensional setup. Figure 7 shows the mean coordination number Z as a function of time for different particle volume fractions. It is observed that the mean coordination number for the 0.20 particle volume fraction case increases at a much faster rate at the initial stage of the simulation, indicating rapid particle aggregation due to dominating van der Waals attraction at small particle separation distances. When equilibrium is achieved, the  s = 0.20 case has a much larger mean coordination number of 3.48 compared to 2.41 for the  s = 0.10 case and 1.59 for the  s = 0.05 case.

Figure 7. Mean coordination number as a function of time for particle volume fraction of 0.05, 0.10 and 0.20 with base fluid pH of 3.5.

23


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 40

Figure 8 shows the longest dimension of the particle clusters as a function of time for different particle volume fractions with a constant base fluid pH of 3.5. It is obvious that the longest dimension for the  s = 0.20 case increases the fastest out of the three test cases due to rapid particle aggregations. Note that the longest dimension for the  s = 0.20 case is capped at 2.07µm, which is approximately the diagonal distance of the square computational domain. Arguably, the longest dimension has lost its physical importance once this value is reached. As shown in Figure 6, the primary particle network spanned across the entire computational domain. Since the computational domain only represents a small window of the real-life system, in reality, the particle network would not be capped, and may grow beyond the size of the computational domain.

Figure 8. Longest dimension as a function of time for particle volume fraction of 0.05, 0.10 and 0.20 with base fluid pH of 3.5. The longest dimension, fractal dimension and coefficient of determination R 2 for different particle volume fractions at t=0.2ms with base fluid pH of 3.5 (refer to Figure 5) are summarised in Table 3. The longest dimension ranges from 0.113µm to 2.062 µm, and the fractal dimension increases from 1.133 to 1.755 as the as the particle volume fraction increases from 0.05 to 0.20, indicating a more complex and compact structure at higher particle volume fractions. This is consistent with the qualitative observations in Figure 5.

24


Page 25 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Table 3. Longest dimension, fractal dimension and regression R2 for different particle volume fractions at t=0.2ms with base fluid pH of 3.5. s 0.05 0.10 0.20

t (ms) 0.2

pH

l (μm)

Df

3.5

0.113 0.306 2.062

1.133 1.283 1.755

R2 0.885 0.943 0.995

In summary, the particle volume fraction largely affects the GNP aggregation behaviour. At high particle volume fractions, GNPs formed a large interconnected particle ‘network’, where van der Waals attraction between GNPs is dominant at smaller particle separation distances, as predicted by DLVO theory. At lower particle volume fractions, chain-like particle aggregates are formed at the initial stage. Further particle aggregation is prevented by dominating electrostatic repulsion forces at larger particle separation distances. Effect of pH on GNP Aggregation It is experimentally proven that the pH of the base fluid has a significant impact on the aggregation of GNPs13. In this study, the effect of base fluid pH on GNP aggregation is examined numerically, and the results are presented and discussed in this section. A set of snapshots of the particle structure at t = 0.2ms with  s = 0.10 at three different base fluid pH levels is included in Figure 9. It can be seen that the particle aggregates formed for the pH = 3.5 case is much thicker and longer compared to the pH = 6.7 and 9.4 cases, where particle aggregates are dominated by short singlestrand particle chains. The change in pH is closely related to the ionic concentration in the base fluid, which in turn affects the surface properties of the GNPs such as zeta/surface potential and Debye length35.

25


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

Page 26 of 40

(c)

Figure 9. Particle structure with particle volume fraction of 0.10 and base fluid (a) pH of 3.5, (b) pH of 6.7 and (c) pH of 9.4, at t=0.2ms. Figure 10 shows the mean coordination number Z as a function of time for different continuous phase pH levels. In consistence with the qualitative observations in Figure 9, GNPs demonstrate a much higher mean coordination number at 3.5 base fluid pH, indicating larger and more interconnected particle aggregates.

Figure 10. Mean coordination number as a function of time for particle volume fraction of 0.10 with base fluid pH of 3.5, 6.7 and 9.4. The longest dimension of the particle clusters as a function of time with  s = 0.10 at different continuous phase pH is shown in Figure 11. It can be seen that in the initial stage of the simulation, the longest dimension 26


Page 27 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

for particle clusters for all three cases are very close, indicating that the variation in base fluid pH level is not the controlling factor at this stage, when particle separation distances are relatively small. After preliminary particle aggregates are formed within the solution, the effect of base fluid pH becomes dominant, and is the main contributing factor in later stages of the particle aggregation process, leading to larger particle clusters for the pH = 3.5 case, compared to the pH = 6.7 and pH = 9.4 cases.

Figure 11. Longest dimension as a function of time for particle volume fraction of 0.10 with base fluid pH of 3.5, 6.7 and 9.4. The longest dimension l , fractal dimension D f and coefficient of determination R 2 at t=0.2ms with  s = 0.10 at different continuous phase pH levels (Refer to Figure 9) are summarised in Table 4.

Table 4. Longest dimension, fractal dimension and regression R2 for different particle volume fractions at t=0.2ms with base fluid pH of 3.5. s 0.10

t (ms)

pH

l (μm)

Df

0.2

3.5 6.7 9.4

0.306 0.193 0.166

1.283 1.173 1.114

R2 0.943 0.945 0.909

27


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 40

The longest dimension at pH = 3.5 reached 0.306μm, and is considerably larger than the other two cases. The fractal dimension for the pH = 3.5 also has a larger value, indicating a more compact structure for the particle aggregates formed. In summary, the base fluid pH has a strong influence on the aggregation behaviour of GNPs. With the three different pH levels investigated, GNPs formed more compact and larger aggregates under acidic conditions (pH = 3.5), while shorter chain-like structures are formed under neutral and alkaline conditions (pH = 6.7 and pH = 9.4). Effect of GNP Aggregation on Plasmonic Resonance Owing to the limitation of measurement techniques and large quantity of GNP aggregations, identifying the structures of GNP aggregations can be rather difficult in experimental studies. In order to validate the predicted particle aggregation structures, the absorption efficiency of various GNP aggregation structures, as shown in Table 5, is calculated using DDA, and the results are compared with the experimental data obtained by Peng et al.40, as depicted in Figure 12. In addition to the transverse surface plasmonic resonance (TSPR) observed at 520nm wavelength, a longitudinal surface plasmonic resonance (LSPR) within the so-called optical window of biological tissues (700nm to 900nm)41 is also observed in both the experiment40 and the numerical prediction for the GNP solution under acidic condition. In the experiment40, the size distribution of GNPs is generally non-uniform and it is useful to characterise the mean diameter of GNPs in the solutions. A mean diameter of 20nm is utilised in the current study. As can be seen for the solution of pH = 12, the GNPs are rather independent; the plasmonic resonance can be compared to the prediction of a single particle. The normalised absorbance of single particle shows a good agreement with the experimental results between the optical wavelength of 450nm and 800nm. However, a minor discrepancy is observed between the experimental result and the numerical predicted value at wavelengths from 400nm to 450nm. The size of the GNPs used in the experiment is 28


Page 29 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

non-uniform with a distribution of 21.3 ± 4.4nm40, whereas the numerical simulation was conducted based on an idealised GNP with a mean diameter of 20nm. Table 5. Typical particle structures characterised in the plasmonic resonance study. Particle Cluster Identifier Number of Particles Particle Structure Monomer

1

Dimer

2

Trimer 3-0

3

Trimer 3-1

3

Tetramer 4-0

4

Tetramer 4-1

4

Tetramer 4-2

4

Octamer 8-0

8

Octamer 8-1

8

With the base fluid pH = 6.5, GNP aggregations are observed during the experiment40, and the experimentally determined absorbance is included in Figure 12. The dash line in Figure 12 is the weighted average of the absorbance of Monomer (29.7%), Dimer (43.9%), Trimer 3-0 (3.7%), Trimer 3-1 (13.6%), Tetramer 4-0 (0.5%), Tetramer 4-1 (3.6%) and Tetramer 4-2 (5.0%) cases, as predicted by the discrete phase simulations with base fluid pH = 6.7. The predicted plasmonic resonance can be seen to be in good agreement with the experimental results.

29


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 40

Figure 12. Normalised absorbance of GNP (diameter = 20 nm) monomer and aggregations compared with experimental data measured from nanoparticle (mean diameter = 20 nm) solutions of pH = 12 and pH = 6.5 Comparing the absorption efficiency spectrum of GNP monomer with GNP dimer and trimer, it can be seen that the chain-like structures have a TSPR at around 520nm as well as a LSPR at nearly 700nm. However, the monomer only has a TSPR which may not be as applicable as the LSPR in biomedical applications because it does not align with the optical window of biological tissues (700nm to 900nm). The wavelength of LSPR shows a red shift from 620nm to 680nm when the number of GNPs in the aggregate increases from 2 to 3, and the length of aggregation expands. However, the wavelength of LSPR depends not only on the number of GNPs, but also on the shape of the structures. Figure 13 (a)-(c) demonstrate the absorption efficiency of various typical GNP aggregate structures (refer to Table 5). Trimer 3-1 shows a TSPR at 520nm as well as a LSPR at 660nm, whereas Trimer 3-0 only has a TSPR at 540nm, although the number of GNPs is identical in these two cases. The cause of the difference is that the corresponding aspect ratio of Trimer 3-1 and Trimer 3-0 is around 3 and 1 respectively. The effect of particle numbers in an aggregate on its plasmon resonance is comparable to the effect of aspect ratio of nanorods42. In Figure 13 (b), the absorption efficiency of Tetramer 4-1 is compared with the spectrum of

30


Page 31 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

nanorod (AR = 3). The effective radius Reff of the nanorod is the same as Tetramer 4-1. Due to the similar length to width ratio, the spectra present comparable LSPRs between these two cases. A similar trend is also shown in Figure 13 (c). Octamer 8-1 shows a TSPR at 520nm and high absorption efficiency of nearly 1.8 at LSPR between 640nm and 660nm. However, the absorption efficiency of Octamer 8-0 is only 1.5 at its LSPR of 660nm. Among all the calculated cases, Tetramer 4-1 reaches the highest absorption efficiency of 1.9 at 680nm, a wavelength which is very close to the range of the optical window. Owing to the high absorption efficiency near 700nm, the tetramer conjugation in a chain-like structure would have a better capability in heat generation when irradiated in potential biomedical applications. In this section, the understanding and controlling of the structure of GNP aggregations provide useful information in the manipulation of the spectrum of absorption efficiency of GNP solutions. For the current investigation, it is demonstrated that the particle aggregation behaviour can be controlled by possibly altering the pH level of the base fluid, as discussed in Section 4.3. As a result, this analysis has predicted that GNP aggregations under acidic condition that can result in an additional LSPR, and thus a higher absorption efficiency. Tetramer 4-1 contributes the highest absorption efficiency in the range of the optical window of biological tissue. It can be inferred that the Tetramer 4-1 structure would have a better performance than the other spherical GNP aggregations in various photothermal applications, such as hyperthermia cancer therapies.

31


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 40

(a)

(b)

(c) Figure 13. DDA calculated absorption efficiency of various GNP aggregation structures (Refer to Table 5): (a) basic short chain-like structures, (b) different structures of trimers and (c) different structures of octamers.

32


Page 33 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Summary and Conclusions In the attempt to understand how aggregation behaviours of gold nanoparticles (GNPs) affect their performance in heat absorbance under plasmonic resonance in biomedical applications, a comprehensive computational model is developed and presented in this paper. Under the present study, the continuous phase (base fluid) is modelled using computational fluid dynamics principles, and the GNPs are modelled using the Discrete Phase Modelling approach with both long-range and contact interactions between particles taken into consideration. Two-way coupling is implemented between the base fluid and GNPs. Discrete Dipole Approximation approach is adopted to examine the absorption and scattering efficiency of GNP aggregates. Numerical simulations are conducted based on a range of different particle volume fractions and base fluid pH levels. It is observed that the particle aggregation process is very rapid, with an equilibrium state being achieved within 0.2ms. At lower particle volume fractions, short chain-like structures are formed at the end of the particle aggregation process, while a more complex interconnected particle ‘network’ structure is formed at higher GNP volume fractions. With the three base fluid pH levels investigated, GNPs formed more compact aggregates with larger fractal dimension and higher mean coordination number at pH = 3.5, whereas a more ‘loose’ structure being formed at pH = 6.7 and 9.4 due to larger electrostatic repulsive forces at a result of changes in GNP zeta potential and Debye length. Furthermore, the structures of GNP aggregates have a significant impact on the optical properties of the GNPs. Among the seven typical GNP aggregate structures investigated, the chain-like tetramer showed the highest absorption efficiency of 1.83 at 700 nm wavelength, within the range of the optical window of biological tissue. DDA predicted GNP optical properties are compared with experimental data40, and a good agreement is achieved. In conclusion, the present approach is capable of characterising the GNP aggregation behaviours and their effects on the absorption and scattering efficiency of GNP clusters. Although the focus of this paper is the 33


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 40

surface application of GNP solutions represented by a 2D computational model, the in-house code is developed to be ‘3D proof’. The present model helps medical practitioners to better understand the effect of particle aggregation on the performance in treatment procedures using GNP solutions, and how to control the particle aggregation by adjusting particle volume fractions and base fluid pH levels. Associated Content Supporting information Mesh grid used for the fluid domain computations; initial particle injection positions; scheme of the twoway coupling approach; Author Information Corresponding author *Email: [email protected] Notes There is no competing financial interest to declare. Acknowledgements This research is supported by the Australian Research Council (ARC Project ID DP150101065 and DP160100021). The support of NVIDIA Corporation with the donation of the Titan Xp and Quadro P6000 GPUs used for this research is also gratefully acknowledged. References 1. Graham, T., X. Liquid diffusion applied to analysis. Philosophical Transactions of the Royal Society of London 1861, 151, 183-224.

34


Page 35 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

2. Masoumi, N.; Sohrabi, N.; Behzadmehr, A., A new model for calculating the effective viscosity of nanofluids. Journal of Physics D: Applied Physics 2009, 42 (5), 055501. 3. Duan, F.; Kwek, D.; Crivoi, A., Viscosity affected by nanoparticle aggregation in Al2O3-water nanofluids. Nanoscale Research Letters 2011, 6 (1), 248. 4. Pastoriza-Gallego, M. J.; Casanova, C.; Páramo, R.; Barbés, B.; Legido, J. L.; Piñeiro, M. M., A study on stability and thermophysical properties (density and viscosity) of Al2O3 in water nanofluid. Journal of Applied Physics 2009, 106 (6), 064301. 5. Xuan, Y.; Li, Q.; Hu, W., Aggregation structure and thermal conductivity of nanofluids. AIChE Journal 2003, 49 (4), 1038-1043. 6. Norman, T. J.; Grant, C. D.; Magana, D.; Zhang, J. Z.; Liu, J.; Cao, D.; Bridges, F.; Van Buuren, A., Near Infrared Optical Absorption of Gold Nanoparticle Aggregates. The Journal of Physical Chemistry B 2002, 106 (28), 7005-7012. 7. Li, D. D.; Yeoh, G. H.; Timchenko, V.; Lam, H. F., Numerical Modeling of Magnetic Nanoparticle and Carrier Fluid Interactions Under Static and Double-Shear Flows. IEEE Transactions on Nanotechnology 2017, 16 (5), 798-805. 8. Li, D. D.; Yeoh, G. H.; Timchenko, V.; Lam, H. F., Numerical modelling of magnetic nanoparticle and carrier fluid interactions. In 2016 IEEE Nanotechnology Materials and Devices Conference (NMDC), Toulouse, France, 2016; pp 1-3. 9. Lee, A.; Li, D. D.; Lau, G. E.; Yeoh, G. H., Thermal Performance of Nanofluids in Microchannel Equipped with a Synthetic Jet Actuator. In International Heat Transfer Conference, International Heat Transfer Conference Digital Library: Kyoto, Japan, 2014; pp 3949-3963. 35


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 40

10. Kelly, M.; Yeoh, G. H.; Timchenko, V., On Computational Fluid Dynamics Study of Magnetic Drug Targeting. The Journal of Computational Multiphase Flows 2015, 7 (1), 43-56. 11. Kennedy, L. C.; Bickford, L. R.; Lewinski, N. A.; Coughlin, A. J.; Hu, Y.; Day, E. S.; West, J. L.; Drezek, R. A., A New Era for Cancer Treatment: Gold‐Nanoparticle‐Mediated Thermal Therapies. Small 2011, 7 (2), 169-183. 12. Watanabe, K.; Tanaka, E.; Ishii, H.; Nagao, D., The plasmonic properties of gold nanoparticle clusters formed via applying an AC electric field. Soft Matter 2018. 13. Nam, J.; Won, N.; Jin, H.; Chung, H.; Kim, S., pH-Induced Aggregation of Gold Nanoparticles for Photothermal Cancer Therapy. Journal of the American Chemical Society 2009, 131 (38), 13639-13645. 14. Albanese, A.; Chan, W. C. W., Effect of Gold Nanoparticle Aggregation on Cell Uptake and Toxicity. ACS Nano 2011, 5 (7), 5478-5489. 15. Mocanu, A.; Cernica, I.; Tomoaia, G.; Bobos, L.-D.; Horovitz, O.; Tomoaia-Cotisel, M., Self-assembly characteristics of gold nanoparticles in the presence of cysteine. Colloids and Surfaces A: Physicochemical and Engineering Aspects 2009, 338 (1), 93-101. 16. Agbangla, G. C.; Bacchin, P.; Climent, E., Collective dynamics of flowing colloids during pore clogging. Soft Matter 2014, 10 (33), 6303-6315. 17. Jezewski, W., Effect of long-range interactions on nanoparticle-induced aggregation. Physical Chemistry Chemical Physics 2016, 18 (33), 22929-22936. 18. Kroupa, M.; Vonka, M.; Kosek, J., Modeling the Mechanism of Coagulum Formation in Dispersions. Langmuir 2014, 30 (10), 2693-2702. 36


Page 37 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

19. Draine, B. T.; Flatau, P. J., Discrete-Dipole Approximation For Scattering Calculations. J. Opt. Soc. Am. A 1994, 11 (4), 1491-1499. 20. Marshall, J. S., Discrete-element modeling of particulate aerosol flows. Journal of Computational Physics 2009, 228 (5), 1541-1561. 21. Ounis, H.; Ahmadi, G.; McLaughlin, J. B., Brownian diffusion of submicrometer particles in the viscous sublayer. Journal of Colloid and Interface Science 1991, 143 (1), 266-277. 22. Xiao, H.; Sun, J., Algorithms in a Robust Hybrid CFD-DEM Solver for Particle-Laden Flows. Communications in Computational Physics 2015, 9 (2), 297-323. 23. Khare, P.; Wang, S.; Yang, V., Modeling of finite-size droplets and particles in multiphase flows. Chinese Journal of Aeronautics 2015, 28 (4), 974-982. 24. Hunter, R. J., Foundations of colloid science. Second edition. ed.; p xii, 806 pages. 25. Min, Y.; Akbulut, M.; Kristiansen, K.; Golan, Y.; Israelachvili, J., The role of interparticle and external forces in nanoparticle assembly. Nature Materials 2008, 7, 527. 26. Kroupa, M.; Vonka, M.; Soos, M.; Kosek, J., Size and Structure of Clusters Formed by Shear Induced Coagulation: Modeling by Discrete Element Method. Langmuir 2015, 31 (28), 7727-7737. 27. Thompson, D. W.; Collins, I. R., Electrolyte-Induced Aggregation of Gold Particles on Solid Surfaces. Journal of Colloid and Interface Science 1994, 163 (2), 347-354. 28. DeVoe, H., Optical Properties of Molecular Aggregates. I. Classical Model of Electronic Absorption and Refraction. The Journal of Chemical Physics 1964, 41 (2), 393-400.

37


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 40

29. McPeak, K. M.; Jayanti, S. V.; Kress, S. J. P.; Meyer, S.; Iotti, S.; Rossinelli, A.; Norris, D. J., Plasmonic Films Can Easily Be Better: Rules and Recipes. ACS Photonics 2015, 2 (3), 326-333. 30. Kreibig, U.; Vollmer, M., Optical Properties of Metal Clusters. Springer Berlin Heidelberg: 1995. 31. Ordal, M. A.; Bell, R. J.; Alexander, R. W.; Long, L. L.; Querry, M. R., Optical properties of fourteen metals in the infrared and far infrared: Al, Co, Cu, Au, Fe, Pb, Mo, Ni, Pd, Pt, Ag, Ti, V, and W. Appl. Opt. 1985, 24 (24), 4493-4499. 32. Novo, C.; Gomez, D.; Perez-Juste, J.; Zhang, Z.; Petrova, H.; Reismann, M.; Mulvaney, P.; Hartland, G. V., Contributions from radiation damping and surface scattering to the linewidth of the longitudinal plasmon band of gold nanorods: a single particle study. Physical Chemistry Chemical Physics 2006, 8 (30), 3540-3546. 33. Handley, D. A., 2 - Methods for Synthesis of Colloidal Gold A2 - Hayat, M.A. In Colloidal Gold, Academic Press: San Diego, 1989; pp 13-32. 34. Wampler, H. P.; Ivanisevic, A., Nanoindentation of gold nanoparticles functionalized with proteins. Micron 2009, 40 (4), 444-448. 35. Kimura, K.; Takashima, S.; Ohshima, H., Molecular Approach to the Surface Potential Estimate of Thiolate-Modified Gold Nanoparticles. The Journal of Physical Chemistry B 2002, 106 (29), 7260-7266. 36. Issa, R. I., Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics 1986, 62 (1), 40-65. 37. Chan, Q. N.; Medwell, P. R.; Nathan, G. J., Algorithm for soot sheet quantification in a piloted turbulent jet non-premixed natural gas flame. Experiments in Fluids 2014, 55 (10), 1827.

38


Page 39 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

38. Qamar, N. H.; Nathan, G. J.; Alwahabi, Z. T.; Chan, Q. N., Soot sheet dimensions in turbulent nonpremixed flames. Combustion and Flame 2011, 158 (12), 2458-2464. 39. Wang, C.; Chan, Q. N.; Zhang, R.; Kook, S.; Hawkes, E. R.; Yeoh, G. H.; Medwell, P. R., Automated determination of size and morphology information from soot transmission electron microscope (TEM)generated images. Journal of Nanoparticle Research 2016, 18 (5), 127. 40. Peng, Z.; Walther, T.; Kleinermanns, K., Influence of Intense Pulsed Laser Irradiation on Optical and Morphological Properties of Gold Nanoparticle Aggregates Produced by Surface Acid−Base Reactions. Langmuir 2005, 21 (10), 4249-4253. 41. Bashkatov, A. N.; Genina, E. A.; Kochubey, V. I.; Tuchin, V. V., Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm. Journal of Physics D: Applied Physics 2005, 38 (15), 2543. 42. Gu, X.; Timchenko, V.; Yeoh, G. H.; Dombrovsky, L. A.; Taylor, R. A. In Heat Generation in Gold Nanorods Solutions due to Absorption of Near-Infrared Radiation, ICHMT International Symposium on Advances in Computational Heat Transfer, Napoli, Italy, Begel House Inc.: Napoli, Italy, 2017.

39


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 40

Table of Contents Graphic

40


Study of Morphology and Optical Properties of Gold Nanoparticle

Recommend Documents