# Quantum Algorithm for Estimating Largest Eigenvalues

Nhat A. Nghiem  Department of Physics and Astronomy, State University of New
York at Stony Brook, Stony Brook, NY 11794-3800, USA  Tzu-Chieh Wei
Department of Physics and Astronomy, State University of New York at Stony
Brook, Stony Brook, NY 11794-3800, USA  C. N. Yang Institute for Theoretical
Physics, State University of New York at Stony Brook, Stony Brook, NY
11794-3840, USA

###### Abstract

A majority of numerical scientific computation relies heavily on handling and
manipulating matrices, such as solving linear equations, finding eigenvalues
and eigenvectors, and so on. Many quantum algorithms have been developed to
advance these computational tasks, and in some cases, such as solving linear
equations, can be shown to yield exponential speedup. Here, employing the
techniques in the HHL algorithm and the ideas of the classical power method,
we provide a simple quantum algorithm for estimating the largest eigenvalue in
magnitude of a given Hermitian matrix. As in the case of the HHL algorithm,
our quantum procedure can also yield exponential speedup compared to classical
algorithms that solve the same problem. We also discuss a few possible
extensions and applications of our quantum algorithm, such as a version of a
hybrid quantum-classical Lanczos algorithm.

##  I Introduction

Quantum computation and information has undergone a long way since its
earliest conceptions feynman2018simulating ; deutsch1985quantum . Many quantum
algorithms have been proposed to tackle computational problems, ranging from
searching grover1996fast , factoring shor1999polynomial , simulating quantum
systems childs2010quantum ; feynman2018simulating ; berry2007efficient ,
solving linear systems harrow2009quantum ; childs2017quantum , machine
learning and data science lloyd2013quantum ; lloyd2014quantum ;
lloyd2016quantum ; lloyd2020quantum ; schuld2018supervised ;
mitarai2018quantum , and to problems in computational topology and geometry
lloyd2013quantum ; nghiem2022constant . Many quantum algorithms and
applications remain to be discovered. One of the most versatile results in the
above is the quantum algorithm–the so-called HHL algorithm–for solving linear
systems, proposed Harrow, Hassidim and Lloyd harrow2009quantum  (as well as
its improvement clader2013preconditioned ; childs2017quantum ;
wossnig2018quantum ), as linear systems play a central role in many
engineering and science domains. Aside from its exponential speedup, another
important result from Ref. harrow2009quantum  is the BQP-completeness of the
matrix inversion, which is a more general statement of the complexity
regarding solving linear systems. The HHL algorithm has also been improved in
Ref. childs2017quantum , which yields an exponential improvement in the error
dependence.

We recall that in the original HHL work harrow2009quantum , the aim is to
solve the following linear equation

$\displaystyle Ax=b,$ (1)

where $A$, for simplicity, is a Hermitian matrix of size $N\times N$, and
without loss of generality, $b$ is given and assumed to be normalized. It is
easy to see that if the unique solution exists, it should be

$\displaystyle x=A^{-1}b.$ (2)

The main tool of Ref. harrow2009quantum  is the ability to exponentiate a
sparse matrix berry2007efficient , i.e, to perform $\exp(-iAt)$, where $A$ is
the sparse matrix given at hand. As emphasized by HHL harrow2009quantum ,
their algorithm essentially applies $A^{-1}$ to some given initial state
$\ket{b}$, via quantum phase estimation of $A$’s eigenvalues and controlled
rotation to achieve the inversion by multiplication of the inverse eigenvalues
to the coefficients in the decomposition of $b$ into the eigenbasis, thereby
resulting in the (normalized) solution $\ket{x}$. The running time of their
algorithm harrow2009quantum  is

$\displaystyle\mathcal{O}\left(\log(N)\frac{s^{4}\kappa^{2}}{\epsilon}\right),$
(3)

where $s$ is the sparsity of an $N\times N$ matrix $A$, $\kappa$ is the
conditional number of $A$, and $\epsilon$ is the error.

In this article, motivated by the results in Refs. harrow2009quantum ;
childs2017quantum , we consider the following problem:

Given a sparse Hermitian matrix $A$ of size $n\times n$ with bounded
eigenvalues, estimate its largest absolute eigenvalue.

For simplicity, we assume that the eigenvalues of $A$ satisfy the following

$\displaystyle 0<|\lambda_{1}|\leq|\lambda_{2}|\leq\cdots\leq|\lambda_{n}|<D,$
(4)

where $D$ is some constant (note that since $A$ is Hermitian, all of
eigenvalues are real).

Below, we construct a quantum algorithm that, given access to entries of a
sparse Hermitian matrix $A$ in a similar manner to Ref. harrow2009quantum ,
can estimate its largest eigenvalue in magnitude with a running time that is
logarithmic with respect to the size of $A$. The main tool of our algorithm is
built upon the quantum linear system harrow2009quantum ; wiebe2012quantum  and
the Hadamard test nielsen2002quantum . As with the HHL algorithm, therefore,
our algorithm provides an exponential speedup against, e.g., the classical
power method.

The remaining structure of this paper is as follows. First, in Sec. II, we
review the classical algorithm by the power method, which can estimate the
largest absolute eigenvalue of some matrix $A$, which in general, does not
have to be sparse or well-conditioned. In Sec. III.1, we present a
corresponding quantum algorithm that solves the largest-eigenvalue (in
magnitude) problem mentioned above, and it is a direct translation of the
classical power method into the quantum setting. In Sec. IV, we then discuss a
few possible extensions of our quantum algorithm, including the non-sparse
case and how we can modify our algorithm for the minimum eigenvalue, as well
as extension to a few largest eigenvalues in magnitude. In the Appendix A, we
provide an analysis of the power method, showing that the number of iteration
is efficiently low in order to reach the convergence of power method. A more
careful elaboration of the running time of our quantum algorithm with respect
to the number of iterations is provided in Appendix B. In Appendix C, we
explain how the original HHL approach is adapted to our case, and in Appendix
D, we show a simplified version of the multiplication application of the
matrix, reducing the time complexity.

##  II Classical Algorithm

The classical approach to tackle the above problem is called the power method
chu1993multivariate , which is a very standard approach parlett1982estimating
; o1979estimating , and literally can be found in any textbook on numerical
methods. The main procedure is as follows:

$\bullet$ First initialize some random vector $x_{0}\in C^{n}$.

$\bullet$ Iterate the matrix application process $k$ times and form the
following sequence:

$\displaystyle x_{1}=Ax_{0},$ $\displaystyle x_{2}=Ax_{1}=A^{2}x_{0},$
$\displaystyle x_{3}=Ax_{2}=A^{3}x_{0},$ $\displaystyle\vdots$ $\displaystyle
x_{k}=A^{k}x_{0}.$

$\bullet$ Compute the approximated largest eigenvalue $\bar{\lambda}_{n}$:

$\displaystyle\bar{\lambda}_{n}\equiv\frac{Ax_{k}^{*}\cdot
x_{k}}{x_{k}^{*}\cdot x_{k}},$

where $(\cdot)$ denotes the inner product, and $*$ denotes the complex
conjugate. We remind that, since $A$ is Hermitian, its eigenvalues are always
real.

The following lemma shows that $x_{k}$ can approximate the eigenvector $E_{n}$
with eigenvalue $\lambda_{n}$.

###### Lemma 1

Let $x_{k}$ be defined as above. Assuming the largest absolute eigenvalue
$\lambda_{n}$ (with corresponding eigenvector $E_{n}$) is unique, i.e,
$|\lambda_{i}|<|\lambda_{n}|$ for all $1<i<n$. Then

$\displaystyle x_{k}\approx\alpha\cdot E_{n},$

where $\alpha$ is a some number that depends on $k$–the number of iterations.

Proof: Since $A$ is Hermitian, its eigenvectors $\\{E_{n}\\}_{1}^{N}$ can be
made orthonormal if degenerate. We can choose them to be the basis of $A$.
Therefore, the initial random vector $x_{0}$ can be expressed as:

$\displaystyle x_{0}=c_{1}E_{1}+c_{2}E_{2}+\cdots+c_{n}E_{n}.$ (5)

After the first iteration, we have

$\displaystyle x_{1}$ $\displaystyle=Ax_{0}$ (6)
$\displaystyle=A\big{(}c_{1}E_{1}+c_{2}E_{2}+\cdots+c_{n}E_{n}\big{)}$ (7)
$\displaystyle=c_{1}\lambda_{1}E_{1}+c_{2}\lambda_{2}E_{2}+\cdots+c_{n}\lambda_{n}E_{n}.$
(8)

Proceeding similarly, after $k$ iterations, we have

$\displaystyle
x_{k}=\lambda_{n}^{k}\Big{(}c_{1}\frac{\lambda_{1}^{k}}{\lambda_{n}^{k}}E_{1}+c_{2}\frac{\lambda_{2}^{k}}{\lambda_{n}^{k}}E_{2}+\cdots+c_{n}E_{n}\Big{)}.$
(9)

We have assumed that the largest eigenvalue in magnitude is non-degenerate,

$0\leq|\lambda_{1}|\leq|\lambda_{2}|\leq\cdots<|\lambda_{n}|.$

Therefore,

$\displaystyle\frac{\lambda_{i}^{k}}{\lambda_{n}^{k}}\rightarrow 0,\
\mbox{as}\,k\rightarrow\infty.$ (10)

Hence, when $k$ is sufficiently large, we have:

$x_{k}\approx\lambda_{n}^{k}c_{n}E_{n}\sim E_{n}.$

The proof is then completed. $\blacksquare$

From here, we can see that the factor $\alpha$ in Lemma. 1 is

$\alpha=\lambda_{n}^{k}c_{n}.$

Now we try to compute:

$\frac{Ax_{k}^{*}\cdot x_{k}}{x_{k}^{*}\cdot x_{k}}.$

Since $x_{k}\approx\alpha E_{n}$, then $Ax_{k}\approx\lambda_{n}x_{k}$. The
inner product is then:

$Ax^{*}_{k}\cdot x_{k}\approx\lambda_{n}(x^{*}_{k}\cdot x_{k}),$

which directly leads to

$\frac{Ax^{*}_{k}\cdot x_{k}}{x^{*}_{k}\cdot x_{k}}\approx\lambda_{n}.$

In Appendix A, we provide a careful analysis of the above approximation. It
turns out that the number of iteration $k$ grows logarithmically with the
inverse of the error, thus showing the efficiency of the power method.

The power method is improved by the Lanczos algorithm chu1993multivariate ,
which essentially builds up the Krylov subspace spanned by
$[x_{0},Ax_{0},A^{2}x_{0},\dots,A^{k}x_{0}]$. By reducing $A$ to this
subspace, one can solve the largest few eigenvalues. We note that in this way
the non-degeneracy constraint above (in Lemma 1) can be lifted.

##  III Quantum Algorithm

In this section, we first outline a general procedure to obtain the desired
quantum states that corresponds to the k-times application of matrix A. We
then show how to process these states further with application of CNOT, plus
Hadamard test to estimate the maximum eigenvalue
$\lambda_{max}\equiv\lambda_{n}$. We then provide the running time analysis
with a particular remark on how to improve the running time further by a
factor of k.

###  III.1 Main Procedure

In this section we provide the quantum algorithm for estimating $\lambda_{n}$,
based on the consecutive application of matrix multiplication. The strategy is
similar to that used in Ref. wiebe2012quantum , which is based upon the result
of Ref. harrow2009quantum . We first recall the matrix multiplication
algorithm.

Suppose we have a normalized vector $x_{0}$ and we want to perform $Ax_{0}$.
Assume that we can prepare the corresponding (normalized) state $\ket{x_{0}}$
with some unitary $U_{0}$ circuit. We note that $x_{0}$ can be arbitrary as
long as it is not a zero vector. Therefore, $U_{0}$ can be some random
circuit. The encoding can be nontrivial. Following Ref. wiebe2012quantum , we
obtain a unitary circuit $U_{A}$ that acts as follows:

$\displaystyle
U_{A}U_{0}\ket{00...0}=\ket{\Phi_{1}}=C|Ax_{0}|\ket{Ax_{0}}\ket{0}+\ket{G_{1}}\ket{1},$
(11)

where $\ket{Ax_{0}}$ is the normalized state corresponding to the vector
$Ax_{0}$, $\ket{G_{1}}$ is some state (not normalized) entangled to $\ket{1}$;
$C$ is a rotational factor chosen to guarantee proper normalization so that
$C|Ax_{0}|<1$. Next, we would like to perform $A^{2}x_{0}=A(Ax_{0})$. The
corresponding unitary $U_{A^{2}}$ is simply $U_{A}^{2}U_{0}$ plus the action
on a second ancilla in the controlled rotation, and this results in:

$\displaystyle\ket{\Phi_{2}}$
$\displaystyle=U_{A}\ket{\Phi_{1}}\ket{0}=U_{A^{2}}\ket{00...0}\ket{0}$ (12)
$\displaystyle=C^{2}|A^{2}x_{0}|\ket{A^{2}x_{0}}\ket{00}+\ket{G_{1}}\ket{10}+\ket{G_{2}}\ket{01},$

where, $|A^{2}x_{0}\rangle$ is the normalized state corresponding to the
vector $A^{2}x_{0}$, $\ket{G_{2}}$ is another irrelevant state (not
normalized). We note that the same constant $C$ is used again, as it is chosen
to be larger than the largest eigenvalue in magnitude. There is an important
detail regarding the above state: recall that the HHL algorithm
harrow2009quantum  makes use of quantum phase estimation to extract the
eigenvalues of matrix $A$, and then rotates an ancilla initialized in
$\ket{0}$ controlled by the phase register, followed by the uncomputation of
the phase register. We do not explicitly include the phase register here, but
we provide further explanation in Appendix C and a further simplification in
Appendix D. In the second application of $U_{A}$ (Eqn. (27)), the rotation
step needs to be controlled further by the state of the first ancilla with
$\ket{0}$ (which is entangled to $\ket{Ax_{0}}$). Repeating the above
procedure $k$ times, we obtain the following: (noting that
$A|Ax_{0}||Ax_{0}\rangle=|A^{2}x_{0}||A^{2}x_{0}\rangle$; see also Appendix C)

$\displaystyle\ket{\Phi_{k}}=$ $\displaystyle U_{A^{k}}\ket{00..0}\ket{0}$
(13) $\displaystyle=$ $\displaystyle
C^{k}|A^{k}x_{0}|\ket{A^{k}x_{0}}\ket{0}^{\otimes
k}+\sum_{i=1}^{k}\ket{G_{i}}\ket{\\{1\\}^{i}},$ (14)

where $\ket{\\{1\\}^{i}}$ is the $k$-bit string that has 0 everywhere but 1 at
$i$-th position.
We remind that our goal is to compute the approximated eigenvalue
$\bar{\lambda}_{n}$:

$\displaystyle\bar{\lambda}_{n}=\frac{Ax_{k}^{*}\cdot x_{k}}{x_{k}^{*}\cdot
x_{k}}=\frac{(A^{k+1}x_{0})^{*}\cdot(A^{k}x_{0})}{(A^{k}x_{0})^{*}\cdot(A^{k}x_{0})}.$
(15)

Now we provide a procedure to estimate $\bar{\lambda}_{n}$.

$\bullet$ Apply another round of matrix multiplication to
$\ket{\Phi_{k}}\ket{0}$ to obtain

$\displaystyle\ket{\Phi_{k+1}}=C^{k+1}|A^{k+1}x_{0}|\ket{A^{k+1}x_{0}}\ket{0}^{\otimes
k}\ket{0}+\sum_{i=1}^{k+1}\ket{G_{i}}\ket{1^{i}}.$ (16)

$\bullet$ Add an ancilla $\ket{00}$, then we have the state
$\ket{\Phi_{k+1}}\ket{00}$. Apply $c^{k+1}-X$ gates, which is the controlled
gate that flips the second bit of the two extra ancillas conditioned on one of
the $(k+1)$ controlled bits being 1 (or equivalently by $(k+1)$ two-qubit
controlled-NOT gates with $(k+1)$ respective controls),

$\displaystyle\ket{\Phi_{k+1}^{{}^{\prime}}}=$ $\displaystyle
C^{k+1}|A^{k+1}x_{0}|\ket{A^{k+1}x_{0}}\ket{0}^{\otimes k}\ket{0}\ket{00}$
$\displaystyle+\sum_{i=1}^{k+1}\ket{G_{i}}\ket{1^{i}}\ket{01}.$ (17)

We denote the whole process leading to $\ket{\Phi^{{}^{\prime}}_{k+1}}$ as
${\cal U}_{\Phi^{{}^{\prime}}_{k+1}}$.

$\bullet$ At an earlier step we obtain $\ket{\Phi_{k}}\ket{0}$. If we add two
ancillas initialized in $\ket{00}$ to $\ket{\Phi_{k}}\ket{0}$, then apply
$c^{k+1}-X$ gates as the previous step, but this time the gate flips the first
bit of the two extra ancillas, we arrive at

$\displaystyle\ket{\Phi_{k}^{{}^{\prime}}}=C^{k}|A^{k}x_{0}|\ket{A^{k}x_{0}}\ket{0}^{\otimes
k}\ket{0}\ket{00}+\sum_{i=1}^{k}\ket{G_{i}}\ket{1^{i}}\ket{0}\ket{10}.$ (18)

We denote the whole process leading to $\ket{\Phi^{{}^{\prime}}_{k}}$ as
${\cal U}_{\Phi^{{}^{\prime}}_{k}}$.

$\bullet$ We note that the overlap between $\ket{\Phi^{{}^{\prime}}_{k+1}}$
and $\ket{\Phi^{{}^{\prime}}_{k}}$ is proportional to the inner product of the
two vectors $x_{k+1}=Ax_{k}=A^{k+1}x_{0}$ and $x_{k}=A^{k}x_{0}$,

$\displaystyle\braket{\Phi^{{}^{\prime}}_{k+1},\Phi^{{}^{\prime}}_{k}}=C^{k}C^{k+1}(x^{*}_{k+1}\cdot
x_{k}).$ (19)

Therefore, we can use the SWAP test nielsen2002quantum  to estimate the above
overlap. To do this, we construct a controlled circuit using with the above
${\cal U}_{k+1}$ and ${\cal U}_{k}$ as well as one additional ancilla qubit,
$|0\rangle\langle 0|\otimes{\cal U}_{k+1}+|1\rangle\langle 1|\otimes{\cal
U}_{k}$, so that it creates a superposition of
$(|0\rangle\otimes\ket{\Phi_{k+1}^{{}^{\prime}}}+|1\rangle\otimes\ket{\Phi_{k}^{{}^{\prime}}})/\sqrt{2}$
from an initial state $|+\rangle\otimes|0\dots 0\rangle$. Measuring the
ancilla in the $|\\!+\\!/\\!-\rangle$ basis (Pauli X) gives the real part of
the overlap. In addition, measuring ancilla in the $|\\!+\\!i/\\!-\\!i\rangle$
(Pauli Y) basis gives the imaginary part of the overlap. In order to estimate
such an overlap up to an additive error $\Delta$, it takes
$\mathcal{O}(1/\Delta^{2})$ repetitions.

$\bullet$ Now we repeat the second step (the step that yields Eqn. 17), but
use the state $\ket{\Phi_{k}}\ket{0}$ instead. We then obtain

$\displaystyle\ket{\Phi_{k}^{{}^{\prime\prime}}}=C^{k}|A^{k}x_{0}|\ket{A^{k}x_{0}}\ket{0}^{\otimes
k}\ket{0}\ket{00}+\sum_{i=1}^{k}\ket{G_{i}}\ket{1^{i}}\ket{0}\ket{01}.$ (20)

The overlap

$\displaystyle\braket{\Phi^{{}^{\prime}}_{k},\Phi_{k}^{{}^{\prime\prime}}}=C^{2k}(x^{*}_{k}\cdot
x_{k})$ (21)

can be estimated by the Hadamard test similar to the previous step.

$\bullet$ Now we compute the ratio (by a classical computer):

$\displaystyle\frac{\braket{\Phi_{k+1}^{{}^{\prime}},\Phi_{k}^{{}^{\prime}}}}{\braket{\Phi_{k}^{{}^{\prime}},\Phi_{k}^{{}^{\prime\prime}}}}=C\frac{x^{*}_{k+1}\cdot
x_{k}}{x_{k}^{*}\cdot x_{k}}=C\bar{\lambda}_{n}.$ (22)

Since $C$ is known, then we can estimate $\bar{\lambda}_{n}$ by dividing the
above ratio by $C$.

###  III.2 Time Complexity

The time required to execute the above quantum algorithm involves the time it
takes to apply $A$ (i.e, running the HHL-like procedure without measuring the
ancilla associated with the controlled rotation) supposedly $k$ times, plus
the time it takes to repeat the Hadamard test algorithm to compute the desired
overlaps. Therefore, at first, the total time is supposedly

$\mathcal{O}\Big{(}\frac{\log(n)s^{4}\kappa k^{2}}{\epsilon\Delta^{2}}\Big{)}$
(23)

where $n$ is the dimension of the matrix $A$, $s$ is the sparsity of $A$,
$\kappa$ is the conditional number of $A$, $\epsilon$ is the tolerance error
in the matrix application, and $\Delta$ is the additive error of the
estimation of $\bar{\lambda}_{n}$. We will now go deeper into some details of
running time. We will see that it can be optimized further.

Remark:

$\bullet$ Note that we repeat multiplication of $A$ by $k$ times to extract
the eigenvalue, which means that the running time grows linearly with $k$. The
reason why we have $k^{2}$ in the above expression instead of $k$ is because
we need to take into account of error accumulation. In each of the iteration,
we aim to apply $A$ to the resulting state of the previous application of $A$.
However, each iteration induces an error $\epsilon$. Then after $k$
iterations, the error accumulates as $\mathcal{O}(k\epsilon)$. In order to
obtain an error of $\mathcal{O}(\epsilon)$ at the end, we require an error
$\epsilon/k$ in each individual iteration. Fortunately, as we have provided
the analysis in appendix A, and found that if

$k\sim\log(1/\delta)$

then $\bar{\lambda}_{n}$ is a $\delta$-multiplicative approximation of the
real eigenvalue $\lambda_{n}$. Therefore, the running time is still very
efficient with respect to the error tolerance.

$\bullet$ The running time provided in Eqn. 23 can in fact be reduced by a
factor of $k$, due to the fact that we do not need to repeat the most time-
consuming step in the algorithm, which is the simulation of $\exp(-iAt)$, as
well the phase estimation. As explicitly elaborated in Appendix D, for the
purpose of extracting eigenvalues of $A$, it suffices to conditionally rotate
multiple ancilla systems (after the phase estimation for eigenvalues), and
hence we arrive at our desired state Eqn. 14, up to the uncomputation of the
phase estimation. Therefore, the total running time of our algorithm is
reduced to

$\displaystyle\mathcal{O}\Big{(}\frac{\log(n)s^{4}\kappa
k}{\epsilon\Delta^{2}}\Big{)}.$ (24)

The classical running time of the power method is

$\mathcal{O}({\rm poly}(n)\cdot k),$

as matrix multiplication classically proceeds by direct computation, hence
requires time at least as much as the square of the dimension $n$ of the
matrix. Eventually, the inner product is evaluated, and hence it requires
further additional linear time. Therefore, with respect to the size $n$ of the
matrix, the quantum algorithm yields an exponential speedup if $s$ and
$\kappa$ grow as fast as only a polynomial of logarithm of $n$.

As a final remark, in addition to the maximum eigenvalue, our procedure also
yields a corresponding approximate eigenvector, which can be used for further
processing, such as computing expectation values of observables.

##  IV Some Extensions and Relevant Works

The problem that we solve in this paper is part of the matrix and eigenvalue
problem. instead. Despite seeming simple, such a problem has far reaching
consequences, in both theory and application. In pure and applied mathematics,
one of the very difficult problems concerns the distribution of eigenvalues of
random matrix tao2011random . In physics, a physical system is characterized
by the Hamiltonian, which is represented by a Hermitian matrix $H$. The ground
state is the particular state that corresponds to the eigenstate of $H$ which
has the smallest eigenvalue. We have assumed the matrix $A$ is Hermitian, but
the extension to non-Hermitian matrix is straightforward as in Ref.
harrow2009quantum  by placing $A$ and $A^{\dagger}$ in the two off-diagonal
blocks,

$\tilde{A}=\begin{pmatrix}0&A\\\ A^{\dagger}&0\end{pmatrix},$ (25)

and proceed with the Hermitian matrix $\tilde{A}$. Below, we will discuss
extension of our work and point out some overlaps/connection to other works.

###  IV.1 Non-Sparse Matrix

Remind that our work assumes certain sparsity of the matrix $A$. At the core
of our algorithm is the application of $A$, which requires the ability to
simulate $\exp(-iAt)$ to be used by the quantum phase estimation. As in
harrow2009quantum , such simulation is doable if $A$ is sparse and row-
computable, thanks to the method in berry2007efficient . However, even if the
Hermitian matrix $A$ is non-sparse, it can be decomposed as

$A=B^{\dagger}B,$

and simulating the time evolution is still possible. If there exists a way to
prepare the following state

$C\sum_{i}|B_{i}|\ket{i}\ket{B_{i}},$

where $C$ is an overall constant, $B_{i}$ refers to the vector columns of $B$,
and $|B_{i}|$ is the $l_{2}$-norm of such vector, we can simulate $\exp(-iAt)$
by employing the method proposed in Ref. lloyd2014quantum , namely, the
density matrix exponentiation. With such an ability to simulate $\exp(-iAt)$,
quantum phase estimation of eigenvalues of $A$ is thus possible and hence we
can execute our quantum algorithm as above to find the maximum eigenvalue of
$A$ (possibly with some scaling/or normalization factor).

In fact, a similar problem was encountered in Ref. lloyd2014quantum , and the
authors proposed to use density matrix exponentiation combined with quantum
phase estimation to solve the problem called principle component analysis.
Roughly speaking, they aimed to find a few largest eigenvalues (and
corresponding eigenvectors) of the so-called covariant matrix, represented by
some density matrix $\sum$ and there exists a way to prepare $\sum$. Once the
exponent $\exp(-i\sum t)$ is obtained, running quantum phase estimation
algorithm (QPE) kitaev1995quantum ; brassard2002quantum  with the density
matrix $\sum$ itself also prepared as initial state allows us to reveal
highest values of the spectra by sampling. The method works well in the case
where a few top eigenvalues have much higher values than the remaining ones
(which is why they are called principle components). We remark that, in our
problem specifically, once $\exp(-iAt)$ is obtained through a known procedure,
then QPE can be used directly to reveal the maximum eigenvalue. We can, for
example, simply run QPE with $\exp(-iA2\pi)$ as the main unitary and the
maximally mixed state $I/n$ as the input state. We would approximately obtain
the following:

$\displaystyle\frac{1}{n}\Big{(}\sum_{i=1}^{n}\ket{\tilde{\lambda_{i}}}\bra{\tilde{\lambda_{i}}}\otimes\ket{E_{i}}\bra{E_{i}}\Big{)}.$
(26)

Sampling from the above state would in fact, reveal full spectra, including
the maximum eigenvalue. However, as simple as it may seem, this is not
efficient as the probability of obtain different eigenvalues is the same
($=1/n)$ and hence, would require as much as $\mathcal{O}(n)$ time to obtain
the maximum one. If we apply the controlled rotation to ancilla in order to
multiple the eigenvalues to the corresponding components, then this
essentially reduces to the quantum principle component analysis
lloyd2014quantum , discussed above.

An alternative solution to non-sparse $A$ is a direct result of
wossnig2018quantum , where they consider the problem of solving a dense linear
system. Given the particular data structure of $A$ and an oracle access to $A$
(see Lemma 1 in wossnig2018quantum ; kerenidis2016quantum ), they outline a
somewhat similar strategy to Ref. harrow2009quantum  in order to solve the
linear system. The goal is to coherently extract the singular values of $A$
(which are eigenvalues in the Hermitian case) and rotate the ancilla (to
multiply the inverse of the eigenvalue). Similarly using an adaptive method
introduced in wiebe2012quantum , we can multiply the eigenvalues instead, so
as to apply the matrix $A$ instead of $A^{-1}$. A potential caveat of this
approach is that the exponential speedup (with respect to size of matrix)
might not be obtained, but only a quadratic speedup, as discussed concretely
in wossnig2018quantum .

###  IV.2 Finding The Minimum Eigenvalue

For a Hermitian matrix, the conditional number of a matrix $\kappa$ is defined
as:

$\kappa=\frac{\lambda_{max}}{\lambda_{min}}.$

Note that this formula is only exact for a Hermitian matrix. In general,
singular values is taken into account instead of eigenvalues.

In principle, if the conditional number of a matrix is known, once the maximum
eigenvalue of a matrix is revealed by the above algorithm, then the minimum
eigenvalue could be estimated. However, in cases where the conditional number
is not exactly known (we might, for example, know its upper bound), then it is
hard to have a good estimate of the minimum eigenvalue. Here, we aim to tackle
such a problem using the algorithm that we have developed above. Let
$\\{\lambda_{i}\\}_{1}^{n}$ be eigenvalues of $A$. Since $A$ is Hermitian, so
if $A$ is invertible, i.e, does not contain zero eigenvalue (or one can use a
shifted matrix $A-cI$, more on this below), the eigenvalues of $A^{-1}$ is

$\Big{\\{}\frac{1}{\lambda_{i}}\Big{\\}}_{i=1}^{n}$

If $\lambda_{max}$ is the maximum eigenvalue of $A$, then $1/\lambda_{\max}$
is the minimum eigenvalue of $A^{-1}$. Similarly, if $\lambda_{\min}$ is the
minimum eigenvalue of $A$, then $1/\lambda_{\min}$ is the maximum eigenvalue
of $A^{-1}$. So instead of $A$, we simply apply $A^{-1}$ as in the HHL, then
the above procedure can find the maximum eigenvalue of $A^{-1}$, from which we
can find the minimum eigenvalue of $A$.

###  IV.3 Shifting by an Identity Matrix to Compute Ground-State Energy

As we have mentioned above, physical systems are generally characterized by a
Hamiltonian $H$, which is Hermitian. Thus, analyzing its spectrum, including
finding the ground state (and possibly a few low lying excited states) and its
corresponding energy, is one of the most important problems in physics. In our
algorithm, we can also consider a shift by an identity matrix $A-cI$ that
shifts the spectrum by a constant $c$. Thus one can compute the
$\max_{i}|\lambda_{i}-c|$. Shifting has been a useful technique, for example,
in finding the excited states and their spectra in many-body localization
yu2017finding . For example, the energy level statistics in potential many-
body localized systems can be gained from applying such a technique and be
sped up in principle by our quantum algorithm.

If one applies a large enough shift so that all the eigenvalues are negative,
then the one corresponding to the largest magnitude gives the lowest
eigenvalue of $A$ and the procedure also yields its eigenvector, i.e., the
ground state. In a similar way, if one adds a sufficient large positive shift,
so all the eigenvalues are positive, then one can find the largest algebraic
eigenvalues.

###  IV.4 Largest few eigenvalues

In Ref. motta2020determining , a quantum Lanczos algorithm was proposed to
compute ground and excited states of a Hamiltonian $H$, based on analogously
quantum imaginary time evolution $e^{-H\Delta\tau}$ (which is in turn
approximated by a unitary evolution under some other Hamiltonian) and the
classical Lanczos algorithm. In the classical Lanczos algorithm, the Krylov
space is built in powers of $H$, whereas in the particular quantum
implementation by motta2020determining , the Krylov space is built in powers
of $\exp(-\Delta\tau H)$. The two approaches (of $H$ or of $e^{-\Delta\tau
H}$) become identical in the limit $\Delta\tau\rightarrow 0$. Moreover, in the
Lanczos algorithm, a crucial aspect is to compute the following inner product
$v_{k}^{*T}Hv_{k}$ (and possibly with some power of $H$), where $v_{k}$ is
some vector with unit norm. In this regard, the procedure outlined in our work
is building up the Krylov subspace in terms of series of quantum states. It
can then be employed to compute the largest few eigenvalues in magnitude, by
using quantum computers to compute the $H$ matrix in the Krylov subspace and
the overlaps among the vectors $H^{m}x_{0}$ and $H^{n}x_{0}$. We then
classically diagonalize the resultant eigenvalue problem in the Krylov
subspace, yielding the largest few eigenvalues (in magnitude). Moreover, by
combining the shifting and the choice of applying $H$ or $H^{-1}$, we can
generalize our algorithm to obtaining a few largest or smallest eigenvalues
either in magnitude or in real value. Such capability can potentially lead to
more efficient access to, e.g., the ground-state energy and those of a few low
lying excited states, as many physics models possess Hamiltonians whose
matrices are spare.

##  V Conclusion

Inspired by many developments in quantum algorithm, we have outlined an
efficient quantum algorithm for estimating largest absolute eigenvalue of a
sparse Hermitian matrix. Possible extension to non-sparse and solving related
eigenvalues problems are discussed. The core routines of our algorithm are the
consecutive application of adaptive version of HHL algorithm harrow2009quantum
; wiebe2012quantum , plus the Hadamard test, both of which have proved to be
useful in many contexts, such as lloyd2013quantum ; motta2020determining ;
wiebe2012quantum , etc. Our work thus suggests an efficient quantum approach
to a particular computational problem, i.e., solving the eigenvalue problem,
for which classical algorithms requires polynomial time with respect to the
size of input matrix. Identifying difficult computational problems and finding
efficient solutions in both classical and quantum contexts have always been
one of the major drivers to improve computation capability and extend its
applications. Our work has contributed another small step to enrich the
existing quantum algorithm pool, and can be used as a subroutine (just as the
HHL algorithm) in other quantum algorithms or future applications.

Acknowledgements. This work was supported in part by the U. S. Department of
Energy, Office of Science, National Quantum Information Science Research
Centers, Co-design Center for Quantum Advantage (C2QA) under contract number
DE-SC0012704, in particular, for the design of the main algorithm, and by the
National Science Foundation under Grant No. PHY 1915165, in particular, for
the extension of the algorithm and its application to physical Hamiltonians.

## References

  * (1) Richard P Feynman.  Simulating physics with computers.  In Feynman and computation, pages 133–153. CRC Press, 2018.
  * (2) David Deutsch.  Quantum theory, the church–turing principle and the universal quantum computer.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 400(1818):97–117, 1985.
  * (3) Lov K Grover.  A fast quantum mechanical algorithm for database search.  In Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, pages 212–219, 1996.
  * (4) Peter W Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM review, 41(2):303–332, 1999.
  * (5) Andrew M Childs and Wim Van Dam.  Quantum algorithms for algebraic problems.  Reviews of Modern Physics, 82(1):1, 2010.
  * (6) Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders.  Efficient quantum algorithms for simulating sparse hamiltonians.  Communications in Mathematical Physics, 270(2):359–371, 2007.
  * (7) Aram W Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical review letters, 103(15):150502, 2009.
  * (8) Andrew M Childs, Robin Kothari, and Rolando D Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.
  * (9) Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum algorithms for supervised and unsupervised machine learning.  arXiv preprint arXiv:1307.0411, 2013.
  * (10) Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum principal component analysis.  Nature Physics, 10(9):631–633, 2014.
  * (11) Seth Lloyd, Silvano Garnerone, and Paolo Zanardi.  Quantum algorithms for topological and geometric analysis of data.  Nature communications, 7(1):1–7, 2016.
  * (12) Seth Lloyd, Maria Schuld, Aroosa Ijaz, Josh Izaac, and Nathan Killoran.  Quantum embeddings for machine learning.  arXiv preprint arXiv:2001.03622, 2020.
  * (13) Maria Schuld and Francesco Petruccione.  Supervised learning with quantum computers, volume 17.  Springer, 2018.
  * (14) Kosuke Mitarai, Makoto Negoro, Masahiro Kitagawa, and Keisuke Fujii.  Quantum circuit learning.  Physical Review A, 98(3):032309, 2018.
  * (15) Nhat A Nghiem, Xianfeng David Gu, and Tzu-Chieh Wei.  Constant-time quantum algorithm for homology detection in closed curves.  arXiv preprint arXiv:2209.12298, 2022.
  * (16) B. D. Clader, B. C. Jacobs, and C. R. Sprouse.  Preconditioned quantum linear system algorithm.  Phys. Rev. Lett., 110:250504, Jun 2013.
  * (17) Leonard Wossnig, Zhikuan Zhao, and Anupam Prakash.  Quantum linear system algorithm for dense matrices.  Physical review letters, 120(5):050502, 2018.
  * (18) Nathan Wiebe, Daniel Braun, and Seth Lloyd.  Quantum algorithm for data fitting.  Physical review letters, 109(5):050505, 2012.
  * (19) Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * (20) Moody T Chu and J Loren Watterson.  On a multivariate eigenvalue problem, part i: Algebraic theory and a power method.  SIAM Journal on scientific computing, 14(5):1089–1106, 1993.
  * (21) Beresford N Parlett, H Simon, and LM Stringer.  On estimating the largest eigenvalue with the lanczos algorithm.  Mathematics of computation, 38(157):153–165, 1982.
  * (22) Dianne P O’Leary, GW Stewart, and James S Vandergraft.  Estimating the largest eigenvalue of a positive definite matrix.  Mathematics of Computation, 33(148):1289–1292, 1979.
  * (23) Terence Tao and Van Vu.  Random matrices: universality of local eigenvalue statistics.  Acta mathematica, 206(1):127–204, 2011.
  * (24) A Yu Kitaev.  Quantum measurements and the abelian stabilizer problem.  arXiv preprint quant-ph/9511026, 1995.
  * (25) Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  Contemporary Mathematics, 305:53–74, 2002.
  * (26) Iordanis Kerenidis and Anupam Prakash.  Quantum recommendation systems.  arXiv preprint arXiv:1603.08675, 2016.
  * (27) Xiongjie Yu, David Pekker, and Bryan K. Clark.  Finding matrix product state representations of highly excited eigenstates of many-body localized hamiltonians.  Phys. Rev. Lett., 118:017201, Jan 2017.
  * (28) Mario Motta, Chong Sun, Adrian TK Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando GSL Brandão, and Garnet Kin Chan.  Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution.  Nature Physics, 16(2):205–210, 2020.

##  Appendix A Analysis of Error

We remark that the original classical algorithm described in Sec. II only
produces an approximation to $\lambda_{n}$. Here we provide a more detailed
analysis of the error. It is clear that the error vanishes in the limit
$k\rightarrow\infty$.

For simplicity of the presentation, we assume $A$ to be real and symmetric,
but the extension to the Hermitian case is straightforward. Let
$x_{k}=A^{k}x_{0}$ and $\bar{x}_{k}$ be:

$\bar{x}_{k}\equiv\lambda_{n}^{k}c_{n}E_{n}.$

We note that, the exact eigenvalue $\lambda_{n}$ is:

$\displaystyle\lambda_{n}=\frac{\bar{x}_{k+1}\cdot\bar{x}_{k}}{\bar{x}_{k}\cdot\bar{x}_{k}}.$
(27)

The approximated eigenvalue is:

$\displaystyle\bar{\lambda}_{n}=\frac{Ax_{k}\cdot x_{k}}{x_{k}\cdot x_{k}}.$
(28)

From Eqn. (9), we have:

$\displaystyle
Ax_{k}=\lambda_{n}^{k+1}\Big{(}c_{1}\frac{\lambda_{1}^{k+1}}{\lambda_{n}^{k+1}}E_{1}+c_{2}\frac{\lambda_{2}^{k+1}}{\lambda_{n}^{k+1}}E_{2}+\cdots+c_{n}E_{n}\Big{)}.$
(29)

By using the property that $\\{E_{i}\\}_{1}^{n}$ are orthonormal, we have:

$\displaystyle Ax_{k}\cdot
x_{k}=\lambda_{n}^{k+1}\lambda_{n}^{k}\Big{(}c_{1}^{2}\frac{\lambda_{1}^{k+1}}{\lambda_{n}^{k+1}}\frac{\lambda_{1}^{k}}{\lambda_{n}^{k}}+c_{2}^{2}\frac{\lambda_{2}^{k+1}}{\lambda_{n}^{k+1}}\frac{\lambda_{2}^{k}}{\lambda_{n}^{k}}+\cdots+c_{n}^{2}\Big{)}.$
(30)

We also have:

$\displaystyle x_{k}\cdot
x_{k}=\lambda_{n}^{2k}\Big{(}c_{1}^{2}\frac{\lambda_{1}^{2k}}{\lambda_{n}^{2k}}+c_{2}^{2}\frac{\lambda_{2}^{2k}}{\lambda_{n}^{2k}}\cdots+c_{n}^{2}\Big{)}.$
(31)

Straightforward computations yield:

$\displaystyle\bar{x}_{k}\cdot\bar{x}_{k}=\lambda_{n}^{2k}c_{n}^{2},$ (32)
$\displaystyle\bar{x}_{k+1}\cdot\bar{x}_{k}=\lambda_{n}^{k+1}\lambda_{n}^{k}c_{n}^{2}.$
(33)

Now we aim to compute the absolute ratio between the approximated eigenvalue
and exact eigenvalue:

$\displaystyle\Big{|}\frac{\bar{\lambda}_{n}}{\lambda_{n}}\Big{|}$
$\displaystyle=\Big{|}\frac{Ax_{k}\cdot
x_{k}}{\bar{x}_{k+1}\cdot\bar{x}_{k}}\cdot\frac{\bar{x}_{k}\cdot\bar{x}_{k}}{x_{k}\cdot
x_{k}}\Big{|}.$ (34)

A simple substitution yields

$\displaystyle\frac{\bar{x}_{k}\cdot\bar{x}_{k}}{x_{k}\cdot x_{k}}$
$\displaystyle=\frac{1}{1+\sum_{i}(\frac{c_{i}}{c_{n}})^{2}\cdot(\frac{\lambda_{i}}{\lambda_{n}})^{2k}}<1.$
(35)

Therefore, we arrive at

$\Big{|}\frac{\bar{\lambda}_{n}}{\lambda_{n}}\Big{|}<\Big{|}\frac{Ax_{k}\cdot
x_{k}}{\bar{x}_{k+1}\cdot\bar{x}_{k}}\Big{|}.$

Another round of substitution yields:

$\displaystyle\frac{Ax_{k}\cdot x_{k}}{\bar{x}_{k+1}\cdot\bar{x}_{k}}$
$\displaystyle=1+\sum_{i=1}^{n-1}\frac{c_{i}^{2}}{c_{n}^{2}}\frac{\lambda_{i}^{k}}{\lambda_{n}^{k}}\frac{\lambda_{i}^{k+1}}{\lambda_{n}^{k+1}}.$
(36)

We then make use of the following conditions:

$\bullet$ Since $\lambda_{i}/\lambda_{n}<$1 for all $i<n$ then

$\frac{\lambda_{i}^{k}}{\lambda_{n}^{k}}\frac{\lambda_{i}^{k+1}}{\lambda_{n}^{k+1}}<\frac{\lambda_{i}^{2k}}{\lambda_{n}^{2k}}.$

$\bullet$ $0\leq\lambda_{1}\leq\lambda_{2}\leq\cdots<\lambda_{n}$ then for all
$i<n$:

$\frac{\lambda_{i}}{\lambda_{n}}\leq\frac{\lambda_{n-1}}{\lambda_{n}}=p<1.$

$\bullet$ All the components $c_{1},c_{2},...,c_{n}$ are bounded above,
therefore, the ratio $c_{i}/c_{n}$ is bounded above by some value $K$ for all
$i$.

Therefore, the ratio in Eqn. 36 is upper bounded,

$\displaystyle\frac{Ax_{k}\cdot
x_{k}}{\bar{x}_{k+1}\cdot\bar{x}_{k}}<1+(n-1)\cdot K\cdot p^{2k}.$ (37)

Altogether, we have:

$\displaystyle\Big{|}\frac{\bar{\lambda}_{n}}{\lambda_{n}}\Big{|}<1+(n-1)\cdot
K\cdot p^{2k}.$ (38)

In order to estimate the eigenvalue $\lambda_{n}$ to some multiplicative error
$\delta$, we expect:

$\displaystyle\Big{|}\frac{\bar{\lambda}_{n}}{\lambda_{n}}\Big{|}<1+\delta.$
(39)

We thus require that

$(n-1)\cdot K\cdot p^{2k}<\delta,$

which is satisfied if

$k\in\Omega\big{(}\log(1/\delta)/\log(1/p)\big{)}.$

Therefore, the running time of the algorithm also depends on the factor $p$,
which is the ratio between two largest (in absolute value) eigenvalues. For
$p$ close to 1, the above analysis shows that the algorithm is highly
efficient, as the running time scales as the logarithmic of the inverse error.

##  Appendix B Further Elaboration on Time Complexity

In Sec. III.1, we have outlined a quantum algorithm for estimating largest
absolute eigenvalue $\bar{\lambda_{n}}$. The algorithm works by iterating, for
example, $k$ times the application of matrix $A$ to some initial state
$\ket{x_{0}}$, which results in the running time linearly in $k$. However, we
have remarked about the square factor $k^{2}$ in the running time (Eqn. 23),
which is due to the accumulated error in each iteration. Here we elaborate
further on such subtlety.

Let the initial state be $x_{0}$ (we ignore the ket notation, for simplicity).
Ideally, after first iteration, we expect to obtain: $Ax_{0}$. However, as
also noted in [7, 18], the quantum application of $A$ results in a vector that
is $\epsilon$-close to $Ax_{0}$, denoted as $\overline{Ax_{0}}$. After another
round of application, we obtain an $\epsilon$-close to $A\overline{Ax_{0}}$,
denoted as $\overline{A^{2}x_{0}}$. Since $\overline{Ax_{0}}$ is
$\epsilon$-close to $Ax_{0}$, i.e,

$\displaystyle|\overline{Ax_{0}}-Ax_{0}|\leq\mathcal{\epsilon},$ (40)

which means that

$\displaystyle|A\cdot\overline{Ax_{0}}-A\cdot
Ax_{0}|\equiv|A\overline{Ax_{0}}-A^{2}x_{0}|\leq\epsilon$ (41)

Using triangle inequality, we obtain:

$\displaystyle|\overline{A^{2}x_{0}}-A^{2}x_{0}|$ (42)
$\displaystyle\leq|\overline{A^{2}x_{0}}-A\overline{Ax_{0}}|+|A\overline{Ax_{0}}-A^{2}x_{0}|$
(43) $\displaystyle\leq\epsilon+\epsilon=2\epsilon.$ (44)

Proceeding similarly, we can deduce that after $k$ iterations:

$\displaystyle|\overline{A^{k}x_{0}}-A^{k}x_{0}|\leq k\epsilon.$ (45)

The error accumulates at least linearly. If we desire an overall error
tolerance $\mathcal{O}(\epsilon)$, in each iteration we must have error
$\sim\mathcal{O}(\epsilon/k)$. Recall that each iteration is an application of
HHL-like algorithm [18, 7], which requires time depending on error as
$\mathcal{O}(1/\epsilon)$. Replacing

$\displaystyle\epsilon\rightarrow\epsilon/k$

We have the running time in each iteration grows as $\mathcal{O}(k/\epsilon)$.
The number of iterations in total is $k$ so we need to multiply the above time
by $k$, which eventually yields $\mathcal{O}(k^{2}/\epsilon)$.

##  Appendix C Details of the HHL-like approach for multiplication

In this section, we explicitly show how the state in Eqn. (14) is obtained by
following the original HHL algorithm [7], except for multiplication and
particularly its adaptive version [18], which form the core subroutine used in
our algorithm to consecutively apply a matrix $A$.

Let the initial state $\ket{x_{0}}$ be given, and its decomposition in the
basis consisting of eigenvectors of A, $\\{\ket{u_{i}}\\}_{i=1}^{N}$ (with
corresponding eigenvalues $\\{\lambda_{i}\\}_{i=1}^{N}$ ) as follows:

$\ket{x_{0}}=\sum_{i=1}^{N}\beta_{i}\ket{u_{i}}$

Following HHL, we first prepare state:

$\displaystyle
U_{0}\ket{0}\otimes\ket{x_{0}}=\ket{\Phi_{0}}=\sqrt{\frac{2}{T}}\sum_{\tau=0}^{T-1}\sin(\frac{\pi(\tau+1/2)}{T})\ket{\tau}\otimes\ket{x_{0}}$
(46)

We then apply the following controlled unitary

$\displaystyle
U_{T}=\sum_{\tau=0}^{T}\ket{\tau}\bra{\tau}\otimes\exp(\frac{-iA\tau
t_{0}}{T})$ (47)

to the above state, followed by a Fourier transform over the first register.
We approximately arrive at (we will elaborate on the approximation later)

$\displaystyle\ket{\phi}=\sum_{i=1}^{N}\beta_{i}\ket{\lambda_{i}}\ket{u_{i}}.$
(48)

Next, we append an ancilla initialized in $\ket{0}$ and perform the following
rotation on it controlled on the phase register:

$\displaystyle\ket{\lambda_{i}}\ket{0}\rightarrow\ket{\lambda_{i}}\big{(}C\lambda_{i}\ket{0}+\sqrt{1-C^{2}\lambda_{i}^{2}}\ket{1}\big{)},$

where the constant $C$ secures the normalization condition and is chosen so
that

$\displaystyle C\lambda_{max}\leq 1.$

Since we have the condition: $1/\kappa\leq\lambda_{i}\leq 1$ for all i,
therefore, the condition for $C$ can be further simplified as

$\displaystyle C\leq\kappa.$

After the above step, the state $\ket{\phi}$ is transformed to:

$\displaystyle\ket{\phi_{1}}=\sum_{i=1}^{N}\beta_{i}\ket{\lambda_{i}}\ket{u_{i}}\big{(}C\lambda_{i}\ket{0}+\sqrt{1-C^{2}\lambda_{i}^{2}}\ket{1}\big{)}$

We note that, our desired vector $A\ket{x_{0}}$ is

$\displaystyle A\ket{x_{0}}=\sum_{i=1}^{N}\lambda_{i}\beta_{i}\ket{u_{i}}.$

Therefore, the state $\ket{\phi_{1}}$, after we uncompute the phase register,
which could be done by running $U_{T}^{\dagger}$ and $U_{0}^{\dagger}$, would
approximately become

$\displaystyle\ket{\phi_{1}}=C\|A\ket{x_{0}}\\!\|\ket{A\ket{x_{0}}}\ket{0}+\ket{G_{1}}\ket{1},$
(49)

where $\ket{G_{1}}$ (not normalized properly), as mentioned earlier, is some
unimportant state; and $\|A\ket{x_{0}}\|$ is the length ($l_{2}$-norm) of the
vector $A\ket{x_{0}}$. We have emphasized previously that we only obtained an
approximation to Eqn. (48) (and hence, the state Eqn. (49) as well), which is
due to the finite-bit encoding of the phase value [7]. A thorough analysis has
been provided in [7], showing that if $t_{0}=\mathcal{O}(\kappa/\epsilon)$,
then the total final error is $\epsilon$. By the final error, we mean that the
state $\ket{A\ket{x_{0}}}$ in Eqn. (49) is $\epsilon$-close to the actual
solution. Throughout the discussions below, we will keep using the notation
$\ket{A\ket{x_{0}}}$ to denote the approximated state with the understanding
that it is just $\epsilon$-close to the real solution (in the main text, we
actually distinguish the two).

Instead of measuring the ancilla and trying to project to the state
$|Ax_{0}\rangle$, as did in the HHL algorithm, our next step is to apply $A$
again to the state $\ket{A\ket{x_{0}}}$, or more precisely, the vector
$Ax_{0}$ itself. (We note that a simpler alternative way is to continue to
apply the controlled rotation on an additional ancilla without undoing the
phase estimation; see Appendix D below.) Again, let the decomposition of
$\ket{A\ket{x_{0}}}$ in the basis $\\{\ket{u_{i}}\\}_{i=1}^{N}$ as,

$\displaystyle\ket{A\ket{x_{0}}}=\sum_{i=1}^{N}\beta^{\prime}_{i}\ket{u_{i}}.$

It is straightforward to see that since
$A\ket{x_{0}}=\sum_{i=1}^{N}\lambda_{i}\beta_{i}\ket{u_{i}}$, then the
coefficient $\beta^{\prime}_{i}$ in the above state is just

$\displaystyle\beta^{\prime}_{i}=\frac{\lambda_{i}\beta_{i}}{\|A\ket{x_{0}}\|}=\frac{\lambda_{i}\beta_{i}}{\sqrt{\sum_{i=1}^{N}|\lambda_{i}\beta_{i}|^{2}}}.$

Following the same procedure as we did for the application of $A$ to
$\ket{x_{0}}$, we first prepare the state $\ket{\Phi_{0}}$ and append another
ancilla $\ket{0}$, then we obtain approximately the following state

$\displaystyle\ket{\phi_{2}}$
$\displaystyle=C\|A\ket{x_{0}}\|\Big{(}\sum_{i=1}^{N}\beta^{\prime}_{i}\ket{\lambda_{i}}\ket{u_{i}}\Big{)}\ket{0}\ket{0}+U_{T}\ket{\Phi_{0}}\ket{G_{1}}\ket{1}\ket{0}.$
(50)

We take a specific look at the first part of the above state (we ignore
everything entangled with $\ket{G_{1}}$ as they are not important). At this
step, we would rotate the ancilla $\ket{0}$ controlled by the phase register.
However, a subtle point here is that this time, the rotation needs to be
controlled further by the qubit $\ket{0}$ that is resulted from the previous
application of $A$ (note that there are two ancillary qubits $\ket{0}\ket{0}$
in the above state). With such a detail in mind, we will obtain the following
state:

$\displaystyle
C\|A\ket{x_{0}}\|\Big{(}\sum_{i=1}^{N}\beta^{\prime}_{i}\ket{\lambda_{i}}\ket{u_{i}}\ket{0}(C\lambda_{i}\ket{0}+\sqrt{1-C^{2}\lambda^{2}_{i}}\ket{1})\Big{)}+U_{T}\ket{\Phi_{0}}\ket{G_{1}}\ket{1,0}.$
(51)

Now we uncompute the phase register (and omit it afterwards) to obtain:

$\displaystyle
C^{2}\|A\ket{x_{0}}\|\sum_{i=1}^{N}\beta^{\prime}_{i}\lambda_{i}\ket{u_{i}}\ket{00}+\ket{G_{2}}\ket{01}+\ket{G_{1}}\ket{1}\ket{0}.$
(52)

We remark that since $\\{\beta^{\prime}_{i}\\}_{i=1}^{N}$ are the coefficients
of the decomposition of the normalized state $\ket{A\ket{x_{0}}}$ in the basis
$\\{\ket{u_{i}}\\}_{i=1}^{N}$, then multiplying with the factor
$\|A\ket{x_{0}}\|$ (which is the actual length of the original vector
$A\ket{x_{0}}$) would yield the coefficient of the original vector
$A\ket{x_{0}}$ in the same basis. From Eqn. C, it directly leads to:

$\displaystyle\|A\ket{x_{0}}\|\sum_{i=1}^{N}\beta^{\prime}_{i}\lambda_{i}\ket{u_{i}}=A(A\ket{x_{0}})=A^{2}\ket{x_{0}}.$
(53)

Note that this is a vector that is not normalized. Hence, the Eqn. (52) can be
written simpler as

$\displaystyle
C^{2}\|A^{2}\ket{x_{0}}\|\ket{A^{2}\ket{x_{0}}}\ket{00}+\ket{G_{2}}\ket{01}+\ket{G_{1}}\ket{10}.$
(54)

Iterating the procedure $k$ times, we eventually obtain what we got in Eqn.
(14). Such a state allows us to estimate $\lambda_{max}$ as we have outlined
an efficient procedure in Sec. III.1.

As a final remark, the property that we have used in Eqn. (53) is indeed very
useful. After each iteration, there is a length factor (such as
$\|A\ket{x_{0}}\|$ in the above) but eventually it will get ‘absorbed’ in the
next iteration, which is very convenient as we only care about the final
vector $A^{k}\ket{x_{0}}$ to estimate $\lambda_{max}$. Otherwise, we might
have to estimate such a length after each iteration and hence would incur a
further linear time $\mathcal{O}(k)$ processing for $k$ iterations.

##  Appendix D Improved Approach

It is quite surprising that the above procedure can be improved by a factor of
$k$–the number of iterations. The application of $A$ on $\ket{x_{0}}$ can be
written explicitly as:

$\displaystyle A\ket{x_{0}}=\sum_{i=1}^{N}\lambda_{i}\beta_{i}\ket{u_{i}}.$
(55)

Likewise, another application of $A$ yields:

$\displaystyle
A^{2}\ket{x_{0}}=\sum_{i=1}^{N}\lambda^{2}_{i}\beta_{i}\ket{u_{i}}.$ (56)

After $k$ times, we have:

$\displaystyle
A^{k}\ket{x_{0}}=\sum_{i=1}^{N}\lambda^{k}_{i}\beta_{i}\ket{u_{i}},$ (57)

as we have described earlier.

As we have emphasized, the key step in our quantum algorithm is the execution
of $\exp(-iAt)$ in QPE to extract the eigenvalues of $A$, and perform the
controlled rotation. As being shown in [7, 18] (also it can be found
explicitly in our elaboration Appendix C), after running QPE and controlled
rotation in the first iteration, we would have the following state:

$\displaystyle\sum_{i=1}^{N}\beta_{i}\ket{\lambda_{i}}\ket{u_{i}}\Big{(}C\lambda_{i}\ket{0}+\sqrt{1-C^{2}\lambda^{2}_{i}}\ket{1}\Big{)}$
$\displaystyle=C\sum_{i=1}^{N}\beta_{i}\lambda_{i}\ket{\lambda_{i}}\ket{u_{i}}\ket{0}+\ket{G_{1}}\ket{1}.$
(58)

Now, instead of uncompute the phase and repeat the process, we simply append
the ancilla $\ket{0}$ and rotate the ancilla controlled by the phase register
and the state $\ket{0}$ resulted from the previous rotation. In other words,
we have the following:

$\displaystyle
C\sum_{i=1}^{N}\beta_{i}\lambda_{i}\ket{u_{i}}\ket{0}\Big{(}C\lambda_{i}\ket{0}+\sqrt{1-C^{2}\lambda_{i}^{2}}\ket{1}\Big{)}+\ket{G_{1}}\ket{1}\ket{0}$
(59)
$\displaystyle=C^{2}\sum_{i=1}^{N}\beta_{i}\lambda_{i}^{2}\ket{00}+\ket{G_{2}}\ket{01}+\ket{G_{1}}\ket{10}.$
(60)

Keep appending the ancilla initialized in $\ket{0}$ and doing controlled
rotation (note that at the $r$-th iteration, the controlled qubits consists of
the phase register $\ket{\lambda_{i}}$ and $(r-1)$ qubits $\ket{0}$ resulted
from previous rotations) $k$ times, we have our desired state as in Eqn. (14).

The original HHL matrix inversion algorithm [7, 18] uncomputes the phase
register (returning it $|0\dots 0\rangle$) after each iteration. Hence if we
repeat the same process $k$ times, we need to have $\mathcal{O}(k)$ times
simulating $\exp(-iAt)$. By keeping the phase register, we showed that the
running time can be reduced by a factor of $k$, which leads to our actual
running time (24).


# Photothermal effects in small gold nanorod aggregates for therapeutic
applications

Dheeraj Pratapa,, Ram Krishna Shahb, Sameer Khandekarb, Sanjeev Sonia

aBiomedical Applications Group, CSIR-Central Scientific Instruments
Organization, Sector-30C, Chandigarh-160030
bDepartment of Mechanical Engineering, Indian Institute of Technology Kanpur,
Kanpur-208016, India
Corresponding author’s email<EMAIL_ADDRESS>
###### Abstract

Proximity of the metal nanoparticles enhance the plasmonic coupling and shifts
the resonance. This article presents a numerical study of the photothermal
effect in aggregates of small gold nanorods considering the ordered as well as
random aggregates. In the ordered aggregates, there is lateral coupling which
causes blueshifts in the plasmonic resonance, while in the random aggregates
there are redshifts in the plasmonic resonance. The plasmon response of latter
could be tailored up to the second infrared biological therapeutic window. It
has been observed that the aggregates show higher absorption power and
therefore, higher temperature rise compared to the single gold nanorod or
monodispersive nanorods. The absorption resonance peak position of the random
aggregate depends on the incident and polarization angles of the incident
light. The aggregation of the nanoparticles often inherently occurs in the
biological medium which affects the photothermal process. This study helps to
understand the photothermal heating of nanoparticle aggregates and the use of
the optimal light source concerning the absorption peak of the aggregates
suspension for therapeutic uses.

_Keywords_ Nanorods $\cdot$ Photothermal $\cdot$ Optical properties $\cdot$
Aggregation $\cdot$ Scattering $\cdot$ Absorption

##  1 Introduction

Metal nanoparticles have applications in several areas like plasmonic
photothermal treatment of tumor Nam et al. (2009); Choi et al. (2011); Pattani
et al. (2015); Li et al. (2008); Wang et al. (2013), drug delivery Adeyemi and
Sulaiman (2015); Janib et al. (2010); Babu et al. (2013), imaging Lin et al.
(2010); Wiederrecht (2004); Syed et al. (2017); Wang et al. (2020),
nanofluidics Donner et al. (2011); Gómez-Villarejo et al. (2017),
nanophotonics Dionne and Atwater (2012); Fang and Sun (2015), sensors Tang et
al. (2014); Chen et al. (2019); Lee et al. (2014), catalysis Roldan Cuenya
(2013); Scholten et al. (2012); Yan et al. (2010), to name a few. Verity of
nanoparticles like nanosphere Soni et al. (2015), nanorod Chen et al. (2017),
nanoring Larsson et al. (2007), nanocage Skrabalak et al. (2008), nanocube Fan
et al. (2008), triangular nanoplate Apyari et al. (2018), nanostar Dutta et
al. (2019), nanodumble Wang et al. (2019) etc. are reported in literature.
These nanoparticles are prepared from several plasmonic metals like gold,
silver, copper and bimetallic Soni et al. (2015); Apyari et al. (2018);
Ruparelia et al. (2008); Aioub et al. (2017). For all such metallic
candidates, gold is highly preferable due to its chemical inertness and ease
of surface functionalization with a verity of surface ligands, such as
proteins and antibodies, by forming stable bonding pairs. In the case of
biological applications, the nanoparticles get aggregated within the medium.
In some cases, an induced aggregation also can be done inside the biological
medium. Sun et al. showed the salt-induced nanoparticle aggregation could
enhance the photothermal efficiency and photoacoustic imaging Sun et al.
(2016). Using the nuclear-targeted gold nanorods extremely low intensity near-
infrared (NIR) activated photothermal therapy can be realised Pan et al.
(2017). The gold nanorods naturally accumulate in the lysosomes of the cells
after reaching inside it and create clusters of the nanorods where a strong
plasmonic coupling between gold nanorods exist Ungureanu et al. (2011); Shen
et al. (2013). Chen et al. reported the photothermal therapy of a single cell
where gold nanorods had been accumulated within the cell Chen et al. (2017).
In this case, the plasmonic resonance of the monodispersive nanorods was at
800 nm and they used a pulsed laser of 800 nm wavelength for photothermal
effect. In the aggregate of the gold nanorods, there should be a resonance
shift due to coupling. The use of the corresponding resonance wavelength laser
could be better option but there was no mention of laser which have peak at
the shifted wavelength. Xie et al. reported a numerical investigation on
absorption and plasmonic coupling in the near-field limit of the spherical
nanoparticle clusters numerically Xie et al. (2019). Gu et al. used
monodispersive gold nanorods for plasmonic photothermal therapy and showed
that if the concentration of nanorods was increased then there was an increase
of plasmonic coupling Gu et al. (2018).

Upon irradiation, the nanoparticles generate heat due to photothermal effect.
Almost all applications work on the concept of plasmonic resonance. When the
nanoparticle absorbs the light then it gets heated and act as a nano source of
heat. Heat generation in such nanoparticles is affected by several parameters
like constituent materials, particles size and shape, surrounding medium,
irradiation intensity as well as wavelength, angle of incidence etc. Plasmonic
heating of two gold nanospheres have been reported by Baffou et al. Baffou et
al. (2010). The optical properties for spherical nanoparticles have been
extensively reported in literature Lermé et al. (1998); Cottancin et al.
(2006); Kreibig and Vollmer (2013). The shape of the nanoparticles will also
affect the plasmonic coupling, and hence the resonance shift in the aggregate.
In literature, heat generation characteristics of plasmonic nanoparticle
aggregates is not reported. For achieving the plasmonic resonance peak at a
higher wavelength in the NIR range, the gold nanorod is a good candidate
because of ease of synthesis and tunability of plasmon response. One can tune
the plasmonic resonance peak by varying the aspect ratio of the nanorod;
however increasing the aspect ratio increases the length of the nanorod. Very
large size nanoparticles are not desirable for biomedical in vivo
applications. For biological applications, there is a high possibility that
the nanoparticles turn into aggregates because of interaction with the complex
body fluids. Therefore, we consider the aggregates of small gold nanorods for
our current study. Most of the synthesised gold nanorods have the geometry of
a cylinder capped with hemispheres at both ends. When a particle is
geometrically symmetric, like a sphere, then its optical properties and effect
of heating could be computed analytically. However, for the case of random
shape particle aggregate, analytical solutions are difficult to obtain.
Further, asymmetric particle’s interaction with light also depends on the
polarization of the light.

In this article, we consider small sized gold nanorods and their ordered and
random aggregates in water to study the optical properties and photothermal
heating numerically. The organization of the article is as follows: Section-2
describes the theoretical background of optical scattering and heating by the
nanorods aggregate; Section-3 presents the numerical simulation of
photothermal effect/interaction the gold nanorods and their aggregates and the
results. In the final Section-4, we conclude our results.

##  2 Scattering by nanoparticle aggregates

When light interacts with a metallic nanoparticle then the free electrons
start oscillating along with the electric field of the light. This oscillation
of the free electrons is called localised surface plasmon (LSP) Maier (2007).
When the frequency of light matches with the natural frequency of the
nanoparticles then a resonance occurs which is known as localised surface
plasmon resonance (LSPR) Maier (2007). The oscillation of electrons inside the
nanoparticle due to the rapid changing electric field behaves as a varying
local electron current which produces the resistive (Joule) heating in the
nanoparticle Baffou and Quidant (2013); Baffou (2017). This localised heat
gets diffused into the surrounding medium and the nanoparticle acts as a
nanosource of heat.

Consider light wave of frequency $\omega$ and free space wavelength
$\lambda_{0}$ having electric and magnetic field components E0(r,t) and
B0(r,t) interacts with the an assembly or aggregate of nanoparticles situated
in a surrounding medium (water) of permittivity $\varepsilon_{s}$. After
interacting with assembly, some parts of light gets extinct via absorption and
scattering. Geometrical factors which quantify the absorption, scattering and
extinction are called absorption cross-section ($\sigma_{\mathrm{abs}}$),
scattering cross-section ($\sigma_{\mathrm{sca}}$) and extinction cross-
section ($\sigma_{\mathrm{ext}}$) respectively. These quantities are defined
as Mishchenko et al. (2002),

$\sigma_{\mathrm{abs}}=\frac{W_{\mathrm{abs}}}{P_{0}},\quad\sigma_{\mathrm{sca}}=\frac{W_{\mathrm{sca}}}{P_{0}},\quad\sigma_{\mathrm{ext}}=\sigma_{\mathrm{abs}}+\sigma_{\mathrm{sca}},$
(1)

where $W_{\mathrm{abs}}$ and $W_{\mathrm{sca}}$ are the absorbed and scattered
energy by the nanoparticle aggregate and $P_{0}$ is the irradiance (power per
unit area) of the incident light. The absorbed energy by the nanoparticle
aggregate is given by the following expression,

$W_{\mathrm{abs}}=\int_{V}w(\textbf{r})d\textbf{r}$ (2)

where
$w(\textbf{r})=(\varepsilon_{s}\omega/2)\mathrm{Im}(\varepsilon_{\omega})|\textbf{E}(\textbf{r})|^{2}$
is the dissipative power density within the nanoparticle aggregate and
$\varepsilon_{\omega}$ is the permittivity of the material of the
nanoparticle. Due to the absorbed power/dissipative power density, the
nanoparticle acts as a heat source and dissipates/diffuses the heat to the
surrounding medium (e. g. water or tissue). The main mode of heat transfer in
plasmonics is diffusion. To get the temperature field, $T$, of the
nanoparticle and its surrounding medium we have to solve the heat diffusion
equation,

$\rho_{\omega}c_{p,\omega}\partial_{t}T-\kappa_{\omega}\nabla^{2}T=w,$ (3)

inside the nanoparticle, where $\rho_{\omega}$, $c_{p,\omega}$,
$\kappa_{\omega}$ and $T$ are density, specific heat capacity at constant
pressure and thermal conductivity of the material of the nanoparticle, and
local temperature respectively; and

$\rho_{s}c_{p,s}\partial_{t}T-\kappa_{s}\nabla^{2}T=0,$ (4)

outside the nanoparticle, where $\rho_{s}$, $c_{p,s}$ and $\kappa_{s}$ are
density, specific heat capacity at constant pressure and thermal conductivity
of the surrounding medium (water) Baffou and Quidant (2013); Baffou (2017).

For the symmetric geometries, like a sphere, the analytical solutions for the
optical cross-sections, dissipative power and temperature field can be
obtained under the imposed boundary condition. However, analytical solutions
are difficult for irregular shapes or geometries. For random shaped
nanoparticles, numerical modelling tools like boundary element method (BEM) De
Abajo and Howie (2002), Green’s dyadic method (GDM) Martin et al. (1995);
Girard et al. (2008), discrete dipole approximation (DDA) Girard (2005);
Yurkin and Hoekstra (2007), multiple elastic scattering of multipole
expansions (MESME) Myroshnychenko et al. (2008); de Abajo (1999), finite
difference time domain (FDTD) Radhakrishnan and Murugesan (2014); Mohsin and
Salim (2018), finite element method (FEM) Bachelier et al. (2008); Alizadeh et
al. (2019) etc. are used for solution. In the present work, we have used
COMSOL Multiphysics® Multiphysics (1998), a FEM based commercial tool to model
the scattering, absorption and heating of nanorod aggregates.

Figure 1: Schematic of (a) a gold nanorod of length 25 nm and diameter 5 nm.
Default three-dimensional views of geometrical structures of (b-d) ordered
aggregates made of nine nanorods C1, C2 and C3; and (e-g) random aggregates
made of nine nanorods C4, C5 and C6. The projected views of the aggregates in
the $xy$, $yz$ and $zx$ planes are shown in the Appendix Fig. A.1.

Figure 1 delineates the geometries of a gold nanorod and its aggregates. The
central part of the nanorod has a cylindrical shape of length 20 nm and
diameter 5 nm. At both ends of the 20 nm cylinder, there are hemispheres of
diameter 5 nm. The total length of the gold nanorod is 25 nm and the aspect
ratio is 5. The considered nanorod dimension and aspect ratio is not
restricted, one could take nanorods of other size and aspect ratio too.
Choosing the larger single particle may increase the scattering and reduce the
heating. Here, our main objective is to study the plasmonic coupling and
heating in aggregates of small nanorods, as nanoparticles of size below 50 nm
are preferred for therapeutic applications. There are only nine nanorods
assumed in each aggregate. In the first three aggregates the nanorods are
arranged in an ordered manner (C1, C2 and C3), and in the other three
aggregates the nanorods are oriented randomly (C4, C5 and C6). The projected
views of all six considered aggregates in the $xy$, $yz$ and $zx$ planes are
provided in the Appendix Fig. A.1. We choose only nine nanorods to construct
the ordered and random aggregates design so that aggregates should not become
very large and nanorods could be arranged symmetrically also for comparison
with the random assembly configurations. We did not consider a linear chain of
nine nanorods arranged in end-to-end (longitudinal coupling) because in such a
long chain large redshift is obvious, and practically such a singular
arrangement will occur rarely in reality. All nanorods (nine) are identical in
each assembly so that their optical properties and photothermal effects could
be compared. Choosing a larger size of single particle may increase the
scattering and reduce the heating. For instance, a nanosphere of 100 nm
diameter gives SPR in the mid-visible range and at this scale; the scattering
phenomena dominates over the absorption. If we want to tune the LSPR of the
aggregates of nanospheres then there would be too much scattering because of
the large size. For the biomedical in vivo application, smaller nanoparticles
are preferred. The small nanorod of aspect ratio 5 itself gives the LSPR at
near-infrared (NIR) wavelength range which helps to tune the LSPR of the
aggregates for higher wavelength in the second biological window easily. We
have chosen three ordered assemblies (C1, C2 and C3) for comparison. For
practical cases, such ordered nature of assembly would not normally occur.
Naturally, there could be countless type of assemblies possible, and even in a
single assembly, the nanorods could have an infinite number of orientations.
Such a huge possibility of the random system is difficult to model in a single
study, therefore, we have considered here only three aggregates of the random
orientation of nine nanorods to understand the photothermal behaviour of the
same.

To model the nanorod and their aggregates to study the photothermal effects in
the COMSOL®, we first carried out the scattering of a gold nanosphere of 100
nm diameter in water, heated by the light of intensity 1 mW/$\mu$m2 at a
wavelength 530 nm as a benchmark study, which was numerically reported by
Baffou et al. using the BEM method Baffou et al. (2010), and later reproduced
by Chen et. al. using COMSOL Multiphysics® Chen et al. (2012). Results of the
benchmark study are given in the Appendix Fig. B.1. Our results are very
closely matching with the results of Chen Chen et al. (2012) and have slight
deviation compared with the results of Baffou Baffou et al. (2010). A slight
variation of results is due to the fact that Baffou used the BEM and Chen used
the FEM based COMSOL® and both approaches have different working principle.

Figure 2: (a) Absorption cross section and (b) absorption power of the gold
nanorod aggregates C1, C2, C3, C4, C5 and C6 dispersed in water. For
comparison, the absorption cross section and absorption power of a single gold
nanorod is also plotted along with the aggregates. The nanorod has length 25
nm and diameter 5 nm.

The optical scattering of the nanorod and their aggregates were modelled in
COMSOL® by considering that the nanorods assembly is situated in a spherical
domain of water. At the outer side of the water domain, a perfectly matched
layer (PML) was used to avoid reflection from the water boundary. The incident
unpolarised electromagnetic field was assumed in background with propagation
vector
$\textbf{k}=k_{0}(\mathrm{sin}(\theta)\mathrm{cos}(\phi),\mathrm{sin}(\theta)\mathrm{sin}(\phi),\mathrm{cos}(\theta))$
at $\theta=45^{\circ}$ and $\phi=45^{\circ}$. The standard dispersive
permittivity of the gold was taken from Johnson and Christy Johnson and
Christy (1972), and the permittivity of the surrounding medium, i.e., water,
was kept constant = 1.777. We computed the optical cross-sections of the gold
nanorod and their ordered and random aggregates. Figure 2 shows the computed
optical cross-sections and absorption power of the gold nanorod and their
aggregates. In all ordered aggregates the inter-particle distance is fixed to
be 5 nm and in the random aggregates the inter-particle distance is variable
and the nanorods are not penetrating each other. The strongest LSPR peaks are
listed in Table 1. The longitudinal LSPR peak of a single gold nanorod is at
898 nm. When the nanoparticles are ordered in an aggregate (C1, C2 and C3)
then there is a blueshift in the LSPR due to the lateral coupling within the
nanorods. There is the longitudinal coupling in aggregate C3 but still, there
is a blueshift in LSPR because the lateral coupling dominates over
longitudinal coupling. As we had considered only a fixed inter-particle
distance in the ordered aggregates because it was well known that decreasing
(or increasing) inter-particle distance would increase the coupling,
therefore, there would be more blueshift (or redshift) in the LSPR peaks. In
the case of the random orientation of the nanorods in aggregates (C4, C5 and
C6) there are redshifts in the LSPR peaks. Practically, there would be random
distribution in the aggregates, therefore, it is good that we can tune the
LSPR in the second biological therapeutic window, as shown in Fig. 2. The
total absorbed power by an aggregate is also higher than a single nanorod.

Table 1: The strongest LSPR wavelength, total absorption power and computed maximum temperature change of the gold nanorod aggregates. Sample | Strongest LSPR peak (nm) | Max. abs. power ($\mu$W) | Max. $\Delta$T (K)
---|---|---|---
NR | 898 | 0.86 | 18.4
C1 | 767 | 2.59 | 23.3
C2 | 772 | 5.21 | 35.4
C3 | 883 | 6.53 | 44.4
C4 | 967 | 3.71 | 27.3
C5 | 984 | 3.70 | 26.2
C6 | 1130 | 3.38 | 22.5

To study the resistive heating in the gold nanorod and aggregates we had
considered a cubical domain of water where nanorod or aggregate was situated.
At the top face of the cubical water domain, facing along the z-direction, we
made the optical port for input light. We considered a general wave vector at
$\theta=45^{\circ}$ and $\phi=45^{\circ}$ propagating towards negative
z-direction and the unpolarised electric field was perpendicular to the
incident wave propagation vector. A PML layer at the outer surface of the
water domain was taken. A constant temperature boundary condition ambient
temperature, ($T_{\mathrm{amb}}=300.15$K) was used at the outermost surface of
the computational domain. We used the nanorod or aggregate as a heat source
where the already calculated absorbed power was used as total dissipated heat
energy inside the nanorod or aggregate. We solve the stationary state
condition to get the maximum raised temperature of the system. Density,
specific heat capacity and thermal conductivity of gold were $19.3\times
10^{3}$ kg/m3, 129 J/(kg.K) and 318 W/(m.K), respectively. For the water
density, specific heat capacity and thermal conductivity were 1000 kg/m3,
$4.2\times 10^{3}$ J/(kg.K) and 0.6 W/(m.K) respectively. Figure 3 shows the
surface plot of total power dissipation density and the increased temperature
field ($\Delta T$) at the $x=0$, $y=0$ and $z=0$ planes of a single gold
nanorod in the water domain. We can see that the maximum temperature increase
was 18.4 K for a single gold nanorod.

Figure 3: (a) Surface plot of the total absorption power density and (b)
increases temperature field at the $x=0$, $y=0$ and $z=0$ planes of a single
gold nanorod in water domain.

Figure 4 shows the surface plot of the total absorption power density of the
ordered and random aggregates of the gold nanorods. Here, due to the three-
dimensional distribution of the nanorods, the slice plot (at $x=0$, $y=0$ and
$z=0$ planes) was not visible properly, therefore, we used a surface plot of
the absorbed power density on the nanorods’ assembly. In the surface plot,
there were finite minimum and maximum value of the absorption power density.
Since the direction of the propagation vector was at $\theta=45^{\circ}$ and
$\phi=45^{\circ}$, therefore, in the ordered aggregates a relative gradient
(shadow) of the absorption power could be seen from one nanorod to another. If
propagation would have at $\theta=0^{\circ}$ and $\phi=0^{\circ}$, then there
would have no gradient (shadow) of the absorption power density from particle
to particle. The case of random aggregate is comparatively complex than the
ordered aggregates. So we can’t see proper shadowing effects on the random
aggregates. In the case of aggregates, the absorption power density is higher
than the single gold nanorod (Fig. 3 (a)).

Figure 4: Surface plots of the absorption power density of (a-c) the ordered
aggregates C1, C2 and C3; and (b) random aggregates C4, C5 and C6 of gold
nanorods. Figure 5: Temperature field profile ($\Delta T$) at the $x=0$, $y=0$
and $z=0$ planes of the (a-c) ordered (C1, C2 and C3) and (d-f) random (C4, C5
and C6) gold nanorod aggregates.

The temperature rise in the ordered and random aggregates are shown in Fig. 5.
We could see that aggregates got higher temperature compared to the single
nanorod (monodispersive). In all six aggregates, ordered aggregates have
higher temperature compared to the random aggregates, since each entity in the
ordered assembly participates in the plasmonic coupling, leading to more
temperature rise than the random case. However, in the case of random
assembly, every particle has different coupling, some higher and some lower.
One should not worry about the higher temperature rise in ordered assembly
compared to the random assembly because practically ordered assembly would not
be ideally possible in a solution dispersive form. At the nanoscale,
aggregation of the nanorods is random, which have a higher temperature rise
compared to a single nanorod (monodispersive). Since the LSPR peak for the
random aggregates occurs at a larger wavelength compared to a single
(monodispersive) nanorod, therefore, incident light of higher wavelength would
allow deeper penetration and thus higher temperature rise in the tissue.
Because of LSPR at higher wavelength, random aggregates, therapeutically have
more advantages over monodispersive nanoparticles. The size of aggregation is
not fixed by a few numbers of nanorods or nanoparticles, it could be anything.
However, through proper control over the experimental parameters, the
aggregate size and hence its LSPR peak could be controlled. If we could
control the size of the aggregates then we could decide the light source of
desired higher wavelength for therapeutic use. The current photothermal case
of small gold nanorod aggregate is more oriented toward the in vivo
therapeutic use in tumor domain thermal therapy etc., but it is not limited to
other uses also. For example, the case of plasmonic photothermal application
in antimicrobial uses and acne treatments are not much limited over the size
of aggregation and these applications have more advantages to increase the
higher surface temperature and destroying the bacteria.

Figure 6: Variation of absorption characteristic and temperature rise of gold
nanorod aggregate C6 for various incident and polarization angles of the
light. (a) Absorption power of the aggregate with incident wavelength at
different incident angle $\theta$ for polarization angle $\phi=45^{\circ}$,
(b) the strongest LSPR peaks, (c) maximum absorption power at the strongest
LSPR peak and (d) maximum temperature change $\Delta T$ (K) at illumination
with strongest LSPR wavelength with respect to the indecent angle at different
polarization angle.

For most of the practical purpose, the incident light would be randomly
polarized and aggregates would be randomly oriented. The exact experimental
condition could not be obtained in the numerical studies. However, we carried
out the angle-dependent study of the light absorption and heat generation by
the random aggregate system C6 only. The angle $\theta$ is the inclination or
incident angle from the $z$-axis of the light wavevector and $\phi$ is the
polarization angle from the $x$-axis of the electric field. We fixed the
polarization angle $0^{\circ}$, $45^{\circ}$ and $90^{\circ}$ and varied the
incident angle from $0^{\circ}$ to $75^{\circ}$ in the interval of
$15^{\circ}$ and computed the light absorption and heat generation. Figure 6
shows the results of the incident light’s angle-dependent study for the random
aggregate C6. Figure 6(a) presents the total absorption power of the aggregate
with respect to the incident wavelength at a fixed polarization angle
$45^{\circ}$ and at different incident angles. We can see that as the
wavevector inclination angle $\theta$ increases, there is a slight monotonic
blueshift in the strongest LSPR peak. This happens due to the reason that the
nanorods in the aggregates are randomly oriented and at a different angle
there would be variation in the plasmonic coupling within the particles. For
the monodispersive nanorods, the blueshift could not happen by varying the
incident and polarization angles but only there would be a change in the
absorption amplitude. Variation in the absorption peaks for aggregates with
incident and polarization angle implies that in the real experimental system
there would be a broadening in the LSPR peaks. Figure 6(b) shows the strongest
LSPR peaks with respect to the incident angle at fixed polarization
$0^{\circ}$, $45^{\circ}$ and $90^{\circ}$. Here we see how the strongest LSPR
peak position varies, indicating the broadening of the peaks in experimental
conditions. Figure 6(c) and Figure 6(d) show the maximum absorption power and
corresponding maximum temperature change at the strongest LSPR peaks with
respect to the incident angle at different polarization angle $0^{\circ}$,
$45^{\circ}$ and $90^{\circ}$ respectively. We could notice that the
polarization did not change absorption power and temperature rise much due to
the random nature of the aggregate. One should note that all these computed
observations are just for one particular aggregate configuration, however,
these results would naturally vary for other configurations of the aggregates.
By fixing the experimental conditions one could get the well defined LSPR peak
for aggregates suspension in the solution.

##  3 Conclusion

We conclude that aggregation of the small gold nanorods causes a blueshift of
LSPR wavelength in ordered aggregates and a redshift of the LSPR wavelength in
the random aggregates. Even if there was either a blueshift or redshift in the
LSPR peak of any aggregate, yet, there was a higher temperature rise in that
aggregate as compared to the single nanorod (monodispersive). The angle-
dependent study in the random aggregates shows that there is a slight
monotonic blueshift in the LSPR peak which would cause broadening of the LSPR
peak of aggregates suspension in an actual experiment. The size and structure
of the aggregation would depend on the experimental or local environmental or
chemical conditions with optimization of the experimental conditions, the
desired type of a stable aggregate could be obtained. We could say that the
aggregation of the gold nanorods or any shape of metallic nanoparticles would
have clear advantages also over monodispersive form, at least wherever there
would be a limitation of the size of the nanoparticle and need of higher
wavelength plasmonic coupling and greater temperature rise. For example, in
the case of plasmonic photothermal therapy of tumor tissue to get higher light
penetration depth, antimicrobial applications, and plasmonic photothermal
treatment of acne, such aggregates will produce more heat than their
monodispersive counterparts.

## Acknowledgements

D.P. acknowledges the CSIR-India for the Nehru Science Postdoctoral Research
Fellowship number HRDG/CSIR-Nehru PDF/EN, ES & PS/EMR-I/04/2019. D.P. and S.S.
acknowledge the support of the CSIR-CSIO Chandigarh for hosting the research.

##  Appendix A Aggregate Projections

All considered six aggregates, ordered as well as random, are three-
dimensional entities. The geometric structure projections in the $xy$, $yz$
and $zx$ planes are shown in the Fig. A.1.

Figure A.1: Projected views of aggregates C1, C2, C3, C4, C5 and C6 in the
$xy$, $yz$ and $zx$ planes used for computation. Each nanorod has length 25 nm
and diameter 5 nm.

##  Appendix B Model validation

Heat generation study of a 100 nm gold nanosphere situated in water was
carried out as a benchmark problem. The permittivity of the gold was from
Johnson and Christy Johnson and Christy (1972). All values of other parameters
had been kept identical as used by G. Baffou et al. Baffou et al. (2010).
Figure B.1 shows the computation of the normalised electric field intensity,
energy density, temperature evolved in the nanosphere at incident wavelength
530 nm.. The total absorption power of the nanosphere from 400 nm to 800 nm is
also shown. The irradiance of the light was 1 mW/$\mu$m2. In our case, the
maximum raised temperature was 52.2 ∘C which was 55 ∘C in Baffou case and 52.5
∘C in X. Chen et al. Chen et al. (2012). The maximum absorbed power was
obtained 28.27 $\mu$W at 543 nm while the maximum absorbed power reported by
Baffou was 20.5 $\mu$W at 530 nm and by Chen 21.45 $\mu$W at 540 nm.

Figure B.1: Comparison validation of algorithm (a) normalised electric field
intensity, (b) absorbed power density, (c) temperature field and (d) total
absorbed power of a gold nanosphere of 100 nm diameter at incident wavelength
530 nm and irradiance 1 mW/$\mu$m2.

## References

  * Nam et al. [2009] Jutaek Nam, Nayoun Won, Ho Jin, Hyokyun Chung, and Sungjee Kim.  ph-induced aggregation of gold nanoparticles for photothermal cancer therapy.  _Journal of the American Chemical Society_ , 131(38):13639–13645, 2009.
  * Choi et al. [2011] Won Il Choi, Ja-Young Kim, Chul Kang, Clare C Byeon, Young Ha Kim, and Giyoong Tae.  Tumor regression in vivo by photothermal therapy based on gold-nanorod-loaded, functional nanocarriers.  _ACS nano_ , 5(3):1995–2003, 2011.
  * Pattani et al. [2015] Varun P Pattani, Jay Shah, Alexandra Atalis, Anirudh Sharma, and James W Tunnell.  Role of apoptosis and necrosis in cell death induced by nanoparticle-mediated photothermal therapy.  _Journal of Nanoparticle Research_ , 17(1):1–11, 2015.
  * Li et al. [2008] Jing Liang Li, Daniel Day, and Min Gu.  Ultra-low energy threshold for cancer photothermal therapy using transferrin-conjugated gold nanorods.  _Advanced materials_ , 20(20):3866–3871, 2008\.
  * Wang et al. [2013] Yucai Wang, Kvar CL Black, Hannah Luehmann, Weiyang Li, Yu Zhang, Xin Cai, Dehui Wan, Si-Yun Liu, Max Li, Paul Kim, et al.  Comparison study of gold nanohexapods, nanorods, and nanocages for photothermal cancer treatment.  _ACS nano_ , 7(3):2068–2077, 2013.
  * Adeyemi and Sulaiman [2015] Oluyomi S Adeyemi and Faoziyat A Sulaiman.  Evaluation of metal nanoparticles for drug delivery systems.  _Journal of biomedical research_ , 29(2):145, 2015\.
  * Janib et al. [2010] Siti M Janib, Ara S Moses, and J Andrew MacKay.  Imaging and drug delivery using theranostic nanoparticles.  _Advanced drug delivery reviews_ , 62(11):1052–1063, 2010.
  * Babu et al. [2013] Anish Babu, Amanda K Templeton, Anupama Munshi, and Rajagopal Ramesh.  Nanoparticle-based drug delivery for therapy of lung cancer: progress and challenges.  _Journal of Nanomaterials_ , 2013, 2013.
  * Lin et al. [2010] Hsing-Ying Lin, Chen-Han Huang, Chih-Han Chang, Yun-Chiang Lan, and Hsiang-Chen Chui.  Direct near-field optical imaging of plasmonic resonances in metal nanoparticle pairs.  _Optics express_ , 18(1):165–172, 2010.
  * Wiederrecht [2004] GP Wiederrecht.  Near-field optical imaging of noble metal nanoparticles.  _The European Physical Journal Applied Physics_ , 28(1):3–18, 2004.
  * Syed et al. [2017] Abdullah Muhammad Syed, Shrey Sindhwani, Stefan Wilhelm, Benjamin R Kingston, Dennis SW Lee, Jennifer L Gommerman, and Warren CW Chan.  Three-dimensional imaging of transparent tissues via metal nanoparticle labeling.  _Journal of the American Chemical Society_ , 139(29):9961–9971, 2017.
  * Wang et al. [2020] Hui Wang, Hui Yu, Yan Wang, Xiaonan Shan, Hong-Yuan Chen, and Nongjian Tao.  Phase imaging of transition from classical to quantum plasmonic couplings between a metal nanoparticle and a metal surface.  _Proceedings of the National Academy of Sciences_ , 117(30):17564–17570, 2020.
  * Donner et al. [2011] Jon S Donner, Guillaume Baffou, David McCloskey, and Romain Quidant.  Plasmon-assisted optofluidics.  _Acs nano_ , 5(7):5457–5462, 2011.
  * Gómez-Villarejo et al. [2017] Roberto Gómez-Villarejo, Elisa I Martín, Javier Navas, Antonio Sánchez-Coronilla, Teresa Aguilar, Juan Jesús Gallardo, Rodrigo Alcántara, Desiré De los Santos, Iván Carrillo-Berdugo, and Concha Fernández-Lorenzo.  Ag-based nanofluidic system to enhance heat transfer fluids for concentrating solar power: Nano-level insights.  _Applied Energy_ , 194:19–29, 2017.
  * Dionne and Atwater [2012] Jennifer A Dionne and Harry A Atwater.  Plasmonics: metal-worthy methods and materials in nanophotonics.  _Mrs Bulletin_ , 37(8):717–724, 2012.
  * Fang and Sun [2015] Yurui Fang and Mengtao Sun.  Nanoplasmonic waveguides: towards applications in integrated nanophotonic circuits.  _Light: Science & Applications_, 4(6):e294–e294, 2015.
  * Tang et al. [2014] Shurong Tang, Ping Tong, Wei Lu, Jinfeng Chen, Zhiming Yan, and Lan Zhang.  A novel label-free electrochemical sensor for hg2+ based on the catalytic formation of metal nanoparticle.  _Biosensors and Bioelectronics_ , 59:1–5, 2014.
  * Chen et al. [2019] Minrui Chen, Weifeng Luo, Zhongqi Xu, Xueping Zhang, Bo Xie, Guanghou Wang, and Min Han.  An ultrahigh resolution pressure sensor based on percolative metal nanoparticle arrays.  _Nature communications_ , 10(1):1–9, 2019.
  * Lee et al. [2014] Jaehwan Lee, Sanghyeok Kim, Jinjae Lee, Daejong Yang, Byong Chon Park, Seunghwa Ryu, and Inkyu Park.  A stretchable strain sensor based on a metal nanoparticle thin film for human motion detection.  _Nanoscale_ , 6(20):11932–11939, 2014.
  * Roldan Cuenya [2013] Beatriz Roldan Cuenya.  Metal nanoparticle catalysts beginning to shape-up.  _Accounts of chemical research_ , 46(8):1682–1691, 2013.
  * Scholten et al. [2012] Jackson D Scholten, Barbara Caroline Leal, and Jairton Dupont.  Transition metal nanoparticle catalysis in ionic liquids.  _Acs Catalysis_ , 2(1):184–200, 2012.
  * Yan et al. [2010] Ning Yan, Chaoxian Xiao, and Yuan Kou.  Transition metal nanoparticle catalysis in green solvents.  _Coordination Chemistry Reviews_ , 254(9-10):1179–1218, 2010.
  * Soni et al. [2015] Sanjeev Soni, Himanshu Tyagi, Robert A Taylor, and Amod Kumar.  Experimental and numerical investigation of heat confinement during nanoparticle-assisted thermal therapy.  _International Communications in Heat and Mass Transfer_ , 69:11–17, 2015.
  * Chen et al. [2017] Zhimin Chen, Haihua Fan, Jinxiang Li, Shaolong Tie, and Sheng Lan.  Photothermal therapy of single cancer cells mediated by naturally created gold nanorod clusters.  _Optics express_ , 25(13):15093–15107, 2017.
  * Larsson et al. [2007] Elin M Larsson, Joan Alegret, Mikael Käll, and Duncan S Sutherland.  Sensing characteristics of nir localized surface plasmon resonances in gold nanorings for application as ultrasensitive biosensors.  _Nano letters_ , 7(5):1256–1263, 2007.
  * Skrabalak et al. [2008] Sara E Skrabalak, Jingyi Chen, Yugang Sun, Xianmao Lu, Leslie Au, Claire M Cobley, and Younan Xia.  Gold nanocages: synthesis, properties, and applications.  _Accounts of chemical research_ , 41(12):1587–1595, 2008.
  * Fan et al. [2008] Feng-Ru Fan, De-Yu Liu, Yuan-Fei Wu, Sai Duan, Zhao-Xiong Xie, Zhi-Yuan Jiang, and Zhong-Qun Tian.  Epitaxial growth of heterogeneous metal nanocrystals: from gold nano-octahedra to palladium and silver nanocubes.  _Journal of the American Chemical Society_ , 130(22):6949–6951, 2008.
  * Apyari et al. [2018] Vladimir V Apyari, Marina O Gorbunova, Anastasiya V Shevchenko, Aleksei A Furletov, Pavel A Volkov, Alexey V Garshev, Stanislava G Dmitrienko, and Yury A Zolotov.  Towards highly selective detection using metal nanoparticles: a case of silver triangular nanoplates and chlorine.  _Talanta_ , 176:406–411, 2018.
  * Dutta et al. [2019] Susom Dutta, Guinevere Strack, and Pradeep Kurup.  Gold nanostar electrodes for heavy metal detection.  _Sensors and Actuators B: Chemical_ , 281:383–391, 2019\.
  * Wang et al. [2019] Hongzhi Wang, Yuying Gao, Jia Liu, Xinyuan Li, Muwei Ji, Erhuan Zhang, Xiaoyan Cheng, Meng Xu, Jiajia Liu, Hongpan Rong, et al.  Efficient plasmonic au/cdse nanodumbbell for photoelectrochemical hydrogen generation beyond visible region.  _Advanced Energy Materials_ , 9(15):1803889, 2019\.
  * Ruparelia et al. [2008] Jayesh P Ruparelia, Arup Kumar Chatterjee, Siddhartha P Duttagupta, and Suparna Mukherji.  Strain specificity in antimicrobial activity of silver and copper nanoparticles.  _Acta biomaterialia_ , 4(3):707–716, 2008.
  * Aioub et al. [2017] Mena Aioub, Sajanlal R Panikkanvalappil, and Mostafa A El-Sayed.  Platinum-coated gold nanorods: efficient reactive oxygen scavengers that prevent oxidative damage toward healthy, untreated cells during plasmonic photothermal therapy.  _ACS nano_ , 11(1):579–586, 2017.
  * Sun et al. [2016] Mengmeng Sun, Fei Liu, Yukun Zhu, Wansheng Wang, Jin Hu, Jing Liu, Zhifei Dai, Kun Wang, Yen Wei, Jing Bai, et al.  Salt-induced aggregation of gold nanoparticles for photoacoustic imaging and photothermal therapy of cancer.  _Nanoscale_ , 8(8):4452–4457, 2016.
  * Pan et al. [2017] Limin Pan, Jianan Liu, and Jianlin Shi.  Nuclear-targeting gold nanorods for extremely low nir activated photothermal therapy.  _ACS applied materials & interfaces_, 9(19):15952–15961, 2017.
  * Ungureanu et al. [2011] Constantin Ungureanu, Rene Kroes, Wilma Petersen, Tom AM Groothuis, Felicia Ungureanu, Hans Janssen, Fijs WB van Leeuwen, Rob PH Kooyman, Srirang Manohar, and Ton G van Leeuwen.  Light interactions with gold nanorods and cells: implications for photothermal nanotherapeutics.  _Nano letters_ , 11(5):1887–1894, 2011.
  * Shen et al. [2013] Shun Shen, Hongyan Tang, Xiaotong Zhang, Jinfeng Ren, Zhiqing Pang, Dangge Wang, Huile Gao, Yong Qian, Xinguo Jiang, and Wuli Yang.  Targeting mesoporous silica-encapsulated gold nanorods for chemo-photothermal therapy with near-infrared radiation.  _Biomaterials_ , 34(12):3150–3158, 2013.
  * Xie et al. [2019] Bowei Xie, Lanxin Ma, Junming Zhao, and Linhua Liu.  Dependent absorption property of nanoparticle clusters: an investigation of the competing effects in the near field.  _Optics express_ , 27(8):A280–A291, 2019.
  * Gu et al. [2018] Xi Gu, Victoria Timchenko, Guan Heng Yeoh, Leonid Dombrovsky, and Robert Taylor.  The effect of gold nanorods clustering on near-infrared radiation absorption.  _Applied Sciences_ , 8(7):1132, 2018.
  * Baffou et al. [2010] Guillaume Baffou, Romain Quidant, and F Javier García de Abajo.  Nanoscale control of optical heating in complex plasmonic systems.  _ACS nano_ , 4(2):709–716, 2010.
  * Lermé et al. [1998] J Lermé, B Palpant, B Prével, E Cottancin, M Pellarin, M Treilleux, JL Vialle, A Perez, and M Broyer.  Optical properties of gold metal clusters: A time-dependent local-density-approximation investigation.  _The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics_ , 4(1):95–108, 1998.
  * Cottancin et al. [2006] E Cottancin, G Celep, J Lermé, M Pellarin, JR Huntzinger, JL Vialle, and M Broyer.  Optical properties of noble metal clusters as a function of the size: comparison between experiments and a semi-quantal theory.  _Theoretical Chemistry Accounts_ , 116(4):514–523, 2006.
  * Kreibig and Vollmer [2013] Uwe Kreibig and Michael Vollmer.  _Optical properties of metal clusters_ , volume 25.  Springer Science & Business Media, 2013.
  * Maier [2007] Stefan Alexander Maier.  _Plasmonics: fundamentals and applications_.  Springer Science & Business Media, 2007.
  * Baffou and Quidant [2013] Guillaume Baffou and Romain Quidant.  Thermo-plasmonics: using metallic nanostructures as nano-sources of heat.  _Laser & Photonics Reviews_, 7(2):171–187, 2013\.
  * Baffou [2017] Guillaume Baffou.  Thermoplasmonics.  _Heating Metal Nanoparticles Using Light_ , 2017.
  * Mishchenko et al. [2002] Michael I Mishchenko, Larry D Travis, and Andrew A Lacis.  _Scattering, absorption, and emission of light by small particles_.  Cambridge university press, 2002.
  * De Abajo and Howie [2002] FJ Garcia De Abajo and A Howie.  Retarded field calculation of electron energy loss in inhomogeneous dielectrics.  _Physical Review B_ , 65(11):115418, 2002.
  * Martin et al. [1995] Olivier JF Martin, Christian Girard, and Alain Dereux.  Generalized field propagator for electromagnetic scattering and light confinement.  _Physical Review Letters_ , 74(4):526, 1995.
  * Girard et al. [2008] Ch Girard, E Dujardin, G Baffou, and R Quidant.  Shaping and manipulation of light fields with bottom-up plasmonic structures.  _New Journal of Physics_ , 10(10):105016, 2008\.
  * Girard [2005] Christian Girard.  Near fields in nanostructures.  _Reports on progress in physics_ , 68(8):1883, 2005.
  * Yurkin and Hoekstra [2007] Maxim A Yurkin and Alfons G Hoekstra.  The discrete dipole approximation: an overview and recent developments.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ , 106(1-3):558–589, 2007.
  * Myroshnychenko et al. [2008] Viktor Myroshnychenko, Jessica Rodríguez-Fernández, Isabel Pastoriza-Santos, Alison M Funston, Carolina Novo, Paul Mulvaney, Luis M Liz-Marzán, and F Javier García De Abajo.  Modelling the optical response of gold nanoparticles.  _Chemical Society Reviews_ , 37(9):1792–1805, 2008.
  * de Abajo [1999] FJ García de Abajo.  Interaction of radiation and fast electrons with clusters of dielectrics: a multiple scattering approach.  _Physical review letters_ , 82(13):2776, 1999\.
  * Radhakrishnan and Murugesan [2014] Archana Radhakrishnan and Dr V Murugesan.  Calculation of the extinction cross section and lifetime of a gold nanoparticle using fdtd simulations.  In _AIP Conference Proceedings_ , volume 1620, pages 52–57. American Institute of Physics, 2014.
  * Mohsin and Salim [2018] Abu SM Mohsin and Mariam B Salim.  Probing the plasmon coupling, quantum yield, and effects of tip geometry of gold nanoparticle using analytical models and fdtd simulation.  _IEEE Photonics Journal_ , 10(3):1–10, 2018.
  * Bachelier et al. [2008] G Bachelier, I Russier-Antoine, E Benichou, C Jonin, N Del Fatti, F Vallée, and P-F Brevet.  Fano profiles induced by near-field coupling in heterogeneous dimers of gold and silver nanoparticles.  _Physical review letters_ , 101(19):197401, 2008\.
  * Alizadeh et al. [2019] M Alizadeh, Kh Hosseinzadeh, and DD Ganji.  Investigating the effects of hybrid nanoparticles on solid-liquid phase change process in a y-shaped fin-assisted lhtess by means of fem.  _Journal of Molecular Liquids_ , 287:110931, 2019.
  * Multiphysics [1998] COMSOL Multiphysics.  Introduction to comsol multiphysics®.  _COMSOL Multiphysics, Burlington, MA, accessed Feb_ , 9:2018, 1998.
  * Chen et al. [2012] Xi Chen, Yiting Chen, Min Yan, and Min Qiu.  Nanosecond photothermal effects in plasmonic nanostructures.  _ACS nano_ , 6(3):2550–2557, 2012.
  * Johnson and Christy [1972] Peter B Johnson and R-WJPrB Christy.  Optical constants of the noble metals.  _Physical review B_ , 6(12):4370, 1972.