Femtosecond Plasmonic Laser Nanosurgery (fs-PLN) mediated by molecularly targeted gold nanospheres at ultra-low pulse fluences

Plasmonic Laser Nanosurgery (PLN) is a novel photomodification technique that exploits the near-field enhancement of femtosecond (fs) laser pulses in the vicinity of gold nanoparticles. While prior studies have shown the advantages of fs-PLN to modify cells, further reduction in the pulse fluence needed to initiate photomodification is crucial to facilitate deep–tissue treatments. This work presents an in-depth study of fs-PLN at ultra-low pulse fluences using 47 nm gold nanoparticles, conjugated to antibodies that target the epithelial growth factor receptor and excited off-resonance using 760 nm, 270 fs laser pulses at 80 MHz repetition rate. We find that fs-PLN can optoporate cellular membranes with pulse fluences as low as 1.3 mJ/cm2, up to two orders of magnitude lower than those used at lower repetition rates. Our results, corroborated by simulations of free-electron generation by particle photoemission and photoionization of the surrounding water, shed light on the off-resonance fs-PLN mechanism. We suggest that photo-chemical pathways likely drive cellular optoporation and cell damage at these off-resonance, low fluence, and high repetition rate fs-laser pulses, with clusters acting as local concentrators of ROS generation. We believe that the low fluence and highly localized ROS-mediated fs-PLN approach will enable targeted therapeutics and cancer treatment.

). Dashed lines represent aspect ratios (AR) of 1 and 2. Figure 2: Cell death after PLN with circularly polarized light. At higher powers (>2 mW), we find that the percentage of cells that remained viable (~2%) when exposed to circular-polarized laser pulses was significantly lower than the 5% that remained viable after linearly polarized irradiation.

Supplementary Figure 3: Uptake kinetics of PLN-induced transient pore formation.
The signal intensity of 10 kDa FITC-Dextran and siRNA mimics were monitored inside and outside single cells after irradiation with 30 scans of 2.6 mJ/cm 2 laser fluence. Exponential fits to the intensity rise were used to determine the hole-closure time, and are represented as solid lines. Dotted lines show prediction intervals at 95% confidence. Dextran diffuses slower than siRNA due to its larger size.

Supplementary Figure 4: Simulation of the absorption cross-sections (Q abs ) and
Poynting vector enhancements (Q Poyn ) in and around a 47 nm diameter single particle and a 4-particle cluster (d = 47 nm) in a water medium for various s/d values, excited with horizontally polarized 760 nm laser light. (a) Poynting vector enhancement (|S|/|S 0 |) around a single particle in water, showing over a 7-fold of maximum enhancement at the particle surface. (b,c,d) Poynting vector enhancement around a 4-particle cluster with packing factor of (b) s/d = 0.6 with a maximum enhancement of 10.0, (c) s/d = 0.8 with a maximum enhancement of 9, and (d) s/d = 1.0 with a maximum enhancement of 8.2, only varying by 10% between each case. These values result in an increase enhancement ranging from 1.43 (s/d = 0.6) to 1.2 (s/d = 1) over the 7-fold enhancement that a single particle experiences at the surface. A packing factor of 0.8 is taken for our simulations in the manuscript. For the Poynting vector simulations, we used a finite element model to solve Maxwell's equations in frequency domain. The complex permittivity of the gold nanoparticle is calculated using the Brendel-Bormann model 1 . Perfectly matched layer outside of the computational region is applied to absorb the scattered thermal radiation at the far-field in all directions. Also, non-uniform meshes are implemented to save computational time effectively and to produce convergence. The results showed a maximum enhancement of ~9, which is 1.29× over the 7 times enhancement that a single particle experiences at the surface. Figure 5: Thermal simulations for our maximum operating average input fluence of 9.1 mJ/cm 2 (7 mw). The temperature dynamics of a cluster with 4 gold nanoparticles (d = 47 nm) and surrounding water at an average input fluence of 9.1 mJ/cm 2 (pulse energy of 87.5 pJ) over the total exposure time of a spot in the field of view during a single line scan, which comprises of 378 consecutive pulses. For conservative calculations, we assumed the particles are located at the laser's focal center and experience the highest local fluences, namely the peak fluences, which are twice the average input fluence. The black solid line represents gold lattice temperature. The black dashed line represents the temperature evolution of water near the particle interface where the water temperature reaches the highest values. The blue dashed line represents the temperature at the cluster-center and shows that there is no lingering thermal build-up between pulses at this highest fluence. Figure 6: Schematic describing the various regimes of PLN driven by a 760 nm, 270 fs laser pulse interacting with a 47 nm gold nanoparticle pair for a packing factor of s/d = 0.6. We ran the simulations for s/d = 0.6 as the lower bound from our literature analysis, with lower s/d values expected increasing the enhancement effects. Fluences used in the simulations assume the particles are located at the focal center, and therefore experience twice the average input fluence locally. The increased field enhancement effects obtained at an s/d of 0.6 did not change the overall photomodification processes involved in our model where we used s/d = 0.8 ( Figure  6).

Supplementary Tables
Supplementary Table 1: Operating parameters for the literature surveyed in Figure 7 with femtosecond (fs) laser pulse widths. AuNS stands for gold nanospheres, ROS stands for reactive oxygen species.

Supplementary Note 1: Particle characterization
Here we provide a full characterization of our gold nanoparticles. We measured the linear absorption of each particle set using UV-Vis Spectroscopy (Beckman-Coulter) and cuvettes having path length 0.5 cm. The normalized linear absorbance spectrum is presented in Supplementary  Figure 1a. The particles exhibit a narrow absorption peak centered at 531 nm with 70 nm FWHM, which is characteristic of dipolar resonance. Broadening of the plasmon band results from the larger size distribution and ellipticity.
A gold nanoparticle concentration of 3.4 x 10 10 particles/mL was estimated from the UV-Vis spectra with less than a 6% error using where Q ext is the theoretically modeled extinction efficiency, A is the peak absorbance as found from the UV-Vis spectrum, r s and r i are the STEM measured radii in the short and long axis of the spheroid, respectively, and d 0 is the cuvette path length. Due to the spheroidal shape of the generated particles, Discrete Dipole Approximation (DDA) performed by Dr. R.K. Harrison was used to model the extinction efficiency [23][24][25][26][27] . DDA is a well-established technique for simulating the optical properties of non-spherical colloidal nanoparticles, because it can be utilized to account for particle size distribution and orientation. Non-spherical particle geometries can be simulated because the overall structure is organized via the superposition of small cube sections; across each cube, the induced dipole polarization from the incident light is assumed to be constant.
DDA model parameters to determine Q ext : The spheroidal shape was organized such that there were 128 dipoles across the short axis, which is equivalent to each dipole representing a 0.04 nm 3 volume. The complex gold refractive index at the irradiation wavelength was obtained from experimentally determined refractive index values for bulk gold and then corrected for sizerelated surface damping according to the Drude equation with a modified damping constant. Optical properties were simulated for isolated particles in a water environment at intervals of 20 nm across a 400 to 900 nm wavelength range.
Accounting for orientation effects with respect to polarization a and size distribution b of the sample, an average extinction coefficient of Q ext = 3.56 was found. If the spheroidal shape were neglected, the concentration would have been over-estimated. The physical characteristics of our nanoparticles are summarized in Supplementary Table 3. The physical dimensions of the generated nanoparticles were characterized with high-resolution field emission scanning transmission electron microscopy (STEM; Hitachi S-5500); R. K. Harrison obtained particle images. A 2μL aliquot from the stock solution was dried onto a copper grid (Cu-400CN; GridTech) at 120 °C, washed via a series of water droplets, and allowed to air dry. The inset of Supplementary Figure 1a provides a representative image of the generated particles at 150,000 x magnification. Image analysis of the particles were performed using ImageJ. We found that our batch of nanoparticles measured 54 ± 6 nm × 44 ± 4 nm (length × width, equivalent volume sphere d = 47 ± 5 nm). Supplementary Figure 1c provides a size distribution scatter plot for a sample set of 200 particles. The nanoparticles have an aspect ratio between 1 and 2 with an ellipticity of 1.24 ± 0.13, which is consistent with the formation of spheroidal particles, and a size distribution of approximately 10%.
The surface charge of our gold sol was determined through zeta potential measurements with a ZetaPlus® (Brookhaven Instruments Corporation). The nanoparticles were re-suspended at a concentration of 28 pM in a 1 mM KCl solution (pH 7.5). Ten zeta potential readings were obtained and averaged, revealing the GNS530 gold sol is stable, having an average surface charge of -27.6 ± 2.2 mV. Zeta potential indicates the degree of repulsion between adjacent particles, with potentials ~ -11 to -20 mV being the threshold below which agglomeration occurs 28,29 . Particle stability, i.e. monodispersity, is supported by both the high zeta potential and relatively narrow spectral bandwidth.

Supplementary Note 2: Characterizing Transient Pore Formation
To quantify PLN-induced pore formation, we monitored the uptake of 10-kDa Dextran and a siRNA mimic (14.9 kDa) after PLN using two-photon microscopy in separate experiments. A predetermined concentration of fluorescent probe (25 mM Dextran; 20 pM siRNA mimic) was added to the extracellular solution prior to irradiation. The cells were irradiated with 30 scans of 2.6 mJ/cm 2 incident fluence (2 mW power), chosen for its high rate of fluorescent probe retention. Supplementary Figure 3 presents the time dependent variation of extracellular and intracellular fluorescence signal levels as measured from the two-photon images of irradiated cells. Within the whole FOV, close to 75% of cells experienced material influx. Cells that were permeabilized showed an exponential increase in fluorescence intensity immediately after targeting. The fluorescence intensity then plateaued, indicating pore closure, or continued to increase towards the intensity values of the surrounding media, indicating permanent membrane dysfunction.
A simple diffusion-based model was used to estimate the PLN-produced hole radius 2 . Briefly, we assume that hole with a constant radius r is produced in the cell membrane, and is of course, much smaller than the cell radius. The increase in dye concentration within the cell,  can then be represented as an exponential rise using Fick's law of diffusion, where 0  is the final dye concentration within the cell and pore where cell V is the volume of the cell. We fit the measured rise of fluorescence within the optoporated cells to Eq. 2 to determine pore  , and estimate the hole radius r from Eq. 3. We found that the fluorescence intensity within optoporated cells increased with a rise time ( pore  ) of 136.7 ± 5.8 s for Dextran and 54.2 ± 2.0 s for siRNA, with comparable hole radii of 17.8 ± 4.9 nm and 17.9 ± 4.8 nm, respectively. The similar hole-radii estimates obtained using different transfection agents (with different D) points to the validity of the 1D diffusion model. The pore radius should allow us to transfect Dextran up to 2,000 kDa at 2.6 mJ/cm 2 , with a theoretical diameter of 12 nm 30 , albeit at slower diffusion rates.
In a different experiment, to verify if the intensity plateau corresponded with pore closure, i.e., membrane healing, the FITC-Dextran probe was added to the extracellular solution 180 ( influx at the 180 s mark, no influx was found beyond 240 s post irradiation. Comparable closure times and pore radii were observed in optoporation experiments using tightly focused NIR, femtosecond laser pulses at high repetition rates in cells sans nanoparticles 2,31 .

Supplementary Note 3: First order model for photoemission from particles and free electron generation in water
To explain our experimental observations and provide a comprehensive and quantitative description of electron emission and nanoparticle ablation using ultrashort laser pulses, we apply the generalized Fowler-DuBridge theory 32 to plasmonic nanospheres. We find particle Poynting vector enhancements using Finite Difference Frequency Domain solution of Maxwell's equations for a 4-particle cluster of gold spheres with d = 47 nm and s/d =0.8 (Supplementary Figure 4) while electron temperature evolutions are determined using the two-temperature model for laser-metal heating.
Conduction band electrons may be ejected to the surrounding medium from a material at elevated temperatures when their kinetic energies exceed the work function. Electron emission due to thermal kinetic energy, referred to as thermionic emission, exists at all temperatures down to zero Kelvin. In addition to thermionic emission, electrons may be released from surfaces by exposure to high intensity fields at the surface 33,34 . High intensity field strengths at a material surface may result in sufficient photon density so that one or more photons are absorbed by conduction band electrons at the material surface, causing their surface-normal kinetic energy components to exceed the barrier energy necessary for electron escape. For nanoparticles interacting with ultrafast laser pulses, the high absorption cross sections of plasmonic particles coupled with short pulses create conditions favorable for both thermionic and multi-photon emission processes.
The generalized Fowler-DuBridge theory of photoemission current density (J), as described by Bechtel et al. 34 incorporates all orders of emission processes as a series expansion with activation energy being a function of electron temperature (T e ), material work function (ϕ), and multiphoton process order (n): where A 0 is the Richardson coefficient (120 A/cm 2 K 2 ), a n is an order-specific constant, and F is the Fowler function. Thermionic emission is assumed to be a zeroth order process, and described by A slightly modified version of this theory describes specific J n 35,36 for photoemission current for three-photon process (n=3) Here c describes the three-photon ionization cross-section and is set at 10 -7 Acm 4 /MW 3 and ϕ is set to 3.72 eV 35,36 . The surface reflectivity (R) is not considered since the source term here ( I ) accounts for the local intensity in the near-field through Poyn Q . We can then calculate the total number of emitted electrons (N) by integrating the current density 0 3 J J  over the particle surface at each time step. To pursue this calculation, we first need to define the local intensity, which is modified by the particle near-field enhancement according to Gauss' Law along with the assumption of uniform charge distribution within the particle yields the induced electric field strength as a function of radial position, r, within the particle, where d is the particle diameter, ε is the nanoparticle permittivity, and V np is the volume of the nanoparticle. When the induced field at the particle surface reaches a critical value (E cr ) where the electrostatic repulsion of the lattice exceeds the lattice binding strength, the particle begins to ablate. Bulgakova and co-workers 36 estimated this value for gold to be 2.76 10 10 V/m. The validity of a uniform charge distribution is questionable for pulse durations on the order of electron-electron relaxation timescale of femtoseconds. However, this assumption will produce a conservative estimate of ablation initiation. Fowler-Nordheim type emission may be ruled out by examining the Anisimov parameter for laser intensities and wavelengths used here 37 .
In the case of multi-particle interactions, the Poyn Q associated with the particle will vary with the extent of inter-particle coupling. Our simulations showed a 1.29 × increase in the Poynting vector enhancement for the 4-particle cluster over a single particle (Supplementary Figure 4), with a maximum single particle enhancement of 7 at the particle surface. This increased enhancement per particle was used when simulating the photoemission rates from a single particle. Likewise, the increased near-fields and absorption was incorporated in our free electron generation model for water, discussed below.
The free-electron generation/plasma formation in water is handled similar to Vogel et al. 38 , where water is treated as a dielectric medium with an ionization potential   of 6.5 eV. We expect most of the free electron generation to take place within this shell for a given particle and ignore free electron diffusion into the volume from surrounding media. We do not expect the plasma from neighboring particles to interact over the duration of the pulse since plasma diffusion is quite slow. As such, simulations are performed for a single particle, with field enhancements accounting for multi-particle interactions. The time evolution of electron density  due to laser irradiance is calculated using a modified version of the rate equation of the generic form [38][39][40][41] process rate must account for ionization time . ion n    , determined by the number of photons (n) that need to be absorbed to acquire sufficient energy for impact ionization, and the mean free time τ (5 fs for plasma in gaseous plasma in water 43 ).
Therefore, the cascade ionization rate is determined using the electron density at time assuming a 50% probability of having a start electron in the focal volume, with the first order approximation expressed as where  represents the ionization rate per electron participating in the cascade 42 , where M is the mass of the water molecule.
Free electrons are lost from the interaction volume via diffusion out of the volume V and through recombination. In our case, we look at a small volume at the enhancement 'pole' beside the particle, where the enhancement in the particle is at its maximum. Assuming a spherical voxel of radius 125 nm, the resulting ionization rate due to diffusion is given by 38 where the characteristic diffusion length is Λ. Here we assume 1D diffusion radially across the voxel in one direction and use the diffusion length for a spherical shell given by   (1 / ) where R 1 and R 2 are the inner and outer radii.
The average kinetic energy of the free electrons is calculated as the average of the energy of photo-emitted electrons from the particle 41 (1.5 B e k T ), photoemitted electrons from water ( 0.5  ) and the energy of the electron produced by the collision ionization process ( 2.5  ). We assume that no diffusion from the neighboring voxel in water. We can expect the photoemitted electrons from the particle to still predominantly remain within the 5 nm voxel width over the duration of the pulse, considering that our electron temperatures are less than 0.5 eV. Simulations and experiments have shown plasma fronts moving under 20 nm for high density plasmas with electron temperatures in the 10-100 eV range 45,46 . For the recombination rate, we used an empirical value obtained by Docchio et al. 47 A fourth order Runga-Kutta scheme was used to numerically solve the rate equation for various input laser intensities. The various effects of particle photoemission and free electron generation in water are represented schematically in Figure 6 in the main text.
While our analysis builds on a widely used model for free electron generation in water, recent work has further improved on the model through a detailed rework based on a more accurate estimate of the band gap of water, which was updated to 9.5 eV 48 . The rate equation for strong field ionization was also updated to account for excitation to intermediate solvated states in the conduction band, and direct photoionization. The value of τ was also revised down to 0.9 fs (1,050 nm). In our case, our experimental threshold is portrayed with a bandgap that is too low by using a larger collision time is increased. Our model, in essence, provides a sense of the range of electron densities over which PLN operates for our input fluences.