Studies on optimizing potential energy functions for maximal intrinsic
hyperpolarizability
Juefei Zhou, Urszula B. Szafruga, David S. Watkins* and Mark G. Kuzyk
Department of Physics and Astronomy, Washington State University,
Pullman, Washington 99164-2814; and *Department of Mathematics
We use numerical optimization to study the properties of (1) the class of one-dimensional po-
tential energy functions and (2) systems of point charges in two-dimensions that yield the largest
hyperpolarizabilities, which we find to be within 30% of the fundamental limit. We investigate the
character of the potential energy functions and resulting wavefunctions and find that a broad range
of potentials yield the same intrinsic hyperpolarizability ceiling of 0.709.
I. INTRODUCTION
Materials with large nonlinear-optical suscep-
tibilities are central for optical applications such
as telecommunications,[1] three-dimensional nano-
photolithography,[2, 3] and making new materials[4] for
novel cancer therapies.[5] The fact that quantum calcu-
lations show that there is a limit to the nonlinear-optical
response[6, 7, 8, 9, 10, 11] is both interesting from the
basic science perspective; and, provides a target for
making optimized materials. In this work, we focus
on the second-order susceptibility and the underlying
molecular hyperpolarizability, which is the basis of
electro-optic switches and frequency doublers.
The fundamental limit of the off-resonance hyperpo-
larizability is given by,[8]
βMAX =
, (1)
where N is the number of electrons and E10 the energy
difference between the first excited state and the ground
state, E10 = E1 − E0. Using Equation 1, we can de-
fine the off-resonant intrinsic hyperpolarizability, βint, as
the ratio of the actual hyperpolarizability (measured or
calculated), β, to the fundamental limit,
βint = β/βMAX . (2)
We note that since the dispersion of the fundamental
limit of β is also known, [12] it is possible to calculate
the intrinsic hyperpolarizability at any wavelength. In
the present work, we treat only the zero-frequency limit.
Until recently, the largest nonlinear susceptibilities of
the best molecules fell short of the fundamental limit by a
factor of 103/2, [10, 13, 14] so the very best molecules had
a value of βint = 0.03. Since a Sum-Over-States (SOS)
calculation of the hyperpolarizability[15] using the ana-
lytical wavefunctions of the clipped harmonic oscillator
yields a value βint = 0.57,[14] the factor-of-thirty gap is
not of a fundamental nature. Indeed, recently, it was re-
ported that a new molecule with asymmetric conjugation
modulation has a measured value of βint = 0.048.[16]
To investigate how one might make molecules with
a larger intrinsic hyperpolarizability, Zhou and cowork-
ers used a numerical optimization process where a trial
potential energy function is entered as an input, and
the code iteratively deforms the potential energy func-
tion until the intrinsic hyperpolarizability, calculated
from the resulting wavefunctions, converges to a local
maximum.[17] In this work, a hyperbolic tangent func-
tion was used as the starting potential due to the fact that
it is both asymmetric yet relatively flat away from the ori-
gin. This calculation was one-dimensional and included
only one electron, so electron correlation effects were
ignored. Furthermore, the intrinsic hyperpolarizability
was calculated using the new dipole-free sum-over-states
expression[18] and only 15 excited states were included.
The resulting optimized potential energy function showed
strong oscillations, which served to separate the spatial
overlap between the energy eigenfunctions. This led Zhou
and coworkers to propose that modulated conjugation
in the bridge between donor and acceptor ends of such
molecules may be a new paradigm for making molecules
with higher intrinsic hyperpolarizability.[17]
Based on this paradigm, Pérez Moreno reported mea-
surements of a class of chromophores with varying de-
gree of modulated conjugation.[16] The best measured
intrinsic hyperpolarizability was βint = 0.048, about 50%
larger than the best previously-reported. Given the mod-
est degree of conjugation modulation for this molecule,
this new paradigm shows promise for further improve-
ments.
In the present work, we extend Zhou’s calculations to
a larger set of starting potentials. To circumvent trun-
cation problems associated with sum-over-states calcula-
tions, we instead determine the hyperpolarizability using
a finite difference technique. The optimization procedure
is then applied to this non-perturbative hyperpolarizabil-
To study the effects of geometry on the hyperpo-
larizability, Kuzyk and Watkins calculated the hyper-
polarizability of various arrangements of point charges,
representing nuclei, in two-dimensions using a two-
dimensional Coulomb potential.[19] In the present con-
tribution, we apply our numerical optimization technique
to determine the arrangement and charges of the nuclei
in a planar molecule that maximize the intrinsic hyper-
polarizability.
http://arxiv.org/abs/0704.1687v1
II. THEORY
In our previous work, we used a finite-state SOS model
of the hyperpolarizability that derives from perturbation
theory (we used both the standard Orr andWard SOS ex-
pression, βSOS ,[15] and the newer dipole free expression,
βDF [18]). The use of a finite number of states in lieu of
the full infinite sums can result in inaccuracies, so, in the
present work, we use the non-perturbative approach, as
follows. We begin by solving the 1-d Schrodinger Equa-
tion on the interval a < x < b for the ground state wave-
function ψ(x,E) of an electron in a potential well defined
by V (x) and in the presence of a static electric field, E,
that adds to the potential δV = −exE. From this, the
off-resonant hyperpolarizability is calculated with numer-
ical differentiation, i.e. using finite differences, yielding
βNP =
|ψ(x,E)|2 ex dx
. (3)
Equation 3 is evaluated using the standard second-order
approximation to the second derivative:
f ′′(z) ≈ f(z + h)− 2f(z) + f(z − h)
with several h values h0, h0/5, h0/25, . . . . We then refine
these values by Richardson extrapolation [20] and obtain
our estimate from the two closest extrapolated values.
Our computational mesh consists of 200 quadratic fi-
nite elements with a total of 399 degrees of freedom. The
potential energy function is a cubic spline with 40 degrees
of freedom. Thus the numerical calculations in regions
where the potential function is represented by 3 points
in the spline are covered by 15 elements with a total of
about 30 degrees of freedom.
Calculating βint from Equations 3, 2 and 1 for a specific
potential, we use an optimization algorithm that contin-
uously varies the potential in a way that maximizes βint.
We also compute the matrix[17, 26]
τ (N)mp = δm,p −
xmax10
xmax10
, (4)
where xmax10 is the magnitude of the fundamental limit of
the position matrix element x10 for a one electron system,
and is given by,
xmax10 =
2mE10
. (5)
Each matrix element of τ (N), indexed by m and p, is
a measure of how well the (m, p) sum rule is obeyed
when truncated to N states. If the sum rules are ex-
actly obeyed, τ
mp = 0 for all m and p. We note that if
the sum rules are truncated to an N-state model, the sum
rules indexed by a large value of m or p (i.e. m, p ∼ N)
are disobeyed even when the position matrix elements
and energies are exact. We have found that the values
mp are small for exact wavefunctions when m < N/2
and p < N/2. So, when evaluating the τ matrix to
test our calculations, we consider only the components
m≤N/2,p≤N/2
We observe that when using more than about 40 states
in SOS calculations of the hyperpolarizability only a
marginal increase of accuracy results when the poten-
tial energy function is parameterized with 400 degrees of
freedom. So, to ensure overkill, we use 80 states when
calculating the τ matrix or the hyperpolarizability with
an SOS expression so that truncation errors are kept to
a minimum. Since the hyperpolarizability depends crit-
ically on the transition dipole moment from the ground
state to the excited states, we use the value of τ
00 as
one important test of the accuracy of the calculated wave-
functions. Additionally, we use the standard deviation of
τ (N),
∆τ (N) =
, (6)
which quantifies, on average, how well the sum rules are
obeyed in aggregate, making ∆τ (N) a broader test of the
accuracy of a large set of wavefunctions.
Our code is written in MATLAB. For each trial po-
tential we use a quadratic finite element method [21] to
approximate the Schrödinger eigenvalue problem and the
implicitly restarted Arnoldi method [22] to compute the
wave functions and energy levels. To optimize β we use
the Nelder-Mead simplex algorithm [23].
As described in our previous work,[17] we perform op-
timization, but this time using the exact intrinsic hyper-
polarizability β = βNP /βMAX , where βMAX is the fun-
damental limit of the hyperpolarizability, which is pro-
portional to E
10 . To determine E10 ≡ E1 −E0, we also
calculate the first excited state energy E1.
III. RESULTS AND DISCUSSIONS
Figure 1 shows an example of the optimized poten-
tial energy function after 7,000 iterations when starting
with the potential V (x) = 0 and optimizing the non-
perturbative intrinsic hyperpolarizability βNP /βMAX as
calculated with Equation 3. Also shown are the eigen-
functions of the first 15 states computed from the opti-
mized potential. First, we note that the potential energy
function shows the same kinds of wiggles as in our original
paper,[17] though not of sufficient amplitude to localize
the wavefunctions.
For the starting potentials we have investigated, our
results fall into two broad classes. In the first, three
common features are: (1) The best intrinsic hyperpolar-
izabilities are near βint = 0.71; (2) the best potentials
have a series of wiggles; and (3) the systems behave as a
0 5 10 15 20
FIG. 1: Optimized potential energy function and first 15
wavefunctions after 7,000 iterations. Starting potential is
V (x) = 0, using the non-perturbative hyperpolarizability for
optimization.
limited-state model. In the second class of starting po-
tentials, (2) the wiggles are much less pronounced and (3)
more states contribute evenly. Figure 1 is an example of
a Class II potential. However, in both classes, the max-
imum calculated intrinsic hyperpolarizability appears to
be bounded by βint = 0.71. Using the set of potentials
from both classes that lead to optimized βNP /βMAX , we
calculate the lowest 80 eigenfunctions and eigenvalues,
from which we calculate βDF and βSOS . In most cases,
we find that the three different formulas for β converge
to the same value when only the first 20 excited states
are used (using 80 states, the three are often the same to
within at least 4 decimal places) and τ00 ≈ 10−4, showing
that the ground state sum rules are well obeyed. Further-
more, the rms deviation of the τ matrix when including
40 states leads to τ (80) < 0.001.
Figure 2 shows an example of the optimized potential
energy function when starting with the potential V (x) =
tanhx and optimizing the exact (non-perturbative) in-
trinsic hyperpolarizability. Also shown are the eigenfunc-
tions of the first 15 states computed with the optimized
potential. First, we note that the potential energy func-
tion shows the same kinds of wiggles as in our original
paper;[17] and only 2 excited state wavefunctions and the
ground state are localized in the first deep well - placing
this system in Class I.
The observation that such potentials lead to hyper-
polarizabilities that are near the fundamental limit mo-
tivated Zhou and coworkers to suggest that molecules
with modulated conjugation may have enhanced intrin-
sic hyperpolarizabilities.[17] A molecule with a modu-
lated conjugations bridge between the donor and accep-
tor end was later shown to have record-high intrinsic
hyperpolarizability.[16] As such, this result warrants a
more careful analysis.
It is worthwhile to compare our present results charac-
0 5 10 15 20
FIG. 2: Optimized potential energy function and first 15
wavefunctions after 8,000 iterations. Starting potential is
V (x) = tanh(x), using the non-perturbative hyperpolarizabil-
ity for optimization.
terized by Figure 2 with our past work,[17] particularly
for the purpose of examining the impact of the approx-
imations used in the previous work.[17] Figure 3 shows
the optimized potential and wavefunctions obtained by
Zhou and coworkers using a 15-state model and opti-
mizing the dipole-free intrinsic hyperpolarizability. Since
only 15 states were used, the SOS expression for β did not
fully converge; making the result inaccurate as suggested
by the fact that βSOS and βDF did not agree. How-
ever, since the code focused on optimizing the dipole-free
form of β, and τ00 was small when βint was optimized,
the dipole-free expression may have converged to a rea-
sonably accurate value while the commonly-used SOS
expression was inaccurate. Indeed, it was found that
βDF ≈ 0.72 - in contrast to our more precise present
calculations using the non-perturbative approach, which
yields βNP < 0.71. So, the fact that our more precise
calculations, which do not rely on a sum-over states ex-
pression, agree so well with the 15-state model suggests
that in both cases, the limit for a one-dimensional single
electron molecule is just over β ≈ 0.7. This brute force
calculation serves as a numerical illustration of the obser-
vation that the limiting value of β is the same for an exact
non-perturbation calculation and for a calculation that
truncates the SOS expression, which presumedly should
lead to large inaccuracies.[24, 25] At minimum, this result
supports the existence of fundamental limits of nonlinear
susceptibilities that are in line with past calculations.
To state Zhou’s approach more precisely,[17] the cal-
culations optimized the very special case of the intrin-
sic hyperpolarizability for a 15 state model for a poten-
tial energy function that is parameterized with 20 spline
points. As such, the potential energy function can at
most develop about 20 wiggles. As a consequence, there
are enough degrees of freedom in the potential energy
function to force the 15 states to be spatially well sepa-
0 5 10 15 20
 V(x)
FIG. 3: Optimized potential energy function using βDF and
first 15 wavefunctions after 7,000 iterations. Starting po-
tential is the tanh(x) potential. The final potential (shown
above) we refer to as the Zhou potential.
TABLE I: Evolution of Zhou’s Potential. βs is the hyperpo-
larizability of the starting potential using 80 states while the
other ones are after optimization of βNP .
Number of βS βSOS βDF βNP τ
00 ∆τ
Iterations (×10−5) (×10−4)
0 0.5612 0.5612 0.5607 0.5612 11.2 15
1000 0.5612 0.7087 0.6682 0.7083 1810 40
rated. Interestingly, after optimization, only two excited
states overlap with ground state, allowing only these two
states to have nonzero transition dipole moments with
each other and the ground state – forcing the system
into a three-level SOS model for βDF . This behavior is
interesting in light of the three-level ansatz, which asserts
that only three states determine the nonlinear response
of a system when it is near the fundamental limits.
It is interesting to compare the exact non-perturbation
calculation, which does not depend on the excited state
wavefunctions (Figure 2) and Zhou’s contrived system of
15 states (Figure 3). Both cases have wiggles and the
wavefunctions appear to be mostly non-overlapping. So,
for the first 15 states, the wavefunctions appear simi-
larly localized. The situation becomes more interesting
when 80 states are included in calculating the hyperpo-
larizability for Zhou’s potential or when the exact non-
perturbative approach is used. The first line in Table
I summarizes the results with Zhou’s potential and 80
states.
First, let’s focus on the sum-over-states results.
Clearly, when 80 states are used in the calculation, it
is impossible for the excited state wavefunctions to not
overlap with each other, so the three-level approximation
to β breaks down. According to the three-level ansatz,
we would expect the hyperpolarizability to get smaller.
Indeed, the additional excited states result in a smaller
0 5 10 15 20
FIG. 4: Optimized potential energy function and first 15
wavefunctions after 1,000 iterations. Starting potential is
Zhou’s potential, using the non-perturbative hyperpolariz-
ability for optimization.
hyperpolarizability (≈ 0.56). Note that the exact and
SOS expressions agree with each other and that τ
and ∆τ (80) are small.
Figure 4 shows the result after 1000 iterations, us-
ing Zhou’s potential as the starting potential and us-
ing the non-perturbative hyperpolarizability for opti-
mization. First, the non-perturbative hyperpolarizabil-
ity reaches just under 0.71, but, the SOS and dipole-
free expressions do not agree with each other. Further-
more, both convergence metrics (τ
00 and ∆τ
(80)) are
larger than before optimization. It would appear that
for Zhou’s potential, even 80 states are not sufficient
to characterize the nonlinear susceptibility when a sum-
over-states expression is used (either dipole free or tradi-
tional SOS expression - though the SOS expression agrees
better with the non-perturbative approach).
Interestingly, the optimized potential energy function
still retains the wiggles and the wave functions are still
well separated. This result is consistent with the sug-
gestion of Zhou and coworkers that modulation of conju-
gation may be a good design strategy for making large-
hyperpolarizability molecules. We note that wiggles in
the potential energy function are not required to get
a large nonlinear-optical response; but, appears to be
one way that Mother Nature optimizes the hyperpo-
larizability. Since this idea has been used to identify
molecules with experimentally measured record intrin-
sic hyperpolarizability,[16] the concept of modulation of
conjugation warrants further experimental studies.
As a case in point that non-wiggly potentials can lead
to a large nonlinear susceptibility is the clipped harmonic
oscillator, which we calculated to have an intrinsic hyper-
polarizability of about 0.57.[14] Figure 5 shows the opti-
mized non-perturbative hyperpolarizability when using a
clipped harmonic oscillator as the starting potential. The
properties of all of the optimized potentials are summa-
TABLE II: Summary of calculations with different starting
potentials. βs is the hyperpolarizability of the starting po-
tential while the other ones are after optimization.
Function βS βSOS βDF βNP τ
00 ∆τ
V (x) (×10−5) (×10−4)
0 0 0.7089 0.7089 0.7089 37.8 5.33
30 tanh(x) 0.67 0.7084 0.6918 0.7083 779 11.8
x 0.66 0.7088 0.7072 0.7088 78.7 8.79
x2 0.57 0.7089 0.7085 0.7088 18.6 703
x1/2 0.68 0.7087 0.7049 0.7087 190 9.76
x+ sin(x) 0.67 0.7088 0.7073 0.7088 75.0 8.46
x+ 10 sin(x) 0.04 0.7085 0.7085 0.7085 1.65 7.78
0 5 10 15 20
FIG. 5: Optimized potential energy function and first 15
wavefunctions after 8,000 iterations. Starting potential is
V (x) = x2, using the non-perturbative hyperpolarizability
for optimization.
rized in Table II. The clipped square root function also
has a large hyperpolarizability (0.69). The optimized po-
tential is shown in Figure 6. In these cases, the amplitude
of the wiggles are small and all the wavefunctions overlap.
So, these fall into Class II. Note that the lack of wiggles
shows that they are not an inevitable consequence of our
numerical calculations.
We may question whether small wiggles in the poten-
tial energy function lead to large amplitude wiggles as an
artifact of our numerical optimization technique. To test
this hypothesis, we used the trial potential energy func-
tion x + sin(x), where the wiggle amplitude is not large
enough to cause the wavefunctions to localize at the min-
ima. The optimized potential energy function retains an
approximately linear from with only small fluctuation.
In fact, the results are very similar to what we found for
the linear starting potential and the wiggles do not affect
the final result. The similarity between these cases can
be seen in Table II.
Next, we test a starting potential with large wiggles
as shown in the upper portion of Figure 7. The lower
0 5 10 15 20
FIG. 6: Optimized potential energy function and first 15
wavefunctions after 8,000 iterations. Starting potential is
V (x) =
x, using the non-perturbative hyperpolarizability
for optimization.
energy eigenfunctions are found to be localized mostly in
the first two wells. In fact, the lowest four energy eigen-
functions are well approximated by harmonic oscillator
wavefunctions, which are centrosymmetric. As a result,
the first excited state holds most of the oscillator strength
and the value of the intrinsic hyperpolarizability is only
0.04.
After 3000 interactions, this Class I potential energy
function has high amplitude wiggles at a wavelength that
is significantly shorter than the wavelength of the ini-
tial sine function (bottom portion of Figure 7). In com-
mon with the optimized tanh(x) function, the wiggles
are of large but almost chaotically varying amplitude.
This leads to wavefunctions that are spatially separated.
While the wavefunctions are not as well separated as
we find for the tanh(x) starting potential, the optimized
potential yields only two dominant transition from the
ground state; so, this system is well approximated by
a three-level model. As is apparent from Table II, the
ground state sum rule (characterized by τ
00 ) is better
obeyed in this optimized potential than in any others.
So, the wavefunctions are accurate and all of the values
of β have converged to the same value, suggesting that
this calculation may be the most accurate of the set
Our results bring up several interesting questions.
First, all of our extensive numerical calculations, inde-
pendent of the starting potential, yield an optimized in-
trinsic hyperpolarizability with an upper bound of 0.71,
which is about thirty percent lower than what the sum
rules allow. Since numerical optimization can settle in to
a local maximum, it is possible that all of the starting
potentials are far from the global maximum of βint = 1.
Indeed, since most potentials lead to systems that require
more than three dominant states to express the hyper-
polarizability, this may in itself be an indicator that we
are not at the fundamental limit precisely because these
0 5 10 15 20
0 5 10 15 20
FIG. 7: Potential energy function and first 15 wavefunctions
before (top) and after (bottom) 3,000 iterations. Starting
potential is of the form V (x) = x+ 10 sin(x), using the non-
perturbative hyperpolarizability for optimization.
systems have more than three states. Indeed, the orig-
inal results of Zhou and coworkers frames the problem
in a way (i.e. a 15-level model in a potential limited to
about 20 wiggles) that allows a solution to the optimiza-
tion problem to lead to three dominant states. So, while
it may be argued that this system is contrived and un-
physical, we have found value in trying such toy models
when testing various hypotheses. This toy model
• leads to a three-level system as the three-level
ansatz proposes
• has the same qualitative properties as more precise
methods
• has given insights into making new molecules with
record-breaking intrinsic hyperpolarizability
Given the complexity of calculating nonlinear-
susceptibilities, our semi-quantitative method may
be a good way of generating new ideas.
The three-level ansatz proposes that at the fundamen-
tal limit, all transitions are negligible except between
three dominant states. There appears to be no proof of
the ansatz aside from the fact that it leads to an accurate
prediction of the upper bound of nonlinear susceptibili-
ties, both calculated and measured. To understand the
motivation behind the ansatz, it is useful to understand
how the two-level model optimizes the polarizability, α,
without the need to rely on any assumptions. This is
trivial to show by using the fact that the polarizability
depends only on the positive-definite transition moments,
〈0|x |n〉 〈n|x |0〉, the same parameters that are found in
the ground state sum rules.[26]
For nonlinear susceptibilities, the situation is much
more complicated because the SOS expression depends
on quantities such as 〈0|x |n〉 〈n|x |m〉 〈m|x |0〉, where
these terms can be both positive and negative. Fur-
thermore, the sum rules that relate excited states mo-
ments to each other allow for these moments to be much
larger than transition moments to the ground state. So,
it would seem plausible that one could design a system
with many excited states in a way that all of the tran-
sition moments between excited states would add con-
structively to yield a larger hyperpolarizability than what
we calculate with the three-level ansatz. None of our
numerical calculations, independent of the potential en-
ergy function, yield a value greater than 0.71. Since our
potential energy functions are general 1-dimensional po-
tentials (i.e. the potentials are not limited to Coulomb
potentials, nor are the wavefunctions approximated as is
common in standard quantum chemical computations),
our calculations most likely span a broader range of pos-
sible wavefunctions leading to a larger variety of states
that contribute to the hyperpolarizability.
However, there appear to be local maxima associated
with systems that behave as a three-level system and oth-
ers with many states, and, the maximum values both are
0.71. It is interesting that so may different sets of transi-
tion moments and energies can yield the exact same local
maximum. To gain a deeper appreciation of the under-
lying physics, let’s consider the transition moments and
energies in the sum-over-states expression for the hyper-
polarizability as adjustable parameters. For a system
with N states, there are N − 1 energy parameters of the
form En − E0. The moment matrix xij has N2 com-
ponents. If the matrix is real, there are (N2 − N)/2
unique off-diagonal terms and N diagonal dipole mo-
ments. Since all dipole moments appear as differences
of the form xnn − x00, there are only N − 1 dipole mo-
ment parameters. Therefore, the dipole matrix is char-
acterized by (N2 − N)/2 + N − 1 = (N + 2)(N − 1)/2
parameters. Combining the energy and dipole matrix pa-
rameters, there are a total of (N + 2)(N − 1)/2 +N − 1
parameters.
The N-state sum rules are of the form:
(Em + Ep)
〈m|x |n〉 〈n|x |p〉 (7)
δm,p,
so the sum rules comprise a total of N2 equations (i.e. an
equation for each (m, p)). If the sum rules are truncated
to N states, the sum rule indexed by (m = N, p = N)
is nonsensical because it contradicts the other sum rules.
Furthermore, if the transition moments are real, then
xmp = xpm, so only (N
2 − N)/2 of the equations are
independent. As such, there are a total of (N2 −N)/2+
N − 1 = (N + 2)(N − 1)/2 independent equations.
Since the SOS expression for the nonlinear-
susceptibility has (N + 2)(N − 1)/2 +N − 1 parameters
and the sum rules provide (N + 2)(N − 1)/2 equations,
the SOS expression can be reduced to a form with N − 1
parameters. For example, the three-level model for the
hyperpolarizability, which is expressed in terms of 7
parameters, can be reduced to two parameters using 5
sum rule equations. In practice, however, even fewer
sum rule equations are usually available because some
of them lead to physically unreasonable consequences.
While the (N,N) sum rule is clearly unphysical due
to truncation, sum rule equations that are near equa-
tion (N,N) may also be unphysical. In the case of
the three-level model, it is found that the equation
(2, 1) allows for an infinite hyperpolarizability, so that
equation is ignored on the grounds that it violates the
principle of physical soundness.[12, 25, 26] This leads to
a hyperpolarizability in terms of 3 variables, which are
chosen to be E10, E = E10/E20, and X = x10/x
The expression is then maximized with respect to the
two parameters E and X , leaving the final result a
function of E10.
We conclude that the SOS expression for the hyperpo-
larizability can be expressed in terms of at least N − 1
parameters; so, it would appear that as more levels are
included in the SOS expression, there are more free pa-
rameters that can be varied without violating the sum
rules. As N → ∞, there are an infinite number of ad-
justable parameters. So, it is indeed puzzling that the
three-level ansatz yields a fundamental limit that is con-
sistent with all of our calculations for a wide range of
potentials, many of which have many excited states. It
may be that we are only considering a small subset of
potential energy functions; or, perhaps the expression
for the hyperpolarizability depends on the parameters in
such a way that large matrix elements contribute to the
hyperpolarizability with alternating signs so that the big
terms cancel. This is a puzzle that needs to be solved if
we are to understand what makes β large.
To investigate whether the limiting behavior is
due to our use of 1-dimensional potentials, we have
also optimized the intrinsic hyperpolarizability in two-
dimensions. In this case, we focus on the largest ten-
sor component, βxxx and describe the potential as a
superposition of point charges. As described in the
literature,[19] we solve the two-dimensional Schrödinger
eigenvalue problem,
∇2Ψ+ VΨ = EΨ, (8)
for the lowest ten to 25 energy eigenstates, depending on
the degree of convergence of the resulting intrinsic hyper-
polarizability. We use the two-dimensional logarithmic
Coulomb potential, which for k nuclei with charges q1e,
. . . , qke located at points s
(1), . . . , s(k), is given by
V (s) =
qj log ‖s− s(j)‖, (9)
where L is a characteristic length. With L = 2Å, the
force due to a charge at distance 2Å is the same as it
would be for a 3D Coulomb potential.
We discretize the eigenvalue problem given by Equa-
tion 8 using a quadratic finite element method [21, 27]
and solve the resulting matrix eigenvalue problem for the
ten to 25 smallest energy eigenvalues and corresponding
eigenvectors by the implicitly-restarted Arnoldi method
[22] as implemented in ARPACK [28]. Each eigenvector
yields a wave function Ψn corresponding to energy level
En. The moments
xmn =
s1Ψm(s1, s2)Ψn(s1, s2) ds1ds2
are computed, and these and the energy levels En are
used to compute β
Figure 8 shows the intrinsic hyperpolarizability of a
two-nucleus molecule plotted as a function of the distance
between the two nuclei and nuclear charge q1. The total
nuclear charge is q1 + q2 = +e, and is expressed in units
of the proton charge, e. Three extrema are observed.
The positive peak parameters are βint = 0.649 for q1 =
0.58 and d = 4.36Å. The negative one yields βint =
−0.649 for q1 = 0.42 and d = 4.36Å. The local negative
peak that extends past the graph on the right reaches its
maximum magnitude of βint = −0.405 at q1 = 2.959 and
d = 2.0Å.
Applying numerical optimization to the intrinsic hy-
perpolarizability using the charges and separation be-
tween the nuclei as parameters, we get βint = 0.654 at
d = 4.539Å, and q1 = 0.430 when the starting param-
eters are near the positive peak; and βint = −0.651,
d = 4.443Å, and q1 = 0.572 when optimization gives
the negative peak. The peak parameters are the same
within roundoff errors when optimization or plotting is
used, confirming that the optimization procedure yields
the correct local extrema.
Figure 9 shows the intrinsic hyperpolarizability of an
octupolar-like molecule made of three evenly-spaced nu-
clei on a circle plotted as a function of the circle’s diam-
eter and charge fraction ǫ (q = ǫe). The charge on each
of the other nuclei is e(1 − ǫ)/2. The positive peak at
ǫ = 0.333 and diameter D = 6.9Å has a hyperpolariz-
ability βint = 0.326, while βint = −0.605 for a charge
fraction ǫ = 0.44 and a diameter D = 6.8Å.
Charge (q
Intrinsic Hyperpolarizability 
0 0.5 1 1.5 2
FIG. 8: The intrinsic hyperpolarizability of two nuclei as a
function of the distance between them and the charge of one
nucleus, q1 where q1 + q2 = +e.
Intrinsic Hyperpolarizability
0 0.2 0.4 0.6 0.8 1
FIG. 9: The intrinsic hyperpolarizability of three evenly-
spaced nuclei on a circle as a function of the circle’s diameter
and the charge, ǫ (in units of e), on one of the nuclei. The
charge on each of the other nuclei is e(1− ǫ)/2.
When the positions and magnitudes of the three
charges are allowed to move freely, the best intrinsic
hyperpolarizability obtained using numerical optimiza-
tion is βint = 0.685 for charges located at ~r1 = (0, 0).
~r2 = (−4.87Å, 0.33Å), and ~r3 = (−9.57Å,−0.16Å); with
charges q1 = 0.43e, q2 = 0.217e, and q3 = 0.351e. There
are only small differences in the optimized values of βint
depending on the starting positions and charges; and the
best results are for a “molecule” that is nearly linear
along the x-direction. This is not surprising given that
the xxx-component of βint is the optimized quantity.
The two-dimensional analysis illustrates that numer-
ical optimization correctly identifies the local maxima
(peaks and valleys) and that the magnitude of maximum
intrinsic hyperpolarizability (0.65 vs 0.68) is close to the
maximum we get for the one-dimensional optimization
of the potential energy function (0.71). All computa-
tions we have tried, including varying the potential en-
ergy function in one dimension or moving around point
charges in a plane all yield an intrinsic hyperpolarizabil-
ity that is less than 0.71.
An open question is the origin of the factor-of-thirty
gap between the best molecules and the fundamental
limit, which had remained firm for decades through the
year 2006. Several of the common proposed explana-
tions, such as vibronic dilution, have been eliminated.[14]
Perhaps it is not possible to make large-enough varia-
tions of the potential energy function without making
the molecule unstable. Or, perhaps there are subtle is-
sues with electron correlation, which prevents electrons
from responding to light with their full potential. The
fact that the idea of modulation of conjugation has lead
to a 50% increase over the long-standing ceiling - reduc-
ing the gap to a factor of twenty - makes it a promising
approach for potential further improvements. Continued
theoretical scrutiny, coupled with experiment, will be re-
quired to confirm the validity of our approach.
IV. CONCLUSIONS
There appear to be many potential energy functions
that lead to an intrinsic hyperpolarizability that is near
the fundamental limit. These separate into two broad
classes: one in which wiggles in the potential energy
function forces the eigenfunctions to be spatially sepa-
rated and a second class of monotonically varying wave-
functions with small or no wiggles that allow for many
strongly overlapping wavefunctions. Interestingly, all
these one-dimensional “molecules” have the same max-
imal intrinsic hyperpolarizability of 0.71. It is puzzling
that the three-level ansatz correctly predicts the funda-
mental limit even when the ansatz does not apply. A
second open question pertains to the origin of the long-
standing factor of 30 gap between the fundamental limit
and the best molecules. The idea of conjugation modu-
lation may be one promising approach for making wiggly
potential energy profiles that lead to molecules that fall
into the gap. Given that there are so many choices of
potential energy functions that lead to maximal intrinsic
hyperpolarizability, it may be possible to engineer many
new classes of exotic molecules with record intrinsic hy-
perpolarizability.
Acknowledgements: MGK thanks the National
Science Foundation (ECS-0354736) and Wright Paterson
Air Force Base for generously supporting this work.
[1] Q. Y. Chen, L. Kuang, Z. Y. Wang, and E. H. Sargent,
Nano. Lett. 4, 1673 (2004).
[2] B. H. Cumpston, S. P. Ananthavel, S. Barlow, D. L. Dyer,
J. E. Ehrlich, L. L. Erskine, A. A. Heikal, S. M. Kuebler,
I.-Y. S. Lee, D. McCord-Maughon, et al., Nature 398, 51
(1999).
[3] S. Kawata, H.-B. Sun, T. Tanaka, and K. Takada, Nature
412, 697 (2001).
[4] A. Karotki, M. Drobizhev, Y. Dzenis, P. N. Taylor, H. L.
Anderson, and A. Rebane, Phys. Chem. Chem. Phys. 6,
7 (2004).
[5] I. Roy, O. T. Y., H. E. Pudavar, E. J. Bergey, A. R.
Oseroff, J. Morgan, T. J. Dougherty, and P. N. Prasad,
J. Am. Chem. Soc. 125, 7860 (2003).
[6] M. G. Kuzyk, Opt. Lett. 25, 1183 (2000).
[7] M. G. Kuzyk, IEEE Journal on Selected Topics in Quan-
tum Electronics 7, 774 (2001).
[8] M. G. Kuzyk, Phys. Rev. Lett. 85, 1218 (2000).
[9] M. G. Kuzyk, Opt. Lett. 28, 135 (2003).
[10] M. G. Kuzyk, Phys. Rev. Lett. 90, 039902 (2003).
[11] M. G. Kuzyk, J. Nonl. Opt. Phys. & Mat. 13, 461 (2004).
[12] M. G. Kuzyk, J. Chem Phys. 125, 154108 (2006).
[13] M. G. Kuzyk, Optics & Photonics News 14, 26 (2003).
[14] K. Tripathi, P. Moreno, M. G. Kuzyk, B. J. Coe,
K. Clays, and A. M. Kelley, J. Chem. Phys. 121, 7932
(2004).
[15] B. J. Orr and J. F. Ward, Molecular Physics 20, 513
(1971).
[16] J. Pérez Moreno, Y. Zhao, K. Clays, and M. G. Kuzyk,
Opt. Lett. 32, 59 (2007).
[17] J. Zhou, M. Kuzyk, and D. S. Watkins, Opt. Lett. 31,
2891 (2006).
[18] M. G. Kuzyk, Phys. Rev. A 72, 053819 (2005).
[19] M. G. Kuzyk and D. S. Watkins, J. Chem Phys. 124,
244104 (2006).
[20] D. Kincaid and E. W. Cheney, Numerical Analy-
sis: Mathematics of Scientific Computing (Brooks-Cole,
2002), 3rd ed.
[21] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The
Finite Element Method: Its Basis and Fundamentals
(Butterworth-Heinemanm, 2005), 6th ed.
[22] D. C. Sorensen, SIAM J. Matrix Anal. Appl. 13, 357
(1992).
[23] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. Wright,
SIAM J. Optim. 9, 112 (1998).
[24] B. Champagne and B. Kirtman, Phys. Rev. Lett. 95,
109401 (2005).
[25] M. G. Kuzyk, Phys. Rev. Lett. 95, 109402 (2005).
[26] M. G. Kuzyk, J. Nonl. Opt. Phys. & Mat. 15, 77 (2006).
[27] K. Atkinson and W. Han, Theoretical Numerical Anal-
ysis, a Functional Analysis Framework (Springer, New
York, 2001).
[28] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK
Users’ Guide: Solution of Large-Scale Eigenvalue Prob-
lems with Implicitly Restarted Arnoldi Methods (SIAM,
Philadelphia, 1998).
0 5 10 15 20
0 5 10 15 20
0 5 10 15 20

  We use numerical optimization to study the properties of (1) the class of
one-dimensional potential energy functions and (2) systems of point charges in
two-dimensions that yield the largest hyperpolarizabilities, which we find to
be within 30% of the fundamental limit. We investigate the character of the
potential energy functions and resulting wavefunctions and find that a broad
range of potentials yield the same intrinsic hyperpolarizability ceiling of
0.709.

Studies on optimizing potential energy functions for maximal intrinsic
hyperpolarizability
Juefei Zhou, Urszula B. Szafruga, David S. Watkins* and Mark G. Kuzyk
Department of Physics and Astronomy, Washington State University,
Pullman, Washington 99164-2814; and *Department of Mathematics
We use numerical optimization to study the properties of (1) the class of one-dimensional po-
tential energy functions and (2) systems of point charges in two-dimensions that yield the largest
hyperpolarizabilities, which we find to be within 30% of the fundamental limit. We investigate the
character of the potential energy functions and resulting wavefunctions and find that a broad range
of potentials yield the same intrinsic hyperpolarizability ceiling of 0.709.
I. INTRODUCTION
Materials with large nonlinear-optical suscep-
tibilities are central for optical applications such
as telecommunications,[1] three-dimensional nano-
photolithography,[2, 3] and making new materials[4] for
novel cancer therapies.[5] The fact that quantum calcu-
lations show that there is a limit to the nonlinear-optical
response[6, 7, 8, 9, 10, 11] is both interesting from the
basic science perspective; and, provides a target for
making optimized materials. In this work, we focus
on the second-order susceptibility and the underlying
molecular hyperpolarizability, which is the basis of
electro-optic switches and frequency doublers.
The fundamental limit of the off-resonance hyperpo-
larizability is given by,[8]
βMAX =
, (1)
where N is the number of electrons and E10 the energy
difference between the first excited state and the ground
state, E10 = E1 − E0. Using Equation 1, we can de-
fine the off-resonant intrinsic hyperpolarizability, βint, as
the ratio of the actual hyperpolarizability (measured or
calculated), β, to the fundamental limit,
βint = β/βMAX . (2)
We note that since the dispersion of the fundamental
limit of β is also known, [12] it is possible to calculate
the intrinsic hyperpolarizability at any wavelength. In
the present work, we treat only the zero-frequency limit.
Until recently, the largest nonlinear susceptibilities of
the best molecules fell short of the fundamental limit by a
factor of 103/2, [10, 13, 14] so the very best molecules had
a value of βint = 0.03. Since a Sum-Over-States (SOS)
calculation of the hyperpolarizability[15] using the ana-
lytical wavefunctions of the clipped harmonic oscillator
yields a value βint = 0.57,[14] the factor-of-thirty gap is
not of a fundamental nature. Indeed, recently, it was re-
ported that a new molecule with asymmetric conjugation
modulation has a measured value of βint = 0.048.[16]
To investigate how one might make molecules with
a larger intrinsic hyperpolarizability, Zhou and cowork-
ers used a numerical optimization process where a trial
potential energy function is entered as an input, and
the code iteratively deforms the potential energy func-
tion until the intrinsic hyperpolarizability, calculated
from the resulting wavefunctions, converges to a local
maximum.[17] In this work, a hyperbolic tangent func-
tion was used as the starting potential due to the fact that
it is both asymmetric yet relatively flat away from the ori-
gin. This calculation was one-dimensional and included
only one electron, so electron correlation effects were
ignored. Furthermore, the intrinsic hyperpolarizability
was calculated using the new dipole-free sum-over-states
expression[18] and only 15 excited states were included.
The resulting optimized potential energy function showed
strong oscillations, which served to separate the spatial
overlap between the energy eigenfunctions. This led Zhou
and coworkers to propose that modulated conjugation
in the bridge between donor and acceptor ends of such
molecules may be a new paradigm for making molecules
with higher intrinsic hyperpolarizability.[17]
Based on this paradigm, Pérez Moreno reported mea-
surements of a class of chromophores with varying de-
gree of modulated conjugation.[16] The best measured
intrinsic hyperpolarizability was βint = 0.048, about 50%
larger than the best previously-reported. Given the mod-
est degree of conjugation modulation for this molecule,
this new paradigm shows promise for further improve-
ments.
In the present work, we extend Zhou’s calculations to
a larger set of starting potentials. To circumvent trun-
cation problems associated with sum-over-states calcula-
tions, we instead determine the hyperpolarizability using
a finite difference technique. The optimization procedure
is then applied to this non-perturbative hyperpolarizabil-
To study the effects of geometry on the hyperpo-
larizability, Kuzyk and Watkins calculated the hyper-
polarizability of various arrangements of point charges,
representing nuclei, in two-dimensions using a two-
dimensional Coulomb potential.[19] In the present con-
tribution, we apply our numerical optimization technique
to determine the arrangement and charges of the nuclei
in a planar molecule that maximize the intrinsic hyper-
polarizability.
http://arxiv.org/abs/0704.1687v1
II. THEORY
In our previous work, we used a finite-state SOS model
of the hyperpolarizability that derives from perturbation
theory (we used both the standard Orr andWard SOS ex-
pression, βSOS ,[15] and the newer dipole free expression,
βDF [18]). The use of a finite number of states in lieu of
the full infinite sums can result in inaccuracies, so, in the
present work, we use the non-perturbative approach, as
follows. We begin by solving the 1-d Schrodinger Equa-
tion on the interval a < x < b for the ground state wave-
function ψ(x,E) of an electron in a potential well defined
by V (x) and in the presence of a static electric field, E,
that adds to the potential δV = −exE. From this, the
off-resonant hyperpolarizability is calculated with numer-
ical differentiation, i.e. using finite differences, yielding
βNP =
|ψ(x,E)|2 ex dx
. (3)
Equation 3 is evaluated using the standard second-order
approximation to the second derivative:
f ′′(z) ≈ f(z + h)− 2f(z) + f(z − h)
with several h values h0, h0/5, h0/25, . . . . We then refine
these values by Richardson extrapolation [20] and obtain
our estimate from the two closest extrapolated values.
Our computational mesh consists of 200 quadratic fi-
nite elements with a total of 399 degrees of freedom. The
potential energy function is a cubic spline with 40 degrees
of freedom. Thus the numerical calculations in regions
where the potential function is represented by 3 points
in the spline are covered by 15 elements with a total of
about 30 degrees of freedom.
Calculating βint from Equations 3, 2 and 1 for a specific
potential, we use an optimization algorithm that contin-
uously varies the potential in a way that maximizes βint.
We also compute the matrix[17, 26]
τ (N)mp = δm,p −
xmax10
xmax10
, (4)
where xmax10 is the magnitude of the fundamental limit of
the position matrix element x10 for a one electron system,
and is given by,
xmax10 =
2mE10
. (5)
Each matrix element of τ (N), indexed by m and p, is
a measure of how well the (m, p) sum rule is obeyed
when truncated to N states. If the sum rules are ex-
actly obeyed, τ
mp = 0 for all m and p. We note that if
the sum rules are truncated to an N-state model, the sum
rules indexed by a large value of m or p (i.e. m, p ∼ N)
are disobeyed even when the position matrix elements
and energies are exact. We have found that the values
mp are small for exact wavefunctions when m < N/2
and p < N/2. So, when evaluating the τ matrix to
test our calculations, we consider only the components
m≤N/2,p≤N/2
We observe that when using more than about 40 states
in SOS calculations of the hyperpolarizability only a
marginal increase of accuracy results when the poten-
tial energy function is parameterized with 400 degrees of
freedom. So, to ensure overkill, we use 80 states when
calculating the τ matrix or the hyperpolarizability with
an SOS expression so that truncation errors are kept to
a minimum. Since the hyperpolarizability depends crit-
ically on the transition dipole moment from the ground
state to the excited states, we use the value of τ
00 as
one important test of the accuracy of the calculated wave-
functions. Additionally, we use the standard deviation of
τ (N),
∆τ (N) =
, (6)
which quantifies, on average, how well the sum rules are
obeyed in aggregate, making ∆τ (N) a broader test of the
accuracy of a large set of wavefunctions.
Our code is written in MATLAB. For each trial po-
tential we use a quadratic finite element method [21] to
approximate the Schrödinger eigenvalue problem and the
implicitly restarted Arnoldi method [22] to compute the
wave functions and energy levels. To optimize β we use
the Nelder-Mead simplex algorithm [23].
As described in our previous work,[17] we perform op-
timization, but this time using the exact intrinsic hyper-
polarizability β = βNP /βMAX , where βMAX is the fun-
damental limit of the hyperpolarizability, which is pro-
portional to E
10 . To determine E10 ≡ E1 −E0, we also
calculate the first excited state energy E1.
III. RESULTS AND DISCUSSIONS
Figure 1 shows an example of the optimized poten-
tial energy function after 7,000 iterations when starting
with the potential V (x) = 0 and optimizing the non-
perturbative intrinsic hyperpolarizability βNP /βMAX as
calculated with Equation 3. Also shown are the eigen-
functions of the first 15 states computed from the opti-
mized potential. First, we note that the potential energy
function shows the same kinds of wiggles as in our original
paper,[17] though not of sufficient amplitude to localize
the wavefunctions.
For the starting potentials we have investigated, our
results fall into two broad classes. In the first, three
common features are: (1) The best intrinsic hyperpolar-
izabilities are near βint = 0.71; (2) the best potentials
have a series of wiggles; and (3) the systems behave as a
0 5 10 15 20
FIG. 1: Optimized potential energy function and first 15
wavefunctions after 7,000 iterations. Starting potential is
V (x) = 0, using the non-perturbative hyperpolarizability for
optimization.
limited-state model. In the second class of starting po-
tentials, (2) the wiggles are much less pronounced and (3)
more states contribute evenly. Figure 1 is an example of
a Class II potential. However, in both classes, the max-
imum calculated intrinsic hyperpolarizability appears to
be bounded by βint = 0.71. Using the set of potentials
from both classes that lead to optimized βNP /βMAX , we
calculate the lowest 80 eigenfunctions and eigenvalues,
from which we calculate βDF and βSOS . In most cases,
we find that the three different formulas for β converge
to the same value when only the first 20 excited states
are used (using 80 states, the three are often the same to
within at least 4 decimal places) and τ00 ≈ 10−4, showing
that the ground state sum rules are well obeyed. Further-
more, the rms deviation of the τ matrix when including
40 states leads to τ (80) < 0.001.
Figure 2 shows an example of the optimized potential
energy function when starting with the potential V (x) =
tanhx and optimizing the exact (non-perturbative) in-
trinsic hyperpolarizability. Also shown are the eigenfunc-
tions of the first 15 states computed with the optimized
potential. First, we note that the potential energy func-
tion shows the same kinds of wiggles as in our original
paper;[17] and only 2 excited state wavefunctions and the
ground state are localized in the first deep well - placing
this system in Class I.
The observation that such potentials lead to hyper-
polarizabilities that are near the fundamental limit mo-
tivated Zhou and coworkers to suggest that molecules
with modulated conjugation may have enhanced intrin-
sic hyperpolarizabilities.[17] A molecule with a modu-
lated conjugations bridge between the donor and accep-
tor end was later shown to have record-high intrinsic
hyperpolarizability.[16] As such, this result warrants a
more careful analysis.
It is worthwhile to compare our present results charac-
0 5 10 15 20
FIG. 2: Optimized potential energy function and first 15
wavefunctions after 8,000 iterations. Starting potential is
V (x) = tanh(x), using the non-perturbative hyperpolarizabil-
ity for optimization.
terized by Figure 2 with our past work,[17] particularly
for the purpose of examining the impact of the approx-
imations used in the previous work.[17] Figure 3 shows
the optimized potential and wavefunctions obtained by
Zhou and coworkers using a 15-state model and opti-
mizing the dipole-free intrinsic hyperpolarizability. Since
only 15 states were used, the SOS expression for β did not
fully converge; making the result inaccurate as suggested
by the fact that βSOS and βDF did not agree. How-
ever, since the code focused on optimizing the dipole-free
form of β, and τ00 was small when βint was optimized,
the dipole-free expression may have converged to a rea-
sonably accurate value while the commonly-used SOS
expression was inaccurate. Indeed, it was found that
βDF ≈ 0.72 - in contrast to our more precise present
calculations using the non-perturbative approach, which
yields βNP < 0.71. So, the fact that our more precise
calculations, which do not rely on a sum-over states ex-
pression, agree so well with the 15-state model suggests
that in both cases, the limit for a one-dimensional single
electron molecule is just over β ≈ 0.7. This brute force
calculation serves as a numerical illustration of the obser-
vation that the limiting value of β is the same for an exact
non-perturbation calculation and for a calculation that
truncates the SOS expression, which presumedly should
lead to large inaccuracies.[24, 25] At minimum, this result
supports the existence of fundamental limits of nonlinear
susceptibilities that are in line with past calculations.
To state Zhou’s approach more precisely,[17] the cal-
culations optimized the very special case of the intrin-
sic hyperpolarizability for a 15 state model for a poten-
tial energy function that is parameterized with 20 spline
points. As such, the potential energy function can at
most develop about 20 wiggles. As a consequence, there
are enough degrees of freedom in the potential energy
function to force the 15 states to be spatially well sepa-
0 5 10 15 20
 V(x)
FIG. 3: Optimized potential energy function using βDF and
first 15 wavefunctions after 7,000 iterations. Starting po-
tential is the tanh(x) potential. The final potential (shown
above) we refer to as the Zhou potential.
TABLE I: Evolution of Zhou’s Potential. βs is the hyperpo-
larizability of the starting potential using 80 states while the
other ones are after optimization of βNP .
Number of βS βSOS βDF βNP τ
00 ∆τ
Iterations (×10−5) (×10−4)
0 0.5612 0.5612 0.5607 0.5612 11.2 15
1000 0.5612 0.7087 0.6682 0.7083 1810 40
rated. Interestingly, after optimization, only two excited
states overlap with ground state, allowing only these two
states to have nonzero transition dipole moments with
each other and the ground state – forcing the system
into a three-level SOS model for βDF . This behavior is
interesting in light of the three-level ansatz, which asserts
that only three states determine the nonlinear response
of a system when it is near the fundamental limits.
It is interesting to compare the exact non-perturbation
calculation, which does not depend on the excited state
wavefunctions (Figure 2) and Zhou’s contrived system of
15 states (Figure 3). Both cases have wiggles and the
wavefunctions appear to be mostly non-overlapping. So,
for the first 15 states, the wavefunctions appear simi-
larly localized. The situation becomes more interesting
when 80 states are included in calculating the hyperpo-
larizability for Zhou’s potential or when the exact non-
perturbative approach is used. The first line in Table
I summarizes the results with Zhou’s potential and 80
states.
First, let’s focus on the sum-over-states results.
Clearly, when 80 states are used in the calculation, it
is impossible for the excited state wavefunctions to not
overlap with each other, so the three-level approximation
to β breaks down. According to the three-level ansatz,
we would expect the hyperpolarizability to get smaller.
Indeed, the additional excited states result in a smaller
0 5 10 15 20
FIG. 4: Optimized potential energy function and first 15
wavefunctions after 1,000 iterations. Starting potential is
Zhou’s potential, using the non-perturbative hyperpolariz-
ability for optimization.
hyperpolarizability (≈ 0.56). Note that the exact and
SOS expressions agree with each other and that τ
and ∆τ (80) are small.
Figure 4 shows the result after 1000 iterations, us-
ing Zhou’s potential as the starting potential and us-
ing the non-perturbative hyperpolarizability for opti-
mization. First, the non-perturbative hyperpolarizabil-
ity reaches just under 0.71, but, the SOS and dipole-
free expressions do not agree with each other. Further-
more, both convergence metrics (τ
00 and ∆τ
(80)) are
larger than before optimization. It would appear that
for Zhou’s potential, even 80 states are not sufficient
to characterize the nonlinear susceptibility when a sum-
over-states expression is used (either dipole free or tradi-
tional SOS expression - though the SOS expression agrees
better with the non-perturbative approach).
Interestingly, the optimized potential energy function
still retains the wiggles and the wave functions are still
well separated. This result is consistent with the sug-
gestion of Zhou and coworkers that modulation of conju-
gation may be a good design strategy for making large-
hyperpolarizability molecules. We note that wiggles in
the potential energy function are not required to get
a large nonlinear-optical response; but, appears to be
one way that Mother Nature optimizes the hyperpo-
larizability. Since this idea has been used to identify
molecules with experimentally measured record intrin-
sic hyperpolarizability,[16] the concept of modulation of
conjugation warrants further experimental studies.
As a case in point that non-wiggly potentials can lead
to a large nonlinear susceptibility is the clipped harmonic
oscillator, which we calculated to have an intrinsic hyper-
polarizability of about 0.57.[14] Figure 5 shows the opti-
mized non-perturbative hyperpolarizability when using a
clipped harmonic oscillator as the starting potential. The
properties of all of the optimized potentials are summa-
TABLE II: Summary of calculations with different starting
potentials. βs is the hyperpolarizability of the starting po-
tential while the other ones are after optimization.
Function βS βSOS βDF βNP τ
00 ∆τ
V (x) (×10−5) (×10−4)
0 0 0.7089 0.7089 0.7089 37.8 5.33
30 tanh(x) 0.67 0.7084 0.6918 0.7083 779 11.8
x 0.66 0.7088 0.7072 0.7088 78.7 8.79
x2 0.57 0.7089 0.7085 0.7088 18.6 703
x1/2 0.68 0.7087 0.7049 0.7087 190 9.76
x+ sin(x) 0.67 0.7088 0.7073 0.7088 75.0 8.46
x+ 10 sin(x) 0.04 0.7085 0.7085 0.7085 1.65 7.78
0 5 10 15 20
FIG. 5: Optimized potential energy function and first 15
wavefunctions after 8,000 iterations. Starting potential is
V (x) = x2, using the non-perturbative hyperpolarizability
for optimization.
rized in Table II. The clipped square root function also
has a large hyperpolarizability (0.69). The optimized po-
tential is shown in Figure 6. In these cases, the amplitude
of the wiggles are small and all the wavefunctions overlap.
So, these fall into Class II. Note that the lack of wiggles
shows that they are not an inevitable consequence of our
numerical calculations.
We may question whether small wiggles in the poten-
tial energy function lead to large amplitude wiggles as an
artifact of our numerical optimization technique. To test
this hypothesis, we used the trial potential energy func-
tion x + sin(x), where the wiggle amplitude is not large
enough to cause the wavefunctions to localize at the min-
ima. The optimized potential energy function retains an
approximately linear from with only small fluctuation.
In fact, the results are very similar to what we found for
the linear starting potential and the wiggles do not affect
the final result. The similarity between these cases can
be seen in Table II.
Next, we test a starting potential with large wiggles
as shown in the upper portion of Figure 7. The lower
0 5 10 15 20
FIG. 6: Optimized potential energy function and first 15
wavefunctions after 8,000 iterations. Starting potential is
V (x) =
x, using the non-perturbative hyperpolarizability
for optimization.
energy eigenfunctions are found to be localized mostly in
the first two wells. In fact, the lowest four energy eigen-
functions are well approximated by harmonic oscillator
wavefunctions, which are centrosymmetric. As a result,
the first excited state holds most of the oscillator strength
and the value of the intrinsic hyperpolarizability is only
0.04.
After 3000 interactions, this Class I potential energy
function has high amplitude wiggles at a wavelength that
is significantly shorter than the wavelength of the ini-
tial sine function (bottom portion of Figure 7). In com-
mon with the optimized tanh(x) function, the wiggles
are of large but almost chaotically varying amplitude.
This leads to wavefunctions that are spatially separated.
While the wavefunctions are not as well separated as
we find for the tanh(x) starting potential, the optimized
potential yields only two dominant transition from the
ground state; so, this system is well approximated by
a three-level model. As is apparent from Table II, the
ground state sum rule (characterized by τ
00 ) is better
obeyed in this optimized potential than in any others.
So, the wavefunctions are accurate and all of the values
of β have converged to the same value, suggesting that
this calculation may be the most accurate of the set
Our results bring up several interesting questions.
First, all of our extensive numerical calculations, inde-
pendent of the starting potential, yield an optimized in-
trinsic hyperpolarizability with an upper bound of 0.71,
which is about thirty percent lower than what the sum
rules allow. Since numerical optimization can settle in to
a local maximum, it is possible that all of the starting
potentials are far from the global maximum of βint = 1.
Indeed, since most potentials lead to systems that require
more than three dominant states to express the hyper-
polarizability, this may in itself be an indicator that we
are not at the fundamental limit precisely because these
0 5 10 15 20
0 5 10 15 20
FIG. 7: Potential energy function and first 15 wavefunctions
before (top) and after (bottom) 3,000 iterations. Starting
potential is of the form V (x) = x+ 10 sin(x), using the non-
perturbative hyperpolarizability for optimization.
systems have more than three states. Indeed, the orig-
inal results of Zhou and coworkers frames the problem
in a way (i.e. a 15-level model in a potential limited to
about 20 wiggles) that allows a solution to the optimiza-
tion problem to lead to three dominant states. So, while
it may be argued that this system is contrived and un-
physical, we have found value in trying such toy models
when testing various hypotheses. This toy model
• leads to a three-level system as the three-level
ansatz proposes
• has the same qualitative properties as more precise
methods
• has given insights into making new molecules with
record-breaking intrinsic hyperpolarizability
Given the complexity of calculating nonlinear-
susceptibilities, our semi-quantitative method may
be a good way of generating new ideas.
The three-level ansatz proposes that at the fundamen-
tal limit, all transitions are negligible except between
three dominant states. There appears to be no proof of
the ansatz aside from the fact that it leads to an accurate
prediction of the upper bound of nonlinear susceptibili-
ties, both calculated and measured. To understand the
motivation behind the ansatz, it is useful to understand
how the two-level model optimizes the polarizability, α,
without the need to rely on any assumptions. This is
trivial to show by using the fact that the polarizability
depends only on the positive-definite transition moments,
〈0|x |n〉 〈n|x |0〉, the same parameters that are found in
the ground state sum rules.[26]
For nonlinear susceptibilities, the situation is much
more complicated because the SOS expression depends
on quantities such as 〈0|x |n〉 〈n|x |m〉 〈m|x |0〉, where
these terms can be both positive and negative. Fur-
thermore, the sum rules that relate excited states mo-
ments to each other allow for these moments to be much
larger than transition moments to the ground state. So,
it would seem plausible that one could design a system
with many excited states in a way that all of the tran-
sition moments between excited states would add con-
structively to yield a larger hyperpolarizability than what
we calculate with the three-level ansatz. None of our
numerical calculations, independent of the potential en-
ergy function, yield a value greater than 0.71. Since our
potential energy functions are general 1-dimensional po-
tentials (i.e. the potentials are not limited to Coulomb
potentials, nor are the wavefunctions approximated as is
common in standard quantum chemical computations),
our calculations most likely span a broader range of pos-
sible wavefunctions leading to a larger variety of states
that contribute to the hyperpolarizability.
However, there appear to be local maxima associated
with systems that behave as a three-level system and oth-
ers with many states, and, the maximum values both are
0.71. It is interesting that so may different sets of transi-
tion moments and energies can yield the exact same local
maximum. To gain a deeper appreciation of the under-
lying physics, let’s consider the transition moments and
energies in the sum-over-states expression for the hyper-
polarizability as adjustable parameters. For a system
with N states, there are N − 1 energy parameters of the
form En − E0. The moment matrix xij has N2 com-
ponents. If the matrix is real, there are (N2 − N)/2
unique off-diagonal terms and N diagonal dipole mo-
ments. Since all dipole moments appear as differences
of the form xnn − x00, there are only N − 1 dipole mo-
ment parameters. Therefore, the dipole matrix is char-
acterized by (N2 − N)/2 + N − 1 = (N + 2)(N − 1)/2
parameters. Combining the energy and dipole matrix pa-
rameters, there are a total of (N + 2)(N − 1)/2 +N − 1
parameters.
The N-state sum rules are of the form:
(Em + Ep)
〈m|x |n〉 〈n|x |p〉 (7)
δm,p,
so the sum rules comprise a total of N2 equations (i.e. an
equation for each (m, p)). If the sum rules are truncated
to N states, the sum rule indexed by (m = N, p = N)
is nonsensical because it contradicts the other sum rules.
Furthermore, if the transition moments are real, then
xmp = xpm, so only (N
2 − N)/2 of the equations are
independent. As such, there are a total of (N2 −N)/2+
N − 1 = (N + 2)(N − 1)/2 independent equations.
Since the SOS expression for the nonlinear-
susceptibility has (N + 2)(N − 1)/2 +N − 1 parameters
and the sum rules provide (N + 2)(N − 1)/2 equations,
the SOS expression can be reduced to a form with N − 1
parameters. For example, the three-level model for the
hyperpolarizability, which is expressed in terms of 7
parameters, can be reduced to two parameters using 5
sum rule equations. In practice, however, even fewer
sum rule equations are usually available because some
of them lead to physically unreasonable consequences.
While the (N,N) sum rule is clearly unphysical due
to truncation, sum rule equations that are near equa-
tion (N,N) may also be unphysical. In the case of
the three-level model, it is found that the equation
(2, 1) allows for an infinite hyperpolarizability, so that
equation is ignored on the grounds that it violates the
principle of physical soundness.[12, 25, 26] This leads to
a hyperpolarizability in terms of 3 variables, which are
chosen to be E10, E = E10/E20, and X = x10/x
The expression is then maximized with respect to the
two parameters E and X , leaving the final result a
function of E10.
We conclude that the SOS expression for the hyperpo-
larizability can be expressed in terms of at least N − 1
parameters; so, it would appear that as more levels are
included in the SOS expression, there are more free pa-
rameters that can be varied without violating the sum
rules. As N → ∞, there are an infinite number of ad-
justable parameters. So, it is indeed puzzling that the
three-level ansatz yields a fundamental limit that is con-
sistent with all of our calculations for a wide range of
potentials, many of which have many excited states. It
may be that we are only considering a small subset of
potential energy functions; or, perhaps the expression
for the hyperpolarizability depends on the parameters in
such a way that large matrix elements contribute to the
hyperpolarizability with alternating signs so that the big
terms cancel. This is a puzzle that needs to be solved if
we are to understand what makes β large.
To investigate whether the limiting behavior is
due to our use of 1-dimensional potentials, we have
also optimized the intrinsic hyperpolarizability in two-
dimensions. In this case, we focus on the largest ten-
sor component, βxxx and describe the potential as a
superposition of point charges. As described in the
literature,[19] we solve the two-dimensional Schrödinger
eigenvalue problem,
∇2Ψ+ VΨ = EΨ, (8)
for the lowest ten to 25 energy eigenstates, depending on
the degree of convergence of the resulting intrinsic hyper-
polarizability. We use the two-dimensional logarithmic
Coulomb potential, which for k nuclei with charges q1e,
. . . , qke located at points s
(1), . . . , s(k), is given by
V (s) =
qj log ‖s− s(j)‖, (9)
where L is a characteristic length. With L = 2Å, the
force due to a charge at distance 2Å is the same as it
would be for a 3D Coulomb potential.
We discretize the eigenvalue problem given by Equa-
tion 8 using a quadratic finite element method [21, 27]
and solve the resulting matrix eigenvalue problem for the
ten to 25 smallest energy eigenvalues and corresponding
eigenvectors by the implicitly-restarted Arnoldi method
[22] as implemented in ARPACK [28]. Each eigenvector
yields a wave function Ψn corresponding to energy level
En. The moments
xmn =
s1Ψm(s1, s2)Ψn(s1, s2) ds1ds2
are computed, and these and the energy levels En are
used to compute β
Figure 8 shows the intrinsic hyperpolarizability of a
two-nucleus molecule plotted as a function of the distance
between the two nuclei and nuclear charge q1. The total
nuclear charge is q1 + q2 = +e, and is expressed in units
of the proton charge, e. Three extrema are observed.
The positive peak parameters are βint = 0.649 for q1 =
0.58 and d = 4.36Å. The negative one yields βint =
−0.649 for q1 = 0.42 and d = 4.36Å. The local negative
peak that extends past the graph on the right reaches its
maximum magnitude of βint = −0.405 at q1 = 2.959 and
d = 2.0Å.
Applying numerical optimization to the intrinsic hy-
perpolarizability using the charges and separation be-
tween the nuclei as parameters, we get βint = 0.654 at
d = 4.539Å, and q1 = 0.430 when the starting param-
eters are near the positive peak; and βint = −0.651,
d = 4.443Å, and q1 = 0.572 when optimization gives
the negative peak. The peak parameters are the same
within roundoff errors when optimization or plotting is
used, confirming that the optimization procedure yields
the correct local extrema.
Figure 9 shows the intrinsic hyperpolarizability of an
octupolar-like molecule made of three evenly-spaced nu-
clei on a circle plotted as a function of the circle’s diam-
eter and charge fraction ǫ (q = ǫe). The charge on each
of the other nuclei is e(1 − ǫ)/2. The positive peak at
ǫ = 0.333 and diameter D = 6.9Å has a hyperpolariz-
ability βint = 0.326, while βint = −0.605 for a charge
fraction ǫ = 0.44 and a diameter D = 6.8Å.
Charge (q
Intrinsic Hyperpolarizability 
0 0.5 1 1.5 2
FIG. 8: The intrinsic hyperpolarizability of two nuclei as a
function of the distance between them and the charge of one
nucleus, q1 where q1 + q2 = +e.
Intrinsic Hyperpolarizability
0 0.2 0.4 0.6 0.8 1
FIG. 9: The intrinsic hyperpolarizability of three evenly-
spaced nuclei on a circle as a function of the circle’s diameter
and the charge, ǫ (in units of e), on one of the nuclei. The
charge on each of the other nuclei is e(1− ǫ)/2.
When the positions and magnitudes of the three
charges are allowed to move freely, the best intrinsic
hyperpolarizability obtained using numerical optimiza-
tion is βint = 0.685 for charges located at ~r1 = (0, 0).
~r2 = (−4.87Å, 0.33Å), and ~r3 = (−9.57Å,−0.16Å); with
charges q1 = 0.43e, q2 = 0.217e, and q3 = 0.351e. There
are only small differences in the optimized values of βint
depending on the starting positions and charges; and the
best results are for a “molecule” that is nearly linear
along the x-direction. This is not surprising given that
the xxx-component of βint is the optimized quantity.
The two-dimensional analysis illustrates that numer-
ical optimization correctly identifies the local maxima
(peaks and valleys) and that the magnitude of maximum
intrinsic hyperpolarizability (0.65 vs 0.68) is close to the
maximum we get for the one-dimensional optimization
of the potential energy function (0.71). All computa-
tions we have tried, including varying the potential en-
ergy function in one dimension or moving around point
charges in a plane all yield an intrinsic hyperpolarizabil-
ity that is less than 0.71.
An open question is the origin of the factor-of-thirty
gap between the best molecules and the fundamental
limit, which had remained firm for decades through the
year 2006. Several of the common proposed explana-
tions, such as vibronic dilution, have been eliminated.[14]
Perhaps it is not possible to make large-enough varia-
tions of the potential energy function without making
the molecule unstable. Or, perhaps there are subtle is-
sues with electron correlation, which prevents electrons
from responding to light with their full potential. The
fact that the idea of modulation of conjugation has lead
to a 50% increase over the long-standing ceiling - reduc-
ing the gap to a factor of twenty - makes it a promising
approach for potential further improvements. Continued
theoretical scrutiny, coupled with experiment, will be re-
quired to confirm the validity of our approach.
IV. CONCLUSIONS
There appear to be many potential energy functions
that lead to an intrinsic hyperpolarizability that is near
the fundamental limit. These separate into two broad
classes: one in which wiggles in the potential energy
function forces the eigenfunctions to be spatially sepa-
rated and a second class of monotonically varying wave-
functions with small or no wiggles that allow for many
strongly overlapping wavefunctions. Interestingly, all
these one-dimensional “molecules” have the same max-
imal intrinsic hyperpolarizability of 0.71. It is puzzling
that the three-level ansatz correctly predicts the funda-
mental limit even when the ansatz does not apply. A
second open question pertains to the origin of the long-
standing factor of 30 gap between the fundamental limit
and the best molecules. The idea of conjugation modu-
lation may be one promising approach for making wiggly
potential energy profiles that lead to molecules that fall
into the gap. Given that there are so many choices of
potential energy functions that lead to maximal intrinsic
hyperpolarizability, it may be possible to engineer many
new classes of exotic molecules with record intrinsic hy-
perpolarizability.
Acknowledgements: MGK thanks the National
Science Foundation (ECS-0354736) and Wright Paterson
Air Force Base for generously supporting this work.
[1] Q. Y. Chen, L. Kuang, Z. Y. Wang, and E. H. Sargent,
Nano. Lett. 4, 1673 (2004).
[2] B. H. Cumpston, S. P. Ananthavel, S. Barlow, D. L. Dyer,
J. E. Ehrlich, L. L. Erskine, A. A. Heikal, S. M. Kuebler,
I.-Y. S. Lee, D. McCord-Maughon, et al., Nature 398, 51
(1999).
[3] S. Kawata, H.-B. Sun, T. Tanaka, and K. Takada, Nature
412, 697 (2001).
[4] A. Karotki, M. Drobizhev, Y. Dzenis, P. N. Taylor, H. L.
Anderson, and A. Rebane, Phys. Chem. Chem. Phys. 6,
7 (2004).
[5] I. Roy, O. T. Y., H. E. Pudavar, E. J. Bergey, A. R.
Oseroff, J. Morgan, T. J. Dougherty, and P. N. Prasad,
J. Am. Chem. Soc. 125, 7860 (2003).
[6] M. G. Kuzyk, Opt. Lett. 25, 1183 (2000).
[7] M. G. Kuzyk, IEEE Journal on Selected Topics in Quan-
tum Electronics 7, 774 (2001).
[8] M. G. Kuzyk, Phys. Rev. Lett. 85, 1218 (2000).
[9] M. G. Kuzyk, Opt. Lett. 28, 135 (2003).
[10] M. G. Kuzyk, Phys. Rev. Lett. 90, 039902 (2003).
[11] M. G. Kuzyk, J. Nonl. Opt. Phys. & Mat. 13, 461 (2004).
[12] M. G. Kuzyk, J. Chem Phys. 125, 154108 (2006).
[13] M. G. Kuzyk, Optics & Photonics News 14, 26 (2003).
[14] K. Tripathi, P. Moreno, M. G. Kuzyk, B. J. Coe,
K. Clays, and A. M. Kelley, J. Chem. Phys. 121, 7932
(2004).
[15] B. J. Orr and J. F. Ward, Molecular Physics 20, 513
(1971).
[16] J. Pérez Moreno, Y. Zhao, K. Clays, and M. G. Kuzyk,
Opt. Lett. 32, 59 (2007).
[17] J. Zhou, M. Kuzyk, and D. S. Watkins, Opt. Lett. 31,
2891 (2006).
[18] M. G. Kuzyk, Phys. Rev. A 72, 053819 (2005).
[19] M. G. Kuzyk and D. S. Watkins, J. Chem Phys. 124,
244104 (2006).
[20] D. Kincaid and E. W. Cheney, Numerical Analy-
sis: Mathematics of Scientific Computing (Brooks-Cole,
2002), 3rd ed.
[21] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The
Finite Element Method: Its Basis and Fundamentals
(Butterworth-Heinemanm, 2005), 6th ed.
[22] D. C. Sorensen, SIAM J. Matrix Anal. Appl. 13, 357
(1992).
[23] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. Wright,
SIAM J. Optim. 9, 112 (1998).
[24] B. Champagne and B. Kirtman, Phys. Rev. Lett. 95,
109401 (2005).
[25] M. G. Kuzyk, Phys. Rev. Lett. 95, 109402 (2005).
[26] M. G. Kuzyk, J. Nonl. Opt. Phys. & Mat. 15, 77 (2006).
[27] K. Atkinson and W. Han, Theoretical Numerical Anal-
ysis, a Functional Analysis Framework (Springer, New
York, 2001).
[28] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK
Users’ Guide: Solution of Large-Scale Eigenvalue Prob-
lems with Implicitly Restarted Arnoldi Methods (SIAM,
Philadelphia, 1998).
0 5 10 15 20
0 5 10 15 20
0 5 10 15 20
