# Enabling self-induced back-action trapping of gold nanoparticles in
metamaterial plasmonic tweezers

Theodoros D. Bouloumis  Okinawa Institute of Science and Technology Graduate
University, Onna, Okinawa 904-0495, Japan<EMAIL_ADDRESS>Domna
G. Kotsifaki  Okinawa Institute of Science and Technology Graduate University,
Onna, Okinawa 904-0495, Japan  Síle Nic Chormaic  Okinawa Institute of Science
and Technology Graduate University, Onna, Okinawa 904-0495, Japan

###### Abstract

The pursuit for efficient nanoparticle trapping with low powers has led to
optical tweezers technology moving from the conventional free-space
configuration to advanced plasmonic tweezers systems. However, trapping
nanoparticles smaller than 10 nm still remains a challenge even for plasmonic
tweezers. Proper nanocavity design and excitation has given rise to the self-
induced back-action (SIBA) effect offering enhanced trapping stiffness with
decreased laser power. In this work, we investigate the SIBA effect in
metamaterial tweezers and its synergy with the exhibited Fano resonance. We
demonstrate stable trapping of 20 nm gold particles for on-resonant and off-
resonant conditions with experimental trap stiffnesses as high as 4.18 $\pm$
0.2 fN/(nm$\cdot$ mW/$\mu$m2) and very low excitation intensity of about 1
mW/$\mu$m2. Simulations reveal the existence of two different groups of
hotspots per unit cell of the metamaterial array. The two hotspots exhibit
tunable trap stiffnesses, a unique feature that can allow for sorting of
particles and biological molecules based on their characteristics.

###### keywords:

American Chemical Society, LaTeX

Natural and Applied Sciences, Duke Kunshan University, No. 8 Duke Avenue,
Kunshan, Jiangsu Province, 215316, China LSP, FIB, SEM, MSP

Keywords: Metamaterial tweezers, self-induced back-action, plasmonic tweezers,
gold nano-particle trapping, Fano resonance

##  1 Introduction

Since the invention of optical tweezers by Arthur Ashkin1, 2 numerous trapping
platforms have been extensively studied and implemented for a variety of
applications3, 4 with trapping specimens ranging from metallic particles5, 6,
7 through biological samples8, 9, 10, to quantum dots11, 12, 13. One major
direction in optical trapping is to develop a platform that can efficiently
trap nanoparticles that are tens of nanometers or smaller in size, using as
low a light intensity as possible, thereby minimizing heating of the
specimen14, 15, 16. The most common configurations that meet these
requirements are based on metallic nanostructures in which one can make use of
the strong near-field forces that arise due to the excitation of surface
plasmons 17, 18, 19, 20.

However, the design of the nanostructure alone cannot improve the trapping
stiffness beyond a certain point. This realization led to the implementation
of the self-induced back-action (SIBA) trapping21. In SIBA trapping, the
plasmonic cavity resonance is tuned in such a way that the particle itself,
while trapped in the cavity, plays a dynamic role in the reconfiguration of
the intracavity intensity based on the spatial overlap of the particle and the
cavity22. This optomechanical coupling facilitates trapping of the particle
with a lower light intensity, thereby relaxing the intensity requirements
while providing a stable trap with less heating 23, 24, 25, 26.

Trapping gold nanoparticles (AuNPs) is of particular interest, due to the
versatility in their synthesis methods 27, 28 and their contribution in
various fields 29. In 1908, rigorous theoretical work was done by Gustav Mie
to mathematically describe the absorption of light by a subwavelength gold
sphere, based on Maxwell’s equations 30. Since then, owing to the great
tunability of plasmon resonances and their unique properties, AuNPs have
contributed to the development of fields within physics and material science
such as the surface-enhanced Raman scattering (SERS) 31 technique, enhanced
photovoltaic conversion efficiency via embedded particles in active layers of
solar cells 32, and quantum nanophotonic systems coupled to single photon
emitters resulting in ultrabright sources of nonclassical light 33. There are
numerous applications in other fields, too, such as in chemistry for catalysis
purposes 34 and gas sensing 35, in computational methods with the newly
established field of neuromorphic computing via AuNPs cluster assemblies 36,
and in biomedical research with applications in drug and gene delivery 37,
cancer treatment with photothermal therapy 38, vaccine development 39, and
many others 40, 41. All these promising applications have triggered interest
in trapping AuNPs of various sizes and investigations into their
optomechanical properties during manipulation and nanopositioning 5, 42, 24.

In our previous work 43, we demonstrated single trapping of 20 nm diameter
polystyrene particles using an array of metamolecules. The resulting trapping
stiffness was measured to be 8.65 fN/(nm$\cdot$mW) and is the highest reported
normalised trapping stiffness 43. Implementing metamaterial structures in a
plasmonic trapping configuration offers the advantage of being able to excite
and tune a very sharp Fano resonance, which is a result of destructive
interference between a broadband dipole resonance and a sharp quadrupole
resonance 44. This destructive interference leads to a nanocavity with very
small mode volume, offering very strong confinement of particles in the trap
44, 45, 46. A Fano resonator is also very sensitive to changes of the local
refractive index 47, making this platform an ideal candidate for combining
with SIBA trapping. Proper structure design could foster a synergy between the
Fano-resonance and the SIBA effect, resulting in fast optomechanical responses
that could facilitate long duration trapping events with low laser power.

In this work, we experimentally trapped AuNPs with a diameter of 20 nm using
the metamaterial tweezers introduced earlier 43. Theoretical analysis of the
excited plasmonic hotspots revealed a wavelength-dependent trapping stiffness
per hotspot. This offers the ability to tune the hotspots’ strength, with
applications in sorting nanoparticles of different sizes and refractive
indices. Additionally, we investigated the synergy between the metamaterial
tweezers and the SIBA effect and experimentally observed SIBA-assisted
trapping. To the best of our knowledge, there is no prior reporting of this in
the literature. A very high trapping stiffness was obtained for on-resonant
excitation (i.e., SIBA-assisted trapping), as well as for blue-detuned
excitation, with the highest trapping value being 4.18 $\pm$ 0.2
fN/(nm$\cdot$mW/$\mu$m${}^{2})$. This value is comparable to trapping of 20 nm
polystyrene particles with metamaterial tweezers where the SIBA effect could
not be utilized 43.

##  2 Results and discussion

A plasmonic metamaterial nanostructured device was used for the trapping
experiments. The device was fabricated using focused ion beam milling (FIB) on
a gold film, with the geometrical characteristics determined in advance
through finite element method simulations. AuNPs with a diameter of 20 nm were
dispersed in a heavy water/buffer mixture solution and a small amount of
surfactant was also used. The solution was placed onto the device, and the
sample was mounted on a custom-made microscope setup. The excitation laser was
incident from the glass slide of the sample. Detailed description of the
methods can be found in S1.

Figure 1a shows the simulated reflection, transmission and absorption spectra
of the fabricated metamaterial device (Figures 1b and 1(a)). The absorption
peak at 928 nm indicates the resonance of the device, hence the incident
wavelength required to excite the plasmonic field, when the device is immersed
in water solution. Microspectrophotometry measurements revealed that the
experimental resonance peak was in perfect agreement with the simulations -
since we cannot directly measure the absorption resonance experimentally, we
used a formula to calculate the experimental absorption resonance peak from
the transmission spectrum (see Figure S1b). Figure 1c(i) shows the simulated
empty-cavity electric field for 928 nm excitation light, propagating in the
$z$-axis, and linearly polarized along the $y$-axis. Each unit cell has three
hotspots where particles can be trapped; two symmetric and identical ones
along the $y$-axis, from now on referred to as Hotspot 1 (we make no
distinction between them) and one hotspot along the $x$-axis, from now on
referred to as Hotspot 2. We assume that a particle moving under Brownian
motion in the solution has equal probabilities of being in the vicinity of any
of the two hotspots and subsequently trapped. In Figures 1(c)(ii),(iii) the
electric field is plotted at the $yz$ and $xz$ planes while a AuNP 20 nm in
diameter is trapped at the equilibrium position of the Hotspot 1 (see S4 for
details on the equilibrium position). Due to its metallic nature, the particle
also has a plasmonic resonance and behaves like a trapped nanoantenna, leading
to significant field enhancement within the nanocavity (Figure 1c(iii)).

(a)

(b)

(c)

Figure 1: (a) The simulated reflection (R), transmission (T), and absorption
(A) spectra of the metamaterial structure indicating the empty cavity
resonance at 928 nm. (b) Scanning electron microscopy (SEM) image of the
fabricated metamaterial array. (c) The electric field generated at a unit cell
of the metamaterial structure, for the case of (i) an empty cavity and (ii,
iii) when a particle is trapped. The fields are for 928 nm excitation light at
the planes (i) $z~{}=~{}13$ nm (ii) $x~{}=~{}149$ nm, and (iii) $y~{}=~{}75$
nm. The particle is at the equilibrium position of Hotspot 1.

The existence of a trapped particle inside the nanocavity not only creates an
enhancement of the E-field, but also results in a frequency shift of the
cavity resonance given by 24

$\delta\omega(r_{p})=\omega_{c}\frac{\alpha_{d}}{2V_{m}\epsilon_{0}}f(r_{p})$
(1)

where $\omega_{c}$ and $V_{m}$ are the resonant frequency and the mode volume
of the empty cavity, respectively, $\alpha_{d}$ is the polarizability of the
nanoparticle, $\epsilon_{0}$ is the vacuum permittivity, and $f(r_{p})$ is the
cavity intensity profile. Using Equation 1, we calculated the expected cavity
shift in our configuration to be $\sim$4.06 nm (see S2). This is of the same
order of magnitude obtained through simulations, as shown in Figure 2. For a
particle trapped in either Hotspot 1 or Hotspot 2 the resulting cavity
redshift was 6 nm and 5 nm, respectively, and for consecutive trapping events
the cavity shifts by an additional step of 6 nm or 5 nm depending on where the
particle was trapped. This small discrepancy between the analytical and
numerical shift value is attributed to the approximate formula that was used
to calculate the mode volume 46, $V_{m}$. Based on the simulated value of the
cavity shift, we estimated that our system provides an optomechanical coupling
constant of approximately $G=2\pi\cdot 71$ GHz/nm (see S2), which is of the
same order of magnitude as other works 48, 24. The optomechanical constant
defines the response of the cavity (intracavity photon flux) due to the motion
of the particle in the trap and subsequent changes of momentum. It is thus
important for SIBA trapping to have as high an optomechanical coupling
constant as possible to achieve strong confinement of the particle. It is
expected that smaller particles would have a lower optomechanical coupling
constant due to a smaller change of momentum; however, in our case, the
optomechanical coupling value was very high considering the small size of the
particle. We attribute this to the Fano resonator being very sensitive to the
motion of the particle (changes of the local refractive index).

The particle-enabled cavity shift is the principle upon which the SIBA effect
is based. In the simplest case of one trapped particle, to enable the SIBA
effect one should excite the metamaterial with light at 930 nm. This is
approximately the wavelength around which the cavity shift is symmetrical to
the empty cavity resonance (Fig. 2 black vertical line). Under this condition,
the particle would be trapped efficiently with a lower laser intensity and
with its Brownian motion restricted by the intracavity field increasing or
reducing in a dynamical way according to its motion 21.The additional fast
optomechanical response of the Fano nanocavity makes this synergy between the
metamaterial tweezers and SIBA ideal for trapping nanoparticles.

Figure 2: Simulated absorption spectra for different configurations of trapped
particles. Empty cavity absorption resonance: red line with dots. Cavity
shifts created by one (blue and orange lines), two (magenta line) and three
(green line) particles trapped at the same unit cell in various combinations
of the three hotspots. The cavity shifts from the empty cavity resonance (928
nm) are indicated by the arrows in corresponding colors.The solid black line
at around 930 nm indicates the wavelength around which there is a symmetrical
cavity shift between the empty cavity resonance and the shifted cavity
resonance resulting from a single particle trapped.

To further map the optical forces exerted on the nanoparticles around the
hotspots and identify the equilibrium positions of the trapped particles,
extensive simulations were performed. In Figure 3(a) the trapping stiffnesses
for the two hotspots are plotted as derived from the simulated optical forces
(see S4 for the simulated optical forces and potentials). An interesting
behavior is observed, where the trapping stiffnesses around Hotspot 1 are
increasing as we scan the excitation wavelength from 920-930 nm, but an
opposite, decreasing trend is observed for the stiffnesses around Hotspot 2
(same trends observed for the optical forces and potential depths,
respectively). Additionally, particles trapped at Hotspot 2 experience higher
trapping stiffness, $k_{2}$, than particles trapped at Hotspot 1, $k_{1}$, for
every wavelength. However, the two hotspots exhibit similar trapping
stiffnesses when excited with 930 nm light.

To interpret this behavior, we simulated the E-field magnitudes at the two
hotspots for both an empty cavity and a trapped particle (Fig. 3(b)). Despite
the simulated absorption resonance being at 928 nm (Fig. 2), Hotspot 1
exhibits a maximum E-field intensity at 930 nm, while Hotspot 2 has a maximum
at 919 nm when the cavity is empty. Additionally, the maximum E-field
intensity for both hotspots has a similar magnitude. However, the situation is
different when a particle is trapped. The presence of the metallic particle in
the plasmonic nanocavity creates an enhancement of the E-field by
approximately 6 times (Fig. 3(b)). The total E-field in the cavity is stronger
for Hotspot 2 than for Hotspot 1. This is because the slot width at Hotspot 2
($w_{2}~{}\approx~{}33$ nm) is smaller than the slot width at Hotspot 1
($w_{1}~{}\approx~{}39$ nm), see Figure 1(a). This leads to a stronger E-field
and better confinement of the nanoparticle (Fig. 4(a)). This also explains the
higher trapping stiffness at Hotspot 2.

The increasing and decreasing trend of the trapping stiffnesses at Hotspots 1
and 2 (Fig. 3(a)) also follows the generated E-field trend. It is related to
the different E-field resonance values for each hotspot, as indicated by the
dashed lines in Figure 3(b). This offers the ability to excite the hotspots
with a different ratio of resulting trapping stiffnesses, hence the
possibility of sorting particles based on their size, shape, and refractive
index at different hotspots, e.g. bigger particles could be trapped in the
hotspot with the lowest trapping force and larger slot width, while smaller
particles that need higher optical forces could be trapped in the other
hotspot. The optical forces as a function of excitation wavelength are shown
in Figure 4(b) and exhibit the same resonance features as the E-field.

(a)

(b)

Figure 3: (a) Simulated trapping stiffnesses per wavelength for particles
trapped in Hotspot 1 (k1 \- blue points) or Hotspot 2 (k2 \- orange points)
calculated from the associated optical forces (see S4). (b) Electric field as
a function of the excitation wavelength for the empty cavity case at the
center of Hotspot 1 (E0,1) and Hotspot 2 (E0,2), and when a particle is
trapped Ep,1 and Ep,2, respectively.

In experiments, we performed trapping with the same excitation wavelengths
used for the simulations, with incident powers ranging from 1.7 - 4.3 mW. For
each wavelength, up to 10 experiments were performed for increasing incident
laser power with a step of 0.5 mW with a total number of about 100 trapping
events per wavelength. The trapping events were recorded and analyzed based on
the trapping transient analysis 49, in order to derive the experimental
trapping stiffness, $k_{exp}$. In Figure 4 we present the probability
distributions of the trapping stiffnesses normalized to the incident laser
intensity. Data were fitted with a normal distribution function and, for the
wavelengths 920 nm - 928 nm, two distribution functions were required in order
to fit the experimental results. This is in agreement with the simulations
where particles trapped at Hotspot 1 (blue curves) experience a weaker
trapping stiffness compared with particles trapped at Hotspot 2 (orange
curves). However, for trapping with 930 nm wavelength (Fig. 4(d)), a single
distribution function is sufficient since, at this wavelength, the two
hotspots have a similar trapping stiffness and are indistinguishable. This is
confirmed by Figure 3(a). It is also worth noting that the probability of
trapping at Hotspot 1 is significantly higher than at Hotspot 2; this is
expected since each metamolecule is comprised of twice as many hotspots of
type 1.

The averaged normalized experimental trapping stiffnesses, $k_{1,norm}$ and
$k_{2,norm}$ (Fig. 4(e)), although following the same trend as the theoretical
ones (Fig. 3(a)), are approximately 5 - 8 times lower in magnitude. The
maximum normalized trapping stiffness obtained was 4.18 fN/(nm $\cdot$ mW/
$\mu$m2) for excitation at 920 nm which is of the same order of magnitude as
that previously reported for polystyrene particle trapping 43. However, the
deviation from the simulated values indicates the existence of additional
destabilisation forces that were not considered in the simulations and may
affect the trapping process.

(a)

(b)

(c)

(d)

(e)

Figure 4: Probability distributions of the normalized experimental trapping
stiffnesses for the excitation wavelengths (a) 920 nm, (b) 925 nm,(c) 928 nm,
and (d) 930 nm, with the Gaussian fits corresponding to trapping events at the
Hotspot 1 (blue curves) and Hotspot 2 (orange curves). (e) The average
normalized experimental trapping stiffnesses per excitation wavelength for the
two hotspots. The values were obtained from the Gaussian fittings of the
histograms above.

The situation is different when trapping with the 930 nm excitation laser; the
average trapping stiffness magnitude, $k_{930,norm}$, significantly deviates
from the simulated behavior, exhibiting the lowest value amongst all
wavelengths (Fig. 4(e)). To explain this and get a clearer view on the
experimental findings, in Figure 5 we plotted the initial non-normalized
trapping stiffnesses as a function of the incident laser power and trapping
wavelength, based on the statistical analysis for the two hotspots. Although
the behavior is quite complex, an increasing trend of the trapping stiffness
with increasing incident laser power is observed for the wavelengths 920 nm -
928 nm for both hotspots (Fig. 5a,b). This trend is clearly not linear, as
theoretically predicted, and it is attributed to the collective effect of
different forces acting on the nanoparticles. Interestingly, the behavior of
the 930 nm wavelength is less complicated and in agreement with previous works
where the SIBA effect was investigated 24.

As mentioned earlier, exciting the structure with 930 nm light is strongly
desired because at this wavelength the SIBA contribution is enabled, resulting
in a high trapping stiffness with a low laser intensity. Indeed, as shown in
Figure 5a,b, trapping with 930 nm light results in two distinct regimes of
trapping stiffness values. These are the result of a combination of the
optical forces along with the SIBA effect acting on the nanoparticle. The
highest normalized trapping stiffness when trapping with 930 nm was measured
to be 2.92 fN/(nm$\cdot$mW/$\mu$m2) (Fig. 5c,d) for an incident intensity as
low as 0.61 mW/$\mu$m2. In the low power regime, the SIBA effect is the
dominant mechanism keeping the particle trapped since the optical forces are
not strong enough 22. As the incident power increases, the optical forces play
a more dominant role, and the trapping stiffness decreases because the scaling
between laser power and particle confinement, upon where the SIBA is based,
becomes less efficient 22. After a certain power threshold, which in our
experiment lies at 2.82 mW (1 mW/$\mu$m2 intensity), the SIBA effect vanishes
and the optical forces linearly scale up with the power. In principle, based
on Figures 5c,d, we observe that low optical intensities result in higher
trapping stiffnesses when exciting on-resonant for the SIBA (930 nm); however,
for higher intensities the blue-detuned excitation of the cavity provides
stiffer traps.

(a)

(b)

(c)

(d)

Figure 5: Experimental trapping stiffness as a function of incident laser
power for particles trapped at (a) Hotspot 1 and (b) Hotspot 2, and normalized
trapping stiffness as a function of the excitation wavelength, respectively
(c) and (d).

At Hotspot 2, due to better confinement of the E-field, the resulting
plasmonic field was stronger than that at Hotspot 1. Even for low powers and
off-resonant conditions, high optical forces were exerted on the nanoparticles
(Fig. 5b,d). Despite the slightly higher trapping stiffnesses that were
obtained when exciting at off-resonant conditions (920 - 928 nm), trapping
with SIBA was still favorable due to the dynamic reconfiguration of the trap
that allowed the particle to be trapped in low intra-cavity intensity.

An interesting observation in Figures 5a,b is the fluctuations in the trapping
stiffness with increasing incident power and the deviation from the linear
increase that is theoretically predicted. This is attributed to the effect of
additional forces, such as thermal and fluid-related forces, that may play a
significant role in the trapping mechanism.

#### Thermal forces:

The height of the microwell containing the gold particles was 120 $\mu$m.
Photo-induced heating resulting from intrinsic losses of the metallic
nanostructures generates thermal convection currents and thermophoresis that
could influence the trapping process 16. Specifically, as the laser trapping
intensity increased, the optically generated temperature gradients could
induce fluid flows that quickly delivered the nanoparticles above the
nanostructures. To maximize particle-trapping probabilities, the nanoparticle
should be delivered with low fluid velocity to the hotspot where it would be
trapped under the short-range plasmonic fields. However, when the height of
the solution’s microwell is larger than 20 $\mu$m the fluid velocity is large,
thereby preventing stable trapping 11, 50. In addition, thermophoresis induced
by the temperature gradients could alter thermophoretic migration of particles
50. For temperatures higher than a critical temperature, the Soret coefficient
sign changes 51 and thermophoresis repels particles away from a high-
temperature location 16, hindering trapping stability.

#### Fluid-related forces:

An important factor that may impact the trapping performance is the role of
the Triton X-100 in the liquid solution. It has been reported that a small
amount of non-ionic surfactant generates thermodiffusion with opposite sign
from the optical force leading to fluctuations of the particle’s position
inside the trap and a lower trap stiffness magnitude 15.

#### Electrostatic and scattering forces:

An additional effect may be present due to electrostatic repulsive forces as
the gold particle (-24 mV zeta potential) approaches the gold surface, with
both being negatively charged 52. This can result in a decreased measured
trapping stiffness since our evaluation method is based on the motion of the
particle as it enters the trap. Finally, upon illumination of the gold
structure, light scattering may contribute to the destabilisation of the
particle and the deviation from the linear increase of the trapping stiffness.

As the combination of the effects above strongly depend on the particle and
solution properties, as well as local temperature effects, the trapping
landscape can get quite complicated and an optimized trapping scheme is
required that would take advantage of all the positive phenomena and minimise
the negative ones. However, this is beyond the scope of this paper and it will
constitute a future work.

##  3 Conclusions

Trapping of gold nanoparticles was performed via metamaterial plasmonic
tweezers. The investigated trapping conditions were for on-resonant and blue-
detuned wavelengths to enable the SIBA effect. Very high experimental trapping
stiffness values were observed, with on-resonant trapping exhibiting the
highest value at 2.92 $\pm$ 0.42 fN/(nm$\cdot$mW/$\mu$m2) for an excitation
laser intensity as low as 0.61 mW/$\mu$m2. For the off-resonant case, the
highest value was 4.18 $\pm$ 0.2 fN/(nm$\cdot$mW/$\mu$m2) but for an
excitation intensity about 62% higher compared with the on-resonant case. The
experimental results had the characteristic feature of SIBA-assisted trapping,
where the maximum trapping stiffness is exhibited for the lowest excitation
intensity. For the SIBA-resonant wavelength of 930 nm, two distinct regimes
were observed; the low-intensity regime where the SIBA was the main mechanism
contributing to the trap and the high-intensity regime where the SIBA vanished
and the optical forces were dominating in the trapping mechanism, in agreement
with theory. Blue-detuned trapping exhibited a more complex behavior that
deviates from the linear increase of the trapping stiffness with increased
excitation intensity. Additional forces, such as thermal, fluid, and
electrostatic forces, may have contributed to this deviation. Extensive
simulations were also performed and revealed the existence of two groups of
hotspots in each unit cell of the metamaterial array. With varying excitation
wavelength a different ratio of E-field intensities could be excited at each
hotspot, thus resulting in traps of different strengths enabling the sorting
of particles based on their size, shape, and refractive index. This could
prove extremely useful in biological applications for sorting, for example,
healthy from damaged distorted cells, and investigating their morphological,
kinetic, and other characteristics.

###  3.1 Associated Content

#### Supporting Information

Details on the sample fabrication, solution preparation and simulations can be
found in the Methods section of the Supplementary Information. Additional
simulation figures and analytical calculations of the cavity shift and the
optomechanical coupling constant are also included.
S1: Methods
S2: Experimental and simulated structure
S3: Cavity shift and optomechanical coupling constant
S4: Simulated optical forces
S5: Simulated forces and potentials at the two hotspots

###  3.2 Authors Information

Theodoros D. Bouloumis ID :0000-0002-5264-7338,
Light – Matter Interactions for Quantum Technologies Unit, Okinawa Institute
of Science and Technology Graduate University, 1919-1 Tancha, Onna-son,
Okinawa, 904-0495, Japan<EMAIL_ADDRESS>

Domna G. Kotsifaki ID :0000-0002-2023-8345,
Natural and Applied Sciences, Duke Kunshan University, No. 8 Duke Avenue,
Kunshan, Jiangsu Province, 215316, China,
Light – Matter Interactions for Quantum Technologies Unit, Okinawa Institute
of Science and Technology Graduate University, 1919-1 Tancha, Onna-son,
Okinawa, 904-0495, Japan,
<EMAIL_ADDRESS>

Síle Nic Chormaic ID :0000-0003-4276-2014,
Light – Matter Interactions for Quantum Technologies Unit, Okinawa Institute
of Science and Technology Graduate University, 1919-1 Tancha, Onna-son,
Okinawa, 904-0495, Japan<EMAIL_ADDRESS>
###  3.3 Notes

TDB and DGK conceived the idea of this work, DGK prepared the solutions, TDB
performed the experiments, simulations, and analyzed the data. SNC supervised
all stages of the work. All authors contributed to writing the paper.

The authors declare no competing financial interests.

The authors thank M. Ozer for technical support, J. Keloth for useful
discussions, and N. Kokkinidis for his input regarding simulations. They also
acknowledge P. Puchenkov and J. Moren from the Scientific Computing and Data
Analysis Section, and the Engineering Section of the Research Support Division
at OIST. This work was partially supported by the Okinawa Institute of Science
and Technology Graduate University. TDB acknowledges funding from the Japan
Society for the Promotion of Science (JSPS) KAKENHI Grants-in-Aid (Grant No.
22J10196) and from the International Society for Optics and Photonics (SPIE)
for the Optics and Photonics Education Scholarship 2021.

###  S0.4 S1 Methods

The simulations were performed on the software COMSOL Multuphysics with the
Wave Optics module, which is based on the finite element method (FEM). The
geometry of the simulated structures was built entirely on COMSOL with the
dimensions matching the experimental dimensions measured via SEM (Fig. S1a).
The mesh was set to finer - custom with minimum size 1.2 nm, maximum size 44.3
nm and resolution of narrow regions at 0.7 nm. The Floquet periodicity was
used with a periodic port from which we illuminated with a plane wave.
Initially, we performed a frequency domain study in order to simulate the
reflection, transmission, and absorption spectra which gave us the theoretical
resonance wavelength of the plasmonic metamaterial cavity. Next, we introduced
a gold nanoparticle with 20 nm diameter, surrounded by a 2.5 nm thick water
shell and scanned its position along the three directions to calculate the
optical forces exerted on the nanoparticle. For the force calculation, the
Maxwell stress tensor 53 was used with surface integration around the water
shell.

The metamaterial structures were fabricated on a 50 nm gold thin film (PHASIS,
Geneva, BioNano) using focused ion beam (FIB) milling. The FIB system was a
FIB-FEI Helios Nanolab G3UC and the structures were etched at 80 nm depth to
ensure that the structures were cut all the way through the gold film
uniformly. The voltage and beam current used for etching were 30 kV and 2 pA,
respectively and the volume per dose was 0.27 $\mu$m3/nC. The scanning
electron microscopy (SEM) images were acquired with the same system. Each
metamaterial array consisted of 15 $\times$ 15 units with a nominal period of
400 nm. The geometrical characteristics of the metamaterial are shown in
Figure S1a. After fabrication, the device underwent oxygen plasma treatment
for 3 minutes to remove gallium residue that was deposited on it during the
milling process 54. The transmission spectra were obtained with a
microspectrophotometer (CRAIC, 20/30 PV) with the device immersed in a miliQ
water environment.

For the particle solution, we prepared an initial mixture of heavy water (D2O,
Sigma-Aldrich) and phosphate buffered saline (PBS) 10x buffer (pH = 7.4,
Invitrogen) with ration 9/1, and dispersed 20 nm diameter gold particles
(Sigma-Aldrich No:753610) in it with a concentration of 0.3% v/v. Surfactant
Triton X-100 (Sigma-Aldrich) with a concentration of 1.2% v/v was also added
to avoid the aggregation of particles. The solution was sonicated for 10
minutes and then 8.7 $\mu$L of the final solution was placed inside a seal
spacer with 9 mm diameter and 120 $\mu$m thickness (Invitrogen) on the gold
film where the nanostructures were fabricated. Finally, the sample was sealed
from the top with a cover glass, forming a microwell.

The sample was mounted on a custom-made inverted microscope setup. The
excitation laser was incident from the bottom side of the sample (SiO2 glass
substrate) and focused on the structures through a high numerical aperture (NA
=1.25) 100x oil immersion objective lens (Nikon MRP71900, E Plan). The
transmitted light was collected from the other side through a second objective
lens (Nikon CF Plan BD DIC 20x) and detected on an avalanche photodetector
(APD430A/M Thorlabs) for real-time monitoring of the signal. The APD was
connected to a data acquisition module (DAQ USB-6363, NI) and interfaced using
LabVIEW code to acquire measurements at a 100 kHz rate.

###  S0.5 S2 Experimental and simulated structure

(a)

(b)

(c)

(d)

Figure S1: (a) SEM image of a unit cell of the metamaterial structure. The
dimensions are: $w_{1}$ = 39.83 $\pm$ 0.77 nm, $w_{2}$ = 33.33 $\pm$ 0.66 nm,
$\alpha$ = 307.55 $\pm$ 0.55 nm, $t$ = 157.23 $\pm$ 0.54 nm, $g$ = 108.80
$\pm$ 0.91 nm, $D$ = 416.14 $\pm$ 1.00 nm. (b) Experimental transmission
spectrum. (c) The simulated design based on the experimental sample. The
etching of the sample goes about 20-30 nm into the glass substrate to ensure
uniform cut throughout the gold film. (d) Simulated potentials along the
z-axis for the two hotspots and for an excitation wavelength of 930 nm.

The experimental dimensions of the unit array are shown in Figure 1(a) and
defined as $w$, the slot width, $\alpha$, the length of the long arm, $t$, the
length of the short arm, $g$, the separation between the "C" and "I" parts of
the structure, and $D$ the period of the array.

From the experimental transmission spectrum (Fig. 1(b)), we can calculate the
absorption peak wavelength, which indicates the Fano resonance, from 55

$\lambda_{0}=\frac{2\lambda_{peak}\lambda_{dip}}{\lambda_{peak}+\lambda_{dip}}=928.3~{}\text{nm},$
(S1)

where $\lambda_{peak}$ and $\lambda_{dip}$ correspond to the wavelengths of
the Fano peak indicated in Figure 1(b).

###  S0.6 S3 Cavity shift and optomechanical coupling constant

#### Cavity shift:

Following work by Mestres et al. 24, a trapped nanoparticle in a plasmonic
cavity induces a plasmonic resonance shift given by

$\delta\omega(r_{p})=\omega_{c}\frac{\alpha_{d}}{2V_{m}\epsilon_{0}}f(r_{p}),$
(S2)

where $\omega_{c}$ is the empty cavity resonance, $\alpha_{d}$ the effective
polarizability of the particle, $V_{m}$ mode volume of the cavity, and
$f(r_{p})$ the cavity intensity profile.

The cavity resonance according to simulations is at $\lambda_{c}~{}=~{}928$nm
or $\omega_{c}~{}=~{}2~{}\cdot~{}10^{15}$$~{}rad/sec$ and the mode volume of
the cavity is given with a rough approximation by 46

$V_{m}=2(\alpha+t-w_{2})\frac{w_{1}+w_{2}}{2}\cdot h$ (S3)

where $\alpha,t,w_{1},w_{2}$ are the geometrical characteristics noted in
Figure 1(a) and $h$ is the thickness of the gold film, in our case 50 nm. The
average value of $w_{1}$ and $w_{2}$ is assumed as the general slot width for
the equation. The calculated value is
$V_{m}~{}=~{}1.58~{}\cdot~{}10^{-3}~{}\mu m^{3}$.

For the polarizability of the particle, the radiative reaction correction to
the Clausius - Mossotti relation was taken into account. For a particle with
radius $r~{}\ll~{}\lambda$ the effective polarizability, $\alpha_{d}$, is
given by 56

$\alpha_{d}=\dfrac{\alpha_{0}}{1-\frac{\epsilon_{r}-1}{\epsilon_{r}+2}[(k_{0}r)^{2}+\frac{2i}{3}(k_{0}r)^{3}]},$
(S4)

where $\alpha_{0}$ is the polarizability of the particle given by the Clausius
- Mossoti relation

$\alpha_{0}=4\pi r^{3}\epsilon_{0}\dfrac{\epsilon_{r}-1}{\epsilon_{r}+2},$
(S5)

$k_{0}$ is the wavenumber, $\epsilon_{0}$ the vacuum dielectric permittivity,
and $\epsilon_{r}$ is the particle’s dielectric permittivity taken from 57 to
be $\epsilon_{r}~{}=~{}-35.618+1.7121i$.

The intensity profile of the cavity, $f(r_{p})$, was taken equal to 1, which
is its maximum normalized value for maximum intensity obtained when exciting
the cavity with the absorption resonant wavelength of 928 nm 24.

Substituting the values of Equations S3 and S4 in Equation S2, the resulting
cavity shift expressed as a wavelength shift is $\delta\lambda~{}=~{}4.06$ nm.

#### Optomechanical coupling constant:

The optomechanical coupling constant is given approximately by 24

$G~{}\approx~{}\frac{\delta\omega}{\delta x},$ (S6)

where $\delta\omega$ is the cavity shift calculated above equal to
2$\pi\cdot$1.41 THz, and $\delta x$ is the displacement of the particle from
the equilibrium position assumed to be equal to the particle’s diameter as a
convention. Thus, G = 2$\pi\cdot$ 71 GHz/nm.

###  S0.7 S4 Simulated optical forces and potentials

(a)

(b)

(c)

(d)

Figure S2: Simulated optical forces (a), (c) and their corresponding
potentials (b),(d) experienced by a gold nanoparticle with 20 nm diameter
around Hotspot 1. The calculations were performed along the (a),(b) $z$\- and
(c),(d) $y$-axis for excitation wavelengths 920 nm, 925 nm, 928 nm, and 930
nm.

In Figure S2, the optical forces and potentials are plotted along the $z$ and
$y$ axes, for a particle around Hotspot 1 and for excitation wavelengths 920
nm, 925 nm, 928 nm, and 930 nm. The potential depth and optical forces exerted
on the particle increase as we scan the excitation wavelength from 920 nm to
930 nm, with 930 nm creating the strongest force. The maximum force along the
$z$-axis for the case of 930 nm is around 1.5 pN, normalised for 1 mW
excitation power. For all the wavelengths, the maximum force is exerted at the
position $z=-3$ nm which is inside the water domain 3 nm away from the surface
of the gold film, whereas the equilibrium position is at $z=12.6$ nm, inside
the gold domain which is etched and filled with water (Fig. 1(c)). The force
along the $x$-axis, although not presented here, was calculated and it was
always pushing the particle toward the external wall of the cavity.

(a)

(b)

(c)

(d)

Figure S3: Simulated optical forces (a), (c) and their corresponding
potentials (b),(d) experienced by a gold nanoparticle with 20 nm diameter
around Hotspot 2. The calculations were performed along the (a),(b) $z$\- and
(c),(d) $y$-axis for excitation wavelengths 920 nm, 925 nm, 928 nm, and 930
nm.

Similarly, the forces and potentials are plotted for Hotspot 2 in Figure S3;
however, in this case, the potential depths and optical forces are decreasing
as we scan from 920 nm to 930 nm which is the opposite situation compared to
Hotspot 1. Additionally, for all wavelengths, the forces at Hotspot 2 are
slightly higher than those at Hotspot 1, with the maximum force and
equilibrium positions at $z=-2.5$ nm and $z=13.5$ nm, respectively. Note that
the scale bar in all the E-field insets is the same as the one in Figure 4(a).

###  S0.8 S5 Simulated forces and potentials at the two hotspots

(a)

(b)

Figure S4: (a) Simulated electric field at a unit cell with particles trapped
at the hotspots. The narrower slit creates better particle confinement and
higher electric field enhancement. (b) Simulated forces acting on a gold
nanoparticle at the equilibrium position of Hotspot 1 (blue line) and Hotspot
2 (orange line).

## References

  * Ashkin 1970 Ashkin, A. Acceleration and trapping of particles by radiation pressure. _Phys. Rev. Lett._ 1970, _24_ , 156–159
  * Ashkin et al. 1986 Ashkin, A.; Dziedzic, J. M.; Bjorkholm, J. E.; Chu, S. Observation of a single-beam gradient force optical trap for dielectric particles. _Opt. Lett._ 1986, _11_ , 288–290
  * Polimeno et al. 2018 Polimeno, P.; Magazzù, A.; Iatì, M. A.; Patti, F.; Saija, R.; Esposti Boschi, C. D.; Donato, M. G.; Gucciardi, P. G.; Jones, P. H.; Volpe, G.; Maragò, O. M. Optical tweezers and their applications. _Journal of Quantitative Spectroscopy and Radiative Transfer_ 2018, _218_ , 131–150
  * Bouloumis and Nic Chormaic 2020 Bouloumis, T. D.; Nic Chormaic, S. From Far-Field to Near-Field Micro- and Nanoparticle Optical Trapping. _Applied Sciences_ 2020, _10_
  * Wang et al. 2009 Wang, K.; Schonbrun, E.; Crozier, K. B. Propulsion of Gold Nanoparticles with Surface Plasmon Polaritons: Evidence of Enhanced Optical Force from Near-Field Coupling between Gold Particle and Gold Film. _Nano Letters_ 2009, _9_ , 2623–2629
  * Zhang et al. 2010 Zhang, W.; Huang, L.; Santschi, C.; Martin, O. J. F. Trapping and Sensing 10 nm Metal Nanoparticles Using Plasmonic Dipole Antennas. _Nano Letters_ 2010, _10_ , 1006–1011
  * Fränzl and Cichos 2022 Fränzl, M.; Cichos, F. Hydrodynamic manipulation of nano-objects by optically induced thermo-osmotic flows. _Nature Communications_ 2022, _13_ , 656
  * Righini et al. 2009 Righini, M.; Ghenuche, P.; Cherukulappurath, S.; Myroshnychenko, V.; de abajo, F. J. G.; Quidant, R. Nano-optical trapping of Rayleigh particles and Escherichia coli bacteria with resonant optical antennas. _Nano letters_ 2009, _9 10_ , 3387–91
  * Verschueren et al. 2019 Verschueren, D.; Shi, X.; Dekker, C. Nano-Optical Tweezing of Single Proteins in Plasmonic Nanopores. _Small Methods_ 2019, _3_ , 1800465
  * Peng et al. 2021 Peng, X.; Kotnala, A.; Rajeeva, B. B.; Wang, M.; Yao, K.; Bhatt, N.; Penley, D.; Zheng, Y. Plasmonic Nanotweezers and Nanosensors for Point-of-Care Applications. _Adv Opt Mater_ 2021, _9_
  * Xu and Crozier 2019 Xu, Z.; Crozier, K. B. All-dielectric nanotweezers for trapping and observation of a single quantum dot. _Optics express_ 2019, _27 4_ , 4034–4045
  * Hong et al. 2021 Hong, C.; Yang, S.; Kravchenko, I. I.; Ndukaife, J. C. Electrothermoplasmonic Trapping and Dynamic Manipulation of Single Colloidal Nanodiamond. _Nano Letters_ 2021, _21_ , 4921–4927
  * Jiang et al. 2021 Jiang, Q.; Roy, P.; Claude, J.-B.; Wenger, J. Single Photon Source from a Nanoantenna-Trapped Single Quantum Dot. _Nano Letters_ 2021, _21_ , 7030–7036
  * Roxworthy et al. 2012 Roxworthy, B. J.; Ko, K. D.; Kumar, A.; Fung, K. H.; Chow, E. K. C.; Liu, G. L.; Fang, N. X.; Toussaint, K. C. Application of plasmonic bowtie nanoantenna arrays for optical trapping, stacking, and sorting. _Nano letters_ 2012, _12 2_ , 796–801
  * Jiang et al. 2020 Jiang, Q.; Rogez, B.; Claude, J.-B.; Baffou, G.; Wenger, J. Quantifying the Role of the Surfactant and the Thermophoretic Force in Plasmonic Nano-optical Trapping. _Nano letters_ 2020,
  * Kotsifaki and Nic Chormaic 2022 Kotsifaki, D. G.; Nic Chormaic, S. The role of temperature-induced effects generated by plasmonic nanostructures on particle delivery and manipulation: a review. _Nanophotonics_ 2022, _11_ , 2199 – 2218
  * Juan et al. 2011 Juan, M. L.; Righini, M.; Quidant, R. Plasmon nano-optical tweezers. _Nature Photonics_ 2011, _5_ , 349–356
  * Saleh and Dionne 2012 Saleh, A. A. E.; Dionne, J. A. Toward efficient optical trapping of sub-10-nm particles with coaxial plasmonic apertures. _Nano letters_ 2012, _12 11_ , 5581–6
  * Kotsifaki and Nic Chormaic 2019 Kotsifaki, D. G.; Nic Chormaic, S. Plasmonic optical tweezers based on nanostructures: fundamentals, advances and prospects. _Nanophotonics_ 2019, _8_ , 1227 – 1245
  * Bouloumis et al. 2020 Bouloumis, T.; Kotsifaki, D. G.; Han, X.; Nic Chormaic, S.; Truong, V. G. Fast and efficient nanoparticle trapping using plasmonic connected nanoring apertures. _Nanotechnology_ 2020, _32_ , 025507
  * Juan et al. 2009 Juan, M. L.; Gordon, R.; Pang, Y.; Eftekhari, F.; Quidant, R. Self -induced back-action optical trapping of dielectric nanoparticles. _Nature Physics_ 2009, _5_ , 915–919
  * Neumeier et al. 2015 Neumeier, L.; Quidant, R.; Chang, D. E. Self-induced back-action optical trapping in nanophotonic systems. _New Journal of Physics_ 2015, _17_ , 123008
  * Descharmes et al. 2013 Descharmes, N.; Dharanipathy, U. P.; Diao, Z.; Tonin, M.; Houdré, R. Observation of backaction and self-induced trapping in a planar hollow photonic crystal cavity. _Physical review letters_ 2013, _110 12_ , 123601
  * Mestres et al. 2016 Mestres, P.; Berthelot, J.; Acimovic, S. S.; Quidant, R. Unraveling the optomechanical nature of plasmonic trapping. _Light. Sci. Appl._ 2016, _5_ , e16092
  * Zhang et al. 2019 Zhang, C.; Li, J.; Park, J. G.; feng Su, Y.; Goddard, R. E.; Gelfand, R. M. Optimization of metallic nanoapertures at short-wave infrared wavelengths for self-induced back-action trapping. _Applied optics_ 2019, _58 35_ , 9498–9504
  * Zhu et al. 2018 Zhu, T.; Cao, Y.; Wang, L.; Nie, Z.; Cao, T.; Sun, F.; Jiang, Z.; Nieto-Vesperinas, M.; Liu, Y.; Qiu, C.-W.; Ding, W. Self-Induced Backaction Optical Pulling Force. _Phys. Rev. Lett._ 2018, _120_ , 123901
  * Daniel and Astruc 2004 Daniel, M.-C.; Astruc, D. Gold Nanoparticles: Assembly, Supramolecular Chemistry, Quantum-Size-Related Properties, and Applications toward Biology, Catalysis, and Nanotechnology. _Chemical Reviews_ 2004, _104_ , 293–346
  * Grammatikopoulos et al. 2016 Grammatikopoulos, P.; Steinhauer, S.; Vernieres, J.; Singh, V.; Sowwan, M. Nanoparticle design by gas-phase synthesis. _Advances in Physics: X_ 2016, _1_ , 81–100
  * Halas et al. 2011 Halas, N. J.; Lal, S.; Chang, W.-S.; Link, S.; Nordlander, P. Plasmons in Strongly Coupled Metallic Nanostructures. _Chemical Reviews_ 2011, _111_ , 3913–3961
  * Mie 1908 Mie, G. Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen. _Annalen der Physik_ 1908, _330_ , 377–445
  * Jeanmaire and Van Duyne 1977 Jeanmaire, D. L.; Van Duyne, R. P. Surface raman spectroelectrochemistry: Part I. Heterocyclic, aromatic, and aliphatic amines adsorbed on the anodized silver electrode. _Journal of Electroanalytical Chemistry and Interfacial Electrochemistry_ 1977, _84_ , 1–20
  * Notarianni et al. 2014 Notarianni, M.; Vernon, K.; Chou, A.; Aljada, M.; Liu, J.; Motta, N. Plasmonic effect of gold nanoparticles in organic solar cells. _Solar Energy_ 2014, _106_ , 23–37, Third and Fourth Generation Solar Cells
  * Nguyen et al. 2018 Nguyen, M. A. P.; Kim, S.; Tran, T. T.; Xu, Z.; Kianinia, M.; Toth, M.; Aharonovich, I. Nanoassembly of quantum emitters in hexagonal boron nitride and gold nanospheres. _Nanoscale_ 2018, _10 5_ , 2267–2274
  * Liu and Corma 2018 Liu, L.; Corma, A. Metal Catalysts for Heterogeneous Catalysis: From Single Atoms to Nanoclusters and Nanoparticles. _Chemical Reviews_ 2018, _118_ , 4981–5079
  * Wuenschell et al. 2022 Wuenschell, J.; Jee, Y.; Buric, M.; Chorpening, B. Gold nanoparticle incorporated oxide thin films for gas sensing at high temperature. _MRS Communications_ 2022, _12_ , 308–314
  * Mirigliano et al. 2021 Mirigliano, M.; Paroli, B.; Martini, G.; Fedrizzi, M.; Falqui, A.; Casu, A.; Milani, P. A binary classifier based on a reconfigurable dense network of metallic nanojunctions. _Neuromorphic Computing and Engineering_ 2021, _1_ , 024007
  * Ghosh et al. 2008 Ghosh, P. S.; Han, G.; De, M.; Kim, C.; Rotello, V. M. Gold nanoparticles in delivery applications. _Advanced drug delivery reviews_ 2008, _60 11_ , 1307–15
  * Vines et al. 2019 Vines, J. B.; Yoon, J.-H.; Ryu, N.-E.; Lim, D. J.; Park, H. Gold Nanoparticles for Photothermal Cancer Therapy. _Frontiers in Chemistry_ 2019, _7_
  * Carabineiro 2017 Carabineiro, S. A. Applications of Gold Nanoparticles in Nanomedicine: Recent Advances in Vaccines. 2017
  * Daraee et al. 2016 Daraee, H.; Eatemadi, A.; Abbasi, E.; Aval, S. F.; Kouhi, M.; Akbarzadeh, A. Application of gold nanoparticles in biomedical and drug delivery. _Artificial Cells, Nanomedicine, and Biotechnology_ 2016, _44_ , 410 – 422
  * Nejati et al. 2021 Nejati, K.; Dadashpour, M.; Gharibi, T.; Mellatyar, H.; Akbarzadeh, A. Biomedical Applications of Functionalized Gold Nanoparticles: A Review. _Journal of Cluster Science_ 2021, _33_ , 1–16
  * Min et al. 2013 Min, C.; Shen, Z.; Shen, J.; Zhang, Y.; Fang, H.; Yuan, G.; Du, L.; Zhu, S.; Lei, T.; Yuan, X. Focused plasmonic trapping of metallic particles. _Nature Communications_ 2013, _4_ , 2891
  * Kotsifaki et al. 2020 Kotsifaki, D. G.; Truong, V. G.; Nic Chormaic, S. Fano-Resonant, Asymmetric, Metamaterial-Assisted Tweezers for Single Nanoparticle Trapping. _Nano letters_ 2020, _20_ , 3388–3395
  * Fedotov et al. 2007 Fedotov, V. A.; Rose, M.; Prosvirnin, S. L.; Papasimakis, N.; Zheludev, N. I. Sharp Trapped-Mode Resonances in Planar Metamaterials with a Broken Structural Symmetry. _Physical Review Letters_ 2007, _99_ , 147401
  * Papasimakis and Zheludev 2009 Papasimakis, N.; Zheludev, N. I. Metamaterial-Induced Transparency:Sharp Fano Resonances and Slow Light. _Optics and Photonics News_ 2009, _20_ , 22–27
  * Tanaka et al. 2010 Tanaka, K.; Plum, E.; Ou, J. Y.; Uchino, T.; Zheludev, N. I. Multifold Enhancement of Quantum Dot Luminescence in Plasmonic Metamaterials. _Physical Review Letters_ 2010, _105_ , 227403
  * Luk’yanchuk et al. 2010 Luk’yanchuk, B.; Zheludev, N. I.; Maier, S. A.; Halas, N. J.; Nordlander, P.; Giessen, H.; Chong, T. C. The Fano resonance in plasmonic nanostructures and metamaterials. _Nature Materials_ 2010, _9_ , 707–715
  * Thijssen et al. 2015 Thijssen, R.; Kippenberg, T. J.; Polman, A.; Verhagen, E. Plasmomechanical Resonators Based on Dimer Nanoantennas. _Nano Letters_ 2015, _15_ , 3971–3976
  * Simmons et al. 1996 Simmons, R.; Finer, J.; Chu, S.; Spudich, J. Quantitative measurements of force and displacement using an optical trap. _Biophysical Journal_ 1996, _70_ , 1813–1822
  * Braibanti et al. 2008 Braibanti, M.; Vigolo, D.; Piazza, R. Does Thermophoretic Mobility Depend on Particle Size? _Phys. Rev. Lett._ 2008, _100_ , 108303
  * Würger 2010 Würger, A. Thermal non-equilibrium transport in colloids. _Reports on Progress in Physics_ 2010, _73_ , 126601
  * Rodríguez-Sevilla et al. 2018 Rodríguez-Sevilla, P.; Prorok, K.; Bednarkiewicz, A.; Marqués, M. I.; García-Martín, A.; Solé, J. A. G.; Haro-González, P.; Jaque, D. Optical Forces at the Nanoscale: Size and Electrostatic Effects. _Nano letters_ 2018, _18 1_ , 602–609
  * Griffiths 2017 Griffiths, D. J. _Introduction to Electrodynamics_ , 4th ed.; Cambridge University Press, 2017
  * Ko et al. 2007 Ko, D.-S.; Park, Y. M.; Kim, S.-D.; Kim, Y.-W. Effective removal of Ga residue from focused ion beam using a plasma cleaner. _Ultramicroscopy_ 2007, _107_ , 368–373
  * Fan et al. 2018 Fan, Y.; Zhang, F.; Shen, N.; Fu, Q.; Wei, Z.; qiang Li, H.; Soukoulis, C. M. Achieving a high- Q response in metamaterials by manipulating the toroidal excitations. _Physical Review A_ 2018, _97_ , 033816
  * Jones et al. 2015 Jones, P. H.; Maragò, O. M.; Volpe, G. _Optical Tweezers: Principles and Applications_ ; Cambridge University Press, 2015
  * Rosenblatt et al. 2020 Rosenblatt, G.; Simkhovich, B.; Bartal, G.; Orenstein, M. Nonmodal Plasmonics: Controlling the Forced Optical Response of Nanostructures. _Phys. Rev. X_ 2020, _10_ , 011071


# Symmetry Analysis of Initial and Boundary Value Problems for Fractional
Differential Equations in Caputo sense

Dogan Kaya1, Gulistan Iskenderoglu2

###### Abstract

In this work we study Lie symmetry analysis of initial and boundary value
problems for partial differential equations (PDE) with Caputo fractional
derivative. We give generalized definition and theorem for the symmetry method
for PDE with Caputo fractional derivative, according to Bluman’s definition
and theorem for the symmetry analysis of PDE system. We investigate the
symmetry analysis of initial and boundary value problem for fractional
diffusion and the third order fractional partial differential equation (FPDE).
Also we give some solutions.

##  1 Introduction

Theory of differential equations is one of the most important discipline in
mathematics that has extensive application in physics. The descriptions of
fundamental physical, chemical, and biological systems phenomena and processes
in nature are made by using nonlinear equations. Following to Nobel Laureate
Werner Heisenburg’s words ”… the progress of physics will to a large extent
depend on the progress of nonlinear mathematics, of methods to solve nonlinear
equations … and therefore we can learn by comparing different nonlinear
problems”, - we will concentrate our attention on nonlinear PDE with
fractional derivative [1]. Numerous complex nonlinear phenomena and dynamic
processes in nature can be described or modeled by dint of a FPDE or a system
of FPDE which have been used successfully in physics, electromagnetics [2],
acoustics, astrophysics [3], viscoelasticity, chemistry, electrochemistry,
etc. [4, 5, 6, 7, 8]. And there are numerous methods to find the solutions for
the differential equations like: Adomian decomposition method [9, 10],
differential transform method [11], modified simple equation method [12], Lie
symmetry analysis [13, 14] and so on.

In this work we focus on initial and boundary value problems for FPDE with
fractional derivative in the sense of Caputo [15], that has a form

${}^{C}D^{\alpha}_{t}f(t)=\left\\{\begin{array}[]{l}\frac{1}{\Gamma(n-\alpha)}\int_{a}^{t}(t-\tau)^{n-\alpha-1}f^{(n)}(\tau)d\tau,\quad\quad
if\quad n-1<\alpha<n,n\in\mathbb{N},\\\ \mbox{}\\\
\frac{d^{n}}{dt^{n}}f(t),\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
if\quad\alpha=n.\\\ \end{array}\right.$

The FPDE with Caputo derivative is more useful in searching the solution of
boundary value problems, than Riemann-Liouville derivative, which has a form

$\displaystyle
D^{\alpha}_{t}f(t)=\left\\{\begin{array}[]{l}\frac{1}{\Gamma(n-\alpha)}\frac{d^{n}}{dt^{n}}\int_{a}^{t}(t-\tau)^{n-\alpha-1}f(\tau)d\tau,\quad\quad
if\quad n-1<\alpha<n,n\in\mathbb{N},\\\ \mbox{}\\\
\frac{d^{n}}{dt^{n}}f(t),\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
if\quad\alpha=n.\\\ \end{array}\right.$

We can assume that Caputo derivative is more handy since the initial value for
fractional differential equation with Caputo derivative is the same as the
initial value for integer PDE [16]. For instance, the initial value condition
for fractional differential equation ${}^{C}D_{t}^{\alpha}f(t)=u(t,f)$ with
$\alpha\in(0,1)$, $t<0$ is given as $f(0)=f_{0}$. Along with it for the
fractional differential equation $D_{t}^{\alpha}f(t)=u(t,f)$ with
$\alpha\in(0,1)$, $t<0$, the initial value condition can contain fractional
integral or derivative and it is given as $[D_{t}^{\alpha-1}f(t)]_{t=0}=f_{0}$
if $\alpha\in(1,2)$, then the initial value conditions are given as
$[D_{t}^{\alpha-2}f(t)]_{t=0}=f_{0}$,
$[D_{t}^{\alpha-1}f(t)]_{t=0}=f_{0}^{(1)}$. So, we consider the boundary value
problems with Caputo fractional derivative.

It is well-known that Lie group method is a powerful and direct approach to
construct exact solutions of nonlinear PDE, by analysing the symmetries of the
nonlinear PDE. In addition, based on Lie symmetry method, many other types of
exact solutions of PDE can be obtained, such as traveling wave solutions,
soliton solutions, power series solutions, and so on [17, 18, 19]. Lie
symmetry method was investigated by Sophus Lie, later Ovsyannikov, Olver and
many other big minds studied it by contributing to the theory of symmetries
[20, 21, 22]. Bluman investigated the symmetry analysis of initial and
boundary value problems for PDE, by given the properties of boundary
conditions and boundary surfaces under symmetries (see [23, 24], and
references therein). And we in our work by using Lie symmetry method and its
analysis search symmetries and obtain solutions for initial and boundary value
problems for PDE, containing Caputo fractional derivative.

##  2 Symmetry analysis for fractional partial differential equations

Consider a time-fractional PDE with two independent variables $x\in\mathbb{R}$
and $t>0$ is given as following:

$F(x,t,u,{{}^{C}D}^{\alpha}_{t}u,\partial_{x}u,\partial_{x}^{2}u,...,\partial_{x}^{s}u)=0,\quad
0<\alpha\leq 1$ (1)

where ${{}^{C}D}^{\alpha}_{t}u$ is Caputo fractional derivative of $u$ and
$\partial_{x}^{i}u=\frac{\partial^{i}u}{\partial x^{i}}$, $i=1,...,s$.

One point Lie symmetry transformation acting on a space of two independent
variables $(t,x)$ and depended variable $u$ is determined as

$\begin{array}[]{l}\bar{t}=t+\varepsilon\tau(t,x,u)+O(\varepsilon^{2}),\\\
\mbox{}\\\ \bar{x}=x+\varepsilon\xi(t,x,u)+O(\varepsilon^{2}),\\\ \mbox{}\\\
\bar{u}=u+\varepsilon\eta(t,x,u)+O(\varepsilon^{2}),\\\ \end{array}$ (2)

where $\varepsilon>0$ is an infinitesimal group parameter and $\tau,\xi,\eta$
are the infinitesimals of the transformation. An infinitesimal generator
associated with the above transformation has a following form

$X=\xi(t,x,u)\frac{\partial}{\partial x}+\tau(t,x,u)\frac{\partial}{\partial
t}+\eta(t,x,u)\frac{\partial}{\partial u}$ (3)

where $\xi=\frac{d\bar{x}}{{d\varepsilon}}\mid_{\varepsilon=0}$,
$\tau=\frac{d\bar{t}}{{d\varepsilon}}\mid_{\varepsilon=0}$ and
$\eta=\frac{d\bar{u}}{{d\varepsilon}}\mid_{\varepsilon=0}$.

###### Definition 2.1

We say $u=\theta(t,x)$ is an invariant solution of the equation (1) resulting
from the invariance of the equation (1) under the symmetry (2) with
infinitesimal generator (3) if only if (see [22])

  * •

$u=\theta(t,x)$ is an invariance surface of $X$,

  * •

$u=\theta(t,x)$ solves the equation (1).

In other words, we can say that $u=\theta(t,x)$ is the invariant solution of
the equation (1) under the infinitesimal generator $X$ with Lie transformation
(2) if only if $u=\theta(t,x)$ satisfies:

  * •

$X(u-\theta(t,x))=0$ for $u=\theta(t,x)$, which gives us

$\xi(t,x,\theta(t,x))\frac{\partial\theta(t,x)}{\partial
x}+\tau(t,x,\theta(t,x))\frac{\partial\theta(t,x)}{\partial
t}=\eta(t,x,\theta(t,x)),$

  * •

$F(x,t,\theta(t,x),{{}^{C}D}^{\alpha}_{t}\theta(t,x),\partial_{x}\theta(t,x),\partial_{x}^{2}\theta(t,x),...,\partial_{x}^{s}\theta(t,x))=0$,
for $u=\theta(t,x)$.

According to the infinitesimal transformation (2) a fractional prolongation
$pr^{(\alpha,n)}X$ of the equation (1)

$pr^{(\alpha,s)}X(E)\mid_{E=0}=0,\quad\quad where\quad
E=F(x,t,u,\partial^{\alpha}_{t}u,\partial_{x}u,\partial_{x}^{2}u,...,\partial_{x}^{s}u)=0,$

is an invariant criteria for equation (1) and has the form

$pr^{(\alpha,n)}X=X+\eta^{t}_{\alpha}\partial_{{{}^{C}D}^{\alpha}_{t}u}+\eta^{x}_{1}\partial_{u_{x}}+\eta^{x}_{2}\partial_{\partial_{x}^{2}u}+...+\eta^{x}_{s}\partial_{\partial_{x}^{s}u},$

here

$\begin{array}[]{c}\eta^{t}_{\alpha}={{}^{C}D}^{\alpha}_{t}(\eta)+\xi{{}^{C}D}^{\alpha}_{t}(u_{x})-{{}^{C}D}^{\alpha}_{t}(\xi
u_{x})+\tau^{C}D^{\alpha}_{t}(u_{t})-{{}^{C}D}^{\alpha}_{t}(\tau u_{t}),\\\
\mbox{}\\\ \eta^{x}_{1}=D_{x}\eta-u_{x}D_{x}\xi_{1}-u_{t}D_{x}\tau,\\\
\mbox{}\\\ \eta^{x}_{2}=D_{x}\eta^{x}_{1}-u_{xx}D_{x}\xi-u_{xt}D_{x}\tau,\\\
\mbox{}\vdots\mbox{}\\\
\eta^{x}_{s}=D_{x}\eta^{x}_{s-1}-\partial_{x}^{s}uD_{x}\xi-\partial_{x}^{s-1}u_{t}D_{x}\tau,\\\
\end{array}$ (4)

with total derivative $D_{i}$ in a form

$D_{i}=\partial_{i}+u_{i}\partial_{u}+u_{it}\partial_{u_{t}}+u_{jt}\partial_{u_{t}}+u_{ii}\partial_{u_{i}}+u_{jj}\partial_{u_{j}}+....$

The expressions for $\eta^{x}_{i}$, $i=1,...,s$ in (4) can be easily obtained
(see [22]), here we focus our attention on $\eta^{t}_{\alpha}$ [22, 25] .

We can define the relation between Caputo fractional derivative and
Riemann–Liouville fractional derivative in a form

$D^{\alpha}_{t}f(t)=^{C}D^{\alpha}_{t}f(t)+\sum\limits_{i=0}^{n-1}\frac{f^{(i)}(0)}{\Gamma(i-\alpha+1)}t^{i-\alpha}.$

And the fractional integral can be defined as

$I^{\alpha}_{t}f(t)=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}f(\tau)d\tau.$

So, we can see the below representations of Caputo and Riemann–Liouville
fractional derivatives (see [6, 15])

${{}^{C}D}^{\alpha}_{t}f(t)=I^{n-\alpha}D^{n}f(t)$ (5)
${D}^{\alpha}_{t}f(t)=D^{n}I^{n-\alpha}f(t)$ (6)

According to the generalized Leibnitz rule in [26]

$D^{\alpha}_{t}(f(t)g(t))=\sum\limits^{\infty}_{n=0}\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)I^{n-\alpha}_{t}f(t)D^{n}_{t}g(t),\quad\quad\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)=\frac{(-1)^{n-1}\alpha\Gamma(n-\alpha)}{\Gamma(1-\alpha)\Gamma(n+1)},$

we have

$\xi{{}^{C}D}^{\alpha}_{t}(u_{x})-{{}^{C}D}^{\alpha}_{t}(\xi
u_{x})=-\sum\limits^{\infty}_{n=1}\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)I^{n-\alpha}_{t}(u_{x})D^{n}_{t}(\xi)+\sum\limits_{i=0}^{n-1}\frac{(\xi
u_{x})^{(i)}(0)}{\Gamma(i-\alpha+1)}t^{i-\alpha},$

and

$\tau{{}^{C}D}^{\alpha}_{t}(u_{t})-{{}^{C}D}^{\alpha}_{t}(\tau
u_{t})=-\sum\limits^{\infty}_{n=1}\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)I^{n-\alpha}_{t}(u_{t})D^{n}_{t}(\tau)+\sum\limits_{i=0}^{n-1}\frac{(\tau
u_{t})^{(i)}(0)}{\Gamma(i-\alpha+1)}t^{i-\alpha}.$

Thus, we get the expression

$\begin{array}[]{l}\eta^{t}_{\alpha}=D^{\alpha}_{t}(\eta)+\eta_{0}-\sum\limits^{\infty}_{n=1}\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)I^{n-\alpha}_{t}(u_{x})D^{n}_{t}(\xi)-\sum\limits_{i=0}^{n-1}\frac{(\xi
u_{x}+\tau u_{t})^{(i)}(0)}{\Gamma(i-\alpha+1)}t^{i-\alpha}-\\\
\sum\limits^{\infty}_{n=1}\left(\begin{matrix}\alpha\\\
n+1\end{matrix}\right)I^{n-\alpha}_{t}(u)D^{n+1}_{t}(\tau)-\alpha
D_{t}(\tau)D_{t}^{\alpha}(u).\end{array}$

Here by using a generalized form of the chain rule (see [27])

$\frac{d^{m}f(g(t))}{dt^{m}}=\sum\limits^{m}_{k=0}\sum\limits^{k}_{r=0}\left(\begin{matrix}k\\\
r\end{matrix}\right)\frac{1}{k!}\left(-g(t)\right)^{r}\frac{d^{m}}{dt^{m}}\left(g(t)^{k-r}\right)\frac{d^{k}f(g)}{dg^{k}},$

which is

$\frac{d^{\alpha}f(t,g(t))}{dt^{\alpha}}=\sum\limits^{\infty}_{l=0}\sum\limits^{l}_{m=0}\sum\limits^{m}_{k=0}\sum\limits^{k-1}_{r=0}\left(\begin{matrix}\alpha\\\
l\end{matrix}\right)\left(\begin{matrix}l\\\
m\end{matrix}\right)\left(\begin{matrix}k\\\
r\end{matrix}\right)\frac{1}{k!}\frac{t^{l-\alpha}}{\Gamma(l+1-\alpha)}\left(-g(t)\right)^{r}\frac{d^{m}}{dt^{m}}\left(g(t)^{k-r}\right)\frac{\partial^{l-m+k}f(g)}{\partial
t^{l-m}\partial g^{k}},$

we have

$\eta_{0}=\sum\limits_{j=0}^{n-1}\frac{t^{n-\alpha}(\eta)^{(j)}(0)}{\Gamma(j+1-\alpha)}=\sum\limits_{j=0}^{n-1}\frac{t^{(n-\alpha)}}{\Gamma(j+1-\alpha)}\sum\limits^{\infty}_{l=0}\sum\limits^{l}_{m=0}\sum\limits^{m}_{k=0}\sum\limits^{k-1}_{r=0}\left(\begin{matrix}\alpha\\\
l\end{matrix}\right)\left(\begin{matrix}l\\\
m\end{matrix}\right)\left(\begin{matrix}k\\\
r\end{matrix}\right)\frac{1}{k!}\times$
$\frac{t^{l-\alpha}}{\Gamma(l+1-\alpha)}\left(\left(-u\right)^{r}\frac{d^{j}}{dt^{j}}\left(u^{k-r}\right)\frac{\partial^{l-j+k}\eta}{\partial
t^{l-j}\partial u^{k}}\right)(0).$

And so, the infinitesimal $\eta^{t}_{\alpha}$ takes a form

$\begin{array}[]{l}\eta^{t}_{\alpha}=\frac{{}^{C}\partial^{\alpha}\eta}{\partial
t^{\alpha}}+(\eta_{u}-\alpha(\tau_{t}+u_{t}\tau_{u}))\frac{{}^{C}\partial^{\alpha}u}{\partial
t^{\alpha}}-u\frac{{}^{C}\partial^{\alpha}\eta_{u}}{\partial
t^{\alpha}}+\mu-\sum\limits^{\infty}_{n=1}\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)D^{n}_{t}(\xi)I^{n-\alpha}_{t}(u_{x})+\\\
\sum\limits^{\infty}_{n=1}\left[\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)\frac{\partial^{n}\eta_{u}}{\partial
t^{n}}-\left(\begin{matrix}\alpha\\\
n+1\end{matrix}\right)D^{n+1}_{t}\tau\right]I^{n-\alpha}_{t}u+\sum\limits_{i=0}^{n-1}\frac{(\xi
u_{x})^{(i)}(0)}{\Gamma(i-\alpha+1)}t^{i-\alpha},\\\ \end{array}$ (7)

where

$\mu=\sum\limits^{\infty}_{n=2}\sum\limits^{n}_{m=2}\sum\limits^{m}_{k=2}\sum\limits^{k-1}_{r=0}\left(\begin{matrix}\alpha\\\
n\end{matrix}\right)\left(\begin{matrix}n\\\
m\end{matrix}\right)\left(\begin{matrix}k\\\
r\end{matrix}\right)\frac{1}{k!}\frac{t^{n-\alpha}(-u)^{r}}{\Gamma(n+1-\alpha)}\frac{\partial^{m}}{\partial
t^{m}}\left(u^{k-r}\right)\frac{\partial^{n-m+k}\eta}{\partial t^{n-m}\partial
u^{k}}.$ (8)

Here we can proof the next lemma.

###### Lemma 2.1

If (1) is invariant under infinitesimal transformation (2) with infinitesimal
generator (3) and the equation (1) has no the second and higher order
derivative of $u$ with respect to $t$, then $\eta=A(t,x)u+B(t,x)$ with
$A(t,x)$ and $B(t,x)$ arbitrary functions.

Proof. By expanding the ${{}^{C}D}^{\alpha}_{t}\eta$ we get the expression
(8). As the equation has not any variations of second and high order
derivative $u$ with respect to $t$, then for $k=2$ we have

$\mu=\frac{1}{2!}\frac{t^{2-\alpha}(-u)}{\Gamma(3-\alpha)}u_{tt}\eta_{uu},$

and $\eta_{uu}=0$. Here by integration we get $\eta=A(t,x)u+B(t,x)$.
Q.E.D.

##  3 Symmetry analysis for time-fractional initial and boundary value
problems

In this section we study the one point Lie symmetry analysis for initial and
boundary value problems. We know that the PDE can describe real processes
according nature and society if there are given initial and boundary
conditions for the PDE. Although Lie symmetry analysis is one of the most
widely-applicable methods of finding exact solutions of differential
equations, but it was not widely used for solving boundary value problems. The
reason is the initial and boundary conditions usually are not invariant under
any obtained Lie symmetry method transformations [23, 24]. Thereby, an
invariant solution for PDE resulting by applying symmetry transformation
solves a given boundary value problem, when the symmetry transformation leaves
invariant all boundary conditions and the domain of the boundary value problem
[23].

In [23] Bluman gives a definition of Lie symmetry invariance for initial and
boundary by means of that there are some classes initial and boundary problems
which can be solved. Bluman studies the PDE with finite values of $x$ and $t$
boundary and initial conditions. Later Cherniha extends the definitions of
Bluman for the PDE with infinite values of $x$ and $t$ boundary and initial
conditions [28, 29]. In this work we investigate the symmetry analysis of
initial and boundary value problems for fractional nonlinear PDE with Caputo
time-fractional derivative.

Let Lie symmetry infinitesimal generator $X$

$X=\xi(x,t,u)\frac{\partial}{\partial x}+\tau(x,t,u)\frac{\partial}{\partial
t}+\eta(x,t,u)\frac{\partial}{\partial u}$ (9)

is admitted by the boundary value problem defined on a domain $\Omega$:

$u_{t}=f\left(x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k}u}{\partial
x^{k}}\right),\quad\quad(x,t)\in\Omega\subset\mathbb{R}^{2},$ (10)
$d_{a}(x,t)=0:B_{a}\left(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}}\right)=0,\quad\quad
a=1,...,p.$ (11)

Here $B_{a}(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}})$ boundary condition on
$d_{a}(x,t)$. Suppose that the above boundary value problem has a unique
solution.

###### Definition 3.1

[23] The symmetry $X$ which has the form (9) is allowed by the boundary value
problem (10, 11) if:

  * •

$X^{(k)}(u_{t}-f(x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k}u}{\partial x^{k}}))=0$ for
$u_{t}=f(x,u,\frac{\partial u}{\partial x},...,\frac{\partial^{k}u}{\partial
x^{k}})$;

  * •

$Xd_{a}(x,t)=0$ for $d_{a}(x,t)=0$, $a=1,...,p$;

  * •

$X^{(k-1)}B_{a}(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}})=0$ for
$B_{a}(t,x,u,\frac{\partial u}{\partial x},...,\frac{\partial^{k-1}u}{\partial
x^{k-1}})=0$ on $d_{a}(x,t)=0$, $a=1,..,p$.

Further we extend the Blusman’s definition for FPDE with Caputo derivative.
Let us consider the boundary value problem for FPDE defined on a domain
$\Omega$

$^{C}D^{\alpha}_{t}u=g\left(x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k}u}{\partial
x^{k}}\right),\quad\quad(x,t)\in\Omega\subset\mathbb{R}^{2},$ (12)
${{}^{C}d}_{a}(x,t)=0:{{}^{C}B}_{a}\left(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}}\right)=0,\quad\quad
a=1,...,p.$ (13)

Here ${{}^{C}B}_{a}(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}})$ is a boundary condition on
${{}^{C}d}_{a}(x,t)$. Suppose the boundary value problem (12-13) has a unique
solution. Then we can give next definition.

###### Definition 3.2

The symmetry $X$ in the form (9) is allowed by the boundary value problem (12,
13) if:

  * •

$X^{(k)}({{}^{C}D}^{\alpha}_{t}u-g(x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k}u}{\partial x^{k}}))=0$ for
${{}^{C}D}^{\alpha}_{t}u=g(x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k}u}{\partial x^{k}})$;

  * •

$X{{}^{C}d}_{a}(x,t)=0$ for ${{}^{C}d}_{a}(x,t)=0$, $a=1,...,p$;

  * •

$X^{(k-1)}{{}^{C}B}_{a}(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}})=0$ for
${{}^{C}B}_{a}(t,x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k-1}u}{\partial x^{k-1}})=0$ on ${{}^{C}d}_{a}(x,t)=0$,
$a=1,..,p$.

###### Theorem 3.1

The solution $u=\nu(t,x)$ for (12) is invariant if only if for infinitesimal
generator $X$ the curve (mapping) $\nu(t,x)$ admits

$\begin{array}[]{c}\eta(t,x,\nu(t,x))-\xi(t,x,\nu(t,x))\\\
-\tau(t,x,\nu(t,x)){{}^{C}D}^{1-\alpha}_{t}g\left(x,\nu(t,x),\frac{\partial\nu(t,x)}{\partial
x},...,\frac{\partial^{k}\nu(t,x)}{\partial x^{k}}\right)=0.\end{array}$ (14)

Proof: As the solution surface $u=\nu(t,x)$ is invariant for (12) if only if
$X(u-\nu(t,x))=0$ which gives

$0=X(u-\nu(t,x))=\left(\xi(x,t,u)\frac{\partial}{\partial
x}+\tau(x,t,u)\frac{\partial}{\partial t}+\eta(x,t,u)\frac{\partial}{\partial
u}\right)(u-\nu(t,x))=$ $=\eta(x,t,u)\frac{\partial u}{\partial
u}-\xi(x,t,u)\frac{\partial\nu(t,x)}{\partial
x}-\tau(x,t,u)\frac{\partial\nu(t,x)}{\partial t}=$
$=\eta(x,t,\nu(t,x))-\xi(x,t,\nu(t,x))\frac{\partial\nu(t,x)}{\partial
x}-\tau(x,t,\nu(t,x))\frac{\partial\nu(t,x)}{\partial t}.$

Further, as $u=\nu(t,x)$, then $u_{t}=\frac{\partial\nu(t,x)}{\partial t}$.
From the property of Caputo derivative (5)
$u_{t}={{}^{C}D}_{t}^{1-\alpha}g(x,u,\frac{\partial u}{\partial
x},...,\frac{\partial^{k}u}{\partial x^{k}}))$, i.e

$\eta(x,t,\nu(t,x))-\xi(x,t,\nu(t,x))\frac{\partial\nu(t,x)}{\partial x}-\\\ $
$\tau(x,t,\nu(t,x)){{}^{C}D}_{t}^{1-\alpha}g(x,\nu(t,x),\frac{\partial\nu(t,x)}{\partial
x},...,\frac{\partial^{k}\nu(t,x)}{\partial x^{k}}))=0.$

Q.E.D.

##  4 Symmetry analysis for some time-fractional initial and boundary value
problems

In this section we will study the symmetries of initial and boundary value
problems with fractional derivatives. Suppose the initial and boundary value
problems for fractional diffusion equation

${{}^{C}D}^{\alpha}_{t}u=(u^{p}u_{x})_{x},\quad\quad\quad p\neq 0,$ (15)

with initial and boundary conditions

$\left\\{\begin{array}[]{l}u(t,0)=a(t),\quad\quad for\quad t>0,\\\
u(0,x)=b(x),\quad\quad for\quad x>0.\end{array}\right.$ (16)

Gazizov and others in [17] found the symmetries for above equation in form

$X=(c_{1}+\frac{\alpha}{2}xc_{2}+pxc_{3})\frac{\partial}{\partial
x}+tc_{2}\frac{\partial}{\partial t}+2uc_{3}\frac{\partial}{\partial u}.$

The invariance of $t=0$ leads to $\tau(0)=0$ and the invariance of $x=0$ leads
to $\xi(0,t)=0$ for any $t>0$. Thus, we get $c_{1}=0$ for $\xi(0,t)=0$ with
$t>0$

$tc_{2}\frac{d}{dt}a(t)=2a(t)c_{3},$

and for $\tau(0)=0$ we have

$\left(\frac{\alpha}{2}xc_{2}+pxc_{3}\right)\frac{d}{dx}b(x)=2b(x)c_{3},$

so, we get

$a(t)=t^{\frac{2c_{3}}{c_{2}}}k_{1},\quad\quad k_{1}\quad is\quad a\quad
constant,$ $b(x)=x^{\frac{2}{c_{2}\alpha/2c_{3}+p}}k_{2},\quad\quad k_{2}\quad
is\quad a\quad constant.$

It means that the initial and boundary value problems can be invariant if the
initial and boundary conditions have above forms.

Therefore the equation (15) with initial and boundary values (16) have two
symmetry operators

$X_{1}=\frac{\alpha}{2}x\frac{\partial}{\partial x}+t\frac{\partial}{\partial
t},\quad and\quad X_{2}=px\frac{\partial}{\partial
x}+2u\frac{\partial}{\partial u},$

where $X_{1}$ gives $u=f(x^{\frac{2}{\alpha}}t^{-1})$, and from $X_{2}$
operator we have $u=x^{\frac{2}{p}}g(t)$. As the problem (15) has only one
solution, then

$f(x^{\frac{2}{\alpha}}t^{-1})=x^{\frac{2}{p}}g(t),$

or
$f(x^{\frac{2}{\alpha}}t^{-1})=(x^{\frac{2}{\alpha}}g(t)^{\frac{p}{\alpha}})^{\frac{\alpha}{p}}$.
Here we see that $f(z)=z^{\frac{\alpha}{p}}$ and
$g(t)^{\frac{p}{\alpha}}=t^{-1}$, which gives $g(t)=t^{-\frac{\alpha}{p}}k$.
Therefore $u=x^{\frac{2}{p}}t^{-\frac{\alpha}{p}}c$, where $c$ is arbitrary
constant, that we will be obtained.

Now let’s consider the boundary and initial conditions. Here for $p>0$,
$u(t,0)=0$, or in our problem $k_{1}=0$ and we have blow-up in $t=0$. Namely
for $p>0$, $\lim\limits_{t\rightarrow 0}u(t,x)=\lim\limits_{t\rightarrow
0}x^{\frac{2}{p}}t^{-\frac{\alpha}{p}}c=+\infty$.

And similarly for $p<0$, $u(0,x)=0$, in our problem $k_{2}=0$ with blow-up in
$x=0$, i.e. $\lim\limits_{x\rightarrow 0}u(t,x)=\lim\limits_{x\rightarrow
0}x^{\frac{2}{p}}t^{-\frac{\alpha}{p}}c=+\infty$.

By checking of the solution for the equation (15) we have
$k=\sqrt[p]{\frac{2(2-p)\Gamma(1-\frac{\alpha}{p})}{p^{2}\Gamma(1-\frac{\alpha}{p}-\alpha)}}$,
thus

$u(t,x)=x^{\frac{2}{p}}t^{-\frac{\alpha}{p}}\sqrt[p]{\frac{2(2-p)\Gamma(1-\frac{\alpha}{p})}{p^{2}\Gamma(1-\frac{\alpha}{p}-\alpha)}}$

is a solution for

$\left\\{\begin{array}[]{l}{{}^{C}D}^{\alpha}_{t}u=(u^{p}u_{x})_{x},\quad\quad\quad
p>0,\\\ u(t,0)=0,\quad\quad\quad\quad\quad\quad t>0,\\\
\lim\limits_{t\rightarrow 0}u(t,x)=+\infty,\quad\quad\quad
x>0.\end{array}\right.$

and

$\left\\{\begin{array}[]{l}{{}^{C}D}^{\alpha}_{t}u=(u^{p}u_{x})_{x},\quad\quad\quad
p<0,\\\ \lim\limits_{x\rightarrow 0}u(t,x)=+\infty,\quad\quad\quad t>0\\\
u(0,x)=0,\quad\quad\quad\quad\quad\quad x>0.\end{array}\right.$

Further consider the initial and boundary value problem for below equation

${{}^{C}D}^{\alpha}_{t}u=(u^{q}u_{xx})_{x},\quad\quad\quad q\neq 0,$ (17)

with initial and boundary conditions

$\left\\{\begin{array}[]{l}u(t,0)=d(t),\quad\quad for\quad t>0,\\\
u(0,x)=e(x),\quad\quad for\quad x>0.\end{array}\right.$ (18)

By applying Lie symmetry method for equation (17) we obtain

$X=\left(c_{1}x+c_{2}\right)\frac{\partial}{\partial
x}+\left(c_{3}t+c_{4}\right)\frac{\partial}{\partial
t}+\left(\frac{3c_{1}-\alpha c_{3}}{q}\right)u\frac{\partial}{\partial u},$
(19)

Thus the invariance of $t=0$ leads to $c_{4}=0$ and the invariance of $x=0$
leads to $c_{2}=0$. And

$c_{3}t\frac{d}{dt}d(t)=\frac{3c_{1}-\alpha c_{3}}{q}d(t),$

which gives $d(t)=t^{\frac{3c_{1}-\alpha c_{3}}{c_{3}q}}k_{3}$, here $k_{3}$
is an arbitrary constant.

So, for $t=0$ we get

$c_{1}x\frac{d}{dx}e(x)=\frac{3c_{1}-\alpha c_{3}}{q}e(x),$

or $e(x)=x^{\frac{3c_{1}-\alpha c_{3}}{c_{1}q}}k_{4}$, here $k_{4}$ is an
arbitrary constant.

That means the initial and boundary problem for equation (17) is invariant
according to infinitesimal generator (19) if the boundary and initial
conditions have above form and there are two infinitesimal operators

$Y_{1}=x\frac{\partial}{\partial x}+\frac{3}{q}u\frac{\partial}{\partial
u},\quad\quad Y_{2}=t\frac{\partial}{\partial
t}-\frac{\alpha}{q}u\frac{\partial}{\partial u}.$

The first operator gives $u=t^{-\frac{\alpha}{q}}h(x)$, and the second
operator gives $u=x^{\frac{3}{q}}l(t)$. Thus, From the uniqueness of the
solution for initial and boundary problem, we have
$u=t^{-\frac{\alpha}{q}}x^{\frac{3}{q}}c^{\prime}$, where $c^{\prime}$ is an
arbitrary constant.

Therefore for $q>0$, $u(t,0)=0$, i.e. $k_{3}=0$ and we have blow-up in $t=0$
as $\lim\limits_{t\rightarrow 0}u(t,x)=\lim\limits_{t\rightarrow
0}t^{-\frac{\alpha}{q}}x^{\frac{3}{q}}c^{\prime}=+\infty$.

In a like manner for $q<0$, $u(0,x)=0$, with $k_{4}=0$ and blow-up in $x=0$,
with $\lim\limits_{x\rightarrow 0}u(t,x)=\lim\limits_{x\rightarrow
0}t^{-\frac{\alpha}{q}}x^{\frac{3}{q}}c^{\prime}=+\infty$.

By checking the solution for equation (17) we obtain

$u(t,x)=t^{-\frac{\alpha}{q}}x^{\frac{3}{q}}\sqrt[q]{\frac{q^{3}\Gamma(1-\frac{\alpha}{q})}{3(3-q)(3+q)\Gamma(1-\frac{\alpha}{q}-\alpha)}},$

the solution for initial and boundary value problem

$\left\\{\begin{array}[]{l}{{}^{C}D}^{\alpha}_{t}u=(u^{q}u_{xx})_{x},\quad\quad\quad
q>0,\\\ u(t,0)=0,\quad\quad\quad\quad\quad\quad t>0,\\\
\lim\limits_{t\rightarrow 0}u(t,x)=+\infty,\quad\quad\quad
x>0.\end{array}\right.$

and

$\left\\{\begin{array}[]{l}{{}^{C}D}^{\alpha}_{t}u=(u^{q}u_{xx})_{x},\quad\quad\quad
q<0,\\\ \lim\limits_{x\rightarrow 0}u(t,x)=+\infty,\quad\quad\quad t>0,\\\
u(0,x)=0,\quad\quad\quad\quad\quad\quad x>0.\end{array}\right.$

##  5 Conclusion

In this work, we present the applications of Lie group analysis to study the
initial and boundary value problems for time-fractional nonlinear PDE given by
Caputo sense. We give the definition of invariance of the initial and boundary
value problems for time-fractional PDE and proved some theorems. By using Lie
symmetry method and its analysis for the initial and boundary value problems
with time-fractional PDE we obtained some exact solutions.

##  6 Acknowledgment

The authors would like to thank Istanbul Commerce University, 22-2018/34 and
Turkish scholarships for supporting this work.

## References

  * [1] Wilhelms, W.: Linear Second-Order PDEs, Institute of Theoretical Physics, Muenster University, Germany, 2013.
  * [2] Zelenyi, L.M., Milovanov, A.V.: Fractal topology and strange kinetics: from percolation theory to problems in cosmic electrodynamics, Physics Uspekhi 47, 749–788 (2004).
  * [3] Tarasov, V.E.: Gravitational field of fractal distribution of particles, Celes. Mech. Dynam. Astron. 19, 1–15 (2006).
  * [4] He, J.H.: Some applications of nonlinear fractional differential equations and their approximations, Bull. Sci. Technol 15, 86–90 (1999).
  * [5] Diethelm, K.: The Analysis of Fractional Differential Equations, Springer, Berlin, 2010.
  * [6] Miller, K.S., Ross, B.: An Introduction to the Fractional Calculus and Fractional Differential Equations Wiley–Interscience, New York, 1993.
  * [7] Wu, G., Lee, E.W.M.: Fractional variational iteration method and its application, Phys. Lett. A374, 2506–2509 (2010).
  * [8] Barley, R., Torvik, P.: A theoretical basis for the application of fractional calculus to viscoelasticity, J. Rheol. 27, 201–210 (1983).
  * [9] Adomian, G.: A review of the decomposition method in applied mathematics, J. Math. Anal. Appl. 135, 44–501 (1988).
  * [10] Kaya, D.: An application for the higher order modified KdV equation by decomposition method, Commun. Nonlinear Sci. Numer. Simul. 10(6), 693–702 (2005).
  * [11] Hereman, W., Banerjee, P.P., Chatterjee, M.R.: On the nonlocal equations and nonlocal charges associated with the Harry-Dym hierarchy Korteweg-de-Vries equation, J. Phys. A 22, 241–252 (1989).
  * [12] Taghizadeh, N., Mirzazadeh, M., Rahimian, M., Akbari, M.: Application of the simplest equation method to some time-fractional partial differential equations, Ain Shams Eng. J. 4(4), 897–902 (2013).
  * [13] Hu, J., Ye, Y., Shen, S., Zhang, J.: Lie symmetry analysis of the time fractional dV-type equation, Appl. Math. and Comp. 223, 439–444 (2014).
  * [14] Sahadevan, R., Bakkyaraj, T.: Invariant analysis of time fractional generalized Burgers and Kortewegde Vries equations, J. Math. Anal. Appl. 393, 341–347 (2012).
  * [15] Podlubny, I.: Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, Some Methods of their Solution and Some of their Applications, Academic Press, San Diego, 1999.
  * [16] Li, C., Qian, D., Chen, Y.Q.: On Riemann–Liouville and Caputo derivatives, Discrete Dynamics in Nature and Society, 15 (2011), Article ID 562494.
  * [17] Gazizov, R.K., Kasatkin, A.A., Lukashchuk, S.Y.: Continuous transformation groups of fractional differential equations, Vestn. USATU, 9, 125–135 (2007).
  * [18] Khalique, C.M., Adem, K.R.: Exact solutions of the $(2+1)$-dimensional Zakharov Kuznetsov modified equal width equation using Lie group analysis, Math. and Comp. Modelling 54, 184–189 (2011).
  * [19] Iskandarova, G., Kaya, D.: Symmetry solution on fractional equation, IJOCTA, 7(3), 255–259 (2017).
  * [20] L.V. Ovsiannikov.: Group analysis of differential equations, Acad. Press, USA, 1982.
  * [21] Ibragimov, N.: Lie group analysis classical heritage, ALGA Publications Blekinge Institute of Technology Karlskrona, Sweden, 2004.
  * [22] Olver, P.: Applications of Lie groups to differential equations Springer Science, Germany, 2012.
  * [23] Bluman, G. W. and Kumei, S.: Symmetries and Differential Equations, Berlin etc., Springer-Verlag, 1989.
  * [24] Bluman, G. W. and Anco, S.C., Symmetry and Integration Methods for Differential Equations, 154 Applied Mathematical Sciences, Springer-Verlag, New York, 2002.
  * [25] Ibragimov, N.: CRC Handbook of Lie Group Analysis of Differential Equations, 1 CRC Press, Boca Raton, 1994.
  * [26] Samko, S., Kilbas, A., Marichev, O.: Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach Science Publishers, USA, 1993.
  * [27] Oldham, K.B. and Spanier, J.: The fractional calculus Academic Press, New York, 1974.
  * [28] Cherniha, R. and Kovalenko, S.: Lie symmetry of a class of nonlinear boundary value problems with free boundaries, Banach Center Publ. 93, 73–82 (2011).
  * [29] Cherniha, R. and Kovalenko, S.: Lie symmetries of nonlinear boundary value problems, Commun. Nonlinear Sci. Numer. Simulat., 17, 71–84 (2012).